140. 修正斐波那契金块(Modified Fibonacci golden nuggets)

无穷级数\(A_G(x) = x G_1 + x^2 G_2 + x^3 G_3 + \dots\),其中\(G_k\)是二阶递归关系\(G_k=G_{k-1}+G_{k-2}\)的第\(k\)项,其中\(G_1=1,G_2=4\),该序列为\(1, 4, 5, 9, 14, 23, …\)

在这个问题中,我们感兴趣的是那些使得\(A_G(x)\)为正整数的\(x\),对应于前五个自然数的\(x\)如下所示:

\(x\) \(A_G(x)\)
\(\frac{\sqrt{5}-1}{4}\) \(1\)
\(\tfrac{2}{5}\) \(2\)
\(\frac{\sqrt{22}-2}{6}\) \(3\)
\(\frac{\sqrt{137}-5}{14}\) \(4\)
\(\tfrac{1}{2}\) \(5\)

\(x\)是有理数时,我们称\(A_G(x)\)是一个修正斐波那契金块,因为这样的数将会变得越来越稀少,例如,第20个修正斐波那契金块将是211345365。求前30个修正斐波那契金块的和。

分析:这道题与第一百三十七题非常相似,实际上可以用相同的思路加以解决。首先我们来看一下题目中给出的数列的生成函数,我们把\(A_G(x)\)左右同时乘以\(x\)\(x^2\),并将相同次数的项对齐并依次相减:

$$ \begin{array}{ c r l c c } A_{G}( x) & = & xG_{1} + & x^{2} G_{2} + & x^{3} G_{3} +\cdots \\ xA_{G}( x) & = & & x^{2} G_{1} + & x^{3} G_{2} +\cdots \\ x^{2} A_{G}( x) & = & & & x^{3} G_{1} +\cdots \\ \left( 1-x-x^{2}\right) A_{G}( x) & = & 1x+ & 3x^{2} + & 0x^{3} +\cdots \end{array} $$

则显然有:

$$ A_{G}( x) =\frac{x+3x^2}{1-x-x^{2}} $$

根据上式我们根据给定的\(x\)计算出无穷级数\(A_G(x)\)的和,如当\(x=2/5\)时有\(A_G(x)=2\),正是题目中给出的数值。根据题意,要求\(A_G(x)\)为正整数时的\(x\),且\(x\)须为有理数,则设:

$$ \frac{x+3x^2}{1-x-x^{2}} =n\Rightarrow x+3x^2=n-nx-nx^{2} \Rightarrow (n+3)x^{2} +( n+1) x-n=0 $$

\(n\)为正整数时,解以上一元二次方程,则有:

$$ x=\frac{-( n+1) +\sqrt{( n+1)^{2} +4n(n+3)}}{2(n+3)} $$

要使得\(x\)为有理数,则根号中的数字必须为完全平方数,设其为\(m^2\),则有:

$$ ( n+1)^{2} +4n(n+3) =m^{2} \Rightarrow 5n^{2} -m^{2} +14n+1=0 $$

则我们要求上述二元二次不定方程的正整数解,则经过适当变换可得:

$$ 25n^{2} +70n+49=5m^{2} +44\Rightarrow ( 5n+7)^{2} -5m^{2} =44 $$

如果令\(x=5n+7,y=m\),则上式变为:

$$ x^2-5y^2=44 $$

我们又一次遇到了一个广义佩尔方程,因为我在第一百三十七题中已经详细介绍过如何来求解这样的方程,所以这里我就不再详细讲了。在第一百三十七题中,我们根据基本解写出了三个递推关系,这里我们换一种形式。假设方程\(x^2-dy^2=n\)的一个基本解为\((x_1,y_1)\),其预解式\(u^2-dv^2=1\)的基本解为\((u_1,v_1)\),则我们可以使用以下递推关系得到原方程的所有整数解:

$$ x_{n} =x_{1} u_{n-1} +dy_{1} v_{n-1} ,\ y_{n} =y_{1} u_{n-1} +x_{1} v_{n-1} $$

我们可以把上述递推关系表示成矩阵的形式:

$$ \begin{bmatrix} x_{n}\\ y_{n} \end{bmatrix} =\begin{bmatrix} x_{1} & dy_{1}\\ y_{1} & x_{1} \end{bmatrix}\begin{bmatrix} u_{n-1}\\ v_{n-1} \end{bmatrix} $$

