Chinaunix首页 | 论坛 | 博客
  • 博客访问: 1676466
  • 博文数量: 585
  • 博客积分: 14610
  • 博客等级: 上将
  • 技术积分: 7402
  • 用 户 组: 普通用户
  • 注册时间: 2008-05-15 10:52
文章存档

2013年(5)

2012年(214)

2011年(56)

2010年(66)

2009年(44)

2008年(200)

分类: C/C++

2012-04-29 05:17:48

 数值计算程序大放送-特殊函数
分类: C/C++ 937人阅读 评论(0) 收藏 举报

选自<<徐世良数值计算程序集(C)>>

所有的函数声明部分

//////////////////////////////////////////////////////////////
//功能:计算伽马(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);

 

#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)   { 
    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);
} //功能:变型第一类整数阶贝塞尔函数
//参数:n-阶数
//调用:
double lfbesl(int n,double x)

int i,m;
    double t,y,p,b0,b1,q;
    static double a[7]={ 1.0,3.5156229,3.0899424,1.2067492,
  0.2659732,0.0360768,0.0045813};
    static double b[7]={ 0.5,0.87890594,0.51498869,
  0.15084934,0.02658773,0.00301532,0.00032411};
    static double c[9]={ 0.39894228,0.01328592,0.00225319,
  -0.00157565,0.00916281,-0.02057706,
  0.02635537,-0.01647633,0.00392377};
    static double d[9]={ 0.39894228,-0.03988024,-0.00362018,
  0.00163801,-0.01031555,0.02282967,
  -0.02895312,0.01787654,-0.00420059};
    if (n<0)
{
  n=-n;
}
    t=fabs(x);
    if (n!=1)

  if (t<3.75)
  { 
   y=(x/3.75)*(x/3.75); 
   p=a[6];
            for (i=5; i>=0; i--)
   {
    p=p*y+a[i];
   }
  }
        else
  { 
   y=3.75/t; 
   p=c[8];
            for (i=7; i>=0; i--)
   {
    p=p*y+c[i];
   }
            p=p*exp(t)/sqrt(t);
  }
}
    if (n==0) 
{
  return(p);
}
    q=p;
    if (t<3.75)

  y=(x/3.75)*(x/3.75); 
  p=b[6];
        for (i=5; i>=0; i--) 
  {
   p=p*y+b[i];
  }
        p=p*t;
}
    else

  y=3.75/t; 
  p=d[8];
        for (i=7; i>=0; i--)
  {
   p=p*y+d[i];
  }
        p=p*exp(t)/sqrt(t);
}
    if (x<0.0) 
{
  p=-p;
  ]
    if (n==1) 
{
  return(p);
}
    if (x==0.0) 
{
  return(0.0);
}
    y=2.0/t; 
t=0.0; 
b1=1.0; 
b0=0.0;
    m=n+(int)sqrt(40.0*n);
    m=2*m;
    for (i=m; i>0; i--)

  p=b0+i*y*b1; 
  b0=b1; b1=p;
        if (fabs(b1)>1.0e+10)
  { 
   t=t*1.0e-10;
   b0=b0*1.0e-10;
            b1=b1*1.0e-10;
  }
        if (i==n) 
  {
   t=b0;
  }
}
    p=t*q/b1;
    if ((x<0.0)&&(n%2==1)) 
{
  p=-p;
}
    return(p);
} //功能:变型第二类整数阶贝塞尔函数
//参数:n-阶数
//调用:lfbesl();
double lgbesl(int n,double x)

