//////////////////////////////////////////////////////////////
//功能：计算伽马(gamma)函数值,gamma函数积分区间为0到正无穷
//描述：gamma[x]=Integrate[exp[-t]*t^(x-1),{t,0,∞}]
//调用：
double lagam(double x);
//////////////////////////////////////////////////////////////
//功能：不完全伽马函数
//描述：gamma[a,x]=P[a,x]/gamma[x]
//描述：P[a,x]=Integrate[exp[-t]*t^(a-1),{t,0,x}]
//参数：a-参数
//调用：lagam(x)函数
double lbgam(double a,double x);
//////////////////////////////////////////////////////////////
//功能：误差函数
//描述：erf[x]=gamma[0.5,x^2]
//描述：erf[x]=2/sqrt[pi]*Integrate[exp[-t^2],{t,0,x}]
//调用：lagam(),lbgam()
double lcerf(double x);
//////////////////////////////////////////////////////////////
//功能：第一类整数阶贝塞尔函数
//参数：n-阶数
//调用：
double ldbesl(int n,double x);
//////////////////////////////////////////////////////////////
//功能：第二类整数阶贝塞尔函数
//参数：n-阶数
//调用：ldbesl()
double lebesl(int n,double x);
//////////////////////////////////////////////////////////////
//功能：变型第一类整数阶贝塞尔函数
//参数：n-阶数
//调用：
double lfbesl(int n,double x);
//////////////////////////////////////////////////////////////
//功能：变型第二类整数阶贝塞尔函数
//参数：n-阶数
//调用：lfbesl();
double lgbesl(int n,double x);
//////////////////////////////////////////////////////////////
//功能：不完全贝塔(beta)函数
//描述：Bx[a,b]=Integrate[t^(a-1)*(1-t)^(b-1),{t,0,x}]/B[a,b]
//描述：B[a,b]=gamma[a]*gamma[b]/gamma[a+b]
//参数：a-参数,b-参数
//调用：lagam();
double lhbeta(double a,double b,double x);
//////////////////////////////////////////////////////////////
//功能：正态分布函数
//参数：a-均值,b-方差
//调用：lcerf(),lagam(),lbgam();
double ligas(double a,double d,double x);
//////////////////////////////////////////////////////////////
//功能：t-分布函数
//参数：n-自由度
//调用：lhbeta(),lagam();
double ljstd(double t,int n);
//////////////////////////////////////////////////////////////
//功能：X^2-分布函数
//参数：n-自由度
//调用：lbgam(),lagam();
double lkchi(double x,int n);
//////////////////////////////////////////////////////////////
//功能：F-分布函数
//参数：n1-自由度,n2-自由度
//调用：lhbeta(),lagam();
double llf(double f,int n1,int n2);
//////////////////////////////////////////////////////////////
//功能：正弦积分
//参数：
//调用：
double lmsi(double x);
//////////////////////////////////////////////////////////////
//功能：余弦积分
//参数：
//调用：
double lnci(double x);
//////////////////////////////////////////////////////////////
//功能：指数积分
//参数：
//调用：
double loei(double x);
//////////////////////////////////////////////////////////////
//功能：第一类椭圆积分
//参数：
//调用：
double lpfk(double k,double f);
//////////////////////////////////////////////////////////////
//功能：第二类椭圆积分
//参数：
//调用：
double lqek(double k,double f);

后面还有，一次发不完，要一点点发

[此贴子已经被作者于2007-2-3 23:15:08编辑过]

