这是一个非常有趣的序列。它几乎是但不完全是 4 阶斐波那契（又名 Tetranacci）数。从它的伴生矩阵中提取了doubling formulas for Tetranacci 之后，我忍不住为了这个非常相似的递归关系再做一次。
在我们进入实际代码之前，一些定义和所用公式的简短推导是有序的。定义一个整数序列A，这样：
A(n) := A(n-1) + A(n-3) + A(n-4)
初始值为A(0), A(1), A(2), A(3) := 1, 1, 1, 2。
对于n >= 0，这是integer compositions 的数量n 从集合{1, 3, 4} 中的部分。这是我们最终希望计算的序列。
为方便起见，定义一个序列T，这样：
T(n) := T(n-1) + T(n-3) + T(n-4)
初始值为T(0), T(1), T(2), T(3) := 0, 0, 0, 1。
请注意，A(n) 和 T(n) 只是相互转换。更准确地说，A(n) = T(n+3) 表示所有整数 n。因此，正如another answer 所阐述的，两个序列的伴随矩阵是：
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]
[1 1 0 1]
调用这个矩阵C，然后让：
a, b, c, d := T(n), T(n+1), T(n+2), T(n+3)
a', b', c', d' := T(2n), T(2n+1), T(2n+2), T(2n+3)
通过归纳，很容易证明：
[0 1 0 0]^n = [d-c-a c-b b-a a]
[0 0 1 0] [ a d-c c-b b]
[0 0 0 1] [ b b+a d-c c]
[1 1 0 1] [ c c+b b+a d]
如上所示，对于任何n，C^n 都可以仅从其最右侧的列中完全确定。此外，将C^n 与其最右边的列相乘会产生C^(2n) 的最右边的列：
[d-c-a c-b b-a a][a] = [a'] = [a(2d - 2c - a) + b(2c - b)]
[ a d-c c-b b][b] [b'] [ a^2 + c^2 + 2b(d - c)]
[ b b+a d-c c][c] [c'] [ b(2a + b) + c(2d - c)]
[ c c+b b+a d][d] [d'] [ b^2 + d^2 + 2c(a + b)]
因此，如果我们希望通过重复平方计算某些n 的C^n，我们只需要在每一步执行矩阵向量乘法，而不是完整的矩阵矩阵乘法。
现在，用 Python 实现：
# O(n) integer additions or subtractions
def A_linearly(n):
a, b, c, d = 0, 0, 0, 1 # T(0), T(1), T(2), T(3)
if n >= 0:
for _ in range(+n):
a, b, c, d = b, c, d, a + b + d
else: # n < 0
for _ in range(-n):
a, b, c, d = d - c - a, a, b, c
return d # because A(n) = T(n+3)
# O(log n) integer multiplications, additions, subtractions.
def A_by_doubling(n):
n += 3 # because A(n) = T(n+3)
if n >= 0:
a, b, c, d = 0, 0, 0, 1 # T(0), T(1), T(2), T(3)
else: # n < 0
a, b, c, d = 1, 0, 0, 0 # T(-1), T(0), T(1), T(2)
# Unroll the final iteration to avoid computing extraneous values
for i in reversed(range(1, abs(n).bit_length())):
w = a*(2*(d - c) - a) + b*(2*c - b)
x = a*a + c*c + 2*b*(d - c)
y = b*(2*a + b) + c*(2*d - c)
z = b*b + d*d + 2*c*(a + b)
if (n >> i) & 1 == 0:
a, b, c, d = w, x, y, z
else: # (n >> i) & 1 == 1
a, b, c, d = x, y, z, w + x + z
if n & 1 == 0:
return a*(2*(d - c) - a) + b*(2*c - b) # w
else: # n & 1 == 1
return a*a + c*c + 2*b*(d - c) # x
print(all(A_linearly(n) == A_by_doubling(n) for n in range(-1000, 1001)))
因为编码相当简单，所以序列以通常的方式扩展到负n。还提供了一个简单的线性实现作为参考点。
对于足够大的n，通过简单（即不严格，并且可能存在缺陷）时序比较，上述对数实现比直接用numpy 对伴随矩阵求幂快10-20 倍。据我估计，计算A(10**12) 仍需要大约 100 年的时间！尽管上面的算法还有改进的空间，但这个数字实在是太大了。另一方面，为某些M 计算A(10**12) mod M 更容易实现。
与卢卡斯数和斐波那契数直接相关
事实证明，T(n) 更接近斐波那契，Lucas numbers 比它更接近 Tetranacci。要看到这一点，请注意T(n) 的特征多项式是x^4 - x^3 - x - 1 = 0，它会影响(x^2 - x - 1)(x^2 + 1) = 0。第一个因素是斐波那契和卢卡斯的特征多项式！ (x^2 - x - 1)(x^2 + 1) = 0的4个根是两个斐波那契根，phi和psi = 1 - phi，以及i和-i——-1的两个平方根。
T(n) 的封闭式表达式或“Binet”公式将具有一般形式：
T(n) = U(n) + V(n)
U(n) = p*(phi^n) + q*(psi^n)
V(n) = r*(i^n) + s*(-i)^n
对于一些常数系数p, q, r, s。
使用T(n) 的初始值，求解系数，应用一些代数，并注意到卢卡斯数具有闭式表达式：L(n) = phi^n + psi^n，我们可以推导出以下关系：
L(n+1) - L(n) L(n-1) F(n) + F(n-2)
U(n) = ------------- = -------- = ------------
5 5 5
其中L(n) 是L(0), L(1) := 2, 1 的第n 个卢卡斯数，F(n) 是F(0), F(1) := 0, 1 的第n 个斐波那契数。我们还有：
V(n) = 1 / 5 if n = 0 (mod 4)
| -2 / 5 if n = 1 (mod 4)
| -1 / 5 if n = 2 (mod 4)
| 2 / 5 if n = 3 (mod 4)
这很丑陋，但对代码来说微不足道。请注意V(n)can also be succinctly expressed 的分子为cos(n*pi/2) - 2sin(n*pi/2) 或(3-(-1)^n) / 2 * (-1)^(n(n+1)/2)，但为了清楚起见，我们使用分段定义。
这是一个更好、更直接的身份：
T(n) + T(n+2) = F(n)
本质上，我们可以使用斐波那契数和卢卡斯数来计算 T(n)（因此是 A(n)）。从理论上讲，这应该比类似 Tetranacci 的方法更有效。
众所周知，卢卡斯数可以比斐波那契数更有效地计算，因此我们将根据卢卡斯数计算A(n)。我所知道的最有效、最简单的卢卡斯数算法是 L.F. Johnson 的算法（请参阅他的 2010 paper：Middle and Ripple，卢卡斯数的快速简单 O(lg n) 算法）。一旦我们有了 Lucas 算法，我们就使用恒等式：T(n) = L(n - 1) / 5 + V(n) 来计算 A(n)。
# O(log n) integer multiplications, additions, subtractions
def A_by_lucas(n):
n += 3 # because A(n) = T(n+3)
offset = (+1, -2, -1, +2)[n % 4]
L = lf_johnson_2010_middle(n - 1)
return (L + offset) // 5
def lf_johnson_2010_middle(n):
"-> n'th Lucas number. See [L.F. Johnson 2010a]."
#: The following Lucas identities are used:
#:
#: L(2n) = L(n)^2 - 2*(-1)^n
#: L(2n+1) = L(2n+2) - L(2n)
#: L(2n+2) = L(n+1)^2 - 2*(-1)^(n+1)
#:
#: The first and last identities are equivalent.
#: For the unrolled iteration, the following is also used:
#:
#: L(2n+1) = L(n)*L(n+1) - (-1)^n
#:
#: Since this approach uses only square multiplications per loop,
#: It turns out to be slightly faster than standard Lucas doubling,
#: which uses 1 square and 1 regular multiplication.
if n >= 0:
a, b, sign = 2, 1, +1 # L(0), L(1), (-1)^0
else: # n < 0
a, b, sign = -1, 2, -1 # L(-1), L(0), (-1)^(-1)
# unroll the last iteration to avoid computing unnecessary values
for i in reversed(range(1, abs(n).bit_length())):
a = a*a - 2*sign # L(2k)
c = b*b + 2*sign # L(2k+2)
b = c - a # L(2k+1)
sign = +1
if (n >> i) & 1:
a, b = b, c
sign = -1
if n & 1:
return a*b - sign
else:
return a*a - 2*sign
您可以验证 A_by_lucas 产生的结果与之前的 A_by_doubling 函数相同，但速度大约快 5 倍。仍然不够快，无法在任何合理的时间内计算 A(10**12)！