#include<stdio.h>
#include<math.h>
#include<alloc.h>
<?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" />

 

zhuiganfa(n,a,b,c,f,x)
    int n;
    double *a,*b,*c,*f,*x;
{
    int i;
    double *B,*y;
    B=malloc(n*sizeof(double));
    y=malloc(n*sizeof(double));
    B[0]=c[0]/b[0];
    for(i=1;i<n-1;i++)
       B[i]=c[i]/(b[i]-a[i]*B[i-1]);
    y[0]=f[0]/b[0];
    for(i=1;i<n;i++)
       y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*B[i-1]);
    x[n-1]=y[n-1];
    for(i=n-2;i>=0;i--)
       x[i]=y[i]-B[i]*x[i+1];
    free(B);
    free(y);
}

 

double SplineIntp1(int n,double *X,double *Y,double ma,double mb,double x)
{
    int i;
    double *a,*b,*c,*g,*M;
    double h0,h1,f0,f1,d0,d1,y;
    a=malloc((n+1)*sizeof(double));
    b=malloc((n+1)*sizeof(double));
    c=malloc((n+1)*sizeof(double));
    g=malloc((n+1)*sizeof(double));
    M=malloc((n+1)*sizeof(double));
    b[0]=b[n]=2;
    a[n]=c[0]=1;
    for(i=1;i<n;i++)
    {
       h0=X[i]-X[i-1]; h1=X[i+1]-X[i];
       f0=(Y[i]-Y[i-1])/h0; f1=(Y[i+1]-Y[i])/h1;
       a[i]=h0/(h0+h1); b[i]=2; c[i]=h1/(h0+h1);
       g[i]=6*(f1-f0)/(h0+h1);
    }
    h0=X[1]-X[0]; g[0]=6*((Y[1]-Y[0])/h0-ma)/h0;
    h1=X[n]-X[n-1]; g[n]=6*(<?xml:namespace prefix = st1 ns = "urn:schemas-microsoft-com:office:smarttags" />mb-(Y[n]-Y[n-1])/h1)/h1;
    zhuiganfa(n+1,a,b,c,g,M);
    i=1;
    while(x>X[i]) i++;
    h0=X[i]-X[i-1]; d0=x-X[i-1]; d1=X[i]-x;
y=M[i-1]*d1*d1*d1/(6*h0)+M[i]*d0*d0*d0/(6*h0)+(Y[i-1]-M[i-1]*h0*h0/6)*d1/h0+(Y[i]-M[i]*h0*h0/6)*d0/h0;
    free(a); free(b); free(c); free(g); free(M);
    return y;
}

 

main()
{
   double X[4]={0,1,2,3},Y[4]={0,0,0,0},ma=1,mb=0,x1=0.5,x2=1.5,x3=2.5,y1,y2,y3;
   y1=SplineIntp1(3,X,Y,ma,mb,x1);
   y2=SplineIntp1(3,X,Y,ma,mb,x2);
   y3=SplineIntp1(3,X,Y,ma,mb,x3);
   printf("y(0.5)=%f\n",y1);
   printf("y(1.5)=%f\n",y2);
   printf("y(2.5)=%f\n",y3);
}