EEOJ/OJ8/main.cpp

68 lines
1.7 KiB
C++

#include<stdio.h>
#include<math.h>
// using newton interpolation
int main(){
int n,m;
double x[100],y[100];
double t[100];
double xt;
scanf("%d",&n);
scanf("%d",&m);
for(int i=0;i<n;i++){
scanf("%lf %lf",&x[i],&y[i]);
for(int j=0;j<i;j++){
// if points are too close, just keep one of them.
if(fabs(x[i]-x[j])<1e-6){
i--;
n--;
continue;
}
}
}
// newton interpolation
for(int i=0;i<n;i++){
for(int j=i-1;j>=0;j--){
double p=1;
for(int k=j-1;k>=0;k--){
p=p*(x[i]-x[k]);
}
y[i]=y[i]-p*t[j];
}
double p=1;
for(int k=i-1;k>=0;k--){
p=p*(x[i]-x[k]);
}
/*
if value is precise enough then higher terms are not necessary (sometimes even
troublesome, when x[k] is too close to x[0], x[1],...,x[k-1] (distance < 1)
then p=(x[k]-x[0])*(x[k]-x[1])*...*(x[k]-x[k-1]) is too small, this will result
in a super large t[k] even when the residual error of a k-1 degree polynomial
is small enough.)
*/
if(fabs(y[i])<1e-9){
y[i]=0;
}
t[i]=y[i]/p;
}
int k=0;
for(int i=0;i<n;i++){
// if a coefficient is too small, then regard it as zero.
if(fabs(t[i])>1e-6){
k=i;
}
}
// degree of the polynomial
printf("%d\n",k);
n=k+1;
for(int i=0;i<m;i++){
scanf("%lf",&xt);
double yt=0;
double p=1;
for(int j=0;j<n;j++){
yt+=t[j]*p;
p*=(xt-x[j]);
}
printf("%lf\n",yt);
}
return 0;
}