CSDN博客

img dcyu

我的算法头文件

发表于2002/9/29 10:51:00  1343人阅读

/*
  FileName : MyLib.h
 
  Author    : dcyu
 
  Date       : 2002.9.9
*/

#include <stdio.h>
#include <conio.h>
#include <time.h>
#include <bios.h>
#include <stdlib.h>
#include <malloc.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <graphics.h>

#define infinity INT_MAX
#define True  1
#define False 0

typedef double Graph ;
typedef int _bool ;
typedef struct complex Complex;


void   Error(char *message);
void   _Pause(void);
void   _delay(double delaytime);
double _time(void);
double _dif(double start,double end);
double **_Alloc2(int r,int c);
int    **_Alloc2_int(int r,int c);
void   _Alloc2Free(double **x);
void   _dtoa(double value, char *string);
void   _putstring(int x, int y, char *msg, int color);
long   _fac(int n);
int    _fac2(int n,int *a);
double _random(void);
int    _rand(int seed);
double _avg(double a,double b);
double _sta(double mu,double sigma);
double _sta2(double mu,double sigma);
double _kai(int n);
double _tdis(int n);
double _F(int n1,int n2);
double _Poisson(int k,double lam);
double _mean(double *a,int start,int end);
double _variance(double *a,int start,int end);
double _Linear(double *x, double *y, double *a, double *b, int n);
double _Floyd(int i,int j,int k,double **d);
double _Dijkstra(Graph **G,int n,int s,int t, int *path, int *count);
double _0618(double (*f)(double x),double start,double end,double eps);
double _integral(double (*f)(double x),double start,double end,int n);
double _integral2(double (*f)(double x),double start,double end,int div);
long   _comb(int n,int m);
long   _rank(int n,int m);
void   _allcomb_dg(int* ps, int* pe, int elems, int *buf, int bufsz,
      int **comb, int* iCount);
void   _allcomb(int *list, int n, int elems, int **comb);
void   _allrank_dg(int m, int k, int s, int *a, int *flag, int *iCount, int **rank);
void   _allrank(int m, int **rank);
int    _bolziman(double de,double t);
int    _isprime(long p);
long   _prime(int n);
void   _swapi(int *a,int *b);
void   _swapl(long *a,long *b);
void   _swapd(double *a,double *b);
void   _inv(int *x,int n);
void   _sort(double *x,int n);
long   _max2(long n,long m);
long   _min2(long n,long m);
double _max(double *array,int n,int *id);
double _min(double *array,int n,int *id);
long   _GCD(long m,long n);
long   _LCM(long m,long n);
void   _initgraph(void);
void   _setPlotdefault(void);
void   _Plot(double (*f)(double x),double start,double end);
double _Line_equation(double x, double a, double b);
void   _Fit(double *x, double *y, int n);
Complex _Cplus(Complex z1,Complex z2);
Complex _Cminus(Complex z1,Complex z2);
Complex _Cmul(Complex z1,Complex z2);
Complex _Cdiv(Complex z1,Complex z2);
double  _Cabs(Complex z);
void    _Croot(Complex z1, int n, Complex *z);
Complex _Cpow(Complex z1, double w);
double _LAG(double *x, double *y, double t, int n);
double _NEWT(double *x, double *y, int n, double t);
int    _Mgauss(double **a, double *b, int n, double ep);
int    _Mgauss2(double **a, double **b, int n, int m, double ep);
int    _Mdjn(double **a, double **c, int n, int m);
double _NOR(double x,int l);
double _Mdet(double **a, int n);
int    _Minv(double t0, double *t, double *tt, int n, int m, double **b);


void Error(char *message)
{
 clrscr();
 fprintf(stderr,"Error: %s/n",message);
 exit(1);

}

void _Pause(void)
{
 while(bioskey(1)==0)  ;

}

void _delay(double delaytime)
{
   long bios_time;
   double start,end;

   bios_time = biostime(0, 0L);
   start = bios_time / CLK_TCK;

   while(1)
   {
    bios_time = biostime(0, 0L);
    end = bios_time / CLK_TCK;
    if(end-start >= delaytime)  return ;

   }

}

double _time(void)
{
 long bios_time;
 double cur;

 bios_time = biostime(0, 0L);
 cur = bios_time / CLK_TCK;
 return cur;

}

double _dif(double start,double end)
{
 return end-start;

}

double **_Alloc2(int r,int c)
{
 double *x,**y;
 int n;

 x=(double *)calloc(r*c,sizeof(double));
 y=(double **)calloc(r,sizeof(double *));
 for(n=0;n<=r-1;n++)
   y[n]=&x[c*n];

 return y;

}

int **_Alloc2_int(int r,int c)
{
 int *x,**y;
 int n;

 x=(int *)calloc(r*c,sizeof(int));
 y=(int **)calloc(r,sizeof(int *));
 for(n=0;n<=r-1;n++)
   y[n]=&x[c*n];

 return y;

}

void _Alloc2Free(double **x)
{
 free(x[0]);
 free(x);

}

void _dtoa(double value, char *string)
{
 unsigned long wn,df;
 char *str1,*str2,*str3,*str4;

 str1=value<0.0?"-":"";
 value=value<0.0?(-value):value;
 wn=(unsigned long)(floor(value));
 df=(unsigned long)(floor((value-wn)*1000));
 str2=(char *)malloc(10*sizeof(char));
 str4=(char *)malloc(10*sizeof(char));
 ultoa(wn,str2,10);
 str3=".";
 ultoa(df,str4,10);
 strcpy(string, str1);
 strcat(string, str2);
 strcat(string, str3);
 strcat(string, str4);

}