#include "stdio.h"
#include "math.h"
//功能：计算伽马(gamma)函数值,gamma函数积分区间为0到正无穷
//描述：Integrate[exp[-t]*t^(x-1),{t,0,∞}]
//调用：
double lagam(double x)
{
int i;
double y,t,s,u;
static double a[11]={ 0.0000677106,-0.0003442342,
0.0015397681,-0.0024467480,0.0109736958,
-0.0002109075,0.0742379071,0.0815782188,
0.4118402518,0.4227843370,1.0};
if (x<=0.0)
{
printf("err**x<=0!\n");
return(-1.0);
}
y=x;
if (y<=1.0)
{
t=1.0/(y*(y+1.0));
y=y+2.0;
}
else if (y<=2.0)
{
t=1.0/y;
y=y+1.0;
}
else if (y<=3.0)
{
t=1.0;
}
else
{
t=1.0;
while (y>3.0)
{
y=y-1.0;
t=t*y;
}
}
s=a[0];
u=y-2.0;
for (i=1; i<=10; i++)
{
s=s*u+a[i];
}
s=s*t;
return(s);
} //功能：不完全伽马函数
//描述：gamma[a,x]=P[a,x]/gamma[x]
//描述：P[a,x]=Integrate[exp[-t]*t^(a-1),{t,0,x}]
//参数：a-参数
//调用：lagam(x)函数
double lbgam(double a,double x)
double a,x;
{
int n;
double p,q,d,s,s1,p0,q0,p1,q1,qq;
if ((a<=0.0)||(x<0.0))
{
if (a<=0.0)
{
printf("err**a<=0!\n");
}
if (x<0.0)
{
printf("err**x<0!\n");
}
return(-1.0);
}
if (x+1.0==1.0)
{
return(0.0);
}
if (x>1.0e+35)
{
return(1.0);
}
q=log(x);
q=a*q;
qq=exp(q);
if (x<1.0+a)
{
p=a;
d=1.0/a;
s=d;
for (n=1; n<=100; n++)
{
p=1.0+p;
d=d*x/p;
s=s+d;
if (fabs(d)<fabs(s)*1.0e-07)
{
s=s*exp(-x)*qq/lagam(a);
return(s);
}
}
}
else
{
s=1.0/x;
p0=0.0;
p1=1.0;
q0=1.0;
q1=x;
for (n=1; n<=100; n++)
{
p0=p1+(n-a)*p0;
q0=q1+(n-a)*q0;
p=x*p0+n*p1;
q=x*q0+n*q1;
if (fabs(q)+1.0!=1.0)
{
s1=p/q;
p1=p;
q1=q;
if (fabs((s1-s)/s1)<1.0e-07)
{
s=s1*exp(-x)*qq/lagam(a);
return(1.0-s);
}
s=s1;
}
p1=p;
q1=q;
}
}
printf("a too large !\n");
s=1.0-s*exp(-x)*qq/lagam(a);
return(s);
} //功能：误差函数
//描述：erf[x]=gamma[0.5,x^2]
//描述：erf[x]=2/sqrt[pi]*Integrate[exp[-t^2],{t,0,x}]
//调用：lagam(),lbgam()
double lcerf(double x)
{
double y;
if (x>=0.0)
{
y=lbgam(0.5,x*x);
}
else
{
y=-lbgam(0.5,x*x);
}
return(y);
} //功能：第一类整数阶贝塞尔函数
//参数：n-阶数
//调用：
double ldbesl(int n,double x)
{
int i,m;
double t,y,z,p,q,s,b0,b1;
static double a[6]={ 57568490574.0,-13362590354.0,
651619640.7,-11214424.18,77392.33017,-184.9052456};
static double b[6]={ 57568490411.0,1029532985.0,
9494680.718,59272.64853,267.8532712,1.0};
static double c[6]={ 72362614232.0,-7895059235.0,
242396853.1,-2972611.439,15704.4826,-30.16036606};
static double d[6]={ 144725228443.0,2300535178.0,
18583304.74,99447.43394,376.9991397,1.0};
static double e[5]={ 1.0,-0.1098628627e-02,0.2734510407e-04,
-0.2073370639e-05,0.2093887211e-06};
static double f[5]={ -0.1562499995e-01,0.1430488765e-03,
-0.6911147651e-05,0.7621095161e-06,-0.934935152e-07};
static double g[5]={ 1.0,0.183105e-02,-0.3516396496e-04,
0.2457520174e-05,-0.240337019e-06};
static double h[5]={ 0.4687499995e-01,-0.2002690873e-03,
0.8449199096e-05,-0.88228987e-06,0.105787412e-06};
t=fabs(x);
if (n<0)
{
n=-n;
}
if (n!=1)
{
if (t<8.0)
{
y=t*t;
p=a[5];
q=b[5];
for (i=4; i>=0; i--)
{
p=p*y+a[i];
q=q*y+b[i];
}
p=p/q;
}
else
{
z=8.0/t;
y=z*z;
p=e[4];
q=f[4];
for (i=3; i>=0; i--)
{
p=p*y+e[i];
q=q*y+f[i];
}
s=t-0.785398164;
p=p*cos(s)-z*q*sin(s);
p=p*sqrt(0.636619772/t);
}
}
if (n==0)
{
return(p);
}
b0=p;
if (t<8.0)
{
y=t*t;
p=c[5];
q=d[5];
for (i=4; i>=0; i--)
{
p=p*y+c[i];
q=q*y+d[i];
}
p=x*p/q;
}
else
{
z=8.0/t;
y=z*z;
p=g[4];
q=h[4];
for (i=3; i>=0; i--)
{
p=p*y+g[i];
q=q*y+h[i];
}
s=t-2.356194491;
p=p*cos(s)-z*q*sin(s);
p=p*x*sqrt(0.636619772/t)/t;
}
if (n==1)
{
return(p);
}
b1=p;
if (x==0.0)
{
return(0.0);
}
s=2.0/t;
if (t>1.0*n)
{
if (x<0.0)
{
b1=-b1;
}
for (i=1; i<=n-1; i++)
{
p=s*i*b1-b0;
b0=b1;
b1=p;
}
}
else
{
m=(n+(int)sqrt(40.0*n))/2;
m=2*m;
p=0.0;
q=0.0;
b0=1.0;
b1=0.0;
for (i=m-1; i>=0; i--)
{
t=s*(i+1)*b0-b1;
b1=b0;
b0=t;
if (fabs(b0)>1.0e+10)
{
b0=b0*1.0e-10;
b1=b1*1.0e-10;
p=p*1.0e-10;
q=q*1.0e-10;
}
if ((i+2)%2==0)
{
q=q+b0;
}
if ((i+1)==n)
{
p=b1;
}
}
q=2.0*q-b0;
p=p/q;
}
if ((x<0.0)&&(n%2==1))
{
p=-p;
}
return(p);
}
//功能：第二类整数阶贝塞尔函数
//参数：n-阶数
//调用：ldbesl()
double lebesl(int n,double x)
{
int i;
double y,z,p,q,s,b0,b1;
extern double ldbesl();
static double a[6]={ -2.957821389e+9,7.062834065e+9,
-5.123598036e+8,1.087988129e+7,-8.632792757e+4,
2.284622733e+2};
static double b[6]={ 4.0076544269e+10,7.452499648e+8,
7.189466438e+6,4.74472647e+4,2.261030244e+2,1.0};
static double c[6]={ -4.900604943e+12,1.27527439e+12,
-5.153438139e+10,7.349264551e+8,-4.237922726e+6,
8.511937935e+3};
static double d[7]={ 2.49958057e+13,4.244419664e+11,
3.733650367e+9,2.245904002e+7,1.02042605e+5,
3.549632885e+2,1.0};
static double e[5]={ 1.0,-0.1098628627e-02,
0.2734510407e-04,-0.2073370639e-05,
0.2093887211e-06};
static double f[5]={ -0.1562499995e-01,
0.1430488765e-03,-0.6911147651e-05,
0.7621095161e-06,-0.934935152e-07};
static double g[5]={ 1.0,0.183105e-02,
-0.3516396496e-04,0.2457520174e-05,
-0.240337019e-06};
static double h[5]={ 0.4687499995e-01,
-0.2002690873e-03,0.8449199096e-05,
-0.88228987e-06,0.105787412e-06};
if (n<0)
{
n=-n;
}
if (x<0.0)
{
x=-x;
}
if (x==0.0)
{
return(-1.0e+70);
}
if (n!=1)
{
if (x<8.0)
{
y=x*x;
p=a[5];
q=b[5];
for (i=4; i>=0; i--)
{
p=p*y+a[i];
q=q*y+b[i];
}
p=p/q+0.636619772*ldbesl(0,x)*log(x);
}
else
{
z=8.0/x;
y=z*z;
p=e[4];
q=f[4];
for (i=3; i>=0; i--)
{
p=p*y+e[i];
q=q*y+f[i];
}
s=x-0.785398164;
p=p*sin(s)+z*q*cos(s);
p=p*sqrt(0.636619772/x);
}
}
if (n==0)
{
return(p);
}
b0=p;
if (x<8.0)
{
y=x*x;
p=c[5];
q=d[6];
for (i=4; i>=0; i--)
{
p=p*y+c[i];
q=q*y+d[i+1];
}
q=q*y+d[0];
p=x*p/q+0.636619772*(ldbesl(1,x)*log(x)-1.0/x);;
}
else
{
z=8.0/x;
y=z*z;
p=g[4];
q=h[4];
for (i=3; i>=0; i--)
{
p=p*y+g[i];
q=q*y+h[i];
}
s=x-2.356194491;
p=p*sin(s)+z*q*cos(s);
p=p*sqrt(0.636619772/x);
}
if (n==1)
{
return(p);
}
b1=p;
s=2.0/x;
for (i=1; i<=n-1; i++)
{
p=s*i*b1-b0;
b0=b1;
b1=p;
}
return(p);
}

