Description
Abacus教授刚刚完成了一个制作数表的计算引擎的设计。它被设计用于同时计算一个多项式在许多点的取值。例如对于多项式 f(x)=x^2+2x+1 ,一种可能的计算结果是 f(0)=1,f(1)=4,f(2)=9.f(3)=16,f(4)=25 。不幸的是,引擎存在一个故障使得计算出的值总有一个是错的,例如对于上述多项式,它可能输出 1,4,12,16,25 而不是 1,4,9,16,25 。请你帮教授找出发生故障的是哪个点值。
Input
输入包含多组测试数据。
每组数据第一行包含一个正整数 d 表示多项式的度数,即多项式最高次项的项数,保证 d≤5 。
接下来 d+3 行,每行一个实数,第 i 行表示输出的 f(i) 的值,保证-100.0≤f(i)≤100.0 。
你可以认为恰好只有一个点值出故障,且与实际值的误差超过 1.0 。
由于不可避免的误差,其他数字与精确值的误差不超过10^(-6) 。
输入以一个零作为结束。
Output
对于每组数据,输出一个非负整数 i 表示 f(i) 的值发生故障。
Solution
高斯消元裸题。
Code
#include<bits/stdc++.h>
using namespace std;
const double eps=1e-6;
double a[110][110];
double b[110];
double ans[110];
int m,n;
bool guass(int n2){
for(int i=0;i<n;i++){
if(fabs(a[i][i])<=eps){
bool ok=0;
for(int j=i+1;j<n;j++){
if(fabs(a[j][i])>eps){
ok=1;
for(int k=0;k<=n;k++) swap(a[j][k],a[i][k]);
break;
}
}
if(!ok) return 0;
}
for(int j=i+1;j<n;j++){
double rt=a[j][i]/a[i][i];
for(int k=0;k<=n;k++)
a[j][k]-=a[i][k]*rt;
}
}
for(int i=n-1;~i;i--){
ans[i]=a[i][n];
for(int j=i+1;j<n;j++) ans[i]-=a[i][j]*ans[j];
ans[i]/=a[i][i];
}
double sum=0,s=1;
for(int i=n-1;i>=0;i--,s*=(n2-1)) sum+=ans[i]*s;//,cout<<ans[i]<<" ";cout<<endl;
if(fabs(sum-b[n2])<=0.001) return 1;
return 0;
}
bool check(int x){
int tot=0,las;
for(int i=1;i<=m+3;i++)
if(i!=x){
if(tot<=m){
double sum=1;
for(int j=n-1;j>=0;j--)
a[tot][j]=sum,sum*=(i-1);
a[tot][n]=b[i];
tot++;
}
else las=i;
}
//cout<<a<<endl;
return guass(las);
}
int main(){
while(114514){
scanf("%d",&m);
if(!m) break;
n=m+1;
for(int i=1;i<=m+3;i++)
scanf("%lf",&b[i]);
for(int i=1;i<=m+3;i++)
if(check(i)){
//cout<<i<<endl;
printf("%d\n",i-1);
break;
}
}
}