Newton均匀插值多项式

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>

#define Q(a,b) Q[(a)*(a+1)/2+b]

ifstream infile("in.txt");
ofstream outfile("out.txt");

int n,i,j;
double *x,*root,*Q,*b;

void input()
{
infile>>n;
x=new double[n];
root=new double[n-1];
b=new double[n];
Q=new double[n*(n+1)/2];

for(i=0;i<n;i++)
{
infile>>x[i]>>Q(i,0);
}
}

void setroot()
{
double t;
for(i=0;i<n;i++)
for(j=i;j<n;j++)
if(x[j]<x[i])
{
t=x[j];
x[j]=x[i];
x[i]=t;
t=Q(j,0);
Q(j,0)=Q(i,0);
Q(i,0)=t;
}
for(i=0;i<n-1;i++)
{
root[i]=(x[i]+x[i+1])/2;
}

}

void Newton()
{
for(i=1;i<n;i++)
for(j=1;j<=i;j++)
Q(i,j)=(Q(i,j-1)-Q(i-1,j-1))/(x[i]-x[i-j]);
b[n-1]=Q(n-1,n-1);
}

void f(double x0)
{
for(i=n-1;i>0;i--)
b[i-1]=Q(i-1,i-1)+b[i]*(x0-x[i-1]);
outfile<<"f("<<x0<<")的近似值为:"<<b[0]<<endl;
}

void main(void)
{
input();
setroot();
Newton();
outfile<<"差分表为："<<endl;
for(i=0;i<n;i++)
{ outfile<<i<<’ ’<<x[i]<<’ ’;
for(j=0;j<=i;j++)
outfile<<Q(i,j)<<’ ’;
outfile<<endl;
}

for(j=0;j<n-1;j++)
f(root[j]);

delete x,b,Q,root;
infile.close();
outfile.close();
}

input===============

7
-3 -4.5
6 2.1
-1 -3.2
5 2.8
-2 -1.8
3 1.7
0 0.5

output==============

差分表为：
0 -3 -4.5
1 -2 -1.8 2.7
2 -1 -3.2 -1.4 -2.05
3 0 0.5 3.7 2.55 1.53333
4 3 1.7 0.4 -0.825 -0.675 -0.368056
5 5 2.8 0.55 0.03 0.1425 0.116786 0.0606052
6 6 2.1 -0.7 -0.416667 -0.0744444 -0.0309921 -0.0184722 -0.00878638
f(-2.5)的近似值为:-1.06517
f(-1.5)的近似值为:-3.0675
f(-0.5)的近似值为:-1.86002
f(1.5)的近似值为:5.90552
f(4)的近似值为:-0.577778
f(5.5)的近似值为:4.59237