int i;
    double y,p,b0,b1;
    static double a[7]={ -0.57721566,0.4227842,0.23069756,
  0.0348859,0.00262698,0.0001075,0.0000074};
    static double b[7]={ 1.0,0.15443144,-0.67278579,
  -0.18156897,-0.01919402,-0.00110404,-0.00004686};
    static double c[7]={ 1.25331414,-0.07832358,0.02189568,
  -0.01062446,0.00587872,-0.0025154,0.00053208};
    static double d[7]={ 1.25331414,0.23498619,-0.0365562,
  0.01504268,-0.00780353,0.00325614,-0.00068245};
    if (n<0) 
{
  n=-n;
}
    if (x<0.0) 
{
  x=-x;
}
    if (x==0.0) 
{
  return(1.0e+70);
}
    if (n!=1)

  if (x<=2.0)
  {
   y=x*x/4.0; 
   p=a[6];
            for (i=5; i>=0; i--) 
   {
    p=p*y+a[i];
   }
            p=p-lfbesl(0,x)*log(x/2.0);
  }
        else
  { 
   y=2.0/x; 
   p=c[6];
            for (i=5; i>=0; i--) 
   {
    p=p*y+c[i];
   }
            p=p*exp(-x)/sqrt(x);
  }
}
    if (n==0) 
{
  return(p);
}
    b0=p;
    if (x<=2.0)

  y=x*x/4.0; 
  p=b[6];
        for (i=5; i>=0; i--)
  {
   p=p*y+b[i];
  }
        p=p/x+lfbesl(1,x)*log(x/2.0);
}
    else

  y=2.0/x; 
  p=d[6];
        for (i=5; i>=0; i--) 
  {
   p=p*y+d[i];
  }
        p=p*exp(-x)/sqrt(x);
}
    if (n==1)
{
  return(p);
}
    b1=p;
    y=2.0/x;
    for (i=1; i{
  p=b0+i*y*b1; 
  b0=b1; 
  b1=p;
}
    return(p);
} //功能:不完全贝塔(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)

double y;
    if (a<=0.0)

  printf("err**a<=0!"); 
  return(-1.0);
}
    if (b<=0.0)

  printf("err**b<=0!"); 
  return(-1.0);
}
    if ((x<0.0)||(x>1.0))

  printf("err**x<0 or x>1 !");
        return(1.0e+70);
}
    if ((x==0.0)||(x==1.0))
{
  y=0.0;
}
    else

  y=a*log(x)+b*log(1.0-x);
        y=exp(y);
        y=y*lagam(a+b)/(lagam(a)*lagam(b));
}
    if (x<(a+1.0)/(a+b+2.0))
{
  y=y*beta(a,b,x)/a;
}
    else
{
  y=1.0-y*beta(b,a,1.0-x)/b;
}
    return(y);
}

static double beta(double a,double b,double x)
{ int k;
    double d,p0,q0,p1,q1,s0,s1;
    p0=0.0; q0=1.0; p1=1.0; q1=1.0;
    for (k=1; k<=100; k++)

  d=(a+k)*(a+b+k)*x;
        d=-d/((a+k+k)*(a+k+k+1.0));
        p0=p1+d*p0; 
  q0=q1+d*q0; 
  s0=p0/q0;
        d=k*(b-k)*x;
        d=d/((a+k+k-1.0)*(a+k+k));
        p1=p0+d*p1; 
  q1=q0+d*q1; 
  s1=p1/q1;
        if (fabs(s1-s0)  {
   return(s1);
  }
}
    printf("a or b too big !");
    return(s1);
}

//功能:正态分布函数
//参数:a-均值,b-方差
//调用:lcerf(),lagam(),lbgam();
double ligas(double a,double d,double x)

double y;
    if (d<=0.0) 
{
  d=1.0e-10;
}
    y=0.5+0.5*lcerf((x-a)/(sqrt(2.0)*d));
    return(y);
} //功能:t-分布函数
//参数:n-自由度
//调用:lhbeta(),lagam();
double ljstd(double t,int n)

double y;
    if (t<0.0)
{
  t=-t;
}
    y=1.0-lhbeta(n/2.0,0.5,n/(n+t*t));
    return(y);
} //功能:X^2-分布函数
//参数:n-自由度
//调用:lbgam(),lagam();
double lkchi(double x,int n)