void _putstring(int x, int y, char *msg, int color)
{
 int len;
 const int sidex=8, sidey=10;

 setcolor(color);
 len=strlen(msg);
 x=x-len*sidex;
 y=y-sidey;
 outtextxy(x,y,msg);

}

long _fac(int n)
{
 if(n<0||n>12) Error("n is not valid in _fac!");
 if(n==0||n==1) return 1;
 return _fac(n-1)*n;

}

int _fac2(int n,int *a)
{
   int m,i,j,c,t;

   a[0]=1;
   m=1;
 for(i=2;i<=n;i++)
 {
   for(c=0,j=0;j<m;j++)
      {
          t=a[j]*i+c;
          a[j]=t%10;
          c=t/10;
       }
      while(c)
      {
         a[m++]=c%10;
         c/=10;
      }
 }

 for(j=0;j<=(m-1)/2;j++) { t=a[j]; a[j]=a[m-1-j]; a[m-1-j]=t; }
 return m;

}

double _random(void)
{
 int a;
 double r;

 a=rand()%32767;
 r=(a+0.00)/32767.00;

 return r;

}

int _rand(int seed)
{
 return rand()%seed;

}

double _avg(double a,double b)
{
 double _random(void);

 return a+_random()*(b-a);

}

double _sta(double mu,double sigma)
{
 double _random(void);

 int i;
 double r,sum=0.0;

 if(sigma<=0.0) Error("Sigma<=0.0 in _sta!");
 for(i=1;i<=12;i++)
  sum = sum + _random();
 r=(sum-6.00)*sigma+mu;

 return r;

}

double _sta2(double mu,double sigma)
{
 double r1,r2;

 r1=_random();
 r2=_random();

 return sqrt(-2*log(r1))*cos(2*M_PI*r2)*sigma+mu ;

}

double _kai(int n)
{
 double _sta(double mu,double sigma);

 int i;
 double sum=0.0;
 for(i=1;i<=n;i++)
     sum=sum + pow(_sta(0,1),2);
 return sum;

}

double _tdis(int n)
{
 double _kai(int n);
 double _sta(double mu,double sigma);

 double x,y;

 x=_sta(0,1);
 y=_kai(n);

 return x/sqrt(y/n);

}

double _F(int n1,int n2)
{
 double _kai(int n);

 double u1,u2;
 u1=_kai(n1);
 u2=_kai(n2);

 return (u1*n2)/(u2*n1);

}

double _Poisson(int k,double lam)
{
 if(k>12||k<0) Error("k is not valid in _Poisson!");

 return pow(lam,k)*exp(-lam)/_fac(k);

}

double _mean(double *a,int start,int end)
{
 int i;
 double sum=0.0;

 if(start<0||end-start<0)
   Error("Start is larger than end , or start is not valid in _mean!");

 for(i=start;i<=end;i++)
   sum=sum+a[i];

 return sum/(end-start+1);

}

double _variance(double *a,int start,int end)
{
 double _mean(double *a,int start,int end);

 int i;
 double mean,sum=0.0;

 mean=_mean(a,start,end);
 for(i=start;i<=end;i++)
     sum=sum+(a[i]-mean)*(a[i]-mean);

 return sum/(end-start);

}

double _Linear(double *x, double *y, double *a, double *b, int n)
{
 double Sxx,Syy,Sxy,meanx,meany,Qe;
 int i;

 if(n<=2) Error("n too small in _Linear!");
 meanx=_mean(x,0,n-1);
 meany=_mean(y,0,n-1);
 Sxx=_variance(x,0,n-1)*(n-1);
 Syy=_variance(y,0,n-1)*(n-1);
 Sxy=0.0;
 for(i=0;i<n;i++)
   Sxy=Sxy+(x[i]-meanx)*(y[i]-meany);
 if(Sxx<=1e-7) Error("Curve too craggedness in _Linear!");
 *b=Sxy/Sxx;
 *a=meany-(*b)*meanx;
 Qe=Syy-(*b)*Sxy;
 return Qe/(n-2);

}

double **_Adj_Form(double (*a)[3], int n, int m)
{
 double **G;
 int i,j;

 G=_Alloc2(n,n);
 for(i=0;i<n;i++)
  for(j=0;j<n;j++)
 G[i][j]=infinity;

 for(i=0;i<m;i++)
 {
  G[(int)a[i][0]][(int)a[i][1]]=a[i][2];
  G[(int)a[i][1]][(int)a[i][0]]=a[i][2];
 }
 return G;

}

double _Floyd(int i,int j,int k,double **d)
{
 double temp_ij,temp_ik,temp_kj;
 if(k>0) {
  temp_ij=_Floyd(i,j,k-1,d);
  temp_ik=_Floyd(i,k,k-1,d);
  temp_kj=_Floyd(k,j,k-1,d);
  return temp_ij<(temp_ik+temp_kj)?temp_ij:(temp_ik+temp_kj);
 }
 if(k==0) return *(*(d+i)+j);
 else { Error("k<0 in _Floyd!"); return 1; }

}

