题目137：考察系数为斐波那契数的无限多项式序列的值

欧拉计划

Fibonacci golden nuggets

Consider the infinite polynomial series AF(x) = xF₁ + x²F₂ + x³F₃ + ..., where Fk is the kth term in the Fibonacci sequence: 1, 1, 2, 3, 5, 8, ... ; that is, F_k = F_k−1 + F_k−2, F₁ = 1 and F₂ = 1.

For this problem we shall be interested in values of x for which AF(x) is a positive integer.

Surprisingly A_F(1/2)         =         (1/2).1 + (1/2)².1 + (1/2)³.2 + (1/2)⁴.3 + (1/2)⁵.5 + ...
                                         =         1/2 + 1/4 + 2/8 + 3/16 + 5/32 + ...
                                         =         2
The corresponding values of x for the first five natural numbers are shown below.

We shall call A_F(x) a golden nugget if x is rational, because they become increasingly rarer; for example, the 10th golden nugget is 74049690.

Find the 15th golden nugget.

题目：

考虑无限多项式序列 A_F(x) = xF₁ + x²F₂ + x³F₃ + ...，其中 F_k 是斐波那契数列中的第 k 项。1，1，2，3，5，8，...；也就是说 F_k = F_k−1 + F_k−2 ，F₁ = 1，F₂ = 1.

在本题中我们考虑使得 A_F(x) 为正整数的 x 值。

有趣的是，A_F(1/2)        =        (1/2).1 + (1/2)².1 + (1/2)³.2 + (1/2)⁴.3 + (1/2)⁵.5 + ...
                                    =        1/2 + 1/4 + 2/8 + 3/16 + 5/32 + ...
                                    =        2

前五个自然数对应的 x 值如下表所示：

如果 x 是有理数，我们称 A_F(x) 为一个金块，因为这样的数非常之少；例如，第 10 个金块是 74049690。

找出第 15 个金块。

jerryxjr1220

这题一开始拿到手感觉是比较难的，没方向，然后去维基百科查了Af（x）的公式：x / (1-x-x**2)，有了公式就好办了。假设 x / (1-x-x**2) = n，nx**2+(n+1)x-n=0，根据b^2-4ac = 5n^2+2n+1必须是完全平方数。

1. from numba import jit
   
2. @jit
   
3. def solve(c):
   
4.         count = 0
   
5.         n = 0
   
6.         while count < c:
   
7.                 n += 1
   
8.                 if int((5*n*n+2*n+1)**0.5)**2 == 5*n*n+2*n+1:
   
9.                         print(n)
   
10.                         count += 1
    
11.         return n
    
12. solve(10)

复制代码

2
15
104
714
4895
33552
229970
1576239
10803704
74049690
前10个都是正确的，其实可以一直算到第15个，但是由于大数开方的误差问题，导致开方后的数不正确。

我们通过观察前10个数，可以发现，其实他们就是Fib数列的连续项的乘积，例如，1*2=2, 3*5=15, 8*13=104...
既然这样，这个问题就又可以转化为：

1. from functools import lru_cache
   
2. @lru_cache(maxsize=None)
   
3. def fib(n):
   
4.         return n if n < 3 else fib(n-1)+fib(n-2)
   
5. 
   
6. print(fib(29)*fib(30))

复制代码

答案就是：
1120149658760

SixPy

jerryxjr1220 发表于 2017-7-20 10:26
这题一开始拿到手感觉是比较难的，没方向，然后去维基百科查了Af（x）的公式：x / (1-x-x**2)，有了公式就 ...

如果你的推理是正确的话
计算就很简单了

1. >>> def fib_gen(n):
   
2.     if n<=0: return None
   
3.     a,b = 0,1
   
4.     if n>=1: yield a
   
5.     if n>=2: yield b
   
6.     for i in range(n-2):
   
7.         a, b = b, a+b
   
8.         yield b
   
9. 
   
10.         
    
11. >>> def fib_gold(n):
    
12.         return [x*y for x,y in zip(*([fib_gen(n*2+2)]*2))]
    
13. 
    
14. >>> fib_gold(15)
    
15. [0, 2, 15, 104, 714, 4895, 33552, 229970, 1576239, 10803704, 74049690, 507544127, 3478759200, 23843770274, 163427632719, 1120149658760]

复制代码

guosl