double y;
    if (x<0.0) 
{
  x=-x;
}
    y=lbgam(n/2.0,x/2.0);
    return(y);
} //功能:F-分布函数
//参数:n1-自由度,n2-自由度
//调用:lhbeta(),lagam();
double llf(double f,int n1,int n2)

double y;
    if (f<0.0) 
{
  f=-f;
}
    y=lhbeta(n2/2.0,n1/2.0,n2/(n2+n1*f));
    return(y);
} //功能:正弦积分
//参数:
//调用:
double lmsi(double x)

int n,k,jt;
    double h,t1,t2,t,s1,s2,p;
    if (x==0.0) 
{
  return(0.0);
}
    h=fabs(x);
    n=1; 
t1=h*(1.0+sin(x)/x)/2.0;
    s1=t1;
jt=1;
    while (jt==1)

  p=0.0;
        for (k=0; k<=n-1; k++)
  {
   t=(k+0.5)*h; 
   p=p+sin(t)/t;
  }
        t2=(t1+h*p)/2.0;
        s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<1.0e-07) 
  {
   jt=0;
  }
        else
  { 
   t1=t2; 
   s1=s2; 
   n=n+n; 
   h=0.5*h;
  }
}
    if (x<0.0) 
{
  s2=-s2;
}
    return(s2);
}

//功能:正弦积分
//参数:
//调用:
double lmsi(double x)

int n,k,jt;
    double h,t1,t2,t,s1,s2,p;
    if (x==0.0) 
{
  return(0.0);
}
    h=fabs(x);
    n=1; 
t1=h*(1.0+sin(x)/x)/2.0;
    s1=t1;
jt=1;
    while (jt==1)

  p=0.0;
        for (k=0; k<=n-1; k++)
  {
   t=(k+0.5)*h; 
   p=p+sin(t)/t;
  }
        t2=(t1+h*p)/2.0;
        s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<1.0e-07) 
  {
   jt=0;
  }
        else
  { 
   t1=t2; 
   s1=s2; 
   n=n+n; 
   h=0.5*h;
  }
}
    if (x<0.0) 
{
  s2=-s2;
}
    return(s2);
} //功能:余弦积分
//参数:
//调用:
double lnci(double x)

int n,k,jt;
    double h,t1,t2,t,s1,s2,p,r,q;
    if (x==0.0) 
{
  x=1.0e-10;
}
    if (x<0.0) 
{
  x=-x;
}
    r=0.57721566490153286060651;
    q=r+log(x);
h=x;
    n=1; 
t1=0.5*h*(1.0-cos(x))/x;
    s1=t1; 
jt=1;
    while (jt==1)

  p=0.0;
        for (k=0; k<=n-1; k++)
  { 
   t=(k+0.5)*h;
   p=p+(1.0-cos(t))/t;
  }
        t2=(t1+h*p)/2.0;
        s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<1.0e-07) 
  {
   jt=0;
  }
        else
  { 
   t1=t2; 
   s1=s2; 
   n=n+n; 
   h=0.5*h;
  }
}
    q=q-s2;
    return(q);
} //功能:指数积分
//参数:
//调用:
double loei(double x)

int n,k,jt;
    double h,t1,t2,t,s1,s2,p,r,q;
    if (x==0.0) 
{
  x=1.0e-10;
}
    if (x<0.0) 
{
  x=-x;
}
    r=0.57721566490153286060651;
    q=r+log(x);
h=x;
    n=1; 
t1=0.5*h*(-1.0+(exp(-x)-1.0)/x);
    s1=t1; 
jt=1;
    while (jt==1)

  p=0.0;
        for (k=0; k<=n-1; k++)
  { 
   t=(k+0.5)*h;
   p=p+(exp(-t)-1.0)/t;
  }
        t2=(t1+h*p)/2.0;
        s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<1.0e-07)
  {
   jt=0;
  }
        else
  { 
   t1=t2;
   s1=s2; 
   n=n+n;
   h=0.5*h;
  }
}
    q=q+s2;
    return(q);
}