double _Dijkstra(Graph **G,int n,int s,int t, int *path, int *count)
{
 int i, j, *mark, *pathd, *pathbk;
 double w,minc,*d;
 int from, to;

 d=(double *)malloc(n*sizeof(double));
 mark=(int *)malloc(n*sizeof(int));
 pathd=(int *)malloc(n*sizeof(int));
 pathbk=(int *)malloc(n*sizeof(int));
 for (i=0; i<n; i++) mark[i]=0;
 for (i=0; i<n; i++) pathd[i]=0;
 for (i=0; i<n; i++)
 {
 d[i]=G[s][i];
 pathd[i]=s;
 }
 mark[s]=1; pathd[s]=0; d[s]=0;
 for(i=1; i<n; i++)
 {
 minc = infinity;
 w = 0;
 for( j = 0; j < n; j++ )
   if( ( mark[j]==0 ) && ( minc >= d[j] ) )
   {  minc=d[j]; w=j; }

 mark[w]=1;
 for(j=0; j<n; j++)
   if( (mark[j]==0) && ( G[w][j] != infinity ) && ( d[j] > d[w]+G[w][j] ) )
   {
    d[j]=d[w]+G[w][j];
    pathd[j]=w;
   }
 }
 from=s; to=t;
  for(i=0;i<n;i++)
 path[i]=pathbk[i]=-1;
 j=0;
 if(d[t]!=infinity)
 {
  while(from!=to)
  {
 pathbk[j]=pathd[to];
 to=pathd[to];
 j++;

  }

 }

  j=0;
  for(i=n-1;i>=0;i--)
 if(pathbk[i]!=-1)
 {
  path[j]=pathbk[i];
  j++;
 }
  path[j]=t;
  *count=j+1;

  return d[t];

}

double _0618(double (*f)(double x),double start,double end,double eps)
{
 double a,b,c,d;

 if(eps<0.0||end-start<eps)
  Error("end-start<eps , or eps<0.0 in _0618!");
 a=start; b=end;
 c=a+(b-a)*0.381966011;
 d=a+(b-a)*0.618033989;
 while(b-a>eps)
 {
  if((*f)(a)>(*f)(c) && (*f)(c)>(*f)(d))  {
   a=c;    c=d;    d=a+b-c;    }

  else if((*f)(b)>(*f)(d) && (*f)(d)>(*f)(c))  {
   b=d;    d=c;    c=a+b-d;    }

  else {
   a=c;    b=d;    c=a+(b-a)*0.381966011;
   d=a+(b-a)*0.618033989;      }

 }
 return (a+b)/2;

}


int _2div(double (*f)(double x), double a, double b, double h, double eps,
    double *x, int n, int *m)
{
 double z,z0,z1,y,y0,y1;
 *m=0;
 z=a;
 y=(*f)(z);
 while(1)
 {
  if((z>b+h/2)||(*m==n)) return 1;
  if(fabs(y)<eps)
  {
   *m+=1;
   x[*m-1]=z;
   z+=h/2;
   y=(*f)(z);
   continue;
  }
  z1=z+h;
  y1=(*f)(z1);
  if(fabs(y1)<eps)
  {
   *m+=1;
   x[*m-1]=z1;
   z=z1+h/2;
   y=(*f)(z);
   continue;
  }
  if(y*y1>0)
  {
   y=y1;
   z=z1;
   continue;
  }
  while(1)
  {
   if(fabs(z1-z)<eps)
   {
    *m+=1;
    x[*m-1]=(z1+z)/2;
    y=(*f)(z);
    break;
   }
   z0=(z1+z)/2;
   y0=(*f)(z0);
   if(fabs(y0)<eps)
   {
    *m+=1;
    x[*m-1]=z0;
    z=z0+h/2;
    y=(*f)(z);
    break;
   }
   if(y*y0<0)
   {
    z1=z0;
    y1=y0;
   }
   else
   {
    z=z0;
    y=y0;
   }
  }
 }
 return 1;

}


double _integral(double (*f)(double x),double start,double end,int n)
{
 double h,al,be,sum;
 int i,m;

 if(n<=0||end-start<=0.0)
  Error("Start is larger than end , or n<=0 in _integral!");
 h=(end-start)/n;
 al=start+h*0.4226497308;
 be=start+h*1.577350269;
 sum=(*f)(al)+(*f)(be);
 m=n/2-1;
 for(i=1;i<=m;i++)
 {
  al=al+2*h;
  be=be+2*h;
  sum=sum+(*f)(al)+(*f)(be);

 }
 sum=sum*sum;
 return sum;

}

double _integral2(double (*f)(double x),double start,double end,int div)
{
 double s,h,y;
 int i;

 if(div<=0||end-start<=0.0)
  Error("Start is larger than end , or div<=0 in _integral2!");
 s=((*f)(start)+(*f)(end))/2.00;
 h=(end-start)/div;
 for(i=1;i<div;i++)
     s=s+(*f)(start+i*h);
 y=s*h;
 return y;

}

long _comb(int n,int m)
{
 long _fac(int n);

 int i;
 long temp=1;

 if(n>15&&m>8||n<=0||m<0||n<m) Error("m,n is not valid in _comb!");

 for(i=n;i>=n-m+1;i--)
  temp=temp*i;
 temp=temp/_fac(m);
 return temp;

}

long _rank(int n,int m)
{
 int i;
 long temp=1;

 if(n>12||n<=0||m<0||n<m)
  Error("m,n is not valid in _rank!");

 for(i=n;i>=n-m+1;i--)
     temp=temp*i;
 return temp;

}

