数学知识——矩阵乘法

矩阵乘法

引入

由于线性递推式可以表示成矩阵乘法的形式，也通常用矩阵快速幂来求线性递推数列的某一项。
利用结合律，矩阵乘法可以利用 快速幂 的思想来优化。

例题

斐波那契前 n 项和

大家都知道 Fibonacci 数列吧，f1=1,f2=1,f3=2,f4=3,…,fn=fn−1+fn−2。

现在问题很简单，输入 n 和 m，求 fn 的前 n 项和 Snmodm。

输入格式
共一行，包含两个整数 n 和 m。

输出格式
输出前 n 项和 Snmodm 的值。

数据范围
1≤n≤2000000000
,
1≤m≤1000000010
输入样例：
5 1000
输出样例：
12

思路

已知递推公式$S_{n+1}=S_n+f_{n+1}$,对于求$S_{n+1}$我们要知道$S_n和f_{n+1}$。
构造矩阵$F_n=[f_n,f_{n+1},S_n]$。
对于$F_{n+1}=[f_{n+1},f_{n+2},S_{n+1}]$等于$F_n$与矩阵$A=[[0,1,0],[1,1,1],[0,0,1]]$相乘，
对于$F_n=F_1\ast A^{n-1}$

代码

from copy import deepcopy
N = 3
F1 = [[0] * N for _ in range(N)]
F1[0] = [1, 1, 1]
A = [[0, 1, 0], [1, 1, 1], [0, 0, 1]]

def mul(a, b) :
	tmp = [[0] * N for _ in range(N)]
	for i in range(N) : # i行
		for j in range(N) : # j列
			for k in range(N) :
				tmp[i][j] = (tmp[i][j] + a[i][k] * b[k][j]) % m
	return tmp

def qmi(a, k) :
	res = deepcopy(F1)
	while k :
		if k & 1 :
			res = mul(res, a)
		k >>= 1
		a = mul(a, a)
	return res
n, m = map(int, input().split())
res = qmi(A, n - 1)
print(res[0][2])

佳佳的斐波那契

佳佳对数学，尤其对数列十分感兴趣。

在研究完 Fibonacci 数列后，他创造出许多稀奇古怪的数列。

例如用 S(n) 表示 Fibonacci 前 n项和 modm
的值，即 S(n)=(F1+F2+…+Fn)modm，其中 F1=F2=1,Fi=Fi−1+Fi−2。

可这对佳佳来说还是小菜一碟。

终于，她找到了一个自己解决不了的问题。

用 T(n)=(F1+2F2+3F3+…+nFn)modm 表示 Fibonacci 数列前 n 项变形后的和 modm 的值。

现在佳佳告诉你了一个 n 和 m，请求出 T(n) 的值。

输入格式
共一行，包含两个整数 n 和 m。

输出格式
共一行，输出 T(n) 的值。

数据范围
1≤n,m≤2^31−1
输入样例：
5 5
输出样例：
1
样例解释
T(5)=(1+2×1+3×2+4×3+5×5)mod5=1

思路

$T_{n+1}=T_n+(n+1)\ast f_{n+1}$，矩阵乘法要使用快速幂推导，等式必须是线性的，也就是保证转移矩阵中只有常数。
对于求取$T_n$，必须构造一个线性的递推公式(必须包含$T_n$)的函数。
对于$P_n=n{S}_n-T_n$有这样的递推公式$P_{n+1}=P_n+S_n$
构造$F_n=[f_n,f_{n+1},S_n,P_n]$，$F_{n+1}=[f_{n+1},f_{n+2},S_{n+1},P_{n+1}]$,
转移矩阵为$A=[[0,1,0,0],[1,1,1,0],[0,0,1,1],[0,0,0,1]]$
$F_n=F_1\ast A^{n-1}$

代码

from copy import deepcopy
N = 4
F1 = [[0] * N for _ in range(N)]
F1[0] = [1, 1, 1, 0]
A = [[0, 1, 0, 0], 
    [1, 1, 1, 0], 
    [0, 0, 1, 1],
    [0, 0, 0, 1]]
def mul(a, b) :
	tmp = [[0] * N for _ in range(N)]
	for i in range(N) :
		for j in range(N) :
			for k in range(N) :
				tmp[i][j] = (tmp[i][j] + a[i][k] * b[k][j]) % m
	return tmp

def qmi(a, k) :
	res = deepcopy(F1)
	while k :
		if k & 1 :
			res = mul(res, a)
		k >>= 1
		a = mul(a, a)
	return res

n, m = map(int, input().split())
res = qmi(A, n - 1)
print((n * res[0][2] - res[0][3]) % m)