刚才出了点问题，没有发完，现在接着发


#include "stdio.h"
#include "math.h"
//功能：计算伽马(gamma)函数值,gamma函数积分区间为0到正无穷
//描述：Integrate[exp[-t]*t^(x-1),{t,0,∞}]
//调用：
double lagam(double x)
{
int i;
double y,t,s,u;
static double a[11]={ 0.0000677106,-0.0003442342,
0.0015397681,-0.0024467480,0.0109736958,
-0.0002109075,0.0742379071,0.0815782188,
0.4118402518,0.4227843370,1.0};
if (x<=0.0)
{
printf("err**x<=0!\n");
return(-1.0);
}
y=x;
if (y<=1.0)
{
t=1.0/(y*(y+1.0));
y=y+2.0;
}
else if (y<=2.0)
{
t=1.0/y;
y=y+1.0;
}
else if (y<=3.0)
{
t=1.0;
}
else
{
t=1.0;
while (y>3.0)
{
y=y-1.0;
t=t*y;
}
}
s=a[0];
u=y-2.0;
for (i=1; i<=10; i++)
{
s=s*u+a[i];
}
s=s*t;
return(s);
} //功能：不完全伽马函数
//描述：gamma[a,x]=P[a,x]/gamma[x]
//描述：P[a,x]=Integrate[exp[-t]*t^(a-1),{t,0,x}]
//参数：a-参数
//调用：lagam(x)函数
double lbgam(double a,double x)
double a,x;
{
int n;
double p,q,d,s,s1,p0,q0,p1,q1,qq;
if ((a<=0.0)||(x<0.0))
{
if (a<=0.0)
{
printf("err**a<=0!\n");
}
if (x<0.0)
{
printf("err**x<0!\n");
}
return(-1.0);
}
if (x+1.0==1.0)
{
return(0.0);
}
if (x>1.0e+35)
{
return(1.0);
}
q=log(x);
q=a*q;
qq=exp(q);
if (x<1.0+a)
{
p=a;
d=1.0/a;
s=d;
for (n=1; n<=100; n++)
{
p=1.0+p;
d=d*x/p;
s=s+d;
if (fabs(d)<fabs(s)*1.0e-07)
{
s=s*exp(-x)*qq/lagam(a);
return(s);
}
}
}
else
{
s=1.0/x;
p0=0.0;
p1=1.0;
q0=1.0;
q1=x;
for (n=1; n<=100; n++)
{
p0=p1+(n-a)*p0;
q0=q1+(n-a)*q0;
p=x*p0+n*p1;
q=x*q0+n*q1;
if (fabs(q)+1.0!=1.0)
{
s1=p/q;
p1=p;
q1=q;
if (fabs((s1-s)/s1)<1.0e-07)
{
s=s1*exp(-x)*qq/lagam(a);
return(1.0-s);
}
s=s1;
}
p1=p;
q1=q;
}
}
printf("a too large !\n");
s=1.0-s*exp(-x)*qq/lagam(a);
return(s);
} //功能：误差函数
//描述：erf[x]=gamma[0.5,x^2]
//描述：erf[x]=2/sqrt[pi]*Integrate[exp[-t^2],{t,0,x}]
//调用：lagam(),lbgam()
double lcerf(double x)
{
double y;
if (x>=0.0)
{
y=lbgam(0.5,x*x);
}
else
{
y=-lbgam(0.5,x*x);
}
return(y);
} //功能：第一类整数阶贝塞尔函数
//参数：n-阶数
//调用：
double ldbesl(int n,double x)
{
int i,m;
double t,y,z,p,q,s,b0,b1;
static double a[6]={ 57568490574.0,-13362590354.0,
651619640.7,-11214424.18,77392.33017,-184.9052456};
static double b[6]={ 57568490411.0,1029532985.0,
9494680.718,59272.64853,267.8532712,1.0};
static double c[6]={ 72362614232.0,-7895059235.0,
242396853.1,-2972611.439,15704.4826,-30.16036606};
static double d[6]={ 144725228443.0,2300535178.0,
18583304.74,99447.43394,376.9991397,1.0};
static double e[5]={ 1.0,-0.1098628627e-02,0.2734510407e-04,
-0.2073370639e-05,0.2093887211e-06};
static double f[5]={ -0.1562499995e-01,0.1430488765e-03,
-0.6911147651e-05,0.7621095161e-06,-0.934935152e-07};
static double g[5]={ 1.0,0.183105e-02,-0.3516396496e-04,
0.2457520174e-05,-0.240337019e-06};
static double h[5]={ 0.4687499995e-01,-0.2002690873e-03,
0.8449199096e-05,-0.88228987e-06,0.105787412e-06};
t=fabs(x);
if (n<0)
{
n=-n;
}
if (n!=1)
{
if (t<8.0)
{
y=t*t;
p=a[5];
q=b[5];
for (i=4; i>=0; i--)
{
p=p*y+a[i];
q=q*y+b[i];
}
p=p/q;
}
else
{
z=8.0/t;
y=z*z;
p=e[4];
q=f[4];
for (i=3; i>=0; i--)
{
p=p*y+e[i];
q=q*y+f[i];
}
s=t-0.785398164;
p=p*cos(s)-z*q*sin(s);
p=p*sqrt(0.636619772/t);
}
}
if (n==0)
{
return(p);
}
b0=p;
if (t<8.0)
{
y=t*t;
p=c[5];
q=d[5];
for (i=4; i>=0; i--)
{
p=p*y+c[i];
q=q*y+d[i];
}
p=x*p/q;
}
else
{
z=8.0/t;
y=z*z;
p=g[4];
q=h[4];
for (i=3; i>=0; i--)
{
p=p*y+g[i];
q=q*y+h[i];
}
s=t-2.356194491;
p=p*cos(s)-z*q*sin(s);
p=p*x*sqrt(0.636619772/t)/t;
}
if (n==1)
{
return(p);
}
b1=p;
if (x==0.0)
{
return(0.0);
}
s=2.0/t;
if (t>1.0*n)
{
if (x<0.0)
{
b1=-b1;
}
for (i=1; i<=n-1; i++)
{
p=s*i*b1-b0;
b0=b1;
b1=p;
}
}
else
{
m=(n+(int)sqrt(40.0*n))/2;
m=2*m;
p=0.0;
q=0.0;
b0=1.0;
b1=0.0;
for (i=m-1; i>=0; i--)
{
t=s*(i+1)*b0-b1;
b1=b0;
b0=t;
if (fabs(b0)>1.0e+10)
{
b0=b0*1.0e-10;
b1=b1*1.0e-10;
p=p*1.0e-10;
q=q*1.0e-10;
}
if ((i+2)%2==0)
{
q=q+b0;
}
if ((i+1)==n)
{
p=b1;
}
}
q=2.0*q-b0;
p=p/q;
}
if ((x<0.0)&&(n%2==1))
{
p=-p;
}
return(p);
}
//功能：第二类整数阶贝塞尔函数
//参数：n-阶数
//调用：ldbesl()
double lebesl(int n,double x)
{
int i;
double y,z,p,q,s,b0,b1;
extern double ldbesl();
static double a[6]={ -2.957821389e+9,7.062834065e+9,
-5.123598036e+8,1.087988129e+7,-8.632792757e+4,
2.284622733e+2};
static double b[6]={ 4.0076544269e+10,7.452499648e+8,
7.189466438e+6,4.74472647e+4,2.261030244e+2,1.0};
static double c[6]={ -4.900604943e+12,1.27527439e+12,
-5.153438139e+10,7.349264551e+8,-4.237922726e+6,
8.511937935e+3};
static double d[7]={ 2.49958057e+13,4.244419664e+11,
3.733650367e+9,2.245904002e+7,1.02042605e+5,
3.549632885e+2,1.0};
static double e[5]={ 1.0,-0.1098628627e-02,
0.2734510407e-04,-0.2073370639e-05,
0.2093887211e-06};
static double f[5]={ -0.1562499995e-01,
0.1430488765e-03,-0.6911147651e-05,
0.7621095161e-06,-0.934935152e-07};
static double g[5]={ 1.0,0.183105e-02,
-0.3516396496e-04,0.2457520174e-05,
-0.240337019e-06};
static double h[5]={ 0.4687499995e-01,
-0.2002690873e-03,0.8449199096e-05,
-0.88228987e-06,0.105787412e-06};
if (n<0)
{
n=-n;
}
if (x<0.0)
{
x=-x;
}
if (x==0.0)
{
return(-1.0e+70);
}
if (n!=1)
{
if (x<8.0)
{
y=x*x;
p=a[5];
q=b[5];
for (i=4; i>=0; i--)
{
p=p*y+a[i];
q=q*y+b[i];
}
p=p/q+0.636619772*ldbesl(0,x)*log(x);
}
else
{
z=8.0/x;
y=z*z;
p=e[4];
q=f[4];
for (i=3; i>=0; i--)
{
p=p*y+e[i];
q=q*y+f[i];
}
s=x-0.785398164;
p=p*sin(s)+z*q*cos(s);
p=p*sqrt(0.636619772/x);
}
}
if (n==0)
{
return(p);
}
b0=p;
if (x<8.0)
{
y=x*x;
p=c[5];
q=d[6];
for (i=4; i>=0; i--)
{
p=p*y+c[i];
q=q*y+d[i+1];
}
q=q*y+d[0];
p=x*p/q+0.636619772*(ldbesl(1,x)*log(x)-1.0/x);;
}
else
{
z=8.0/x;
y=z*z;
p=g[4];
q=h[4];
for (i=3; i>=0; i--)
{
p=p*y+g[i];
q=q*y+h[i];
}
s=x-2.356194491;
p=p*sin(s)+z*q*cos(s);
p=p*sqrt(0.636619772/x);
}
if (n==1)
{
return(p);
}
b1=p;
s=2.0/x;
for (i=1; i<=n-1; i++)
{
p=s*i*b1-b0;
b0=b1;
b1=p;
}
return(p);
}

[此贴子已经被作者于2007-2-3 23:14:02编辑过]