技術標籤：矩陣乘法

Description

求數列f[n]=f[n-2]+f[n-1]+n+1的第Ｎ項,其中f[1]=1,f[2]:=1.

Input

N(１＜Ｎ＜2^31-1)

Output

第n項結果 mod 9973

Sample Input

10000

Sample Output

4399

思路：
考慮1×4的矩陣【f[n-2],f[n-1],n,1】，希望求得某4×4的矩陣A，使得此1×4的矩陣乘以A得到矩陣：
【f[n-1],f[n],n+1,1】＝【f[n-1],f[n-1]+f[n-2]+n+1,n+1,1】
容易構造出這個4×4的矩陣A，即：

#include <cstring>
#include <cstdio>
#include <iostream>
#define ll long long
using namespace std;
const ll Mod = 9973;
struct node
{
	ll x, y, a[40][40];
} A, B, C;
ll n;
node operator *(node x, node y)  //矩陣乘法
{
	node ans;
	memset(ans.a, 0, sizeof(ans.a));
	ans.x = x.x, ans.y = y.y;
 

	for(ll k = 1; k <= x.y; k++)
	 for(ll i = 1; i <= x.x; i++)
	  for(ll j = 1; j <= y.y; j++)
	  	ans.a[i][j] = (ans.a[i][j] + x.a[i][k] * y.a[k][j] % Mod) % Mod;
	return ans;

}
void ksm(ll k)  //乘法結合律
{
	if(k == 1) {C = B; return ;}
	ksm(k >> 1);
	C = C * C;
	if(k & 1) C = C * B;
}
 

int main()
{
	scanf("%lld", &n);
	if(n < 3) {printf("1\n"); return 0;}
	A.a[1][2] = A.a[1][1] = A.a[1][4] = 1, A.a[1][3] = 3, A.x = 1; A.y = 4;
	
	B.a[1][1] = B.a[1][3] = B.a[1][4] = 0, B.a[1][2] = 1;  //構造單位矩陣
	B.a[2][1] = B.a[2][2] = 1, B.a[2][3] = B.a[2][4] = 0;
	B.a[3][2] = B.a[3][3] = 1, B.a[3][1] = B.a[3][4] = 0;
	B.a[4][4] = B.a[4][2] = B.a[4][3] = 1, B.a[4][1] = 0;
	B.x = B.y = 4;
	ksm(n - 1), C = A * C;
	printf("%lld", C.a[1][1] % Mod);
	return 0;
}