如何求解：f(n) = f(n-1) + 3*f(n-2) + 3*f(n-3) + f(n-4)

Question

我该如何解决这个问题呢？

f(n) = f(n-1) + 3*f(n-2) + 3*f(n-3) + f(n-4) maximum value of n = 10^18 minimum is 1
initial conditions are
f(1) = 1
f(2) = 3
f(3) = 3
f(4) = 1

当 f(n) 可能很大的时候，我需要 print f(n) modulus 10000007。

我对这个问题的尝试是这样的（可能在使用取模时有误）：

1st test case:
3
2
5
9

output:
3
20
695
(working fine)

2nd tst case:
3
1554894
5959595
2562651

output:
7505501
9551828
6592801
(working fine)

但是对于更大的数字，程序就失败了；这是为什么呢？

#!/usr/bin/python
T = int(raw_input())
def fib_iter(n):
    if n == 1:
        return 1
    if n == 2:
        return 3
    if n == 3:
        return 3
    if n == 4:
        return 1   
    prev1, prev2, prev3, prev4 = 1, 3, 3, 1

    i = 5
    while i < n+1:
        curr = (prev1%10000007 + 3*prev2%10000007 + 3*prev3%10000007 + prev4%10000007)%10000007
        prev1, prev2, prev3, prev4 = prev2%10000007, prev3%10000007, prev4%10000007, curr%10000007
        i = i + 1
    return curr%10000007

for t in xrange(T):
    n = int(raw_input())
    print fib_iter(n)

Answer 1

我在这里写了一个更详细的回答。

1. 首先要知道，1000007是一个质数。这意味着你可以使用欧拉定理。欧拉定理的意思是，对于一个质数p，如果你有两个数字a和x、y，那么只要满足x % (p-1) == y % (p-1)，就可以得出a**x % p == a**y % p（这里用的是Python的写法）。

2. 接下来，试着找到一种矩阵表示法，这样你就可以计算f(n)而不需要先递归地计算f(n-1)等等。然后，使用平方取幂的方法来计算它（按照上面提到的模运算）。

Answer 2

更新 - 尝试把循环中的那两行改成：

    curr = (prev1 + 3*prev2 + 3*prev3 + prev4)%10000007
    prev1, prev2, prev3, prev4 = curr, prev1, prev2, prev3

用10000007这个数，n = 15616511848，你应该得到5370032，就像你提到的。注意，10000007 = 941 x 10627，并不是一个质数。我不太确定为什么会选择这个数。通常这类问题会用1000000007（7 + 10^9），这是一个质数。

加快这个算法的速度其实更多是个数学问题，而不是编程问题。它的矩阵形式是：

       |  1  3  3  1 |
   a = |  1  0  0  0 |
       |  0  1  0  0 |
       |  0  0  1  0 |

       |  1 |
x[4] = |  3 |
       |  3 |
       |  1 |

x[i+1] = a x[i]

x[i+j] = a^j x[i]

可以用重复平方的方法来加速计算a^j。

如果你感兴趣的话，x[0]到x[4]的值是不带模的。对于模运算，负数要加上模（10000007）：

x[0] =  -14,  39, -73, 117
x[1] =    1, -14,  39, -73
x[2] =    3,   1, -14,  39
x[3] =    3,   3,   1, -14
x[4] =    1,   3,   3,   1

Answer 3

先看看这本书，然后自己动手解决问题。

补充说明：有人让我提供更多细节：我所说的“解决”是指推导出一个明确的公式，这个公式可以直接给出系数 f(n) 的值，而不需要用递归的方法。对于普通的斐波那契数列，这个公式就是比内公式。按照链接书中的步骤，你可以这样做：

1. 推导出第 n 个系数的生成函数。简单来说，就是考虑一个函数，这个函数的幂级数展开式是你想要的系数 f(n)，然后把几个项合并成一个函数。

2. 有了生成函数后，推导出在单项式 x^n 前面的展开系数。因为你指定的生成函数最多包含四次多项式，所以这一步甚至可以用解析的方法来完成。

如果这里有数学模式的话，我会提供更多的数学细节，运气好的话也能给出答案。不过就这样试试吧，这很简单，而且书也很好读。

Answer 4

使用矩阵方法后，执行速度变得非常快而且可靠。