使用矩阵找出将n写为1、3和4之和的不同方法的数量？

Question

这是本演示文稿中给出的问题。Dynamic Programming

现在我已经使用递归实现了这个算法，它对小值很有效。但是当n大于30时，它就变成了真的，演示文稿提到，对于较大的n值，人们应该考虑一些类似于斐波那契数的东西，我很难理解如何使用斐波那契数的矩阵形式来得出solution.Can，有人给了我一些提示或伪代码

谢谢

Answer 1

这是一个非常有趣的序列。它几乎但不完全是-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要容易得多。

与Lucas和Fibonacci数有直接关系

事实证明，T(n)比Lucas numbers更接近斐波那契和Tetranacci。要了解这一点，请注意T(n)的特征多项式是x^4 - x^3 - x - 1 = 0，它是(x^2 - x - 1)(x^2 + 1) = 0的因数。第一个因子是Fibonacci & Lucas的特征多项式！(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)的初始值，求解系数，应用一些代数，并注意到Lucas数具有封闭形式的表达式：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，fast simple O(lg ) algorithms for Lucas number )。一旦我们有了Lucas算法，我们就使用identity：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)！

Answer 2

通过添加memoization，您可以轻松地改进当前的递归实现，从而使解决方案再次变快。C#代码：

// Dictionary to store computed values
private static Dictionary<int, long> s_Solutions = new Dictionary<int, long>();

private static long Count134(int value) {
  if (value == 0)
    return 1;
  else if (value <= 0)
    return 0;

  long result;

  // Improvement: Do we have the value computed? 
  if (s_Solutions.TryGetValue(value, out result))
    return result;

  result = Count134(value - 4) + 
           Count134(value - 3) + 
           Count134(value - 1);

  // Improvement: Store the value computed for future use
  s_Solutions.Add(value, result);

  return result;
}

所以你可以很容易的调用

Console.Write(Count134(500));

结果(大约需要2毫秒)是

3350159379832610737