void _allcomb_dg(int* ps, int* pe, int elems, int buf[], int bufsz, int **comb, int* iCount)
{
  int* pi, i;
  if(elems == 1) {
 for(pi = ps; pi < pe; pi++)
 {
   for(i = bufsz - 1; i > 0; i--)
   comb[*iCount][i]=buf[i];
   comb[*iCount][0]=*pi;
   *iCount=*iCount+1;
 }
 return ;
  }
  for(pi = ps; pi < pe - elems; pi++) {
 buf[elems - 1] = *pi;
 _allcomb_dg(pi + 1, pe, elems - 1, buf, bufsz, comb, iCount);
  }
  if(pi == pe - elems) {
 for(i = bufsz - 1; i > elems - 1; i--)
  comb[*iCount][i]=buf[i];
 for(;pi < pe; pi++,i--)
  comb[*iCount][i]=*pi;
 *iCount=*iCount+1;
  }

}

void _allcomb(int *list, int n, int elems, int **comb)
{
  int *buf;
  int i,iCount;

  iCount=0;
  buf=(int *)malloc(elems*sizeof(int));
  _allcomb_dg(&list[0], &list[n], elems, buf, elems, comb, &iCount);
  for(i=0;i<iCount;i++)
  _inv(comb[i], elems);

}

void _allrank_dg(int m, int k, int s, int *a, int *flag, int *iCount, int **rank)
{
  int i;

  if(s>=k)
  {
    for (i = 0; i < k; i++)
      rank[*iCount][i] = a[i] ;

 *iCount=*iCount+1;
  }
  else
  {
 for (i = 1; i <= m; i++)
   if (flag[i]==0)
      {
        flag[i]=1;
    a[s]=i;
    _allrank_dg(m,k,s+1,a,flag,iCount,rank);
    a[s]=0;
        flag[i]=0;
      }
  }

}

void _allrank(int m, int **rank)
{
 int *flag, *a ;
 int i,iCount,n,k;

 n=m; k=0;
 flag=(int *)malloc((m+1)*sizeof(int));
 a=(int *)malloc((m+1)*sizeof(int));
 for(i=1;i<=m;i++) flag[i]=0;
 iCount=0;
 _allrank_dg(m,n,k,a,flag,&iCount,rank);

}

int _bolziman(double de,double t)
{
  double _random(void);

  return de<0.0 || _random()<exp(-de/(t));

}

int _isprime(long p)
{
 long i;

 for(i=2;i<=sqrt(p);i++)
  if(p%i==0) return 0;

 return 1;

}

long _prime(int n)
{
 int _isprime(long p);

 long i=2;
 int num=1;

 while(num<=n) {
  if(_isprime(i++))  num++;
 }
 return i-1;

}

void _swapi(int *a,int *b)
{
 int p;
 p=*a;
 *a=*b;
 *b=p;

}

void _swapl(long *a,long *b)
{
 long p;
 p=*a;
 *a=*b;
 *b=p;

}

void _swapd(double *a,double *b)
{
 double p;
 p=*a;
 *a=*b;
 *b=p;

}

void _inv(int *x,int n)
{
 int *p,t,*i,*j,m;

 m=(n-1)/2;
 i=x; j=x+n-1; p=x+m;
 for(;i<=p;i++,j--)
 { t=*i; *i=*j; *j=t; }

}

void _sort(double *x,int n)
{
 int i,j,k;
 double t;

 for(i=0;i<n-1;i++)
 {
  k=i;
  for(j=i+1;j<n;j++)
  if(*(x+j)>*(x+k))  k=j;
  if(k!=i)
  { t=*(x+i); *(x+i)=*(x+k); *(x+k)=t; }

 }

}

long _max2(long n,long m)
{
 return n>m?n:m;

}

long _min2(long n,long m)
{
 return n<m?n:m;

}

double _max(double *array,int n,int *id)
{
 double *p, *array_end,dmax;

 array_end=array+n;
 dmax=*array;  *id=0;
 for(p=array+1;p<array_end;p++)
  if(*p>dmax)  {
 dmax=*p;  *id=p-array;  }
 return dmax;

}

double _min(double *array,int n,int *id)
{
 double *p, *array_end, dmin;

 array_end=array+n;
 dmin=*array;   *id=0;
 for(p=array+1;p<array_end;p++)
  if(*p<dmin)  {
 dmin=*p;  *id=p-array;  }

 return dmin;

}

long _GCD(long m,long n)
{
 long _max2(long n,long m);
 long _min2(long n,long m);

 long t;

 t=_max2(n,m);
 m=_min2(n,m);
 n=t;
 while(n-m!=0)
 {
  n=n-m;
  t=_max2(n,m);
  m=_min2(n,m);
  n=t;

 }
 return t;

}

long _LCM(long m,long n)
{
 long _GCD(long n,long m);

 return m*n/_GCD(m,n);

}

void _initgraph(void)
{
 int gdriver=VGA, gmode=VGAHI, errorcode;

 initgraph(&gdriver,&gmode,".//");
 errorcode = graphresult();
 if (errorcode != grOk)
 {
 printf("/nGraphics error: %s/n", grapherrormsg(errorcode));
 printf("Press any key to halt!");
 getch();
 exit(1);
 }

}