没有通过严格数学证明的解法为令x=F(i)/F(i+1)。代入函数
f(x)中得f(F(i)/F(i+1)=F(i)*F(i+1)/(F(i-1)*F(i+1)-F(i)^2)
再根据斐波那契数列的特性知：F(i-1)*F(i+1)-F(i)^2=(-1)^i
所以f(F(i)/F(i+1))=(-1)^i F(i)*F(i+1)，由于函数值是正的故i只能取偶数。
所以取i=2得第一个黄金块2，取i=4得第二个黄金块15
...
取i=20得第10个黄金块74049690。
所以我们猜想第i个黄金块的值是F(2*i)*F(2*i+1)
结论是对的，但我无法证明方程的有理数解就只有这种形式！
答案：1120149658760

guosl

另一种解法：
假设 x / (1-x-x**2) = n，nx**2+(n+1)x-n=0，根据b^2-4ac = 5n^2+2n+1必须是完全平方数。令：m^2=5n^2+2n+1，则得二元二次不定方程：5x^2-y^2+2x+1=0，而根据
https://www.alpertron.com.ar/QUAD.HTM中给出的这个方程的解的递推公式，我们可以得到这个方程的所有整数解，而问题要求的解就是这些解中的正整数解。

1. #include <iostream>
   
2. #include <vector>
   
3. #include <algorithm>
   
4. using namespace std;
   
5. 
   
6. int main(void)
   
7. {
   
8.   long long x1 = 0, y1 = 1, x2, y2;
   
9.   vector<long long> vr;
   
10.   for (int i = 1; i <= 10; ++i)
    
11.   {
    
12.     x2 = -9 * x1 + 4 * y1 - 2;
    
13.     y2 = 20 * x1 - 9 * y1 + 4;
    
14.     if (x2 >0)
    
15.       vr.push_back(x2);
    
16.     x1 = x2;
    
17.     y1 = y2;
    
18.   }
    
19.   x1 = 2;
    
20.   y1 = 5;
    
21.   for (int i = 1; i <= 10; ++i)
    
22.   {
    
23.     x2 = -9 * x1 + 4 * y1 - 2;
    
24.     y2 = 20 * x1 - 9 * y1 + 4;
    
25.     if (x2 > 0)
    
26.       vr.push_back(x2);
    
27.     x1 = x2;
    
28.     y1 = y2;
    
29.   }
    
30.   x1 = 0;
    
31.   y1 = -1;
    
32.   for (int i = 1; i <= 10; ++i)
    
33.   {
    
34.     x2 = -9 * x1 + 4 * y1 - 2;
    
35.     y2 = 20 * x1 - 9 * y1 + 4;
    
36.     if (x2 > 0)
    
37.       vr.push_back(x2);
    
38.     x1 = x2;
    
39.     y1 = y2;
    
40.   }
    
41.   x1 = -1;
    
42.   y1 = 2;
    
43.   for (int i = 1; i <= 10; ++i)
    
44.   {
    
45.     x2 = -9 * x1 + 4 * y1 - 2;
    
46.     y2 = 20 * x1 - 9 * y1 + 4;
    
47.     if (x2 > 0)
    
48.       vr.push_back(x2);
    
49.     x1 = x2;
    
50.     y1 = y2;
    
51.   }
    
52.   x1 = -1;
    
53.   y1 = -2;
    
54.   for (int i = 1; i <= 10; ++i)
    
55.   {
    
56.     x2 = -9 * x1 + 4 * y1 - 2;
    
57.     y2 = 20 * x1 - 9 * y1 + 4;
    
58.     if (x2 > 0)
    
59.       vr.push_back(x2);
    
60.     x1 = x2;
    
61.     y1 = y2;
    
62.   }
    
63.   x1 = 0;
    
64.   y1 = 1;
    
65.   for (int i = 1; i <= 10; ++i)
    
66.   {
    
67.     x2 = -9 * x1 - 4 * y1 - 2;
    
68.     y2 = -20 * x1 - 9 * y1 - 4;
    
69.     if (x2 >0)
    
70.       vr.push_back(x2);
    
71.     x1 = x2;
    
72.     y1 = y2;
    
73.   }
    
74.   x1 = 2;
    
75.   y1 = 5;
    
76.   for (int i = 1; i <= 10; ++i)
    
77.   {
    
78.     x2 = -9 * x1 - 4 * y1 - 2;
    
79.     y2 = -20 * x1 - 9 * y1 - 4;
    
80.     if (x2 > 0)
    
81.       vr.push_back(x2);
    
82.     x1 = x2;
    
83.     y1 = y2;
    
84.   }
    
85.   x1 = 0;
    
86.   y1 = -1;
    
87.   for (int i = 1; i <= 10; ++i)
    
88.   {
    
89.     x2 = -9 * x1 - 4 * y1 - 2;
    
90.     y2 = -20 * x1 - 9 * y1 - 4;
    
91.     if (x2 > 0)
    
92.       vr.push_back(x2);
    
93.     x1 = x2;
    
94.     y1 = y2;
    
95.   }
    
96.   x1 = -1;
    
97.   y1 = 2;
    
98.   for (int i = 1; i <= 10; ++i)
    
99.   {
    
100.     x2 = -9 * x1 - 4 * y1 - 2;
     
101.     y2 = -20 * x1 - 9 * y1 - 4;
     
102.     if (x2 > 0)
     
103.       vr.push_back(x2);
     
104.     x1 = x2;
     
105.     y1 = y2;
     
106.   }
     
107.   x1 = -1;
     
108.   y1 = -2;
     
109.   for (int i = 1; i <= 10; ++i)
     
110.   {
     
111.     x2 = -9 * x1 - 4 * y1 - 2;
     
112.     y2 = -20 * x1 - 9 * y1 - 4;
     
113.     if (x2 > 0)
     
114.       vr.push_back(x2);
     
115.     x1 = x2;
     
116.     y1 = y2;
     
117.   }
     
118.   sort(vr.begin(), vr.end());
     
119.   unique(vr.begin(), vr.end());
     
120.   cout << vr[14] << endl;
     
121.   return 0;
     
122. }

复制代码