ch12.h

将C语言的常用程序集移植到VC开发环境的源代码

      }
    p=t*q/b1;
    if ((x<0.0)&&(n%2==1)) p=-p;
    return(p);
}
/////////////////////////////////////////////////////////////
double lbsl4(int n,double x)
{
	int i;
    double y,p,b0,b1;
    //extern double lbsl3();
    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-lbsl3(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+lbsl3(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<n; i++)
      { p=b0+i*y*b1; b0=b1; b1=p;}
    return(p);
}
/////////////////////////////////////////////////////////////
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)<fabs(s1)*1.0e-07)
          return(s1);
      }
    printf("a or b too big !");
    return(s1);
}

double lbeta(double a,double b,double x)
{ 
	double y;
    //double beta();
    //extern double lgam1();
    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*lgam1(a+b)/(lgam1(a)*lgam1(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);
}
/////////////////////////////////////////////////////////////
double lgass(double a,double d,double x)
{ 
	double y;
    //extern double lerrf();
    if (d<=0.0) d=1.0e-10;
    y=0.5+0.5*lerrf((x-a)/(sqrt(2.0)*d));
    return(y);
}
/////////////////////////////////////////////////////////////
double lstdt(double t,int n)
{ 
	double y;
    //extern double lbeta();
    if (t<0.0) t=-t;
    y=1.0-lbeta(n/2.0,0.5,n/(n+t*t));
    return(y);
}
/////////////////////////////////////////////////////////////
double lchii(double x,int n)
{ 
	double y;
    //extern double lgam2();
    if (x<0.0) x=-x;
    y=lgam2(n/2.0,x/2.0);
    return(y);
}
/////////////////////////////////////////////////////////////
double lffff(double f,int n1,int n2)
{
	double y;
    //extern double lbeta();
    if (f<0.0) f=-f;
    y=lbeta(n2/2.0,n1/2.0,n2/(n2+n1*f));
    return(y);
}
/////////////////////////////////////////////////////////////
double lsinn(double x)
{ 
	int m,i,j;
    double s,p,ep,h,aa,bb,w,xx,g;
    static double t[5]={-0.9061798459,-0.5384693101,0.0,
                         0.5384693101,0.9061798459};
    static double c[5]={0.2369268851,0.4786286705,0.5688888889,
                        0.4786286705,0.2369268851};
    m=1;
    if (x==0) return(0.0);
    h=fabs(x);  s=fabs(0.0001*h);
    p=1.0e+35; ep=0.000001; g=0.0;
    while ((ep>=0.0000001)&&(fabs(h)>s))
      { g=0.0;
        for (i=1;i<=m;i++)
          { aa=(i-1.0)*h; bb=i*h;
            w=0.0;
            for (j=0;j<=4;j++)
              { xx=((bb-aa)*t[j]+(bb+aa))/2.0;
                w=w+sin(xx)/xx*c[j];
              }
            g=g+w;
          }
        g=g*h/2.0;
        ep=fabs(g-p)/(1.0+fabs(g));
        p=g; m=m+1; h=fabs(x)/m;
      }
    return(g);
}
/////////////////////////////////////////////////////////////
double lcoss(double x)
{ 
	int m,i,j;
    double s,p,ep,h,aa,bb,w,xx,g,r,q;
    static double t[5]={-0.9061798459,-0.5384693101,0.0,
                         0.5384693101,0.9061798459};
    static double c[5]={0.2369268851,0.4786286705,0.5688888889,
                        0.4786286705,0.2369268851};
    m=1;
    if (x==0) x=1.0e-35;
    if (x<0.0) x=-x;
    r=0.57721566490153286060651;
    q=r+log(x);
    h=x; s=fabs(0.0001*h);
    p=1.0e+35; ep=0.000001; g=0.0;
    while ((ep>=0.0000001)&&(fabs(h)>s))
      { g=0.0;
        for (i=1;i<=m;i++)
          { aa=(i-1.0)*h; bb=i*h;
            w=0.0;
            for (j=0;j<=4;j++)
              { xx=((bb-aa)*t[j]+(bb+aa))/2.0;
                w=w+(1.0-cos(xx))/xx*c[j];
              }
            g=g+w;
          }
        g=g*h/2.0;
        ep=fabs(g-p)/(1.0+fabs(g));
        p=g; m=m+1; h=x/m;
      }
    g=q-g;
    return(g);
}
/////////////////////////////////////////////////////////////
double lexpp(double x)
{ 
	int m,i,j;
    double s,p,ep,h,aa,bb,w,xx,g,r,q;
    static double t[5]={-0.9061798459,-0.5384693101,0.0,
                         0.5384693101,0.9061798459};
    static double c[5]={0.2369268851,0.4786286705,0.5688888889,
                        0.4786286705,0.2369268851};
    m=1;
    if (x==0) x=1.0e-10;
    if (x<0.0) x=-x;
    r=0.57721566490153286060651;
    q=r+log(x);
    h=x; s=fabs(0.0001*h);
    p=1.0e+35; ep=0.000001; g=0.0;
    while ((ep>=0.0000001)&&(fabs(h)>s))
      { g=0.0;
        for (i=1;i<=m;i++)
          { aa=(i-1.0)*h; bb=i*h;
            w=0.0;
            for (j=0;j<=4;j++)
              { xx=((bb-aa)*t[j]+(bb+aa))/2.0;
                w=w+(exp(-xx)-1.0)/xx*c[j];
              }
            g=g+w;
          }
        g=g*h/2.0;
        ep=fabs(g-p)/(1.0+fabs(g));
        p=g; m=m+1; h=x/m;
      }
    g=q+g;
    return(g);
}
/////////////////////////////////////////////////////////////
static double fk(double k,double f)
{
	int m,i,j;
    double s,p,ep,h,aa,bb,w,xx,g,q;
    static double t[5]={-0.9061798459,-0.5384693101,0.0,
                         0.5384693101,0.9061798459};
    static double c[5]={0.2369268851,0.4786286705,0.5688888889,
                        0.4786286705,0.2369268851};
    m=1; g=0.0;
    h=fabs(f); s=fabs(0.0001*h);
    p=1.0e+35; ep=0.000001;
    while ((ep>=0.0000001)&&(fabs(h)>s))
      { g=0.0;
        for (i=1;i<=m;i++)
          { aa=(i-1.0)*h; bb=i*h;
            w=0.0;
            for (j=0;j<=4;j++)
              { xx=((bb-aa)*t[j]+(bb+aa))/2.0;
                q=sqrt(1.0-k*k*sin(xx)*sin(xx));
                w=w+c[j]/q;
              }
            g=g+w;
          }
        g=g*h/2.0;
        ep=fabs(g-p)/(1.0+fabs(g));
        p=g; m=m+m; h=0.5*h;
      }
    return(g);
}
double lelp1(double k,double f)
{
	int n;
    double pi,y,e,ff;
    //double fk();
    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=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 ek(double k,double f)
{ 
	int m,i,j;
    double s,p,ep,h,aa,bb,w,xx,g,q;
    static double t[5]={-0.9061798459,-0.5384693101,0.0,
                         0.5384693101,0.9061798459};
    static double c[5]={0.2369268851,0.4786286705,0.5688888889,
                        0.4786286705,0.2369268851};
    m=1; g=0.0;
    h=fabs(f); s=fabs(0.0001*h);
    p=1.0e+35; ep=0.000001;
    while ((ep>=0.0000001)&&(fabs(h)>s))
      { g=0.0;
        for (i=1;i<=m;i++)
          { aa=(i-1.0)*h; bb=i*h;
            w=0.0;
            for (j=0;j<=4;j++)
              { xx=((bb-aa)*t[j]+(bb+aa))/2.0;
                q=sqrt(1.0-k*k*sin(xx)*sin(xx));
                w=w+q*c[j];
              }
            g=g+w;
          }
        g=g*h/2.0;
        ep=fabs(g-p)/(1.0+fabs(g));
        p=g; m=m+m; h=0.5*h;
      }
    return(g);
}
double lelp2(double k,double f)
{
	int n;
    double pi,y,e,ff;
    //double ek();
    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);
  }
/////////////////////////////////////////////////////////////
#endif