void _setPlotdefault(void)
{
 const int basex=60, basey=420;
 const int basebkcolor=LIGHTBLUE, basecolor=RED;

 setbkcolor(basebkcolor);
 setcolor(basecolor);

 line(0,basey,getmaxx(),basey);
 line(basex,0,basex,getmaxy());

 setfillstyle(SOLID_FILL,basecolor);
 fillellipse(basex,basey,3,3);

}

void _Plot(double (*f)(double x),double start,double end)
{
 const int basex=60, basey=420, grid=20;
 const double eps=1e-5;
 double stepx, stepy, *point;
 double dmax, dmin, di;
 int i,x,y,lx,ly,maxid,minid;

 if(end-start<eps) { closegraph(); Error("end-start too small in_Plot!"); }

 _setPlotdefault();
 stepx=(end-start)/(getmaxx()-basex+0.0);
 point=(double *)malloc((getmaxx()-basex)*sizeof(double *));
 for(i=basex;i<=getmaxx();i++)
   point[i-basex]=(*f)(start+(i-basex)*stepx);
 dmax=_max(point, getmaxx()-basex,&maxid);
 dmin=_min(point, getmaxx()-basex,&minid);
 if(dmax-dmin<eps) { closegraph(); Error("Curve too smoothless in _Plot!"); }
 stepy=(dmax-dmin)/(basey-10.0);
 lx=basex; ly=basey-(point[0]-dmin)/stepy;
 for(i=basex+1;i<=getmaxx();i++)
 {
  x=i; y=basey-(int)((point[i-basex]-dmin)/stepy);
  line(lx,ly,x,y);
  lx=x; ly=y;

 }
 setcolor(BLUE);
 outtextxy(getmaxx()/2+50,getmaxy()-30,"Press any key to start analysis!");
 if(getch()==27) return ;
 closegraph();
 clrscr();
 printf("Curve analysis:/n");
 printf("From %lf to %lf/n",start,end);
 printf("The maximum of the curve is %lf, when x = %lf./n",
   dmax, start+(maxid+0.0)/(639-basex+0.0)*(end-start) );
 printf("The minimum of the curve is %lf, when x = %lf./n",
   dmin, start+(minid+0.0)/(639-basex+0.0)*(end-start) );
 stepx=(end-start)/(grid+0.0);
 for(di=start;di<=end;di=di+stepx)
   printf("f ( %lf ) = %lf/n", di, (*f)(di));

 _Pause();

}

double _Line_equation(double x, double a, double b)
{
 return a+b*x;

}

void _Fit(double *x, double *y, int n)
{
 double _Linear(double *x, double *y, double *a, double *b, int n);

 const int basex=60, basey=420, grid=20;
 const double eps=1e-5;
 double a,b;
 double dmaxx,dminx,dmaxy,dminy,*pdx,*pdy;
 int i,maxidx,minidx,maxidy,minidy,*px,*py;
 double stepx, stepy, *point;
 double dmax, dmin, di, start, end;
 int cx,cy,lx,ly,maxid,minid;
 char *str;

 _setPlotdefault();
 _Linear(x, y, &a, &b, n);

 dmaxx=_max(x,n,&maxidx);
 dminx=_min(x,n,&minidx);
 dmaxy=_max(y,n,&maxidy);
 dminy=_min(y,n,&minidy);
 if(dmaxx-dminx<eps||dmaxy-dminy<eps)
   Error("Curve too smoothless in _Fit!");

 start=dminx; end=dmaxx;
 stepx=(end-start)/(getmaxx()-basex+0.0);
 point=(double *)malloc((getmaxx()-basex)*sizeof(double *));
 for(i=basex;i<=getmaxx();i++)
   point[i-basex]=a+b*(start+(i-basex)*stepx);
 dmax=_max(point, getmaxx()-basex,&maxid);
 dmin=_min(point, getmaxx()-basex,&minid);
 if(dmax-dmin<eps)
 { closegraph(); Error("Curve too smoothless in _Plot!"); }
 stepy=(dmax-dmin)/(basey-10.0);
 lx=basex; ly=basey-(point[0]-dmin)/stepy;
 for(i=basex+1;i<=getmaxx();i++)
 {
  cx=i; cy=basey-(int)((point[i-basex]-dmin)/stepy);
  line(lx,ly,cx,cy);
  lx=cx; ly=cy;

 }
 str=(char *)malloc(10*sizeof(char));
 for(i=0;i<grid;i++)
 {
 _dtoa(point[(getmaxx()-basex)*i/grid],str);
 _putstring(basex,basey-(basey-10.0)*i/grid,str,BLUE);
 line(basex-4,basey-(basey-10.0)*i/grid,basex+4,basey-(basey-10.0)*i/grid);
 if(i%3==0) { _dtoa((start+(getmaxx()-basex)*i/grid*stepx),str);
 _putstring(basex+(getmaxx()-basex)*i/grid+strlen(str)*4,basey+12,str,BLUE);
 line(basex+(getmaxx()-basex)*i/grid,basey-4,basex+(getmaxx()-basex)*i/grid,basey+4);
 }

 }
 setcolor(RED);
 px=(int *)malloc(n*sizeof(int *));
 py=(int *)malloc(n*sizeof(int *));
 pdx=x; pdy=y;
 for(i=0;i<n;i++)
 {
   *px=basex+(*pdx-dminx)/(dmaxx-dminx)*(getmaxx()-basex-0.0);
   *py=basey-(int)((*pdy-dmin)/stepy);
   setfillstyle(SOLID_FILL,RED);
   fillellipse(*px,*py,3,3);
   px++;  py++;
   pdx++; pdy++;

 }

}