我们知道对于预解式\(u^2-dv^2=1\)也有类似的递推关系,表示成矩阵形式为:

$$ \begin{bmatrix} u_{n}\\ v_{n} \end{bmatrix} =\begin{bmatrix} u_{1} & dv_{1}\\ v_{1} & u_{1} \end{bmatrix}\begin{bmatrix} u_{n-1}\\ v_{n-1} \end{bmatrix} =\begin{bmatrix} u_{1} & dv_{1}\\ v_{1} & u_{1} \end{bmatrix}^{n-1}\begin{bmatrix} u_{1}\\ v_{1} \end{bmatrix} $$

综合上面两个递推式,则得:

$$ \begin{bmatrix} x_{n}\\ y_{n} \end{bmatrix} =\begin{bmatrix} x_{1} & dy_{1}\\ y_{1} & x_{1} \end{bmatrix}\begin{bmatrix} u_{1} & dv_{1}\\ v_{1} & u_{1} \end{bmatrix}^{n-2}\begin{bmatrix} u_{1}\\ v_{1} \end{bmatrix} $$

使用如上矩阵形式可以让递推关系的表达更为简洁紧凑,也有利于我们的编程实现,特别是在基本解比较多时,列出所有递推关系会比较繁琐,使用矩阵形式和相应运算,代码实现就会简单的多。

在具体实现时,我调用了\(sympy\)库中专门用来求解佩尔方程的模块,它可以给出广义佩尔方程的所有基本解。题中佩尔方程有六个基本解是\((-7,1),(13,5),(7,1),(17,7),(-8,2),(8,2)\), 且原方程的预解式\(u^2-5v^2=1\)的基本解为\((9,4)\),则选择其中一个基本解\((13,5)\),有:

$$ \begin{bmatrix} x_{n}\\ y_{n} \end{bmatrix} =\begin{bmatrix} 13 & 25\\ 5 & 13 \end{bmatrix}\begin{bmatrix} 9 & 20\\ 4 & 9 \end{bmatrix}^{n-2}\begin{bmatrix} 9\\ 4 \end{bmatrix} $$

则当\(n=2\)时有:

$$ \begin{bmatrix} x_{2}\\ y_{2} \end{bmatrix} =\begin{bmatrix} 13 & 25\\ 5 & 13 \end{bmatrix}\begin{bmatrix} 9 & 20\\ 4 & 9 \end{bmatrix}^{0}\begin{bmatrix} 9\\ 4 \end{bmatrix} =\begin{bmatrix} 13 & 25\\ 5 & 13 \end{bmatrix}\begin{bmatrix} 9\\ 4 \end{bmatrix} =\begin{bmatrix} 217\\ 97 \end{bmatrix} $$

代换回原方程,则有\(n=(217-7)/5=42\),即为一个符合题目要求的解。使用类似的方式,我们可以找到前三十个符合要求的解。因为我们不知道具体在第多少组解时可以得到三十个解,所以我在实现时使用了\(python\)的生成器,使其可以产生无穷多个广义佩尔方程的解,再从其中筛选出符合要求的解。最后需要注意的是,因为佩尔方程的解实际上关于\(x\)\(y\)轴都是对称的,所在每在一个象限找到解,在剩下三个象限中也可以找到解。当我们要求方程的正整数解时,只需要把任一象限的解取绝对值即可。代码如下:

# time cost = 8.72 ms ± 64.4 µs

import numpy as np
from sympy.solvers.diophantine.diophantine import diop_DN

def pelleq_solution_generator(d,n):
    u,v = diop_DN(d,1)[0]
    uv_zero = np.mat([[u],[v]])
    uv_mat = np.mat([[u,d*v],[v,u]])
    xy_mat = [np.mat([[x,d*y],[y,x]]) for x,y in diop_DN(d,n)]
    new_mat = np.identity(2,dtype=int)
    while True:
        yield [np.abs(x @ new_mat @ uv_zero) for x in xy_mat]
        new_mat = new_mat @ uv_mat

def main(N=30):
    arr = [2]
    for sol in pelleq_solution_generator(5,44):
        for ans in sol:
            if ans[0,0] % 5 == 2:
                arr.append((ans[0,0]-7)//5)
                if len(arr) == N:
                    return sum(arr)