mathlab.h

一个在C语言环境中运行的贝塞尔曲线.rar

#if !defined mathlabheadfile 
#define mathlabheadfile 
#include<math.h> 
#include<time.h> 
#define N 20 
double polyvalue(double *a,int n,double x0); 
/*秦九韶多项式的值*/ 
/*a：系数数组；n：多项式次数；x0：取样点*/ 
double polyrate(double *a,int n,double x0); 
/*秦九韶多项式导数值*/ 
int dtoint(double d); 
/*四舍五入*/ 

void initrand(); 
/*初始化随机数函数*/ 
int randto(int s,int b); 
/*产生s-b的随机整数*/ 

void makeg(double y0[],double gm[],int n0); 
/*构造样条函数生成过程中的g数组*/ 
void makem(double y0[],double gm[],int n0); 
/*建立样条插值的m值*/ 
void makepoly(double y[],int ny,double poly[][4]); 
/*将y样条插值得到的多项式复制到poly*/ 
double spline(double poly[][4],int np,double x0); 
/*取得样条函数生成后在x0点的函数值*/ 

void initrand() 
/*初始化随机数函数*/ 
{ 
srand((unsigned)time(0)); 
} 
int randto(int s,int b) 
/*产生s-b的随机整数*/ 
{ 
if(s>b) 
{ 
s+=b; 
b=s-b; 
s-=b; 
} 
if(s==b) b++; 
return rand()%(b+1-s)+s; 
} 

int dtoint(double d) 
/*四舍五入*/ 
{ 
int i; 
i=(int)d+(d>((int)d+0.5)); 
if(d<0) 
i--; 
return i; 
} 
double polyvalue(double*a,int n,double x0) 
{ 
int k; 
double ans=*a; 
for(k=1;k<n+1;k++) 
ans=x0*ans+*(a+k); 
return ans; 
} 
double polyrate(double*a,int n,double x0) 
{ 
int k; 
for(k=0;k<n-1;k++) 
*(a+k+1)=*(a+k+1)+*(a+k)*x0; 
return polyvalue(a,n-1,x0); 
} 
/***************样条函数********************************/ 
void makepoly(double y[],int ny,double poly[][4]) 
/*将y样条插值得到的多项式复制到poly*/ 
{ 
if(ny>N-1)return; 
{ 
int k; 
double y0[N]; 
int j; 
double i; 
double gm[N]; 
for(k=0;k<ny;k++) 
{ 
y0[k]=y[k]; 
/*printf("%lf ",Y[i]);*/ 
} 
y0[ny]=y0[0]; 
makem(y0,gm,ny+1); 
for(j=1;j<=ny;j++) 
{ 
i=(double)j; 
poly[j-1][0]=(gm[j]-gm[j-1])/6; 
poly[j-1][1]=(i*gm[j-1]-(i-1)*gm[j])/2; 
poly[j-1][2]=(3*(gm[j]*(i-1)*(i-1)-i*i*gm[j-1])+gm[j-1]-gm[j])/6+y0[j]-y0[j-1]; 
poly[j-1][3]=((gm[j-1]*i*i*i-gm[j]*(i-1)*(i-1)*(i-1)+gm[j]*(i-1)-gm[j-1]*i)/6+i*y0[j-1]-(i-1)*y0[j]); 
} 
} 
/* { 
int k,j; 
for(j=0;j<ny;j++) 
{ 
printf("\n"); 
for(k=0;k<4;k++) 
printf("%lf\t",poly[j][k]); 
} 
}*/ 
} 
void makem(double y0[],double gm[],int ny) 
/*建立样条插值的m值*/ 
{ 
double u[N-1],l[N-1],p[N-2],t[N-2],y[N-2],d[N-1],sum; 
int i; 
makeg(y0,gm,ny); 
ny--; 
for(i=0;i<ny;i++) 
d[i]=gm[i]; 
u[0]=2; 
p[0]=0.5; 
t[0]=0.5/u[0]; 
for(i=1;i<ny-1;i++) 
{ 
l[i]=0.5/u[i-1]; 
u[i]=2-l[i]*0.5; 
p[i]=-l[i]*p[i-1]; 
t[i]=-t[i-1]*0.5/u[i]; 
} 
l[ny-1]=0.5/u[ny-2]; 
sum=0; 
for(i=0;i<ny-2;i++) 
sum=sum+t[i]*p[i]; 
u[ny-1]=2-(l[ny-1]+t[ny-2])*(0.5+p[ny-2])-sum; 
y[0]=d[0]; 
for(i=1;i<ny-1;i++) 
y[i]=d[i]-l[i]*y[i-1]; 
sum=0; 
for(i=1;i<ny-2;i++) 
sum=sum+t[i]*y[i]; 
y[ny-1]=d[ny-1]-sum-(l[ny-1]+t[ny-2])*y[ny-2]; 
gm[ny]=y[ny-1]/u[ny-1]; 
for(i=ny-1;i>0;i--) 
gm[i]=(y[i-1]-0.5*gm[i+1]-p[i-1]*gm[ny])/u[i-1]; 
gm[0]=gm[ny]; 
/*printf("\nm\n"); 
for(i=0;i<=ny;i++)printf("\n%lf",gm[i]);*/ 
} 
void makeg(double y0[],double gm[],int ny) 
/*构造样条函数生成过程中的g数组*/ 
{ 
/*int i; 
for(i=0;i<n-2;i++) 
GM[i]=3*(Y[i+2]-2*Y[i+1]+Y[i]); 
GM[n-2]=3*(Y[1]-Y[0]-Y[n-1]+Y[n-2]);*/ 
int i; 
for(i=0;i<ny-2;i++) 
gm[i]=3*(y0[i+2]-2*y0[i+1]+y0[i]); 
gm[ny-2]=3*(y0[1]-y0[0]-y0[ny-1]+y0[ny-2]); 
/*printf("g\n"); 
for(i=0;i<ny-1;i++)printf("\n%lf",gm[i]);*/ 
} 
double spline(double poly[][4],int np,double x0) 
/*取得样条函数生成后在x0点的函数值*/ 
{ 
int k=(int)x0; 
x0=fmod(x0,(double)np); 
if(k<0)k=0; 
else if(k>=np) 
{ 
x0=fmod(x0,(double)np); 
k=0; 
} 
return polyvalue(poly[k],3,x0); 
} 
#endif 

/*end of mathlab.h*/