Complex _Cplus(Complex z1,Complex z2)
{
 Complex z;
 z.x=z1.x+z2.x;
 z.y=z1.y+z2.y;
 return z;

}

Complex _Cminus(Complex z1,Complex z2)
{
 Complex z;
 z.x=z1.x-z2.x;
 z.y=z1.y-z2.y;
 return z;

}

Complex _Cmul(Complex z1,Complex z2)
{
 Complex z;
 z.x=z1.x*z2.x-z1.y*z2.y;
 z.y=z1.x*z2.y+z1.y*z2.x;
 return z;

}

Complex _Cdiv(Complex z1,Complex z2)
{
 double e,f;
 Complex z;

 if(fabs(z2.x)>=fabs(z2.y))
 {
  e=z2.y/z2.x;
  f=z2.x+e*z2.y;
  z.x=(z1.x+z1.y*e)/f;
  z.y=(z1.y-z1.x*e)/f;
 }
 else
 {
  e=z2.x/z2.y;
  f=z2.y+e*z2.x;
  z.x=(z1.x*e+z1.y)/f;
  z.y=(z1.y*e-z1.x)/f;
 }
 return z;

}

double _Cabs(Complex z)
{
 return hypot(z.x, z.y);

}

void _Croot(Complex z1, int n, Complex *z)
{
 int j;
 double r, c, c1, ck, cs, s, s1, sk, sc;

 if(z1.y==0.0)
 {
  r=exp(log(fabs(z1.x))/n);
  cs=0;
 }
 else if(z1.x==0.0)
 {
  if(z1.y>0)  cs=M_PI_2;
  else cs=-M_PI_2;
  r=exp(log(fabs(z1.y))/n);
 }
 else
 {
  r=exp(log(z1.x*z1.x+z1.y*z1.y)/n/2);
  cs=atan(z1.y/z1.x);
 }
 if(z1.x<0) cs+=M_PI;
 cs/=n;
 sc=2*M_PI/n;
 c=cos(cs);   s=sin(cs);
 c1=cos(sc);  s1=sin(sc);
 ck=1;        sk=0;
 z[0].x=c*r;  z[0].y=s*r;
 for(j=1;j<n;j++)
 {
  cs=ck*c1-sk*s1;
  sc=sk*c1+ck*s1;
  z[j].x=(c*cs-s*sc)*r;
  z[j].y=(s*cs-c*sc)*r;
  ck=cs;
  sk=sc;
 }

}

Complex _Cpow(Complex z1, double w)
{
 double r,t;
 Complex z;

 if(z1.x==0&&z1.y==0) return z1;
 if(z1.x==0)
 {
  if(z1.y>0) t=M_PI_2;
  else t=-M_PI_2;
 }
 else
 {
  if(z1.x>0)  t=atan(z1.y/z1.x);
  else
  {
 if(z1.y>=0)  t=atan(z1.y/z1.x)+M_PI;
 else  t=atan(z1.y/z1.x)-M_PI;
  }
 }
 r=exp(w*log(sqrt(z1.x*z1.x+z1.y*z1.y)));
 z.x=r*cos(w*t);
 z.y=r*sin(w*t);
 return z;

}

double _LAG(double *x, double *y, double t, int n)
{
 int i,j;
 double p,s=0.0;

 for(i=0;i<n;i++)
 {
  p=1.0;
  for(j=0;j<n;j++)
 if(i!=j)
  p*=(t-x[j])/(x[i]-x[j]);
  s+=p*y[i];

 }
 return s;

}

double _NEWT(double *x, double *y, int n, double t)
{
 int i,j;
 double *s,p;

 s=(double *)calloc(n+1,sizeof(double));
 if(s==NULL) Error("calloc error in _NEWT");
 for(i=1;i<=n;i++)
  s[i]=y[i];
 for(j=1;j<=n-1;j++)
  for(i=n;i>=j+1;i--)
 s[i]=(s[i]-s[i-1])/(x[i]-x[i-j]);
 p=y[n];
 for(i=n-1;i>=1;i--)
  p=p*(t-x[i])+s[i];
 free(s);
 return p;

}

int _Mgauss(double **a, double *b, int n, double ep)
{
 int i,j,k,l;
 double c,t;

 for(k=1;k<=n;k++)
 {
  c=0.0;
  for(i=k;i<=n;i++)
 if(fabs(a[i-1][k-1])>fabs(c))
 {
  c=a[i-1][k-1];
  l=i;
 }
  if(fabs(c)<=ep) return 0;
  if(l!=k)
  {
 for(j=k;j<=n;j++)
 {
  t=a[k-1][j-1];
  a[k-1][j-1]=a[l-1][j-1];
  a[l-1][j-1]=t;
 }
 t=b[k-1];
 b[k-1]=b[l-1];
 b[l-1]=t;
  }
  c=1/c;
  for(j=k+1;j<=n;j++)
  {
 a[k-1][j-1]=a[k-1][j-1]*c;
 for(i=k+1;i<=n;i++)
  a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1];
  }
  b[k-1]*=c;
  for(i=k+1;i<=n;i++)
 b[i-1]-=b[k-1]*a[i-1][k-1];
 }
 for(i=n;i>=1;i--)
  for(j=i+1;j<=n;j++)
 b[i-1]-=a[i-1][j-1]*b[j-1];

  return 1;
}

