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\)，并将相同次数的项对齐并依次相减：

则显然有：

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

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

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

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

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

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

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

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

综合上面两个递推式，则得：

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

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

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

代换回原方程，则有\(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)