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

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 查看源码
SPECFUNC.CPP
资源名称：BCMATH.rar　[点击查看]
上传用户：zengbais
上传日期：2022-08-08
资源大小：49k
文件大小：6k
源码类别：

数值算法/人工智能

开发平台：

C++ Builder

SPECFUNC.CPP：源码内容
  1. #include "func.h"
  3. DOUBLE gamma(DOUBLE x)
  4. {
  5. size_t i;
  6. DOUBLE y,t,s,u;
  7. static DOUBLE a[11]={ 0.0000677106,-0.0003442342,
  8.   0.0015397681,-0.0024467480,0.0109736958,
  9.   -0.0002109075,0.0742379071,0.0815782188,
  10.   0.4118402518,0.4227843370,1.0};
  11. if (x<=0.0)
  12. throw "err**x<=0!n";
  13. y=x;
  14. if (y<=1.0)
  15. { t=1.0/(y*(y+1.0)); y+=2.0;}
  16. else if (y<=2.0)
  17. { t=1.0/y; y+=1.0;}
  18. else if (y<=3.0) t=1.0;
  19. else
  20. { t=1.0;
  21.   while (y>3.0)
  22.  { y-=1.0; t=t*y;}
  23. }
  24.  s=a[0]; u=y-2.0;
  25.  for (i=1; i<=10; i++)
  26. s=s*u+a[i];
  27.  s*=t;
  28.  return(s);
  29. }
  31. DOUBLE gamma2(DOUBLE a, DOUBLE x)
  32. { size_t n;
  33. DOUBLE p,q,d,s,s1,p0,q0,p1,q1,qq;
  34. if ((a<=0.0)||(x<0.0))
  35. { if (a<=0.0) throw "err**a<=0!n";
  36.   if (x<0.0) throw "err**x<0!n";
  37.   return(-1.0);
  38. }
  39. if (x+1.0==1.0) return(0.0);
  40. if (x>1.0e+35) return(1.0);
  41. q=log(x); q=a*q; qq=exp(q);
  42. if (x<1.0+a)
  43. { p=a; d=1.0/a; s=d;
  44.   for (n=1; n<=100; n++)
  45.  { p=1.0+p; d=d*x/p; s=s+d;
  46.  if (fabs(d)<fabs(s)*defaulterr)
  47.   { s=s*exp(-x)*qq/gamma(a);
  48.  return(s);
  49.   }
  50.  }
  51. }
  52.  else
  53. { s=1.0/x; p0=0.0; p1=1.0; q0=1.0; q1=x;
  54.   for (n=1; n<=100; n++)
  55.  { p0=p1+(n-a)*p0; q0=q1+(n-a)*q0;
  56. p=x*p0+n*p1; q=x*q0+n*q1;
  57. if (fabs(q)+1.0!=1.0)
  58.   { s1=p/q; p1=p; q1=q;
  59.  if (fabs((s1-s)/s1)<defaulterr)
  60. { s=s1*exp(-x)*qq/gamma(a);
  61.   return(1.0-s);
  62. }
  63.  s=s1;
  64.   }
  65. p1=p; q1=q;
  66.  }
  67. }
  68. //  printf("a too large !n");
  69.  s=1.0-s*exp(-x)*qq/gamma(a);
  70.  return(s);
  71. }
  73. DOUBLE erf(DOUBLE x)
  74. { double y;
  75. if (x>=0.0)
  76. y=gamma2(0.5,x*x);
  77. else
  78. y=-gamma2(0.5,x*x);
  79. return(y);
  80. }
  83. static DOUBLE bet(DOUBLE a,DOUBLE b,DOUBLE x);
  85. DOUBLE beta(DOUBLE a,DOUBLE b,DOUBLE x) // 计算不完全贝塔分布函数
  86.   { DOUBLE y;
  87.  if (a<=0.0) throw "err**a<=0!";
  88.  if (b<=0.0) throw "err**b<=0!";
  89.  if ((x<0.0)||(x>1.0)) throw "err**x<0 or x>1 !";
  90.  if ((x==0.0)||(x==1.0)) y=0.0;
  91.  else
  92. { y=a*log(x)+b*log(1.0-x);
  93.   y=exp(y);
  94.   y*=gamma(a+b)/(gamma(a)*gamma(b));
  95. }
  96.  if (x<(a+1.0)/(a+b+2.0))
  97. y*=bet(a,b,x)/a;
  98.  else
  99. y=1.0-y*bet(b,a,1.0-x)/b;
  100.  return(y);
  101.   }
  103. static DOUBLE bet(DOUBLE a,DOUBLE b,DOUBLE x)
  104. { size_t k;
  105.  DOUBLE d,p0,q0,p1,q1,s0,s1;
  106.  p0=0.0; q0=1.0; p1=1.0; q1=1.0;
  107.  for (k=1; k<=100; k++)
  108. { d=(a+k)*(a+b+k)*x;
  109.   d=-d/((a+k+k)*(a+k+k+1.0));
  110.   p0=p1+d*p0; q0=q1+d*q0; s0=p0/q0;
  111.   d=k*(b-k)*x;
  112.   d=d/((a+k+k-1.0)*(a+k+k));
  113.   p1=p0+d*p1; q1=q0+d*q1; s1=p1/q1;
  114.   if (fabs(s1-s0)<fabs(s1)*defaulterr)
  115.  return(s1);
  116. }
  117. //  printf("a or b too big !");
  118.  return(s1);
  119. }
  121. DOUBLE gass(DOUBLE a,DOUBLE d,DOUBLE x)
  122. { DOUBLE y;
  123. if (d<=0.0) d=1.0e-10;
  124. y=0.5+0.5*erf((x-a)/(sqrt(2.0)*d));
  125. return(y);
  126. }
  128. DOUBLE student(DOUBLE t, size_t n)
  129.   { DOUBLE y;
  130.  if (t<0.0) t=-t;
  131.  y=1.0-beta(n/2.0,0.5,n/(n+t*t));
  132.  return(y);
  133.   }
  135. DOUBLE chii(DOUBLE x,size_t n)
  136. {   DOUBLE y;
  137.  if (x<0.0) x=-x;
  138.  y=gamma2(n/2.0,x/2.0);
  139.  return(y);
  140. }
  142. DOUBLE fdisp(DOUBLE f,size_t n1,size_t n2)
  143. { DOUBLE y;
  144. if (f<0.0) f=-f;
  145. y=beta(n2/2.0,n1/2.0,n2/(n2+n1*f));
  146. return(y);
  147. }
  149. DOUBLE integral(DOUBLE (*f)(DOUBLE),DOUBLE a, DOUBLE b, DOUBLE eps)
  150.  // 函数f,在(a,b)区间积分,采用勒让德高斯求积法
  151. {
  152.  size_t m,i,j;
  153.  DOUBLE s,p,ep,h,aa,bb,w,x,g;
  154.  static DOUBLE t[5]={-0.9061798459,-0.5384693101,0.0,
  155.  0.5384693101,0.9061798459};
  156.  static DOUBLE c[5]={0.2369268851,0.4786286705,0.5688888889,
  157. 0.4786286705,0.2369268851};
  158.  m=1;
  159.  h=b-a; s=fabs(0.001*h);
  160.  p=1.0e+35; ep=eps+1.0;
  161.  while ((ep>=eps)&&(fabs(h)>s))
  162. { g=0.0;
  163.   for (i=1;i<=m;i++)
  164.  { aa=a+(i-1.0)*h; bb=a+i*h;
  165. w=0.0;
  166. for (j=0;j<=4;j++)
  167.   { x=((bb-aa)*t[j]+(bb+aa))/2.0;
  168.  w=w+f(x)*c[j];  //flrgsf(x)*c[j];
  169.   }
  170. g=g+w;
  171.  }
  172.   g *= h/2.0;
  173.   ep=fabs(g-p)/(1.0+fabs(g));
  174.   p=g; m=m+1; h=(b-a)/m;
  175. }
  176.  return g;
  177. }
  179. static DOUBLE sinc(DOUBLE x)
  180. {
  181. return sin(x)/x;
  182. }
  184. DOUBLE sinn(DOUBLE x)
  185. {
  186. if(x==0.0)return 0;
  187. if(x<0.0) throw "x less than 0n";
  188. return integral(sinc,0,x);
  189. }
  191. static DOUBLE cosc(DOUBLE x)
  192. {
  193. return (1-cos(x))/x;
  194. }
  196. #define EULER 0.57721566490153286060651
  198. DOUBLE coss(DOUBLE x)
  199. {
  200. DOUBLE r=EULER;
  201. if (x<0) throw "x less than 0n";
  202. if (x==0) x=1.0e-35;
  203. return r+log(x)-integral(cosc,0,x);
  204. }
  206. static DOUBLE ex(DOUBLE x)
  207. {
  208. return (exp(-x)-1)/x;
  209. }
  211. DOUBLE expp(DOUBLE x)
  212. {
  213. if(x<0) throw "x less than 0n";
  214. if(x==0) x=1.0e-10;
  215. DOUBLE r=EULER;
  216. return r+log(x)+integral(ex,0,x);
  217. }
  219. DOUBLE getroot(DOUBLE (*f)(DOUBLE), DOUBLE x0, DOUBLE eps)
  220. { int i,j,m,l;// it
  221. double a[10],y[10],z,h,q;
  222. l=10; q=1.0e+35; h=0.0;
  223. while (l!=0)
  224. { l=l-1; j=0; // it=l;
  225.   while (j<=7)
  226.   { if (j<=2) z=x0+0.1*j;
  227.  else z=h;
  228.  y[j]=f(z);
  229.  h=z;
  230.  if (j==0) a[0]=z;
  231.  else
  232. { m=0; i=0;
  233.   while ((m==0)&&(i<=j-1))
  234.  { if (fabs(h-a[i])==0.0) m=1;
  235. else h=(y[j]-y[i])/(h-a[i]);
  236. i=i+1;
  237.  }
  238.   a[j]=h;
  239.   if (m!=0) a[j]=q;
  240.   h=0.0;
  241.   for (i=j-1; i>=0; i--)
  242.  { if (fabs(a[i+1]+h)==0.0) h=q;
  243. else h=-y[i]/(a[i+1]+h);
  244.  }
  245.   h=h+a[0];
  246. }
  247.  if (fabs(y[j])>=eps) j=j+1;
  248.  else { j=10; l=0;}
  249.   }
  250. x0=h;
  251.  }
  252. if(fabs(f(x0)) > eps) throw "no root!";
  253. return x0;
  254. }
  256. DOUBLE getroot(algo & alg, DOUBLE x0, DOUBLE eps)
  257. { int i,j,m,l;
  258. double a[10],y[10],z,h,q;
  259. l=10; q=1.0e+35; h=0.0;
  260. while (l!=0)
  261. { l--; j=0;
  262.   while (j<=7)
  263.   { if (j<=2) z=x0+0.1*j;
  264.  else z=h;
  265.  y[j]=alg.cal(z);
  266.  h=z;
  267.  if (j==0) a[0]=z;
  268.  else
  269. { m=0; i=0;
  270.   while ((m==0)&&(i<=j-1))
  271.  { if (fabs(h-a[i])==0.0) m=1;
  272. else h=(y[j]-y[i])/(h-a[i]);
  273. i=i+1;
  274.  }
  275.   a[j]=h;
  276.   if (m!=0) a[j]=q;
  277.   h=0.0;
  278.   for (i=j-1; i>=0; i--)
  279.  { if (fabs(a[i+1]+h)==0.0) h=q;
  280. else h=-y[i]/(a[i+1]+h);
  281.  }
  282.   h=h+a[0];
  283. }
  284.  if (fabs(y[j])>=eps) j=j+1;
  285.  else { j=10; l=0;}
  286.   }
  287. x0=h;
  288.  }
  289. if(fabs(alg.cal(x0)) > eps) throw "no root!";
  290. return x0;
  291. }