//功能:第一类椭圆积分
//参数:
//调用:
static double fk(double k,double y);
double lpfk(double k,double f)
{ int n;
    double pi,y,e,ff;
    if (k<0.0) 
{
  k=-k;
}
    if (k>1.0) 
{
  k=1.0/k;
}
    pi=3.141592653589793;
y=fabs(f);
    n=0;
    while (y>=pi)

  n=n+1;
  y=y-pi;
}
    e=1.0;
    if (y>=pi/2.0)

  n=n+1;
  e=-e; 
  y=pi-y;
}
    if (n==0)
  ff=fk(k,y);
    else

  ff=fk(k,pi/2.0);
        ff=2.0*n*ff+e*fk(k,y);
}
    if (f<0.0) 
{
  ff=-ff;
}
    return(ff);
}

static double fk(double k,double y)

int n,i,jt;
    double h,t,t1,t2,s1,s2,p,q;
    if ((1.0-k<1.0e-20)&&(fabs(sin(y)-1.0)<1.0e-20))
{
  return(1.0e+10);
}
    n=1;
h=y;
    q=sqrt(1.0-k*k*sin(y)*sin(y));
    if (q<=1.0e-35)
{
  q=1.0e+35;
}
    else 
{
  q=1.0/q;
}
    t1=0.5*h*(1.0+q);
    s1=t1; 
jt=1;
    while (jt==1)

  p=0.0;
        for (i=0; i<=n-1; i++)
  { 
   t=(i+0.5)*h;
            q=sqrt(1.0-k*k*sin(t)*sin(t));
            if (q<=1.0e-35) 
   {
    q=1.0e+35;
   }
            else 
   {
    q=1.0/q;
   }
            p=p+q;
  }
        t2=(t1+h*p)/2.0;
  s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)  {
   jt=0;
  }
        else if (h<1.0e-02) 
  {
   jt=0;
  }
        else
  { 
   t1=t2; 
   s1=s2; 
   n=n+n; 
   h=0.5*h;
  }
}
    return(s2);
} //功能:第二类椭圆积分
//参数:
//调用:
static double ek(double k,double y);
double lqek(double k,double f)

int n;
    double pi,y,e,ff;
    if (k<0.0) 
  k=-k;
    if (k>1.0) 
  k=1.0/k;
    pi=3.1415926; 
y=fabs(f);
    n=0;
    while (y>=pi)

  n=n+1; 
  y=y-pi;
}
    e=1.0;
    if (y>=pi/2.0)

  n=n+1;
  e=-e; 
  y=pi-y;
}
    if (n==0)
  ff=ek(k,y);
    else

  ff=ek(k,pi/2.0);
        ff=2.0*n*ff+e*ek(k,y);
}
    if (f<0.0) 
  ff=-ff;
    return(ff);
}

static double ek(double k,double y)

int n,i,jt;
    double h,t,t1,t2,s1,s2,p,q;
    n=1;
h=y;
    q=sqrt(1.0-k*k*sin(y)*sin(y));
    t1=0.5*h*(1.0+q);
    s1=t1;
jt=1;
    while (jt==1)

  p=0.0;
        for (i=0; i<=n-1; i++)
  { 
   t=(i+0.5)*h;
            q=sqrt(1.0-k*k*sin(t)*sin(t));
            p=p+q;
  }
        t2=(t1+h*p)/2.0;
  s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<1.0e-07) 
   jt=0;
        else if 
   (h<1.0e-03) jt=0;
        else
  { 
   t1=t2; 
   s1=s2;
   n=n+n; 
   h=0.5*h;
  }
}
    return(s2);
}

阅读(748) | 评论(0) | 转发(0) |
给主人留下些什么吧!~~