/*
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);
}
