开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 查看源码
SpecialFunction.inl
资源名称：C++常用算法之06矩阵特征值与特征向量的计算.rar　[点击查看]
上传用户：fxromeo
上传日期：2010-04-08
资源大小：89k
文件大小：18k
源码类别：

数值算法/人工智能

开发平台：

Visual C++

SpecialFunction.inl：源码内容
  1. //SpecialFunction.inl 特殊函数定义头文件
  2. // Ver 1.0.0.0
  3. // 版权所有(C) 何渝, 2002
  4. // 最后修改: 2002.5.31.
  7. #ifndef _SPECIALFUNCTION_INL //避免多次编译
  8. #define _SPECIALFUNCTION_INL
  10. //template <class _Ty = float>
  12. //伽马函数
  13. template <class _Ty>
  14. _Ty GammaFunction(_Ty x)
  15. {
  16.     _Ty t;
  17.     _Ty a[11] = 
  18. {
  19.  0.0000677106, -0.0003442342,  0.0015397681,
  20. -0.0024467480,  0.0109736958, -0.0002109075,
  21.  0.0742379071,  0.0815782188,  0.4118402518,
  22.  0.4227843370,  1.0
  23. };
  25.     if(x<0.0 || FloatEqual(x,0.0))
  26.     {
  27. cout << "Error **x<=0!" << endl; 
  28. return(-1.0);
  29. }
  30.     _Ty y = x;
  31.     if(y<1.0 || FloatEqual(y,1.0))
  32.     {
  33. t=1.0/(y*(y+1.0)); 
  34. y=y+2.0;
  35. }
  36.     else
  37. if(y<2.0 || FloatEqual(y,2.0))
  38. {
  39. t=1.0/y; 
  40. y=y+1.0;
  41. }
  42. else
  43. if(y<3.0 || FloatEqual(y,3.0)) t=1.0;
  44. else
  45. {
  46. t=1.0;
  47. while(y>3.0)
  48. {
  49. y=y-1.0;
  50. t=t*y;
  51. }
  52. }
  53.     _Ty s=a[0]; 
  54. _Ty u = y - 2.0;
  55.     for(int i=1; i<11; i++) s = s * u + a[i];
  56.     s = s * t;
  57.     return(s);
  58. }
  60. //不完全伽马函数
  61. template <class _Ty>
  62. _Ty IncompleteGammaFunction(_Ty a, _Ty x)
  63. {
  64. int n;
  65.     _Ty p,d,s,s1,p0,q0,p1,q1;
  66.     if((a<0.0 || FloatEqual(a,0.0)) || (x<0.0))
  67.     { 
  68. if(a<0.0 || FloatEqual(a,0.0)) cout <<"Error **a<=0!" <<endl;
  69.         if(x<0.0) cout << "Error **x<0!" << endl;
  70.         return(-1.0);
  71.     }
  72.     if(FloatEqual(x,0.0)) return(0.0);
  73.     if(x>1.0e+35) return(1.0);
  74.     _Ty q = log(x); 
  75. q = a * q; 
  76. _Ty qq = exp(q);
  77.     if(x<1.0+a)
  78.     {
  79. p=a; 
  80. d=1.0/a; 
  81. s=d;
  82.         for(n=1; n<101; n++)
  83.         {
  84. p=1.0+p;
  85. d=d*x/p;
  86. s=s+d;
  87. if(Abs(d)<Abs(s)*FLOATERROR)
  88.             {
  89. s=s*exp(-x)*qq/GammaFunction(a);
  90.                 return(s);
  91.             }
  92.         }
  93.     }
  94.     else
  95.     { 
  96. s=1.0/x; 
  97. p0=0.0; 
  98. p1=1.0;
  99. q0=1.0;
  100. q1=x;
  101.         for(n=1; n<101; n++)
  102.         { 
  103. p0=p1+(n-a)*p0;
  104. q0=q1+(n-a)*q0;
  105.             p=x*p0+n*p1;
  106. q=x*q0+n*q1;
  107.             if(FloatNotEqual(q,0.0))
  108.             {
  109. s1=p/q; 
  110. p1=p; 
  111. q1=q;
  112.                 if(Abs((s1-s)/s1)<FLOATERROR)
  113.                 { 
  114. s=s1*exp(-x)*qq/GammaFunction(a);
  115.                     return(1.0-s);
  116.                 }
  117.                 s=s1;
  118.             }
  119.             p1=p; 
  120. q1=q;
  121.         }
  122.     }
  123.     cout << "a too large !" << endl;
  124.     s=1.0-s*exp(-x)*qq/GammaFunction(a);
  125.     return(s);
  126. }
  128. //误差函数
  129. template <class _Ty>
  130. _Ty ErrorFunction(_Ty x)
  131. {
  132. _Ty y;
  134.     if(x>0.0 || FloatEqual(x,0.0))
  135.       y=IncompleteGammaFunction(0.5,x*x);
  136.     else
  137.       y=-IncompleteGammaFunction(0.5,x*x);
  138.     return(y);
  139. }
  141. //第一类整数阶贝塞尔函数
  142. template <class _Ty>
  143. _Ty IntegerBessel1stFunction(int n, _Ty x)
  144. {
  145. int i,m;
  146.     _Ty y,z,p,q,s,b0,b1;
  147.     _Ty a[6]={ 57568490574.0,-13362590354.0,
  148.              651619640.7,-11214424.18,77392.33017,-184.9052456};
  149.     _Ty b[6]={ 57568490411.0,1029532985.0,
  150.              9494680.718,59272.64853,267.8532712,1.0};
  151.     _Ty c[6]={ 72362614232.0,-7895059235.0,
  152.              242396853.1,-2972611.439,15704.4826,-30.16036606};
  153.     _Ty d[6]={ 144725228443.0,2300535178.0,
  154.              18583304.74,99447.43394,376.9991397,1.0};
  155.     _Ty e[5]={ 1.0,-0.1098628627e-02,
  156.              0.2734510407e-04,-0.2073370639e-05,0.2093887211e-06};
  157.     _Ty f[5]={ -0.1562499995e-01,
  158.              0.1430488765e-03,-0.6911147651e-05,
  159.              0.7621095161e-06,-0.934935152e-07};
  160.     _Ty g[5]={ 1.0,0.183105e-02,
  161.              -0.3516396496e-04,0.2457520174e-05, -0.240337019e-06};
  162.     _Ty h[5]={ 0.4687499995e-01,
  163.              -0.2002690873e-03,0.8449199096e-05,
  164.              -0.88228987e-06,0.105787412e-06};
  165.     _Ty t=Abs(x);
  166.     if(n<0) n=-n;
  167.     if(n!=1)
  168.     { 
  169. if(t<8.0)
  170.         {
  171. y=t*t;
  172. p=a[5];
  173. q=b[5];
  174. for(i=4; i>=0; i--)
  175.             { 
  176. p=p*y+a[i];
  177. q=q*y+b[i];
  178. }
  179. p=p/q;
  180.         }
  181.         else
  182.         { 
  183. z=8.0/t;
  184. y=z*z;
  185.             p=e[4];
  186. q=f[4];
  187.             for(i=3; i>=0; i--)
  188.             {
  189. p=p*y+e[i];
  190. q=q*y+f[i];
  191. }
  192.             s=t-0.785398164;
  193.             p=p*cos(s)-z*q*sin(s);
  194.             p=p*sqrt(0.636619772/t);
  195.         }
  196.     }
  197.     if(n==0) return(p);
  198.     b0=p;
  199.     if(t<8.0)
  200.     { 
  201. y=t*t; 
  202. p=c[5];
  203. q=d[5];
  204.         for(i=4; i>=0; i--)
  205.         {
  206. p=p*y+c[i];
  207. q=q*y+d[i];
  208. }
  209.         p=x*p/q;
  210.     }
  211.     else
  212.     {
  213. z=8.0/t; 
  214. y=z*z;
  215.         p=g[4];
  216. q=h[4];
  217.         for(i=3; i>=0; i--)
  218.         {
  219. p=p*y+g[i];
  220. q=q*y+h[i];
  221. }
  222.         s=t-2.356194491;
  223.         p=p*cos(s)-z*q*sin(s);
  224.         p=p*x*sqrt(0.636619772/t)/t;
  225.     }
  226.     if(n==1) return(p);
  227.     b1=p;
  228.     if(x==0.0) return(0.0);
  229.     s=2.0/t;
  230.     if(t>1.0*n)
  231.     {
  232. if(x<0.0) b1=-b1;
  233.         for(i=1; i<n; i++)
  234.         { 
  235. p=s*i*b1-b0; 
  236. b0=b1; 
  237. b1=p;
  238. }
  239.     }
  240.     else
  241.     {
  242. m=(n+(int)sqrt(40.0*n))/2;
  243.         m=2*m;
  244.         p = q = b1 = 0.0;
  245. b0=1.0;
  246.         for(i=m-1; i>=0; i--)
  247.         { 
  248. t=s*(i+1)*b0-b1;
  249.             b1=b0;
  250. b0=t;
  251.             if(Abs(b0)>1.0e+10)
  252.             {
  253. b0=b0*1.0e-10; 
  254. b1=b1*1.0e-10;
  255.                 p=p*1.0e-10; 
  256. q=q*1.0e-10;
  257.             }
  258.             if((i+2)%2==0) q=q+b0;
  259.             if((i+1)==n) p=b1;
  260.         }
  261.         q=2.0*q-b0;
  262. p=p/q;
  263.     }
  264.     if((x<0.0)&&(n%2==1)) p=-p;
  265.     return(p);
  266. }
  268. //第二类整数阶贝塞尔函数
  269. template <class _Ty>
  270. _Ty IntegerBessel2ndFunction(int n, _Ty x)
  271. {
  272. int i;
  273.     _Ty y,z,p,q,s,b0,b1;
  275.     _Ty a[6]={ -2.957821389e+9,7.062834065e+9,
  276.              -5.123598036e+8,1.087988129e+7,-8.632792757e+4,
  277.              2.284622733e+2};
  278.     _Ty b[6]={ 4.0076544269e+10,7.452499648e+8,
  279.            7.189466438e+6,4.74472647e+4,2.261030244e+2,1.0};
  280.     _Ty c[6]={ -4.900604943e+12,1.27527439e+12,
  281.             -5.153438139e+10,7.349264551e+8,-4.237922726e+6,
  282.              8.511937935e+3};
  283.     _Ty d[7]={ 2.49958057e+13,4.244419664e+11,
  284.             3.733650367e+9,2.245904002e+7,1.02042605e+5,
  285.             3.549632885e+2,1.0};
  286.     _Ty e[5]={ 1.0,-0.1098628627e-02,
  287.              0.2734510407e-04,-0.2073370639e-05,
  288.              0.2093887211e-06};
  289.     _Ty f[5]={ -0.1562499995e-01,
  290.              0.1430488765e-03,-0.6911147651e-05,
  291.              0.7621095161e-06,-0.934935152e-07};
  292.     _Ty g[5]={ 1.0,0.183105e-02,
  293.              -0.3516396496e-04,0.2457520174e-05,
  294.              -0.240337019e-06};
  295.     _Ty h[5]={ 0.4687499995e-01,
  296.              -0.2002690873e-03,0.8449199096e-05,
  297.              -0.88228987e-06,0.105787412e-06};
  298.     if(n<0) n=-n;
  299.     if(x<0.0) x=-x;
  300.     if(FloatEqual(x,0.0)) return(-1.0e+70);
  301.     if(n!=1)
  302.     { 
  303. if(x<8.0)
  304.         {
  305. y=x*x; 
  306. p=a[5]; 
  307. q=b[5];
  308. for(i=4; i>=0; i--)
  309.             {
  310. p=p*y+a[i];
  311. q=q*y+b[i];}
  312. p=p/q+0.636619772*IntegerBessel1stFunction(0,x)*log(x);
  313.         }
  314.         else
  315.         { 
  316. z=8.0/x;
  317. y=z*z;
  318.             p=e[4];
  319. q=f[4];
  320.             for(i=3; i>=0; i--)
  321.             {
  322. p=p*y+e[i];
  323. q=q*y+f[i];
  324. }
  325. s=x-0.785398164;
  326.             p=p*sin(s)+z*q*cos(s);
  327.             p=p*sqrt(0.636619772/x);
  328.         }
  329.     }
  330.     if(n==0) return(p);
  331.     b0=p;
  332.     if(x<8.0)
  333.     {
  334. y=x*x; 
  335. p=c[5];
  336. q=d[6];
  337.         for(i=4; i>=0; i--)
  338.         {
  339. p=p*y+c[i];
  340. q=q*y+d[i+1];
  341. }
  342.         q=q*y+d[0];
  343.         p=x*p/q+0.636619772*(IntegerBessel1stFunction(1,x)*log(x)-1.0/x);;
  344.     }
  345.     else
  346.     { 
  347. z=8.0/x;
  348. y=z*z;
  349.         p=g[4]; 
  350. q=h[4];
  351.         for(i=3; i>=0; i--)
  352.         {
  353. p=p*y+g[i];
  354. q=q*y+h[i];
  355. }
  356.         s=x-2.356194491;
  357.         p=p*sin(s)+z*q*cos(s);
  358.         p=p*sqrt(0.636619772/x);
  359.     }
  360.     if(n==1) return(p);
  361.     b1=p;
  362.     s=2.0/x;
  363.     for(i=1; i<n; i++)
  364.     { 
  365. p=s*i*b1-b0; 
  366. b0=b1; 
  367. b1=p;
  368. }
  369.     return(p);
  370. }
  372. //变形第一类整数阶贝塞尔函数
  373. template <class _Ty>
  374. _Ty TransformativeIntegerBessel1stFunction(int n,_Ty x)
  375. {
  376. int i,m;
  377.     _Ty t,y,p,b0,b1,q;
  378.     _Ty a[7]={ 1.0,3.5156229,3.0899424,1.2067492,
  379.                          0.2659732,0.0360768,0.0045813};
  380.     _Ty b[7]={ 0.5,0.87890594,0.51498869,
  381.               0.15084934,0.02658773,0.00301532,0.00032411};
  382.     _Ty c[9]={ 0.39894228,0.01328592,0.00225319,
  383.                         -0.00157565,0.00916281,-0.02057706,
  384.                          0.02635537,-0.01647633,0.00392377};
  385.     _Ty d[9]={ 0.39894228,-0.03988024,-0.00362018,
  386.                         0.00163801,-0.01031555,0.02282967,
  387.                         -0.02895312,0.01787654,-0.00420059};
  388.     if(n<0) n=-n;
  389.     t=Abs(x);
  390.     if(n!=1)
  391.     { 
  392. if(t<3.75)
  393.         {
  394. y=(x/3.75)*(x/3.75);
  395. p=a[6];
  396.             for(i=5; i>=0; i--) p=p*y+a[i];
  397.         }
  398.         else
  399.         { 
  400. y=3.75/t;
  401. p=c[8];
  402.             for(i=7; i>=0; i--) p=p*y+c[i];
  403.             p=p*exp(t)/sqrt(t);
  404.         }
  405.     }
  406.     if(n==0) return(p);
  407.     q=p;
  408.     if(t<3.75)
  409.     {
  410. y=(x/3.75)*(x/3.75); 
  411. p=b[6];
  412.         for(i=5; i>=0; i--) p=p*y+b[i];
  413.         p=p*t;
  414.     }
  415.     else
  416.     { 
  417. y=3.75/t;
  418. p=d[8];
  419.         for(i=7; i>=0; i--) p=p*y+d[i];
  420.         p=p*exp(t)/sqrt(t);
  421.     }
  422.     if(x<0.0) p=-p;
  423.     if(n==1) return(p);
  424.     if(FloatEqual(x,0.0)) return(0.0);
  425.     y=2.0/t;
  426. t = b0 = 0.0;
  427. b1=1.0;
  428.     m=n+(int)sqrt(40.0*n);
  429.     m=2*m;
  430.     for(i=m; i>0; i--)
  431.     { 
  432. p=b0+i*y*b1; 
  433. b0=b1;
  434. b1=p;
  435.         if(Abs(b1)>1.0e+10)
  436.         {
  437. t=t*1.0e-10; 
  438. b0=b0*1.0e-10;
  439.             b1=b1*1.0e-10;
  440.         }
  441.         if(i==n) t=b0;
  442.     }
  443.     p=t*q/b1;
  444.     if((x<0.0)&&(n%2==1)) p=-p;
  445.     return(p);
  446. }
  448. //变形第二类整数阶贝塞尔函数
  449. template <class _Ty>
  450. _Ty TransformativeIntegerBessel2ndFunction(int n, _Ty x)
  451. {
  452. int i;
  453.     _Ty y,p,b0,b1;
  455.     _Ty a[7]={ -0.57721566,0.4227842,0.23069756,
  456.            0.0348859,0.00262698,0.0001075,0.0000074};
  457.     _Ty b[7]={ 1.0,0.15443144,-0.67278579,
  458.            -0.18156897,-0.01919402,-0.00110404,-0.00004686};
  459.     _Ty c[7]={ 1.25331414,-0.07832358,0.02189568,
  460.            -0.01062446,0.00587872,-0.0025154,0.00053208};
  461.     _Ty d[7]={ 1.25331414,0.23498619,-0.0365562,
  462.            0.01504268,-0.00780353,0.00325614,-0.00068245};
  463.     if(n<0) n=-n;
  464.     if(x<0.0) x=-x;
  465.     if(FloatEqual(x,0.0)) return(1.0e+70);
  466.     if(n!=1)
  467.     { 
  468. if(x < 2.0 || FloatEqual(x, 2.0))
  469. { 
  470. y=x*x/4.0; p=a[6];
  471.             for(i=5; i>=0; i--) p=p*y+a[i];
  472.             p=p-TransformativeIntegerBessel1stFunction(0,x)*log(x/2.0);
  473.         }
  474.         else
  475.         { 
  476. y=2.0/x; p=c[6];
  477.             for(i=5; i>=0; i--) p=p*y+c[i];
  478.             p=p*exp(-x)/sqrt(x);
  479.         }
  480.     }
  481.     if(n==0) return(p);
  482.     b0=p;
  483.     if(x < 2.0 || FloatEqual(x, 2.0))
  484.     { 
  485. y=x*x/4.0;
  486. p=b[6];
  487.         for(i=5; i>=0; i--) p=p*y+b[i];
  488.         p=p/x+TransformativeIntegerBessel1stFunction(1,x)*log(x/2.0);
  489.     }
  490.     else
  491.     { 
  492. y=2.0/x;
  493. p=d[6];
  494.         for(i=5; i>=0; i--) p=p*y+d[i];
  495.         p=p*exp(-x)/sqrt(x);
  496.     }
  497.     if(n==1) return(p);
  498.     b1=p;
  499.     y=2.0/x;
  500.     for(i=1; i<n; i++)
  501.     { 
  502. p=b0+i*y*b1;
  503. b0=b1; 
  504. b1=p;
  505. }
  506.     return(p);
  507. }
  509. //不完全贝塔函数
  510. template <class _Ty>
  511. _Ty IncompleteBetaFunction(_Ty a, _Ty b, _Ty x)
  512. {
  513. _Ty y;
  515.     if(a<0.0 || FloatEqual(a,0.0))
  516.     { 
  517. cout << "err**a<=0!";
  518. return(-1.0);
  519. }
  520.     if(b<0.0 || FloatEqual(b,0.0))
  521.     { 
  522. cout << "err**b<=0!";
  523. return(-1.0);
  524. }
  525.     if((x<0.0)||(x>1.0))
  526.     {
  527. cout << "err**x<0 or x>1 !";
  528.         return(1.0e+70);
  529.     }
  530.     if(FloatEqual(x,0.0)||FloatEqual(x,1.0)) y=0.0;
  531.     else
  532.     { 
  533. y=a*log(x)+b*log(1.0-x);
  534.         y=exp(y);
  535.         y=y*GammaFunction(a+b)/(GammaFunction(a)*GammaFunction(b));
  536.     }
  537.     if(x<(a+1.0)/(a+b+2.0))
  538.       y=y*beta(a,b,x)/a;
  539.     else
  540.       y=1.0-y*beta(b,a,1.0-x)/b;
  541.     return(y);
  542. }
  544. //不完全贝塔函数的从属函数
  545. template <class _Ty>
  546. _Ty beta(_Ty a, _Ty b, _Ty x)
  547. {
  548.     _Ty p0(0.0), q0(1.0), p1(1.0), q1(1.0), s1;
  550.     for(int k=1; k<101; k++)
  551.     {
  552. _Ty d=(a+k)*(a+b+k)*x;
  553. d=-d/((a+k+k)*(a+k+k+1.0));
  554.         p0=p1+d*p0;
  555. q0=q1+d*q0;
  556. _Ty s0=p0/q0;
  557.         d=k*(b-k)*x;
  558.         d=d/((a+k+k-1.0)*(a+k+k));
  559.         p1=p0+d*p1;
  560. q1=q0+d*q1;
  561. s1=p1/q1;
  562.         if(Abs(s1-s0)<Abs(s1)*1.0e-07) return(s1);
  563.     }
  564.     cout << "a or b too big !";
  565.     return(s1);
  566. }
  568. //正态分布函数
  569. template <class _Ty>
  570. _Ty NormalDistributionFunction(_Ty a, _Ty d, _Ty x)
  571. {
  572. _Ty y;
  573.     
  574.     if(d<0.0 || FloatEqual(d,0.0)) d=1.0e-10;
  575.     y=0.5+0.5*ErrorFunction((x-a)/(sqrt(2.0)*d));
  576.     return(y);
  577. }
  579. //t-分布函数
  580. template <class _Ty>
  581. _Ty tDistributionFunction(_Ty t, int n)
  582. {
  583. _Ty y;
  584.     
  585.     if(t<0.0) t = -t;
  586.     y=1.0-IncompleteBetaFunction(n/2.0, 0.5, n/(n+t*t));
  587.     return(y);
  588. }
  590. //X^2-分布函数
  591. template <class _Ty>
  592. _Ty X2DistributionFunction(_Ty x, int n)
  593. {
  594. _Ty y;
  595.     
  596.     if(x<0.0) x=-x;
  597.     y=IncompleteGammaFunction(n/2.0,x/2.0);
  598.     return(y);
  599. }
  601. //F-分布函数
  602. template <class _Ty>
  603. _Ty FDistributionFunction(_Ty f, int n1, int n2)
  604. {
  605. _Ty y;
  606.     //extern double lbeta();
  607.     if(f<0.0) f=-f;
  608.     y=IncompleteBetaFunction(n2/2.0,n1/2.0,n2/(n2+n1*f));
  609.     return(y);
  610. }
  612. //正弦积分
  613. template <class _Ty>
  614. _Ty SineIntegralFunction(_Ty x)
  615. {
  616. int m(1), i, j;
  617.     _Ty aa,bb,w,xx;
  618.     _Ty t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459};
  619.     _Ty c[5]={0.2369268851,0.4786286705,0.5688888889,
  620. 0.4786286705,0.2369268851};
  622.     if(FloatEqual(x,0.0)) return(0.0);
  623.     _Ty h=Abs(x); 
  624. _Ty s=Abs(0.0001*h);
  625.     _Ty p(1.0e+35), ep(FLOATERROR), g(0.0);
  626.     while((ep>0.0000001 || FloatEqual(ep,0.0000001)) && (Abs(h)>s))
  627.     { 
  628. g=0.0;
  629.         for(i=1;i<=m;i++)
  630.         { 
  631. aa=(i-1.0)*h; 
  632. bb=i*h;
  633.             w=0.0;
  634.             for(j=0;j<5;j++)
  635.             {
  636. xx=((bb-aa)*t[j]+(bb+aa))/2.0;
  637.                 w=w+sin(xx)/xx*c[j];
  638.             }
  639.             g=g+w;
  640.         }
  641.         g=g*h/2.0;
  642.         ep=Abs(g-p)/(1.0+Abs(g));
  643.         p=g;
  644. m=m+1;
  645. h=Abs(x)/m;
  646.     }
  647.     return(g);
  648. }
  650. //余弦积分
  651. template <class _Ty>
  652. _Ty CosineIntegralFunction(_Ty x)
  653. {
  654. int m(1), i, j;
  655.     _Ty ep(FLOATERROR),aa,bb,w,xx,g(0.0);
  656.     _Ty t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459};
  657.     _Ty c[5]={0.2369268851,0.4786286705,0.5688888889,
  658. 0.4786286705,0.2369268851};
  659.     
  660.     if(FloatEqual(x,0.0)) x = 1.0e-35;
  661.     if(x<0.0) x = -x;
  662.     _Ty r = 0.57721566490153286060651;
  663.     _Ty q = r + log(x);
  664.     _Ty h = x; 
  665. _Ty s = Abs(0.0001*h);
  666.     _Ty p=1.0e+35;
  667.     while((ep>0.0000001||FloatEqual(ep,0.0000001))&&(Abs(h)>s))
  668.     { 
  669. g=0.0;
  670.         for(i=1;i<=m;i++)
  671.         { 
  672. aa=(i-1.0)*h; 
  673. bb=i*h;
  674.             w=0.0;
  675.             for(j=0;j<5;j++)
  676.             {
  677. xx=((bb-aa)*t[j]+(bb+aa))/2.0;
  678.                 w=w+(1.0-cos(xx))/xx*c[j];
  679.             }
  680.             g=g+w;
  681.         }
  682.         g=g*h/2.0;
  683.         ep=Abs(g-p)/(1.0+Abs(g));
  684.         p=g; m=m+1; h=x/m;
  685.     }
  686.     g=q-g;
  687.     return(g);
  688. }
  690. //指数积分
  691. template <class _Ty>
  692. _Ty ExponentIntegralFunction(_Ty x)
  693. {
  694. int m(1), i, j;
  695.     _Ty ep(FLOATERROR),aa,bb,w,xx,g(0.0);
  696.     _Ty t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459};
  697.     _Ty c[5]={0.2369268851,0.4786286705,0.5688888889,
  698. 0.4786286705,0.2369268851};
  699.     
  700.     if(FloatEqual(x,0.0)) x=1.0e-10;
  701.     if(x<0.0) x=-x;
  702.     _Ty r=0.57721566490153286060651;
  703.     _Ty q=r+log(x);
  704.     _Ty h=x;
  705. _Ty s=Abs(0.0001*h);
  706.     _Ty p=1.0e+35;
  707.     while((ep>0.0000001||FloatEqual(ep,0.0000001))&&(Abs(h)>s))
  708.     { 
  709. g=0.0;
  710.         for(i=1;i<=m;i++)
  711.         {
  712. aa=(i-1.0)*h; bb=i*h;
  713.             w=0.0;
  714.             for(j=0;j<5;j++)
  715.             { 
  716. xx=((bb-aa)*t[j]+(bb+aa))/2.0;
  717.                 w=w+(exp(-xx)-1.0)/xx*c[j];
  718.             }
  719.             g=g+w;
  720.         }
  721.         g=g*h/2.0;
  722.         ep=Abs(g-p)/(1.0+Abs(g));
  723.         p=g; m=m+1; h=x/m;
  724.     }
  725.     g=q+g;
  726.     return(g);
  727. }
  729. //第一类椭圆积分
  730. template <class _Ty>
  731. _Ty Ellipse1stIntegral(_Ty k, _Ty f)
  732. {
  733. int n(0);
  734.     _Ty pi(3.1415926), e(1.0), ff;
  735.     //_Ty fk();
  736.     if(k<0.0) k=-k;
  737.     if(k>1.0) k=1.0/k;
  738.     _Ty y = Abs(f);
  739.     
  740.     while(y>pi||FloatEqual(y,pi))
  741.     {
  742. n=n+1;
  743. y=y-pi;
  744. }
  745.     
  746.     if(y>pi/2.0||FloatEqual(y,pi/2.0))
  747.     {
  748. n=n+1;
  749. e=-e; 
  750. y=pi-y;
  751. }
  752.     if(n==0)
  753. ff=fk(k,y);
  754.     else
  755.     { 
  756. ff=fk(k,pi/2.0);
  757.         ff=2.0*n*ff+e*fk(k,y);
  758.     }
  759.     if(f<0.0) ff=-ff;
  760.     return(ff);
  761. }
  763. //第一类椭圆积分辅助函数
  764. template <class _Ty>
  765. _Ty fk(_Ty k, _Ty f)
  766. { 
  767. int m(1), i, j;
  768.     _Ty p(1.0e+35),ep(FLOATERROR),aa,bb,w,xx,g(0.0),q;
  769.     _Ty t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459};
  770.     _Ty c[5]={0.2369268851,0.4786286705,0.5688888889,
  771. 0.4786286705,0.2369268851};
  772.     
  773.     _Ty h=Abs(f); 
  774. _Ty s=Abs(0.0001*h);
  775.     //p=1.0e+35; ep=0.000001;
  776.     while((ep>0.0000001||FloatEqual(ep,0.0000001))&&(Abs(h)>s))
  777.     { 
  778. g=0.0;
  779.         for(i=1;i<=m;i++)
  780.         {
  781. aa=(i-1.0)*h; 
  782. bb=i*h;
  783.             w=0.0;
  784.             for(j=0;j<5;j++)
  785.             { 
  786. xx=((bb-aa)*t[j]+(bb+aa))/2.0;
  787.                 q=sqrt(1.0-k*k*sin(xx)*sin(xx));
  788.                 w=w+c[j]/q;
  789.             }
  790.             g=g+w;
  791.         }
  792.         g=g*h/2.0;
  793.         ep=Abs(g-p)/(1.0+Abs(g));
  794.         p=g; 
  795. m=m+m; 
  796. h=0.5*h;
  797.     }
  798.     return(g);
  799. }
  801. //第二类椭圆积分
  802. template <class _Ty>
  803. _Ty Ellipse2ndIntegral(_Ty k, _Ty f)
  804. {
  805. int n(0);
  806.     _Ty pi(3.1415926), e(1.0), ff;
  807.     //_Ty ek();
  808.     if(k<0.0) k=-k;
  809.     if(k>1.0) k=1.0/k;
  810.     _Ty y = Abs(f);
  811.     
  812.     while(y>pi||FloatEqual(y,pi))
  813.     {
  814. n=n+1; 
  815. y=y-pi;
  816. }
  817.     
  818.     if(y>pi/2.0||FloatEqual(y,pi/2.0))
  819.     {
  820. n=n+1;
  821. e=-e; 
  822. y=pi-y;
  823. }
  824.     if(n==0)
  825. ff=ek(k,y);
  826.     else
  827.     { 
  828. ff=ek(k,pi/2.0);
  829.         ff=2.0*n*ff+e*ek(k,y);
  830.     }
  831.     if(f<0.0) ff=-ff;
  832.     return(ff);
  833. }
  835. //第二类椭圆积分辅助函数
  836. template <class _Ty>
  837. _Ty ek(_Ty k, _Ty f)
  838. { int m(1), i, j;
  839.     _Ty p(1.0e+35),ep(FLOATERROR),aa,bb,w,xx,g(0.0),q;
  840.     _Ty t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459};
  841.     _Ty c[5]={0.2369268851,0.4786286705,0.5688888889,
  842. 0.4786286705,0.2369268851};
  843.     
  844.     _Ty h=Abs(f);
  845. _Ty s=Abs(0.0001*h);
  846.     
  847.     while((ep>0.0000001||FloatEqual(ep,0.0000001))&&(Abs(h)>s))
  848.     { 
  849. g=0.0;
  850.         for(i=1;i<=m;i++)
  851. {
  852. aa=(i-1.0)*h;
  853. bb=i*h;
  854.             w=0.0;
  855.             for(j=0;j<5;j++)
  856.             { 
  857. xx=((bb-aa)*t[j]+(bb+aa))/2.0;
  858.                 q=sqrt(1.0-k*k*sin(xx)*sin(xx));
  859.                 w=w+q*c[j];
  860.             }
  861.             g=g+w;
  862.         }
  863.         g=g*h/2.0;
  864.         ep=Abs(g-p)/(1.0+Abs(g));
  865.         p=g; 
  866. m=m+m; 
  867. h=0.5*h;
  868.     }
  869.     return(g);
  870. }
  872. #endif // _SPECIALFUNCTION_INL