谁有三次样条插值的源程序？(100分)

谁有能拟合分散数据的源程序，比如三次样条插值的？

如有傅立叶变换的原程序，我再加64分，我就这么多分了。

你好

你好，你可不可以把你的关于最小二乘法进行曲线拟合的代码给我一份！我现在刚好要用，谢谢！

FFT:
PROGRAM d12r3(input,output);

LABEL 1,99;
CONST
eps=1.0e-3;
np=32;
np2=34; (* np+2 *)
width=50.0;
pi=3.1415926;
TYPE
gldarray = ARRAY [1..np2] OF real;
VAR
big,per,scal,small : real;
i,j,n,nlim : integer;
data,size : gldarray;

BEGIN
n := np DIV 2;
1: writeln('Period of sinusoid in channels (2-',np:2,')');
readln(per);
IF (per <= 0.0) THEN GOTO 99;
FOR i := 1 to np DO BEGIN
data := cos(2.0*pi*(i-1)/per)
END;
realft(data,n,+1);
big := -1.0e10;
FOR i := 1 to n DO BEGIN
size := sqrt(sqr(data[2*i-1])+sqr(data[2*i]));
IF (size > big) THEN big := size
END;
size[1] := data[1];
IF (size[1] > big) THEN big := size[1];
scal := width/big;
FOR i := 1 to n DO BEGIN
nlim := round(scal*size+eps);
write(i:4,' ');
FOR j := 1 to nlim+1 DO write('*');
writeln
END;
writeln('press RETURN to continue ...');
readln;
realft(data,n,-1);
big := -1.0e10;
small := 1.0e10;
FOR i := 1 to np DO BEGIN
IF (data < small) THEN small := data;
IF (data > big) THEN big := data
END;
scal := width/(big-small);
FOR i := 1 to np DO BEGIN
nlim := round(scal*(data-small)+eps);
write(i:4,' ');
FOR j := 1 to nlim+1 DO write('*');
writeln
END;
GOTO 1;
99:
END.

Cubic Spline Interpolation:
void spline(float x[], float y[], int n, float yp1, float ypn, float y2[])
{
int i,k;
float p,qn,sig,un,*u;
u = new float[n-1];
if (yp1 > 0.99e30)
y2[1]=u[1]=0.0;
else
{
y2[1] = -0.5;
u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
}
for (i = 1;i < n - 1;i++)
{
sig=(x-x[i-1])/(x[i+1]-x[i-1]);
p=sig*y2[i-1]+2.0;
y2=(sig-1.0)/p;
u=(y[i+1]-y)/(x[i+1]-x) - (y-y[i-1])/(x-x[i-1]);
u=(6.0*u/(x[i+1]-x[i-1])-sig*u[i-1])/p;
}
if (ypn > 0.99e30)
qn=un=0.0;
else
{
qn=0.5;
un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
}
y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
for (k=n-1;k>=1;k--)
y2[k]=y2[k]*y2[k+1]+u[k];
delete u;
}

void splint(float xa[], float ya[], float y2a[], int n, float x, float *y)
{
int klo,khi,k;
float h,b,a;
klo=1;
khi=n;
while (khi-klo > 1)
{
k=(khi+klo) >> 1;
if (xa[k] > x) khi=k;
else klo=k;
}
h=xa[khi]-xa[klo];
if (h == 0.0) cout<<"Bad xa input to routine splint";
a=(xa[khi]-x)/h;
b=(x-xa[klo])/h; //Cubic spline polynomial is now evaluated.
*y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
}

接受答案了.