具有两个相互依赖变量的关系的递归关系解答案

【问题标题】：Recurrence relation solution for relations with two interdependent variables具有两个相互依赖变量的关系的递归关系解
【发布时间】：2019-08-26 05:02:09
【问题描述】：

我在一次在线挑战中遇到了以下问题。

考虑以下向量：

x = [3, 4, ...]
y = [2, 3, ...]

这样对于 i >= 2：

x[i] = x[i-1] + 3 * y[i-2]
y[i] = 2 * y[i-1] + 2 * x[i-2]

什么是 x[10^15]？

虽然这个问题有一个非常简单的解决方案，但问题是 10^15 的值无法在短时间内计算出来。我唯一能想到的是我们必须从递归关系中推导出一个多项式——但这并不容易。我错过了什么吗？

【问题讨论】：

这个：math.stackexchange.com/questions/784710/… 可以直接改编。要快速计算矩阵的 n 次方，请使用二进制取幂：cp-algorithms.com/algebra/binary-exp.html。
@dyukha 我无法形成矩阵。如果我假设基数为 [ [xn, yn-1], [xn-1, yn] ]，乘法数为 [ [1, 2], [3, 2] ]。 yn-1 和 xn-1 在转换后不会保持循环格式。你能告诉我矩阵吗？
如果是斐波那契，你的坐标是 [F[n], F[n-1]]，使用矩阵 [[1 1], [1 0]] 你可以从 [F [n-1]，F[n-2]]。在这种情况下，您的坐标是 [x[n], x[n-1], y[n], y[n-1]]，因此您需要一个 4x4 矩阵。
@dyukha 谢谢。知道了。如果我有时间，我会尝试发布正确的解决方案。

标签： recurrence

【解决方案1】：

问题陈述可以表示为矩阵乘法如下：

A= [
    [1, 0, 0, 3],
    [1, 0, 0, 0],
    [0, 2, 2, 0],
    [0, 0, 1, 0]
]

    [xn+1, xn, yn+1, yn] = A*[xn, xn-1, yn, yn-1]
=>  [xn+1, xn, yn+1, yn] = A^(n-1) * [x1, x0, y1, y0]

[x1, x0, y1, y0] = [4, 3, 3, 2]

虽然问题中没有提到，但由于矩阵乘法超过整数限制，因此需要将解表示为某个素数的余数。让素数为 1000000007。但是我们如何在乘法时不超过整数限制呢？考虑以下几点：

(X * Y) mod p = ((X mod p) * (Y mod p)) mod p

Now, X = A^n
Let, A^n mod p = B
Now, B = B mod p

So,
(X * Y) mod p = 
    ((X mod p) * (Y mod p)) mod p
=>  ((A^n mod p) * (Y mod p)) mod p
=>  ( B * (Y mod p)) mod p
=>  ((B mod p) * (Y mod p)) mod p
=> (B * Y) mod p

所以一个简单的python实现应该是：

import numpy as np

p = 1000000007
A= np.array([
    [1, 0, 0, 3],
    [1, 0, 0, 0],
    [0, 2, 2, 0],
    [0, 0, 1, 0]
])
Y = np.array([4, 3, 3, 2])

# We will use binary exponentiation for fast matrix multiplication
# See: https://cp-algorithms.com/algebra/binary-exp.html
# The `power` list is the array of A's powers needed for that
powers = []
powers.append(A % p)
for i in range(1, 50): # Till 50 since 10^15 ~= 2^50
    Ap = powers[i - 1]
    powers.append(Ap.dot(Ap) % p)

def solve(n):
    pow_of_a = n - 3
    index = 0
    prod = np.identity(4)
    while (pow_of_a > 0):
        if (pow_of_a & 1) == 1:
            prod = prod.dot(powers[index])
        pow_of_a >>= 1  
        index += 1
    B = prod % p
    print(B.dot(Y) % p)

【讨论】：