Project Euler 356

题目

Largest roots of cubic polynomials

Let \(a_n\) be the largest real root of a polynomial \(g(x) = x^3 - 2^n·x^2 + n\).

For example, \(a_2 = 3.86619826\dots\)

Find the last eight digits of \(\sum \limits_{i = 1}^{30} {\left \lfloor a_i^{987654321} \right \rfloor}\).

Note: \(\lfloor a \rfloor\) represents the floor function.

解决方案

牛顿恒等式

牛顿恒等式描述了\(k\)次多项式方程\(F(x)=0\)的根的同次幂之和与\(F(x)\)之间的关系。

设 \(m\) 次多项式\(F(x)=x^m+c_1x^{m-1}+c_2x^{m-2}+\cdots+c_m,\)其根为 \(x_1,\dots,x_m\)。定义幂和 \(\displaystyle{P_k=\sum_{i=1}^m x_i^k (k\ge 0), P_0=m.}\)

牛顿恒等式给出 \(P_k\) 与系数的关系：

• 若 \(1\le k\le m\)，则\(P_k+c_1P_{k-1}+c_2P_{k-2}+\cdots+c_{k-1}P_1+k,c_k=0.\)
• 若 \(k>m\)，则\(P_k+c_1P_{k-1}+c_2P_{k-2}+\cdots+c_mP_{k-m}=0.\)

令\(N=30\)，枚举几个函数点，有如下发现：

\(\begin{aligned} &g(-1)=n-1-2^n<0\\ &g(0)=n>0\\ &g(1)=n+1-2^n\le 0 \end{aligned}\)

因此\(g(n)\)的前两个根分别在区间\((-1,0),(0,1]\)中。

最后令 \(t=2^n\)，计算

\[ \begin{aligned} g(t-1) &=(t-1)^3-t(t-1)^2+n \\ &=(t^3-3t^2+3t-1)-(t^3-2t^2+t)+n\\ &=-t^2+2t+(n-1), \end{aligned} \]

当 \(n\ge 2\) 时显然 \(g(t-1)<0\)；而\(g(t)=t^3-t\cdot t^2+n=n>0.\)

因此还存在一根\(r_1\in(2^n-1,2^n).\)

这三根落在互不相交的区间中，所以 \(g_n\) 有三个互异实根，且最大实根就是\(a_n=r_1.\)

本题\(g(x)=x^3-2^n x^2+0\cdot x+n,\)因此 \(m=3\)，且\(c_1=-2^n, c_2=0, c_3=n.\)

令三根为 \(r_1,r_2,r_3\)，并定义\(p_k=r_1^k+r_2^k+r_3^k (k\ge 0), p_0=3.\)

由牛顿恒等式：

• \(k=1\)：\(p_1+c_1=0 \Rightarrow p_1=2^n.\)

• \(k=2\)：\(p_2+c_1p_1+2c_2=0 \Rightarrow p_2=2^n p_1=2^{2n}.\)

• 对 \(k\ge 3\)（对三次式的“长递推”）：\(p_k+c_1p_{k-1}+c_2p_{k-2}+c_3p_{k-3}=0,\) 即\(p_k=2^np_{k-1}-np_{k-3} (k\ge 3).\)

由于初值 \(p_0,p_1,p_2\) 都是整数且递推系数是整数，所以对所有 \(k\) 都有\(p_k\in\mathbb{Z}.\)

已知 \(a_n=r_1\) 为最大实根，而另两根满足\(r_2\in(0,1], r_3\in(-1,0).\)再用韦达定理（或直接由系数）：\(r_1+r_2+r_3=2^n.\)

因为 \(r_1<2^n\)，所以\(r_2+r_3=2^n-r_1>0,\)从而\(r_2>-r_3 \Rightarrow r_2>|r_3|.\)

由于 \(E\) 是奇数，\(r_3^E<0\)，且 \(|r_3|<r_2\le 1\)，于是\(0<r_2^E-|r_3|^E=r_2^E+r_3^E<r_2^E<1.\)也就是说\(0<r_2^E+r_3^E<1.\)

而\(p_E=r_1^E+r_2^E+r_3^E,\)因此\(p_E-1<r_1^E<p_E,\)立刻得到\(\left\lfloor r_1^E\right\rfloor=p_E-1.\)

所以本题目标化为 \[\sum_{n=1}^{N}\left\lfloor a_n^E\right\rfloor=\sum_{n=1}^{N}(p_E-1)=\left(\sum_{n=1}^{N}p_E\right)-N.\]

递推\(p_k=2^np_{k-1}-np_{k-3}\)是 3 阶线性递推。设状态向量（对固定 \(n\)）\(\mathbf{v}_k= \begin{pmatrix} p_k\\ p_{k-1}\\ p_{k-2} \end{pmatrix}.\)则\(\mathbf{v}_{k+1}=A_n,\mathbf{v}_k,\)其中 \(A_n= \begin{pmatrix} 2^n & 0 & -n\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{pmatrix}.\)

初始向量取 \(k=2\)：\(\mathbf{v}_2=\begin{pmatrix}p_2\\p_1\\p_0\end{pmatrix}=\begin{pmatrix}2^{2n}\\2^n\\3\end{pmatrix}.\)

于是\(\mathbf{v}_E=A_n^{,E-2}\mathbf{v}_2,\)使用矩阵快速幂即可完成本题的计算。

代码

N = 30
MOD = 10 ** 8
E = 987654321

def mat_mul(A, B, mod):
    return [[sum(A[i][k]*B[k][j] for k in range(len(A[0]))) % mod for j in range(len(B[0]))] for i in range(len(A))]


def mat_pow(M, e, mod):
    # 3x3 matrix power
    R = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    while e:
        if e & 1:
            R = mat_mul(R, M, mod)
        M = mat_mul(M, M, mod)
        e >>= 1
    return R


total = 0
for n in range(1, 31):
    a = pow(2, n, MOD)  # 2^n mod MOD
    p0 = 3 % MOD
    p1 = a
    p2 = (a * a) % MOD  # 2^(2n) mod MOD

    # A_n = [[2^n, 0, -n], [1, 0, 0], [0, 1, 0]] mod MOD
    A = [
        [a, 0, (-n) % MOD],
        [1, 0, 0],
        [0, 1, 0],
    ]

    P = mat_pow(A, E - 2, MOD)
    pE = (P[0][0] * p2 + P[0][1] * p1 + P[0][2] * p0) % MOD

    total = (total + (pE - 1)) % MOD

print(f"{total:08d}")