"""
If we are presented with the first k terms of a sequence it is impossible to
say with certainty the value of the next term, as there are infinitely many
polynomial functions that can model the sequence.

As an example, let us consider the sequence of cube numbers. This is defined by
the generating function,
    u_n = n3: 1, 8, 27, 64, 125, 216, ...

Suppose we were only given the first two terms of this sequence. Working on the
principle that "simple is best" we should assume a linear relationship and
predict the next term to be 15 (common difference 7). Even if we were presented
with the first three terms, by the same principle of simplicity, a quadratic
relationship should be assumed.

We shall define OP(k, n) to be the nth term of the optimum polynomial
generating function for the first k terms of a sequence. It should be clear
that OP(k, n) will accurately generate the terms of the sequence for n ≤ k, and
potentially the first incorrect term (FIT) will be OP(k, k+1); in which case we
shall call it a bad OP (BOP).

As a basis, if we were only given the first term of sequence, it would be most
sensible to assume constancy; that is, for n ≥ 2, OP(1, n) = u1.

Hence we obtain the following OPs for the cubic sequence:

    OP(1, n) = 1	        1, 1, 1, 1, ...
    OP(2, n) = 7n−6	        1, 8, 15, ...
    OP(3, n) = 6n2−11n+6  	1, 8, 27, 58, ...
    OP(4, n) = n3	        1, 8, 27, 64, 125, ...

Clearly no BOPs exist for k ≥ 4.

By considering the sum of FITs generated by the BOPs (indicated in red above),
we obtain 1 + 15 + 58 = 74.

Consider the following tenth degree polynomial generating function:

un = 1 − n + n2 − n3 + n4 − n5 + n6 − n7 + n8 − n9 + n10

Find the sum of FITs for the BOPs.

Solution comment: Very simple. Just a straight forward implementation
                  of the stuff in the problem statement. Very quick.
"""

import numpy as np
from itertools import count


def OP(f, k):
    n = list(range(1, k + 1))
    u = [f(n_i) for n_i in n]
    g = np.poly1d(np.polyfit(n, u, k - 1))
    return g


def FIT(guess, correct, rtol=1e-5):
    for i, (g, c) in enumerate(zip(guess, correct)):
        if abs((g - c) / g) > rtol:
            return g
    return None


def generator(n):
    return sum((-1) ** p * n ** p for p in range(11))


fit_sum = 0
for k in count(1):
    op = OP(generator, k)
    fit = FIT(*zip(*[(round(op(n)), generator(n)) for n in range(1, k + 2)]))
    if fit is None:
        break
    fit_sum += fit

print('Answer: ', int(fit_sum))
