Polynomials.inl
上传用户:fxromeo
上传日期:2010-04-08
资源大小:89k
文件大小:5k
开发平台:

Visual C++

  1. //polynomials.inl 多项式与连分式头文件
  2. // Ver 1.0.0.0
  3. // 版权所有(C) 2002
  4. // 版权所有(C) 何渝, 2002
  5. // 最后修改: 2002.5.31.
  6. #ifndef _POLYNOMIALS_INL
  7. #define _POLYNOMIALS_INL
  8. //求一维实(复)多项式值
  9. template<class T, class U>
  10. inline U  //内联函数
  11. PolyValueOneDim(valarray<T>& dCoff, size_t stNo, U dX)
  12. { //dCoff为一维多项式系数数组,stNo为多项式次数,dX为自变量x
  13. U dValue = dCoff[stNo-1];
  14. for(int st = stNo - 2; st > -1; st--)
  15. dValue = dValue * dX + dCoff[st];
  16. return(dValue); //返回多项式值
  17. }
  18. //求一维多项式组值
  19. template<class T, class V, class U>
  20. inline void 
  21. PolyValueOneDimGroup(valarray<T>& dCoff, valarray<V>& dX, valarray<U>& dValue)
  22. { //dCoff为一维多项式系数数组,dX存放自变量,dValue存放对应多项式值
  23.    int i,j,mm,nn,ll,t,s,kk,k;
  24.     V y, z;
  25.     int n = dCoff.size(); //多项式项数,其最高次数为n-1
  26. int m = dX.size(); //给定自变量个数
  27. valarray<V> b(2*n);
  28. //b=malloc(2*n*sizeof(double));
  29.     y=dCoff[n-1];
  30.     for(i=0; i<n; i++) b[i] = dCoff[i]/y;
  31.     k = log(n-0.5) / log(2.0) + 1;
  32.     nn = 1;
  33.     for(i=0; i<k; i++) nn=2*nn;
  34.     for(i=n; i<nn-1; i++) b[i]=0.0;
  35.     b[nn-1] = 1.0;
  36.     t = nn;
  37. s = 1;
  38.     for(i=1; i<=k-1; i++)
  39.     {
  40. t = t / 2;
  41. mm = -t;
  42.         for(j=1; j<=s; j++)
  43.         {
  44. mm = mm + t + t;
  45. b[mm-1] -= 1.0; 
  46.             for(kk=2; kk<=t; kk++)
  47.               b[mm-kk] = b[mm-kk] - b[mm-1] * b[mm+t-kk];
  48.         }
  49.         s = s + s;
  50.     }
  51.     for(kk=1; kk<=m; kk++)
  52.       { 
  53. for(i=0; i<=(nn-2)/2; i++)
  54.            dCoff[i] = dX[kk-1] + b[2*i];
  55.         mm = 1;
  56. z = dX[kk-1];
  57.         for(i=1; i<=k-1; i++)
  58.         {
  59. mm = mm + mm; 
  60. ll = mm + mm;
  61. z = z * z;
  62.             for(j=0; j<nn; j=j+ll)
  63.               dCoff[j/2]=dCoff[j/2]+dCoff[(j+mm)/2]*(z+b[j+mm-1]);
  64.         }
  65.         z = z * z / dX[kk-1];
  66.         if(nn!=n) dCoff[0] -= z;
  67.         dValue[kk-1] = dCoff[0] * y;
  68.     }
  69. }
  70. //求二维多项式值
  71. template<class T, class U>
  72. inline U //内联函数
  73. PolyValueTwoDim(matrix<T>& dCoff, U dX, U dY)
  74. { //dCoff为二维多项式系数数组,dX为自变量x,dY为自变量y
  75. U dValue = 0, dTemp = 1;
  76. size_t styNo = dCoff.GetColNum();//二维多项式自变量y的项数,其最高次数为styNo-1
  77. size_t stxNo = dCoff.GetRowNum();//二维多项式自变量x的项数,其最高次数为stxNo-1
  78. for(size_t st = 0; st < stxNo; st++)
  79. {
  80. U dUS = dCoff(st, styNo-1) * dTemp;
  81. for(int sr = styNo-2; sr > -1; sr--)
  82. dUS = dUS * dY + dCoff(st, sr) * dTemp;
  83. dValue = dValue + dUS;
  84. dTemp = dTemp * dX;
  85. }
  86. return(dValue); //返回多项式值
  87. }
  88. //两一维多项式相乘
  89. template<class T>
  90. inline void
  91. PolyMultip(valarray<T>& dCoffP, valarray<T>& dCoffQ, valarray<T>& dCoffS)
  92. { //dCoffP被乘多项式P(x)系数数组,最高次数为stNoP-1
  93. //dCoffQ乘多项式Q(x)系数数组,最高次数为stNoQ-1
  94. //dCoffS积多项式S(x)系数数组,最高次数为stNoP+stNoQ-1
  95. //最后结果得到积多项式S(x)的系数
  96. size_t stNoP = dCoffP.size();//一维多项式P(x)项数,最高次数为stNoP-1
  97. size_t stNoQ = dCoffQ.size();//一维多项式Q(x)项数,最高次数为stNoQ-1
  98. size_t stNoS = dCoffS.size();//一维多项式S(x)项数,最高次数为stNoS-1,stNoS=stNoP+stNoQ-1
  99. for(size_t st = 0; st < stNoS; st++) dCoffS[st]=0.0;
  100. for(st = 0; st < stNoP; st++)
  101. for(size_t sr = 0; sr < stNoQ; sr++)
  102. dCoffS[st+sr]=dCoffS[st+sr]+dCoffP[st]*dCoffQ[sr];
  103. }
  104. //两一维多项式除法
  105. template<class T>
  106. inline int
  107. PolyDiv(valarray<T>& dCoffP, valarray<T>& dCoffQ, valarray<T>& dCoffS, valarray<T>& dCoffR)
  108. { //dCoffP被除多项式P(x)系数数组,最高次数为stNoP-1
  109. //dCoffQ除多项式Q(x)系数数组,最高次数为stNoQ-1
  110. //dCoffS商多项式S(x)系数数组,最高次数为stNoP-stNoQ
  111. //dCoffR余多项式R(x)系数数组,最高次数为stNoR-1
  112. //最后得到的结果是商多项式S(x)和余多项式R(x)的系数
  113. size_t stNoP = dCoffP.size();//一维多项式P(x)项数,最高次数为stNoP-1
  114. size_t stNoQ = dCoffQ.size();//一维多项式Q(x)项数,最高次数为stNoQ-1
  115. size_t stNoS = dCoffS.size();//一维多项式S(x)项数,最高次数为stNoS-stNoQ
  116. size_t stNoR = dCoffR.size();//一维多项式R(x)项数,最高次数为stNoR-1
  117. for(size_t st = 0; st < stNoS; st++) dCoffS[st] = 0.0;
  118. if(FloatEqual(Abs(dCoffQ[stNoQ-1]),0)) return (0);
  119. size_t stk = stNoP - 1;
  120. for(st = stNoS; st > 0; st--)
  121. {
  122. dCoffS[st-1] = dCoffP[stk] / dCoffQ[stNoQ-1];
  123. size_t stm = stk;
  124. for(size_t sr = 1; sr < stNoQ; sr++)
  125. {
  126. dCoffP[stm-1] = dCoffP[stm-1] - dCoffS[st-1] * dCoffQ[stNoQ-sr-1];
  127. stm = stm-1;
  128. }
  129. stk = stk - 1;
  130. }
  131. for(st = 0; st < stNoR; st++) dCoffR[st] = dCoffP[st];
  132. return (1);
  133. }
  134. // 计算连分式函数值
  135. template<class T>
  136. inline T
  137. FractionValue(valarray<T>& dXpara, valarray<T>& dCoff, T dX)
  138. {
  139. size_t stNo = dXpara.size(); //连分式的节数
  140. T dValue = dCoff[stNo-1];
  141. for(int st=stNo-2; st>-1; st--)
  142. {
  143. if (FloatEqual(dValue,0))
  144. dValue = 1.0e+35 * (dX-dXpara[st]) / Abs(dX-dXpara[st]);
  145.         else
  146. dValue = dCoff[st] + (dX-dXpara[st]) / dValue;
  147.     }
  148.     return(dValue); //返回连分式值
  149. }
  150. #endif //_POLYNOMIALS_INL