int _Mgauss2(double **a, double **b, int n, int m, double ep)
{
 int i,j,k,l;
 double c,t;

 for(k=1;k<=n;k++)
 {
  c=0;
  for(j=k;j<=n;j++)
 if(fabs(a[k-1][j-1])>fabs(c))
 {
  c=a[k-1][j-1];
  l=j;

 }
  if(fabs(c)<=ep)  return 0;
  if(l!=k)
  {
 for(i=k;i<=n;i++)
 {
  t=a[i-1][k-1];
  a[i-1][k-1]=a[i-1][l-1];
  a[i-1][l-1]=t;

 }
 for(i=1;i<=m;i++)
 {
  t=b[i-1][k-1];
  b[i-1][k-1]=b[i-1][l-1];
  b[i-1][l-1]=t;

 }
 c=1/c;
 for(i=k+1;i<=n;i++)
 {
  a[i-1][k-1]*=c;
  for(j=k+1;j<=n;j++)
   a[i-1][j-1]-=a[k-1][j-1]*a[i-1][k-1];

 }
 for(i=1;i<=m;i++)
 {
  b[i-1][k-1]*=c;
  for(j=k+1;j<=n;j++)
   b[i-1][j-1]-=a[k-1][j-1]*b[i-1][k-1];

 }

  }
  for(i=n;i>=1;i--)
  {
 for(j=i+1;j<=n;j++)
 for(k=1;k<=m;k++)
  b[k-1][i-1]-=a[j-1][i-1]*b[k-1][j-1];

  }

 }

 return 1;

}

int _Mdjn(double **a, double **c, int n, int m)
{
 int i,j,k,k1,k2,k3;

 if(fabs(a[0][0])==0) return 0;
 for(i=2;i<=n;i++) a[i-1][0]/=a[0][0];

 for(i=2;i<=n-1;i++)
 {
  for(j=2;j<=i;j++)
 a[i-1][i-1]-=a[i-1][j-2]*a[i-1][j-2]*a[j-2][j-2];
  for(k=i+1;k<=n;k++)
  {
 for(j=2;j<=i;j++)
  a[k-1][i-1]-=a[k-1][j-2]*a[i-1][j-2]*a[j-2][j-2];
 if(fabs(a[i-1][i-1])==0) return 0;
 a[k-1][i-1]/=a[i-1][i-1];

  }

 }
 for(j=2;j<=n;j++)
  a[n-1][n-1]-=a[n-1][j-2]*a[n-1][j-2]*a[j-2][j-2];
 for(j=1;j<=m;j++)
  for(i=2;i<=n;i++)
 for(k=2;k<=i;k++)
  c[i-1][j-1]-=a[i-1][k-2]*c[k-2][j-1];

 for(i=2;i<=n;i++)
  for(j=i;j<=n;j++)
 a[i-2][j-1]=a[i-2][i-2]*a[j-1][i-2];

 if(fabs(a[n-1][n-1])==0)  return 0;
 for(j=1;j<=m;j++)
 {
  c[n-1][j-1]/=a[n-1][n-1];
  for(k=2;k<=n;k++)
  {
 k1=n-k+2;
 for(k2=k1;k2<=n;k2++)
 {
  k3=n-k+1;
  c[k3-1][j-1]-=a[k3-1][k2-1]*c[k2-1][j-1];

 }
 c[k3-1][j-1]/=a[k3-1][k3-1];
  }

 }

 return 1;

}

double _NOR(double x,int l)
{
 double rm,rn,rp,rt,pq,rrt,x1,x2,y,s,h;
 double p1,p2,q1,q2;

 if(x==0) return 0.5;
 rp=1;
 if(x*l<0) rp=-1;
 x1=fabs(x);
 x2=x1*x1;
 y=1/sqrt(2*M_PI)*exp(-0.5*x2);
 rn=y/x1;
 if(rp>=0&&fabs(rn)<1e-12)  return 1.0;
 if(rp<=0&&rn==0)   return 0.0;

 rrt=3.5;
 if(rp==-1) rrt=2.32;
 if(x1<=rrt)
 {
  x1*=y;
  s=x1;
  rn=1.0;
  rt=0;
  do
  {
   rn+=2;
   rt=s;
   x1*=x2/rn;
   s+=x1;
  }
  while(s!=rt);
  pq=0.5+s;
  if(rp==-1) pq=0.5-s;
  return pq;
 }
 q1=x1;
 p2=y*x1;
 rn=1;
 p1=y;
 q2=x2+1;
 if(rp!=-1)
 {
  rm=1-p1/q1;
  s=rm;
  rt=1-p2/q2;
 }
 else
 {
  rm=p1/q1;
  s=rm;
  rt=p2/q2;
 }
 while(1)
 {
  rn++;
  if(rm==rt||s==rt) return rt;
  s=x*p2+rn*p1;
  p1=p2;  p2=s;
  s=x*q2+rn*q1;
  q1=q2;  q2=s;
  s=rm;
  rm=rt;
  rt=1-p2/q2;
  if(rp==-1) rt=p2/q2;
 }

 return rt;

}

