【矩陣乘法】【快速冪】遞推

Description

.　　動態規劃的實現形式之一是遞推，因此遞推在oi中十分重要。在某資訊學的分支學科中，LC學會了如何求一階線性遞推數列。由於他現在正在學習主幹學科，因此希望知道求出N階線性遞推數列。為此，他了解到以下內容：
一個N階線性遞推式是這樣的式子：
　 $F_i=a_0\ast F_{i-n}+a_1\ast F+i-(n-1)+...+a_{n-1}\ast F_{i-1}+a_n$
　　也就是說，這個數列的每一項都是由他之前的連續N項加權相加所得。其中還包括一個常數an。
　　例如，當N=2，a0=a1=1，a2=0時，這個式子就是我們熟悉的斐波那契數列。當然，作為邊界條件。f0，f1，…fn-1都是已知的。
　　Lc對如何去求這個式子一籌莫展，因此他想請你幫忙。你的任務是對於一個給定的N階線性遞推式，求出他的第k項。

Input

第一行兩個整數：n，k。其中n表示這是一個N階線性遞推式，k表示你需要球的那一項。
第二行有n+1個整數：a0,a1,…an，表示這個遞推式的係數。
第三行有n個整數：f0,f1,…,fn-1表示數列的初始值。

Output

只有一行，其中只有一個整數，表示這個數列第k項的值。由於資料較大，你只需輸出mod 9973的值。

Sample Input

2 10
1 1 0
0 1

Sample Output

55

---

解題思路

那我真的很失語啊，執行記憶體爆了一直RE，我以為是我陣列開小了: )

首先構造初始矩陣
（為了方便理解所以從$1$開始存，一共$n$個數那麼陣列大小應該是$1\ast (n+1)$）
$A=\left[\begin{array}{ccccc}f[1]& f[2]& ...& f[n]& 1\end{array}\right]$

然後構造一個轉移矩陣
大小設為$(n+1)\ast (n+1)$

1. $F_i=a_0\ast F_{i-n}+a_1\ast F+i-(n-1)+...+a_{n-1}\ast F_{i-1}+a_n$
   算$f[n+1]$時所用的數是$f[1]f[2]...f[n]$
   算$f[n+2]$時所用的數是$f[2]f[3]...f[n+1]$
   所有數向前推一位，第一個數直接推掉不用
   $B=\left[\begin{array}{ccccc}0& 0& 0& & \\ 1& 0& 0& & \\ 0& 1& 0& & \\ 0& 0& ...& & \\ 0& 0& ...& & \end{array}\right]$

2. 最後一個數更新為$f_n$
   $B=\left[\begin{array}{ccccc}0& 0& 0& a_1& \\ 1& 0& 0& a_2& \\ 0& 1& 0& a_3& \\ 0& 0& ...& ...& \\ 0& 0& ...& a_n& \end{array}\right]$
   這個時候$A$矩陣中最後一個$1$的作用就體現出來了，因為要乘上一個$a_n$

3. 最後一列用來傳遞最後的那個$1$
   $B=\left[\begin{array}{ccccc}0& 0& 0& a_1& 0\\ 1& 0& 0& a_2& 0\\ 0& 1& 0& a_3& 0\\ 0& 0& ...& ...& 0\\ 0& 0& ...& a_n& 1\end{array}\right]$

問題轉換成 $ans=A\ast B^k$，然後輸出$ans[1][1]$
但是$A$裡面包含了$f[1]..f[n]$，所以$ans$其實算出了$f[n+1]..f[n+k]$
那麼$ans=A\ast B^{k-n+1}$，輸出$ans[1][n]$

---

Code

#include <iostream>
#include <cstring>
#include <cstdio>
#define ll long long

using namespace std;

const int Mod = 9973;
struct DT{
	int n, m;
	ll aed[20][20];
}A, B, ans;
int n;
ll k;

DT operator *(DT a, DT b){//矩陣乘法
	DT c;
	memset (c.aed, 0, sizeof (c.aed));
	c.n = a.n, c.m = b.m;
	for (int k = 1; k <= a.m; k++)
		for (int i = 1; i <= c.n; i++)
			for (int j = 1; j <= c.m; j++)
				c.aed[i][j] = (c.aed[i][j] + a.aed[i][k] * b.aed[k][j] % Mod) % Mod;
	return c;
}

void power (ll x){
	if (x == 1)
	{
		B = A;
		return;
	}
	power (x / 2);
	B = B * B;
	if (x % 2) B = B * A;
}

int main(){
	scanf ("%d%lld", &n, &k);
	A.n = n + 1, A.m = n + 1;//此處的A不是思路中的A，A和B都是用來算冪的
	for (int i = 1; i <= n + 1; i++)
	{
		scanf ("%lld", &A.aed[i][n]);//a陣列放入轉移矩陣
		A.aed[i][n] %= Mod;
	}
	for (int i = 1; i < n; i++)
		A.aed[i + 1][i] = 1;//轉移f
	A.aed[n + 1][n + 1] = 1;//轉移最後一個1
	power (k - n + 1);//快速冪
	
	ans.n = 1, ans.m = n + 1;//這個ans是思路中的初始矩陣A和答案矩陣ans
	for (int i = 1; i <= n; i++)
	{
		scanf ("%lld", &ans.aed[1][i]);
		ans.aed[1][i] %= Mod; 
	}
	ans.aed[1][n + 1] = 1;
	ans = ans * B;
	printf ("%lld", ans.aed[1][n]);	
}