Problem:
Problem D: Find the Outlier Professor Abacus has just built a new computing engine for making numerical tables. It was designed to calculate the values of a polynomial function in one variable at several points at a time. With the polynomial function f ( x ) = x 2 + 2 x + 1, for instance, a possible expected calculation result is 1 (= f (0)), 4 (= f (1)), 9 (= f (2)), 16 (= f (3)), and 25 (= f (4)). It is a pity, however, the engine seemingly has faulty components and exactly one value among those calculated simultaneously is always wrong. With the same polynomial function as above, it can, for instance, output 1, 4, 12, 16, and 25 instead of 1, 4, 9, 16, and 25. You are requested to help the professor identify the faulty components. As the first step, you should write a program that scans calculation results of the engine and finds the wrong values. Input The input is a sequence of datasets, each representing a calculation result in the following format.  d  v 0  v 1 ...  v d +2  Here, d in the first line is a positive integer that represents the degree of the polynomial, namely, the highest exponent of the variable. For instance, the degree of 4 x 5 + 3 x + 0.5 is five and that of 2.4 x + 3.8 is one. d is at most five. The following d + 3 lines contain the calculation result of f (0), f (1), ... , and f ( d + 2) in this order, where f is the polynomial function. Each of the lines contains a decimal fraction between -100.0 and 100.0, exclusive. You can assume that the wrong value, which is exactly one of f (0), f (1), ... , and f ( d +2), has an error greater than 1.0. Since rounding errors are inevitable, the other values may also have errors but they are small and never exceed 10 -6 . The end of the input is indicated by a line containing a zero. Output For each dataset, output i in a line when v i is wrong. 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 Output for the Sample Input 2
1
1
6

Solution (Python):
from functools import reduce

EPSILON = 1e-5
d = 0
Vs = []

times = lambda x,y: x*y

def interpolate(x, e):
    return sum(
        reduce(
            times, [(x-j)/(i-j) for j in range(len(Vs)) if j not in [i,x,e]]
            )*v for i,v in enumerate(Vs) if i not in [x,e])

def outlier(e):
    for i,v in enumerate(Vs):
        if i==e: continue
        if abs(interpolate(i,e) - v) < EPSILON: return False
    return True

def solve():
    for i in range(d+3):
        if not outlier(i): return i
    return -1

while True:
    d = int(input())
    if d==0: break
    Vs = [float(input()) for i in range(d+3)]
    print(solve())

Solution (Python):
while 1:
    d = input()
    if d==0:
        break
    V = [float(raw_input()) for i in xrange(d+3)]

    a = [0]*(d+3)
    cnts = [0]*(d+3)
    import itertools
    for K in itertools.combinations(range(d+3), d+1):
        for k in K:
            res = 1
            for j in K:
                if k==j: continue
                res *= k-j
            a[k] = V[k] / res
        for i in xrange(d+3):
            if i in K: continue
            res = 0.
            for k in K:
                tmp = 1
                for j in K:
                    if k==j: continue
                    tmp *= i-j
                res += a[k]*tmp
            if abs(V[i]-res) > 0.5:
                cnts[i] += 1
    print cnts.index(max(cnts))

Solution (Python):
from itertools import combinations

def build(X, Y):
    A = []
    for x in X:
        res = 1
        for xi in X:
            if x == xi:
                continue
            res *= x - xi
        A.append(Y[x] / res)
    return A

def calc(X, A, x):
    base = 1
    for xi in X:
        base *= i - xi
    return sum(base / (i - x) * a for x, a in zip(X, A))


while 1:
    d = int(input())
    if d == 0:
        break
    N = d+3
    Y = [float(input()) for i in range(N)]

    cnts = [0]*N
    for X in combinations(range(N), d+1):
        U = [0]*N
        for x in X:
            U[x] = 1
        A = build(X, Y)

        for i in range(N):
            if U[i]:
                continue
            res = calc(X, A, i)

            if abs(Y[i] - res) > 0.5:
                cnts[i] += 1
    print(cnts.index(max(cnts)))

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**13
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    rr = []

    def f(n):
        a = [F() for _ in range(n+3)]
        kd = {}
        for i in range(n+3):
            k = [i**c for c in range(n,-1,-1)]
            kd[(i,)] = [a[i], k]

        for _ in range(n):
            nd = {}
            kl = sorted(list(kd.keys()))
            for i in range(len(kl)):
                k1 = kl[i]
                t1 = kd[k1]
                if t1[-1][-1] == 0:
                    nd[k1] = [t1[0], t1[1][:-1]]
                    continue
                for j in range(i+1,len(kl)):
                    k2 = kl[j]
                    sf = 0
                    for c in k1:
                        if c not in k2:
                            sf += 1
                    if sf != 1:
                        continue
                    t2 = kd[k2]
                    if t2[-1][-1] == 0:
                        continue
                    w = t2[-1][-1] / t1[-1][-1]
                    mt = t2[0] - t1[0] * w
                    ma = [c2 - c1*w for c1,c2 in zip(t1[-1][:-1],t2[-1][:-1])]
                    mk = tuple(sorted(set(k1) | set(k2)))
                    nd[mk] = [mt, ma]
            kd = nd
        msa = inf
        msi = -1
        kl = sorted(list(kd.keys()))
        
            
        for i in range(n+3):
            mi = inf
            ma = -inf
            for k in kl:
                if i in k or abs(kd[k][1][0]) < eps:
                    
                    continue
                t = kd[k][0] / kd[k][1][0]
                
                if mi > t:
                    mi = t
                if ma < t:
                    ma = t
            sa = ma - mi
            
            if msa > sa:
                msa = sa
                msi = i

        return msi


    while 1:
        n = I()
        if n == 0:
            break
        rr.append(f(n))
        

    return '\n'.join(map(str,rr))


print(main())