题目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必须是完全平方数。

from numba import jit
@jit
def solve(c):
        count = 0
        n = 0
        while count < c:
                n += 1
                if int((5*n*n+2*n+1)**0.5)**2 == 5*n*n+2*n+1:
                        print(n)
                        count += 1
        return n
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...
既然这样，这个问题就又可以转化为：

from functools import lru_cache
@lru_cache(maxsize=None)
def fib(n):
        return n if n < 3 else fib(n-1)+fib(n-2)

print(fib(29)*fib(30))

答案就是：
1120149658760

SixPy

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

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

>>> def fib_gen(n):
    if n<=0: return None
    a,b = 0,1
    if n>=1: yield a
    if n>=2: yield b
    for i in range(n-2):
        a, b = b, a+b
        yield b

        
>>> def fib_gold(n):
        return [x*y for x,y in zip(*([fib_gen(n*2+2)]*2))]

>>> fib_gold(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中给出的这个方程的解的递推公式，我们可以得到这个方程的所有整数解，而问题要求的解就是这些解中的正整数解。

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;

int main(void)
{
  long long x1 = 0, y1 = 1, x2, y2;
  vector<long long> vr;
  for (int i = 1; i <= 10; ++i)
  {
    x2 = -9 * x1 + 4 * y1 - 2;
    y2 = 20 * x1 - 9 * y1 + 4;
    if (x2 >0)
      vr.push_back(x2);
    x1 = x2;
    y1 = y2;
  }
  x1 = 2;
  y1 = 5;
  for (int i = 1; i <= 10; ++i)
  {
    x2 = -9 * x1 + 4 * y1 - 2;
    y2 = 20 * x1 - 9 * y1 + 4;
    if (x2 > 0)
      vr.push_back(x2);
    x1 = x2;
    y1 = y2;
  }
  x1 = 0;
  y1 = -1;
  for (int i = 1; i <= 10; ++i)
  {
    x2 = -9 * x1 + 4 * y1 - 2;
    y2 = 20 * x1 - 9 * y1 + 4;
    if (x2 > 0)
      vr.push_back(x2);
    x1 = x2;
    y1 = y2;
  }
  x1 = -1;
  y1 = 2;
  for (int i = 1; i <= 10; ++i)
  {
    x2 = -9 * x1 + 4 * y1 - 2;
    y2 = 20 * x1 - 9 * y1 + 4;
    if (x2 > 0)
      vr.push_back(x2);
    x1 = x2;
    y1 = y2;
  }
  x1 = -1;
  y1 = -2;
  for (int i = 1; i <= 10; ++i)
  {
    x2 = -9 * x1 + 4 * y1 - 2;
    y2 = 20 * x1 - 9 * y1 + 4;
    if (x2 > 0)
      vr.push_back(x2);
    x1 = x2;
    y1 = y2;
  }
  x1 = 0;
  y1 = 1;
  for (int i = 1; i <= 10; ++i)
  {
    x2 = -9 * x1 - 4 * y1 - 2;
    y2 = -20 * x1 - 9 * y1 - 4;
    if (x2 >0)
      vr.push_back(x2);
    x1 = x2;
    y1 = y2;
  }
  x1 = 2;
  y1 = 5;
  for (int i = 1; i <= 10; ++i)
  {
    x2 = -9 * x1 - 4 * y1 - 2;
    y2 = -20 * x1 - 9 * y1 - 4;
    if (x2 > 0)
      vr.push_back(x2);
    x1 = x2;
    y1 = y2;
  }
  x1 = 0;
  y1 = -1;
  for (int i = 1; i <= 10; ++i)
  {
    x2 = -9 * x1 - 4 * y1 - 2;
    y2 = -20 * x1 - 9 * y1 - 4;
    if (x2 > 0)
      vr.push_back(x2);
    x1 = x2;
    y1 = y2;
  }
  x1 = -1;
  y1 = 2;
  for (int i = 1; i <= 10; ++i)
  {
    x2 = -9 * x1 - 4 * y1 - 2;
    y2 = -20 * x1 - 9 * y1 - 4;
    if (x2 > 0)
      vr.push_back(x2);
    x1 = x2;
    y1 = y2;
  }
  x1 = -1;
  y1 = -2;
  for (int i = 1; i <= 10; ++i)
  {
    x2 = -9 * x1 - 4 * y1 - 2;
    y2 = -20 * x1 - 9 * y1 - 4;
    if (x2 > 0)
      vr.push_back(x2);
    x1 = x2;
    y1 = y2;
  }
  sort(vr.begin(), vr.end());
  unique(vr.begin(), vr.end());
  cout << vr[14] << endl;
  return 0;
}