Project Euler 101

题目

Optimum polynomial

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 = n^3: 1, 8, 27, 64, 125, 216, \dots\)

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 \(n^{\text{th}}\) 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 \leq 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 \ge 2\), \(OP(1, n) = u_1\). Hence we obtain the following OPs for the cubic sequence:

\(OP(1, n) = 1\)	\(1, \mathbf{1}, 1, 1, \dots\)
\(OP(2, n) = 7n−6\)	\(1, 8, \mathbf{15}, \dots\)
\(OP(3, n) = 6n^2−11n+6\)	\(1, 8, 27,\mathbf{58}, \dots\)
\(OP(4, n) = n^3\)	\(1, 8, 27, 64, 125, \dots\)

Clearly no BOPs exist for \(k \ge 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:

\[u_n = 1 − n + n^2 − n^3 + n^4 − n^5 + n^6 − n^7 + n^8 − n^9 + n^{10}\]

Find the sum of FITs for the BOPs.

解决方案

对于一元函数上任意\(n\)个不同的点，可以用唯一一个\(n-1\)次多项式确定。求出这个\(n-1\)次的多项式\(p(x)\)最简单的办法是使用待定系数法，设\(p(x)=c_0+c_1x+c_2x^2+\dots+c_{n-1}x^{n-1}\)。

设已知的\(n\)个点的坐标为\((x_1,y_1),(x_2,y_2),\dots,(x_n,y_n)\)。由于已经知道的是数列的前\(n\)项，因此\(x_k=k\)。每一个\(n-1\)次多项式都有\(n\)个系数，所以代入\((1,u_1),(2,u_2),\dots,(n,u_n)\)，列出一个有\(n\)条方程的\(n\)元一次方程组：

\[ \begin{bmatrix} 1^0 & 1^1 & 1^2 &\cdots & 1^{n-1} \\ 2^0 & 2^1 & 2^2 & \cdots & 2^{n-1} \\ 3^0 & 3^1 & 3^2 & \cdots & 3^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ n^0 & n^1 & n^2 &\cdots & n^{n-1} \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{n-1} \end{bmatrix} = \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_n \end{bmatrix} \]

使用高斯消元方法，解得系数\([c_0,c_1,c_2,\dots,c_{n-1}]^T\)后，即可确定\(p(x)\)本身，那么，第\(n\)个多项式对应的FIT为\(p(n+1)\)。当\(p\)的次数等于多项式\(u\)的次数时，就不存在FIT了。

因此将各个FIT相加即可。

另外一种做法则是使用多项式插值。本代码使用的是sympy库的interpolating_poly函数，它是基于拉格朗日插值实现的。给定\(n\)个点\((x_1,y_1),(x_2,y_2),\dots,(x_n,y_n)\)，那么这个函数就会返回一个\(n-1\)次多项式。拉格朗日插值的公式如下：

\[p(x)=\sum_{j=1}^n y_j\prod_{k=1,k\neq j}^n \dfrac{x-x_k}{x_j-x_k}\]

代码

from sympy import interpolating_poly
from sympy.abc import x

N = 10
x_list = list(range(1, N + 1))
y_list = list(
    map(lambda n: 1 - n + n ** 2 - n ** 3 + n ** 4 - n ** 5 + n ** 6 - n ** 7 + n ** 8 - n ** 9 + n ** 10, x_list))
ans = 0
for i in range(1, N + 1):
    poly = interpolating_poly(i, x, x_list, y_list)
    if poly.is_Number:
        ans += poly
    else:
        ans += poly.as_poly()(i + 1)
print(ans)

from fractions import Fraction


def fun(n: int):
    return 1 - n + n ** 2 - n ** 3 + n ** 4 - n ** 5 + n ** 6 - n ** 7 + n ** 8 - n ** 9 + n ** 10


def solve(mat: list):
    n = len(mat)
    for j in range(n):
        for i in range(j, n):
            if mat[i][j]:
                mat[i], mat[j] = mat[j], mat[i]
                break
        if mat[j][j] == 0:
            continue
        for i in range(n):
            if i == j or mat[i][j] == 0:
                continue
            val = -mat[i][j] / mat[j][j]
            for k in range(j, n + 1):
                mat[i][k] += mat[j][k] * val
    b = []
    for j in range(n):
        b.append(mat[j][n] / mat[j][j])
    return b


ans = 0
for m in range(1, 11):
    mat = []
    for i in range(1, m + 1):
        ls = []
        for j in range(m):
            ls.append(Fraction(i ** j))
        ls.append(Fraction(fun(i)))
        mat.append(ls)
    b = [int(x) for x in solve(mat)]
    val = 0
    for i in range(m):
        val += b[i] * ((m + 1) ** i)
    ans += val
print(ans)