一元插值，呵呵

bonZ

//下面贴出来的函数程序在borland c++5.0上都能编译运行



/*                         一元全区间不等距插值
* x----双精度实型一维数组，长度n。存放给定的n个结点的值xi,要求x0<x1<....<xn-1
* y----双精度实型一维数组，长度n。存放n个给定结点上的函数值yi=f(xi)
* n----整型变量。给定结点的个数
* t----双精度实型变量。存放给定插值点的值。
*/

double aalgrg(x,y,n,t)
int n;
double t,x[],y[];
{
        int i,j,k,m;
   double z,s;
   z = 0.0;
   if(n < 1)
           return(z);
   if(n == 1)
   {
           z = y[0];
      return(z);
   }
   if(n == 2)
   {
           z = (y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
      return(z);
   }
   i = 0;
   while((x < t) && (i < n))
           i = i + 1;
   k = i - 4;
   if(k < 0)
           k = 0;
   m = i + 3;
   if(m>n-1)
           m=n-1;
   for(i=k;i<=m;i++)
   {
           s=1.0;
      for(j=k;j<=m;j++)
              if(j != i)
                 s=s*(t-x[j])/(x-x[j]);
      z=z+s*y;
   }
   return(z);
}
/*                       一元全区间等距插值
* x0----双精度实型变量。等距结点中的第一个结点值
× h ----双精度实型变量。等距结点的步长
× n ----整型变量。等距结点的个数
× y ----双精度实型一维数组，长n。存放n个等距结点上的函数值
× t ----双精度实型变量。制定插值点的值
*/


double ablgrg(x0,h,n,y,t)
int n;
double x0,h,t,y[];
{
        int i,j,k,m;
        double z,s,xi,xj;
        float p,q;
        z=0.0;
        if(n<1)
                return(z);
        if(n==1)
        {
                z=y[0];
                return(z);
        }
        if(n==2)
        {
                z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
                return(z);       
        }
        if(t>x0)
        {
                p=(t-x0)/h;
                i=(int)p;
                q=(float)i;
                if(p>q)
                        i=i+1;       
        }
        else i=0;
        k=i-4;
        if(k<0)
                k=0;
        m=i+3;
        if(m>n-1)
                m=n-1;
        for(i=k;i<=m;i++)
        {
                s=1.0;
                xi=x0+i*h;
                for(j=k;j<=m;j++)
                        if(j!=i)
                        {
                                xj=x0+j*h;
                                s=s*(t-xj)/(xi-xj);
                        }
                z=z+s*y;
        }
        return(z);
               
}
/*                                                        一元三点不等距插值
* x----双精度实型一维数组，长n。存放给定n个不等距结点的值
× y----双精度实型一维数组，长n。存放给定n个不等距结点上的函数值
× n----整形变量。给定不等距结点的个数
* t----双精度实型变量。指定插值点的值
*/

#include <math.h>

double aclgrg(x,y,n,t)
int n;
double t,x[],y[];
{
        int i,j,k,m;
        double z,s;
        z=0.0;
        if(n<1)
                return(z);
        if(n==1)
        {
                z=y[0];
                return(z);
        }
        if(n==2)
        {
                z=(y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
                return(z);       
        }       
        if(t<=x[1])
                {k=0;m=2;}
        else if(t<=x[n-2])
                {k=n-3;m=n-1;}
        else
                {
                        k=1;m=n;
                        while(m-k != 1)
                        {
                                i=(k+m)/2;
                                if(t<x[i-1]) m=i;
                                else k=i;       
                        }       
                        k=k-1;m=m-1;
                        if(fabs(t-x[k])<fabs(t-x[m])) k=k-1;
                        else m=m+1;
                }       
        z=0.0;
        for(i=k;i<=m;i++)
        {
                s=1.0;
                for(j=k;j<=m;j++)
                        if(j!=i) s=s*(t-x[j])/(x-x[j]);
                z=z+s*y;       
        }
        return(z);
}
/*                                                一元三点等距插值
* x0----双精度实型变量。给定n个等距结点中的第一个结点值
* h ----双精度实型变量。给定n个等距结点的步长
* n ----int，结点个数
* y ----双精度实型一维数组，长n。存放给定的n个等距结点上的函数值
* t ----双精度实型变量。存放制定插值点的值
*/
#include <math.h>
double adlgrg(x0,h,n,y,t)
int n;
double x0,h,t,y[];
{
        int i,j,k,m;
        double z,s,xi,xj;
        z=0.0;
        if(n<1)
                return(z);
        if(n==1)
        {
                z=y[0];
                return(z);
        }
        if(n==2)
        {
                z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
                return(z);       
        }       
        if(t<=x0+h)
                {k=0;m=2;}
        else if(t>=x0+(n-3)*h)
                {k=n-3;m=n-1;}
        else
                {
                        i=(int)((t-x0)/h)+1;
                        if(fabs(t-x0-i*h)>=fabs(t-x0-(i-1)*h))
                                {k=i-2;m=i;}
                        else
                                {k=i-1;m=i+1;}
                }                                
        z=0.0;
        for(i=k;i<=m;i++)
        {
                s=1.0;
                xi=x0+i*h;
                for(j=k;j<=m;j++)
                        if(j!=i)
                        {
                                xj=x0+j*h;
                                s=s*(t-xj)/(xi-xj);
                        }
                z=z+s*y;
        }
        return(z);       
}

[/B][/B][/I][/B]

[此贴子已经被作者于2004-4-16 2:05:49编辑过]

bonZ

再来
/*                                                        有理不等距插值
* x----双精度实型一维数组，长n。存放给定的n个不等距结点值
* y----双精度实型一维数组，长n。存放给定的n个不等距结点上的函数值
* n----int,结点个数
* t----双精度实型变量，给定插值点的值
*/
#include <math.h>

double aepq(x,y,n,t)
int n;
double t,x[],y[];
{
        int i,j,k,m,l;
        double z,h,b[8];
        z=0.0;
          if(n<1)
                return(z);
        if(n==1)
        {
                z=y[0];
                return(z);
        }               
        if(n<=8)
                {k=0;m=n;}
        else if(t<x[4]) {k=0;m=8;}
        else if(t>x[n-5]) {k=n-8;m=8;}
        else
        {
                k=1;j=n;
                while(j-k != 1)
                {
                        i=(k+j)/2;
                        if(t<x[i-1]) j=i;
                        else k=i;       
                }       
                k=k-4;m=8;
        }
        b[0]=y[k];
        for(i=2;i<=m;i++)
        {
                h=y[i+k-1];l=0;
                for(j=1;j<=i-1;j++)
                {
                        if(fabs(h-b[j-1])+1.0 == 1.0) l=l+1;
                        else h=(x[i+k-1]-x[j+k-1])/(h-b[j-1]);       
                }       
                b[i-1]=h;
                if(l!=0) b[i=1]=1.0e+35;
        }
        z=b[m-1];
        for(i=m-1;i>=1;i--) z=b[i-1]+(t-x[i+k-1])/z;
        return(z);
}
/*                        有理等距插值
* x0----双精度实型变量。n个等距结点中的第一个结点值
* h ----双精度实型变量。n个等距结点的步长
* n ----int,等距结点的个数
* y ----双精度实型一维数组，长n。存放n个等距结点上的函数值
* t ----双精度实型变量。指定插值点的值
*/

#include <math.h>
double afpq(x0,h,n,y,t)
int n;
double x0,h,t,y[];
{
        int i,j,k,m,l;
        double z,hh,xi,xj,b[8];
        z=0.0;
        if(n<1) return(z);
        if(n==1) {z=y[0];return(z);}
        if(n<=8) {k=0;m=n;}
        else if(t<(x0+4.0*h)) {k=0;m=8;}
        else if(t>(x0+(n-5)*h)) {k=n-8;m=8;}
        else {k=(int)((t-x0)/h)-3;m=8;}
        b[0]=y[k];
        for(i=2;i<=m;i++)
        {
                hh=y[i+k-1];l=0;
                for(j=1;j<=i-1;j++)
                        if(l==0)
                        {
                                if(fabs(hh-b[j-1])+1.0==1.0) l=l+1;
                                else
                                {
                                        xi=x0+(i+k-1)*h;
                                        xj=x0+(j+k-1)*h;
                                        hh=(xi-xj)/(hh-b[j-1]);       
                                }       
                        }       
                b[i-1]=hh;
                if(l!=0) b[i-1]=1.0e+35;
        }       
        z=b[m-1];
        for(i=m-1;i>=1;i--)
                z=b[i-1]+(t-(x0+(i+k-1)*h))/z;
        return(z);
}
/*                                                埃尔米特不等距插值
* 给定n个不等距结点xi上的函数值f(xi)和一阶导数值f'(xi),用埃尔米特（Hermite）插值公式计算指定插值点t处的函数近似值f(t)
* x ----双精度实型一维数组，长n。存放给定n个不等距结点的值
* y ----双精度实型一维数组，长n。存放给定n个不等距结点上的函数值
* dy----双精度实型一维数组，长n。存放给定n个不等距结点上的一阶导数
* n ----整形变量。给定不等距结点的个数
* t ----双精度实型变量。指定插值点的值
*/
double aghmt(x,y,dy,n,t)
int n;
double t,x[],y[],dy[];
{
        int i,j;
        double z,p,q,s;
        z=0.0;
        for(i=1;i<=n;i++)
        {
                s=1.0;
                for(j=1;j<=n;j++)
                        if(j!=i) s=s*(t-x[j-1])/(x[i-1]-x[j-1]);
                s=s*s;
                p=0.0;
                for(j=1;j<=n;j++)
                        if(j!=i) p=p+1.0/(x[i-1]-x[j-1]);
                q=y[i-1]+(t-x[i-1])*(dy[i-1]-2.0*y[i-1]*p);
                z=z+q*s;       
        }       
        return(z);
}
/*                                        埃尔米特等距插值
* 给定n个等距结点xi上的函数值f(xi)和一阶导数值f'(xi),用埃尔米特（Hermite）插值公式计算指定插值点t处的函数近似值f(t)
* x0----双精度实型变量。n个等距结点中的第一个结点值
* h ----双精度实型变量。n个等距结点的步长
* n ----int,等距结点的个数
* y ----双精度实型一维数组，长n。存放n个等距结点上的函数值
* dy----dy----双精度实型一维数组，长n。存放给定n个等距结点上的一阶导数
* t ----双精度实型变量。指定插值点的值
*/
double ahhmt(x0,h,n,y,dy,t)
int n;
double x0,h,t,y[],dy[];
{
        int i,j;
        double z,s,p,q;
        z=0.0;
        for(i=1;i<=n;i++)
        {
                s=1.0;q=x0+(i-1)*h;
                for(j=1;j<=n;j++)
                {
                        p=x0+(j-1)*h;
                        if(j!=i) s=s*(t-p)/(q-p);       
                }       
                for(j=1;j<=n;j++)
                        if(j!=i) p=p+1.0/(q-(x0+(j-1)*h));
                q=y[i-1]+(t-q)*(dy[i-1]-2.0*y[i-1]*p);
                z=z+q*s;
        }       
        return(z);
}

was1100

强

bonZ

继续
/*                                                        埃特金不等距逐步插值
* 给定n个不等距结点xi上的函数值f(xi)及精度要求,用埃特金(Aitken)插值公式计算指定插值点t处的函数近似值f(t)
* x ----双精度实型一维数组，长n。存放给定n个不等距结点的值
* y ----双精度实型一维数组，长n。存放给定n个不等距结点上的函数值
* n ----整形变量。给定不等距结点的个数
* t ----双精度实型变量。指定插值点的值
* eps---双精度实型变量。插值的精度要求
*/

#include <math.h>
double aiatkn(x,y,n,t,eps)
int n;
double t,eps,x[],y[];
{
        int i,j,k,m,l;
        double z,xx[10],yy[10];
        z=0.0;
        if(n<1) return(z);
        if(n==1) {z=y[0];return(z);}
        m=10;
        if(m>n) m=n;
        if(t<=x[0]) k=1;
        else if(t>=x[n-1]) k=n;
        else
        {
                k=1;j=n;
                while((k-j!=1) && (k-j!=-1))
                {
                        l=(k+j)/2;
                        if(t<x[l-1]) j=l;
                        else k=l;
                }
                if(fabs(t-x[l-1])>fabs(t-x[j-1])) k=j;
        }
        j=1;l=0;
        for(i=1;i<=m;i++)
        {
                k=k+j*l;
                if((k<1)||(k>n)) {l=l+1;j=-j;k=k+j*l;}
                xx[i-1]=x[k-1];
                yy[i-1]=y[k-1];
                l=l+1;j=-j;
        }
        i=0;
        do
        {
                i=i+1;z=yy;
                for(j=0;j<=i-1;j++)
                        z=yy[j]+(t-xx[j])*(yy[j]-z)/(xx[j]-xx);
                yy=z;       
        }
        while((i!=m-1) && (fabs(yy-yy[i-1])>eps));
        return(z);
       
}

/*                                                        埃特金等距逐步插值
* x0----双精度实型变量。n个等距结点中的第一个结点值
* h ----双精度实型变量。n个等距结点的步长
* n ----int,等距结点的个数
* y ----双精度实型一维数组，长n。存放n个等距结点上的函数值
* t ----双精度实型变量。指定插值点的值
* eps---双精度实型变量。插值的精度要求
*/
#include <math.h>
double ajatkn(x0,h,n,y,t,eps)
int n;
double x0,h,t,eps,y[];
{
        int i,j,k,m,l;
        double z,xx[10],yy[10];
        z=0.0;
        if(n<1) return(z);
        if(n==1) {z=y[0];return(z);}
        m=10;
        if(m>n) m=n;
        if(t<x0) k=1;
        else if(t>=x0+(n-1)*h) k=n;
        else
        {
                k=1;j=n;
                while((k-j!=1) && (k-j!=-1))
                {
                        l=(k+j)/2;
                        if(t<x0+(l-1)*h) j=l;
                        else k=l;       
                }               
                if(fabs(t-(x0+(l-1)*h)) > fabs(t-(x0+(j-1)*h))) k=j;
        }       
        j=1;l=0;
        for(i=1;i<=m;i++)
        {
                k=k+j*l;
                if((k<1) || (k>n))
                        {l=l+1;j=-j;k=k+j*l;}
                xx[i-1]=x0+(k-1)*h;yy[i-1]=y[k-1];
                l=l+1;j=-j;       
        }
        i=0;
        do
        {
                i=i+1;z=yy;
                for(j=0;j<=i-1;j++)
                        z=yy[j]+(t-xx[j])*(yy[j]-z)/(xx[j]-xx);
                yy=z;       
        }
        while((i!=m-1) && (fabs(yy-yy[i-1])>eps));
        return(z);
}

/*                                                光滑不等距插值
* 给定n个不等钜结点xi上的函数值yi=f(xi),利用阿克玛（Akima）方法，计算指定子区间上的三次插值多项式与指定插值点上的函数值
* x ----双精度实型一维数组，长n。存放给定的n个不等距结点值
* y ----双精度实型一维数组，长n。存放给定的n个不等距结点上的函数值
* n ----int,给定不等距节点的个数
* k ----int,若k>=0，表示只计算第k个子区间[xk,xk+1]上的三次多项式的系数s0,s1,s2,s3.若k<0,表示需要计算指定插值点t处的函数近似值f(t)和系数
* t ----双精度实型变量。制定插值点值。当k>=0时，此形参对本函数不起作用，可任意
* s ----双精度实型一维数组，长5。s0到s3返回三次多项式的系数，s4返回制定插值点t处的函数近似值(k<0时)或任意(k>=0时)
*/
#include <math.h>
void akspl(x,y,n,k,t,s)
int n,k;
double t,x[],y[],s[5];
{
        int kk,m,l;
        double u[5],p,q;
        s[4]=0.0;s[0]=0.0;s[1]=0.0;s[2]=0.0;s[3]=0.0;
        if(n<1) return;
        if(n==1) {s[0]=y[0];s[4]=y[0];return;}
        if(n==2)
                {
                        s[0]=y[0];s[1]=(y[1]-y[0])/(x[1]-x[0]);
                        if(k<0)
                                s[4]=(y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
                        return;       
                }       
        if(k<0)
                {
                        if(t<=x[1]) kk=0;
                        else if(t>=x[n-1]) kk=n-2;
                        else
                                {
                                        kk=1;m=n;
                                        while(((kk-m)!=1) && ((kk-m)!=-1))
                                                {
                                                        l=(kk+m)/2;
                                                        if(t<x[l-1]) m=l;
                                                        else kk=l;
                                                }       
                                        kk=kk-1;
                                }       
                }
        else kk=k;
        if(kk>=n-1) kk=n-2;
        u[2]=(y[kk+1]-y[kk])/(x[kk+1]-x[kk]);
        if(n==3)
                {
                        if(kk==0)
                                {
                                        u[3]=(y[2]-y[1])/(x[2]-x[1]);
                                        u[4]=2.0*u[3]-u[2];
                                        u[1]=2.0*u[2]-u[3];
                                        u[0]=2.0*u[1]-u[2];       
                                }       
                        else
                                {
                                        u[1]=(y[1]-y[0])/(x[1]-x[0]);
                                        u[0]=2.0*u[1]-u[2];
                                        u[3]=2.0*u[2]-u[1];
                                        u[4]=2.0*u[3]-u[2];       
                                }       
                }
        else
                {
                        if(kk<=1)
                                {
                                        u[3]=(y[kk+2]-y[kk+1])/(x[kk+2]-x[kk+1]);
                                        if(kk==1)
                                                {
                                                        u[1]=(y[1]-y[0])/(x[1]-x[0]);
                                                        u[0]=2.0*u[1]-u[2];
                                                        if(n==4) u[4]=2.0*u[3]-u[2];
                                                        else u[4]=(y[4]-y[3])/(x[4]-x[3]);       
                                                }       
                                        else
                                                {
                                                        u[1]=2.0*u[2]-u[3];
                                                        u[0]=2.0*u[1]-u[2];
                                                        u[4]=(y[3]-y[2])/(x[3]-x[2]);       
                                                }
                                }       
                        else if(kk>=(n-3))
                                {
                                        u[1]=(y[kk]-y[kk-1])/(x[kk]-x[kk-1]);
                                        if(kk==(n-3))
                                                {
                                                        u[3]=(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
                                                        u[4]=2.0*u[3]-u[2];
                                                        if(n==4) u[0]=2.0*u[1]-u[2];
                                                        else u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);       
                                                }       
                                        else
                                                {
                                                        u[3]=2.0*u[2]-u[1];
                                                        u[4]=2.0*u[3]-u[2];
                                                        u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);       
                                                }
                                }
                        else
                                {
                                        u[1]=(y[kk]-y[kk-1])/(x[kk]-x[kk-1]);
                                        u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);
                                        u[3]=(y[kk+2]-y[kk+1])/(x[kk+2]-x[kk+1]);
                                        u[4]=(y[kk+3]-y[kk+2])/(x[kk+3]-x[kk+2]);       
                                }
                }
        s[0]=fabs(u[3]-u[2]);
        s[1]=fabs(u[0]-u[1]);
        if((s[0]+1.0==1.0) && (s[1]+1.0==1.0))
                p=(u[1]+u[2])/2.0;
        else p=(s[0]*u[1]+s[1]*u[2])/(s[0]+s[1]);
        s[0]=fabs(u[3]-u[4]);
        s[1]=fabs(u[2]-u[1]);
        if((s[0]+1.0==1.0) && (s[1]+1.0==1.0))
                q=(u[2]+u[3])/2.0;
        else q=(s[0]*u[2]+s[1]*u[3])/(s[0]+s[1]);
        s[0]=y[kk];
        s[1]=p;
        s[3]=x[kk+1]-x[kk];
        s[2]=(3.0*u[2]-2.0*p-q)/s[3];
        s[3]=(q+p-2.0*u[2])/(s[3]*s[3]);
        if(k<0)
                {
                        p=t-x[kk];
                        s[4]=s[0]+s[1]*p+s[2]*p*p+s[3]*p*p*p;       
                }
        return;               
}

/*                                                光滑不等距插值
* 给定n个不等钜结点xi上的函数值yi=f(xi),利用阿克玛（Akima）方法，计算指定子区间上的三次插值多项式与指定插值点上的函数值
* x0----双精度实型变量。n个等距结点中的第一个结点值
* h ----双精度实型变量。n个等距结点的步长
* n ----int,等距结点的个数
* y ----双精度实型一维数组，长n。存放n个等距结点上的函数值
* k ----int,若k>=0，表示只计算第k个子区间[xk,xk+1]上的三次多项式的系数s0,s1,s2,s3.若k<0,表示需要计算指定插值点t处的函数近似值f(t)和系数
* t ----双精度实型变量。指定插值点的值
* s ----双精度实型一维数组，长5。s0到s3返回三次多项式的系数，s4返回制定插值点t处的函数近似值(k<0时)或任意(k>=0时)
*/
#include <math.h>

void alspl(x0,h,n,y,k,t,s)
int n,k;
double t,x0,h,y[],s[5];
{
        int kk,m,l;
        double u[5],p,q;
        s[4]=0.0;s[0]=0.0;s[1]=0.0;s[2]=0.0;s[3]=0.0;
        if(n<1) return;
        if(n==1) {s[0]=y[0];s[4]=y[0];return;}
        if(n==2)
                {
                        s[0]=y[0];s[1]=(y[1]-y[0])/h;
                        if(k<0)
                                s[4]=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
                        return;       
                }
        if(k<0)
                {
                        if(t<=x0+h) kk=0;
                        else if(t>=x0+(n-1)*h) kk=n-2;
                        else
                                {
                                        kk=1;m=n;
                                        while(((kk-m)!=1) && ((kk-m)!=-1))
                                                {
                                                        l=(kk+m)/2;
                                                        if(t<x0+(l-1)*h) m=l;
                                                        else kk=l;
                                                }       
                                        kk=kk-1;
                                }                                       
                }
        else kk=k;
        if(kk>=n-1) kk=n-2;
        u[2]=(y[kk+1]-y[kk])/h;
        if(n==3)
                {
                        if(kk==0)
                                {
                                        u[3]=(y[2]-y[1])/h;
                                        u[4]=2.0*u[3]-u[2];
                                        u[1]=2.0*u[2]-u[3];
                                        u[0]=2.0*u[1]-u[2];       
                                }       
                        else
                                {
                                        u[1]=(y[1]-y[0])/h;
                                        u[0]=2.0*u[1]-u[2];
                                        u[3]=2.0*u[2]-u[1];
                                        u[4]=2.0*u[3]-u[2];       
                                }       
                }
        else
                {
                        if(kk<=1)
                                {
                                        u[3]=(y[kk+2]-y[kk+1])/h;
                                        if(kk==1)
                                                {
                                                        u[1]=(y[1]-y[0])/h;
                                                        u[0]=2.0*u[1]-u[2];
                                                        if(n==4) u[4]=2.0*u[3]-u[2];
                                                        else u[4]=(y[4]-y[3])/h;       
                                                }       
                                        else
                                                {
                                                        u[1]=2.0*u[2]-u[3];
                                                        u[0]=2.0*u[1]-u[2];
                                                        u[4]=(y[3]-y[2])/h;       
                                                }
                                }       
                        else if(kk>=(n-3))
                                {
                                        u[1]=(y[kk]-y[kk-1])/h;
                                        if(kk==(n-3))
                                                {
                                                        u[3]=(y[n-1]-y[n-2])/h;
                                                        u[4]=2.0*u[3]-u[2];
                                                        if(n==4) u[0]=2.0*u[1]-u[2];
                                                        else u[0]=(y[kk-1]-y[kk-2])/h;       
                                                }       
                                        else
                                                {
                                                        u[3]=2.0*u[2]-u[1];
                                                        u[4]=2.0*u[3]-u[2];
                                                        u[0]=(y[kk-1]-y[kk-2])/h;       
                                                }
                                }
                        else
                                {
                                        u[1]=(y[kk]-y[kk-1])/h;
                                        u[0]=(y[kk-1]-y[kk-2])/h;
                                        u[3]=(y[kk+2]-y[kk+1])/h;
                                        u[4]=(y[kk+3]-y[kk+2])/h;       
                                }
                }
        s[0]=fabs(u[3]-u[2]);
        s[1]=fabs(u[0]-u[1]);
        if((s[0]+1.0==1.0) && (s[1]+1.0==1.0))
                p=(u[1]+u[2])/2.0;
        else p=(s[0]*u[1]+s[1]*u[2])/(s[0]+s[1]);
        s[0]=fabs(u[3]-u[4]);
        s[1]=fabs(u[2]-u[1]);
        if((s[0]+1.0==1.0) && (s[1]+1.0==1.0))
                q=(u[2]+u[3])/2.0;
        else q=(s[0]*u[2]+s[1]*u[3])/(s[0]+s[1]);
        s[0]=y[kk];
        s[1]=p;
        s[3]=h;
        s[2]=(3.0*u[2]-2.0*p-q)/s[3];
        s[3]=(q+p-2.0*u[2])/(s[3]*s[3]);
        if(k<0)
                {
                        p=t-(x0+kk*h);
                        s[4]=s[0]+s[1]*p+s[2]*p*p+s[3]*p*p*p;       
                }
        return;                                                                                       
}

chenhorse

<>都是你自己写的嘛？</P>
<>强！</P>