gauss.cpp

bpsk 调制

#include "gauss.h"
#include "stdio.h"
#include <math.h>
 
  double lgam1(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);
  }


  double lgam2(double a,double x) 
  { int n;
    double p,q,d,s,s1,p0,q0,p1,q1,qq;
    //extern double lgam1();
    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/lgam1(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/lgam1(a);
                    return(1.0-s);
                  }
                s=s1;
              }
            p1=p; q1=q;
          }
      }
    printf("a too large !\n");
    s=1.0-s*exp(-x)*qq/lgam1(a);
    return(s);
  }


  double lerrf(double x)
  { double y;
  //  extern double lgam2();
    if (x>=0.0)
      y=lgam2(0.5,x*x);
    else
      y=-lgam2(0.5,x*x);
    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 mgrn1(double u,double g,double *r)
 { int i,m;
    double s,w,v,t;
    s=65536.0; w=2053.0; v=13849.0;
    t=0.0;
    for (i=1; i<=12; i++)
      { *r=(*r)*w+v; m=(int)(*r/s);
        *r=*r-m*s; t=t+(*r)/s;
      }
    t=u+g*(t-6.0);
    return(t);
  }