【矩阵乘法】【快速幂】幼儿园数学题II

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)=(n+4)(n-1)/2+f(3)+f(4)+f(5)+...+f(n-3)+f(n-2)$
$f(n)=$ $\frac{n^2}{2}$ +$\frac{3n}{2}-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)=$ $\frac{n^2}{2}+\frac{5n}{2}$$+f(3)+f(4)+f(5)+...f(n-3)$$+f(n-2)+f(n-1)$

$f(n)=>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)=f(n-2)+f(n-1)+n+1$

考虑1×5的矩阵$\left[\begin{array}{ccccc}f[n-2]& f[n-1]& s[n-2]& n& 1\end{array}\right]$,
我们需要找到一个5×5的矩阵A，使得它乘以A得到如下1×5的矩阵：
$\left[\begin{array}{ccccc}f[n-1]& f[n]& s[n-1]& n+1& 1\end{array}\right]=\left[\begin{array}{ccccc}f[n-1]& f[n-1]+f[n-2]+n+1& s[n-2]+f[n-1]& n+1& 1\end{array}\right]$
容易构造出A为：
$\left[\begin{array}{ccccc}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{array}\right]$

---

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]的位置
	}
}