技術標籤：矩陣乘法快速冪

Description

這天,當一頭霧水的LZH同學在考場上痛哭的時候,一旁的YMW早就如切菜一樣cut掉了簡單至極的第一題,風輕雲淡的衝擊著滿分,然而最後一道題著實難道了他,畢竟是幼兒園副園長樹皮和著名毒瘤秋彪為了防止人AK而出的,可是YMW作為ACrush的著名粉絲,向來以AK為目標,永不言敗,而他能不能AK就看你了
題目是醬紫的,f(n)-f(3)-f(4)-f(5)-…-f(n-3)-f(n-2)=(n+4)(n-1)/2,f(1)=1,f(2)=1
求f(n)的前n項和

Input

輸入 一個正整數n（保證0<=n<=2^31-1）

Output

輸出 一個正整數,表示這個圖形的整點個數,需要對1000000007求餘

Sample Input

樣例輸入1

1

樣例輸入2

2

Sample Output

樣例輸出1

1

樣例輸出2

2

解題思路

f ( n ) − f ( 3 ) − f ( 4 ) − f ( 5 ) − . . . − f ( n − 3 ) − f ( n − 2 ) = ( n + 4 ) ( n − 1 ) / 2 f(n)-f(3)-f(4)-f(5)-...-f(n-3)-f(n-2)=(n+4)(n-1)/2 f(n)−f(3)−f(4)−f(5)−...−f(n−3)−f(n−2)=(n+4)(n−1)/2

f ( n ) = ( n + 4 ) ( n − 1 ) / 2 + f ( 3 ) + f ( 4 ) + f ( 5 ) + . . . + f ( n − 3 ) + f ( n − 2 ) f(n) = (n+4)(n-1)/2 +f(3)+f(4)+f(5)+...+f(n-3)+f(n-2)

f ( n + 1 ) = ( n + 5 ) ( n ) / 2 + f ( 4 ) + f ( 5 ) + . . . + f ( n − 3 ) + f ( n − 2 ) + f ( n − 1 ) f(n+1) = (n+5)(n)/2 +f(4)+f(5)+...+f(n-3)+f(n-2)+f(n-1)

f ( n ) = > f ( n + 1 ) f(n)=>f(n+1) f(n)=>f(n+1)
f ( n + 1 ) = f ( n ) + 2 n 2 + 2 + f ( n − 1 ) = f ( n ) + ( n + 1 ) + 1 + f ( n − 1 ) f(n+1) = f(n) + \frac{2n}{2} +2+ f(n-1)=f(n) + (n+1)+1+f(n - 1) f(n+1)=f(n)+22n​+2+f(n−1)=f(n)+(n+1)+1+f(n−1)
f ( n ) = f ( n − 2 ) + f ( n − 1 ) + n + 1 f(n) = f(n-2) + f(n - 1) + n + 1 f(n)=f(n−2)+f(n−1)+n+1

考慮1×5的矩陣 [ f [ n − 2 ] f [ n − 1 ] s [ n − 2 ] n 1 ] \begin{bmatrix} f[n-2]&f[n-1]&s[n-2]&n&1 \end{bmatrix} [f[n−2]​f[n−1]​s[n−2]​n​1​],
我們需要找到一個5×5的矩陣A，使得它乘以A得到如下1×5的矩陣：
[ f [ n − 1 ] f [ n ] s [ n − 1 ] n + 1 1 ] = [ f [ n − 1 ] f [ n − 1 ] + f [ n − 2 ] + n + 1 s [ n − 2 ] + f [ n − 1 ] n + 1 1 ] \begin{bmatrix} f[n-1]&f[n]&s[n-1]&n+1&1 \end{bmatrix}=\begin{bmatrix}f[n-1]&f[n-1]+f[n-2]+n+1&s[n-2]+f[n-1]&n+1&1\end{bmatrix} [f[n−1]​f[n]​s[n−1]​n+1​1​]=[f[n−1]​f[n−1]+f[n−2]+n+1​s[n−2]+f[n−1]​n+1​1​]
容易構造出A為：
[ 0 1 0 0 0 1 1 1 0 0 0 0 1 0 0 0 1 0 1 0 0 1 0 1 1 ] \begin{bmatrix} 0& 1& 0& 0& 0\\ 1& 1& 1& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 1& 0& 1& 0\\ 0& 1& 0& 1& 1 \end{bmatrix} ⎣⎢⎢⎢⎢⎡​01000​11011​01100​00011​00001​⎦⎥⎥⎥⎥⎤​

Code

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

using namespace std;

const ll Mod = 1000000007;
struct DT{
	int n, m;
	ll aed[10][10];
}A, B, Bc;
int n;

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 (int x){//快速冪
	if (x == 1)
	{
		Bc = B;
		return; 
	}
	power (x / 2);
	Bc = Bc * Bc;
	if (x % 2)
		Bc = Bc * B;
}

void init (){//初始化
	A.n = 1, A.m = 5;
	A.aed[1][1] = A.aed[1][2] = A.aed[1][3] = A.aed[1][5] = 1, A.aed[1][4] = 3;//從f(3)開始初始化
	B.n = 5, B.m = 5;
	B.aed[1][2] = B.aed[2][1] = B.aed[2][2] = B.aed[2][3] = 1;
	B.aed[3][3] = B.aed[4][2] = B.aed[4][4] = 1;
	B.aed[5][2] = B.aed[5][4] = B.aed[5][5] = 1;
}

int main(){
	init();
	scanf ("%d", &n);
	if (n == 1)
		printf ("1");
	else
	{
		power (n - 1);
		A = A * Bc;
		printf("%lld", A.aed[1][3]);//A.aed[1][3]是s[n-2]的位置
	}
}