double _Mdet(double **a, int n)
{
 int i,j,k,i0,j0,sgn;
 double q,mp,pv,pr,t,det;

 sgn=1;
 pr=1.0;
 for(k=1;k<=n-1;k++)
 {
  mp=0.0;
  for(i=k;i<=n;i++)
   for(j=k;j<=n;j++)
   {
    pv=a[i-1][j-1];
    if(fabs(pv)>=fabs(mp))
    {
     mp=pv;
     i0=i;
     j0=j;
    }
   }
   if(mp==0.0)
    Error("Failed in _Mdet!");
   if(i0!=k)
   {
    sgn=-sgn;
    for(j=k;j<=n;j++)
    {
     t=a[i0-1][j-1];
     a[i0-1][j-1]=a[k-1][j-1];
     a[k-1][j-1]=t;
    }
   }
   if(j0!=k)
   {
    sgn=-sgn;
    for(i=k;i<=n;i++)
    {
     t=a[i-1][j0-1];
     a[i-1][j0-1]=a[i-1][k-1];
     a[i-1][k-1]=t;
    }
   }
   pr=pr*mp;
   mp=-1/mp;
   for(i=k+1;i<=n;i++)
   {
    q=a[i-1][k-1]*mp;
    for(j=k+1;j<=n;j++)
     a[i-1][j-1]=a[i-1][j-1]+q*a[k-1][j-1];
   }
 }
 det=sgn*pr*a[n-1][n-1];
 return det;

}

int _Minv(double t0, double *t, double *tt, int n, int m, double **b)
{
 int i,j,k;
 double a,*c,*r,*p,s;

 c=(double *)malloc(m*sizeof(double));
 r=(double *)malloc(m*sizeof(double));
 p=(double *)malloc(m*sizeof(double));
 if(c==NULL||r==NULL||r==NULL) Error("calloc error in _Minv!");
 if(fabs(t0)==0)
 {
  free(c);
  free(r);
  free(p);
  return 0;
 }
 a=t0;
 c[0]=tt[0]/t0;
 r[0]=t[0]/t0;
 for(k=0;k<=n-2;k++)
 {
  s=0;
  for(j=1;j<=k+1;j++)
   s=s+c[k+1-j]*tt[j-1];
  s=(s-tt[k+1])/a;
  for(i=1;i<=k+1;i++)
   p[i-1]=c[i-1]+s*r[k+1-i];
  c[k+1]=-s;
  s=0;
  for(j=1;j<=k+1;j++)
   s=s+r[k+1-j]*t[j-1];
  s=(s-t[k+1])/a;
  for(i=1;i<=k+1;i++)
  {
   r[i-1]=r[i-1]+s*c[k+1-i];
   c[k+1-i]=p[k+1-i];
  }
  r[k+1]=-s;
  a=0;
  for(j=1;j<=k+2;j++)
   a=a+t[j-1]*c[j-1];
  a=t0-a;
  if(fabs(a)==0) return 0;

 }
  b[0][0]=1/a;
  for(i=1;i<=n;i++)
  {
   b[0][i]=-r[i-1]/a;
   b[i][0]=-c[i-1]/a;
  }
  for(i=1;i<=n;i++)
   for(j=1;j<=n;j++)
   {
    b[i][j]=b[i-1][j-1]-c[i-1]*b[0][j];
    b[i][j]=b[i][j]+c[n-j]*b[0][n+1-j];
   }
 free(c);
 free(r);
 free(p);
 return 1;

}
 
void _FFT(double *fr, double *fi, int n, int flag)
{
 int mp,arg,q,cntr,p1,p2;
 int i,j,a,b,k,m;
 double sign,pr,pi,harm,x,y,t;
 double *ca,*sa;

 ca=(double *)calloc(n,sizeof(double));
 sa=(double *)calloc(n,sizeof(double));
 if(ca==NULL||sa==NULL) Error("calloc error in _FFT!");
 j=0;
 if(flag!=0)
 {
  sign=1.0;
  for(i=0;i<=n-1;i++)
  {
   fr[i]=fr[i]/n;
   fi[i]=fi[i]/n;
  }
 }
 else
  sign=-1.0;
 for(i=0;i<=n-2;i++)
 {
  if(i<j)
  {
   t=fr[i];
   fr[i]=fr[j];
   fr[j]=t;
   t=fi[i];
   fi[i]=fi[j];
   fi[j]=t;
  }
  k=n/2;
  while(k<=j)
  {
   j-=k;
   k/=2;
  }
  j+=k;
 }
 mp=0;
 i=n;
 while(i!=1)
 {
  mp+=1;
  i/=2;
 }
 harm=2*M_PI/n;
 for(i=0;i<=n-1;i++)
 {
  sa[i]=sign*sin(harm*i);
  ca[i]=cos(harm*i);
 }
 a=2;
 b=1;
 for(cntr=1;cntr<=mp;cntr++)
 {
  p1=n/a;
  p2=0;
  for(k=0;k<=b-1;k++)
  {
   i=k;
   while(i<n)
   {
    arg=i+b;
    if(k==0)
    {
     pr=fr[arg];
     pi=fi[arg];
    }
    else
    {
     pr=fr[arg]*ca[p2]-fi[arg]*sa[p2];
     pi=fr[arg]*sa[p2]+fi[arg]*ca[p2];
    }
    fr[arg]=fr[i]-pr;
    fi[arg]=fi[i]-pi;
    fr[i]+=pr;
    fi[i]+=pi;
    i+=a;
   }
   p2+=p1;
  }
  a*=2;
  b*=2;
 }
 free(ca);
 free(sa);
}

 

0 0

相关博文

我的热门文章

img
取 消
img