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) 的值发生故障。
Sample Input
2
1.0
4.0
12.0
16.0
25.0
1
-30.5893962764
5.76397083962
39.3853798058
74.3727663177
4
42.4715310246
79.5420238202
28.0282396675
-30.3627807522
-49.8363481393
-25.5101480106
7.58575761381
5
-21.9161699038
-48.469304271
-24.3188578417
-2.35085940324
-9.70239202086
-47.2709510623
-93.5066246072
-82.5073836498
0
1.0
4.0
12.0
16.0
25.0
1
-30.5893962764
5.76397083962
39.3853798058
74.3727663177
4
42.4715310246
79.5420238202
28.0282396675
-30.3627807522
-49.8363481393
-25.5101480106
7.58575761381
5
-21.9161699038
-48.469304271
-24.3188578417
-2.35085940324
-9.70239202086
-47.2709510623
-93.5066246072
-82.5073836498
0
Sample Output
2
1
1
6
1
1
6
HINT
……什么鬼。。
D如此之小。。
给出D+3个方程,有一个错掉的,问错了哪个。
那当然是暴力枚举了。。然后剩下(D+2)个方程,
由于只有(D+1)个元,所以会有一个多余。
那么只要任意取(D+1)个方程,解出它们的解,然后判断剩下的方程满不满足即可。
假如满足,根据题意,当前枚举到的不合法方程一定就是答案了。
好久没打高斯消元都快不会了= =
只好手推一遍2333
#include<bits/stdc++.h>
using namespace std;
const double
eps=1e-3;
const int
N=15;
int D,ksm[9][6];
double a[N][N],tG[N];
bool gauss(int tt){
double x[N];
for (int k=1;k<=D+1;k++){
for (int i=1;i<=D+1;i++) x[i]=a[i][k]/a[k][k];
x[k]=0.0;
for (int i=k+1;i<=D+1;i++){
tG[i]-=tG[k]*x[i];
for (int j=1;j<=D+1;j++) a[i][j]-=a[k][j]*x[i];
}
}
double ANS[N];
ANS[D+1]=tG[D+1]/a[D+1][D+1];
for (int i=D;i;i--){
double tmp=0.0;
for (int j=i+1;j<=D+1;j++)
tmp+=ANS[j]*a[i][j];
ANS[i]=(tG[i]-tmp)/a[i][i];
}
double tmp=0.0;
for (int i=1;i<=D+1;i++) tmp+=ANS[i]*a[D+2][i];
if (fabs(tmp-tG[D+2])<=eps) return 1;
return 0;
}
void solve(){
double G[N];
for (int i=1;i<=D+3;i++) scanf("%lf",&G[i]);
for (int i=1;i<=D+3;i++){
int n=0;
for (int j=1;j<=D+3;j++)
if (i!=j){
tG[++n]=G[j];
for (int k=D;~k;k--) a[n][k+1]=1.0*ksm[j][k];
}
if (gauss(i)){printf("%d\n",i-1);return;}
}
}
int main(){
for (int i=1;i<9;i++){
ksm[i][0]=1;
for (int j=1;j<6;j++) ksm[i][j]=ksm[i][j-1]*i;
}
while (scanf("%d",&D),D) solve();
return 0;
}