一本通-P1799-数列

题目链接

(f_0)	(f_1)	(f_2)	(f_3)	(f_4)	(f_5)	(f_6)	(f_7)	(f_8)	(f_9)	(f_{10})
0	0	0	2	10	32	84	198	438	932	1936

我们发现

[f_i=2f_{i-1}+(i-1)* (i-2) ]

我们可以设

[g_i=i * (i-1) ]

则

[g_i=i*(i-1)=(i-1) * (i-2)+2(i-1)=g_{i-1}+2(i-1) ]

则
我们可以推出初始矩阵

[egin{bmatrix}& f_0&g_0&0&1 & end{bmatrix}=egin{bmatrix}& 0&0&0&1 & end{bmatrix} ]

第一个数是(f_0)，第二个数是(g_0)，第三个数是i，第四个数是1。
以及转移矩阵

[egin{bmatrix} & 2&0&0&0 & \\ & 1&1&0&0 & \\ & 0&2&1&0 & \\ & 0&0&1&1 & end{bmatrix}]

第一行对应的是(f_0)，第二行对应的是(g_0)，第三行对应的是i，第四行对应的是1。

对于题目给出的(10^{5000})，本蒟蒻因为太弱，不想打高精，想到了使用类似于快速幂的方法，因为$a^n=(a{x})^{10}* a^y (10* x+y=n) $

可以直接避免高精

于是，我们就可以写出代码

#include<cstdio>
#include<cctype>
#include<cstring>
#define re register
#define ll long long
using namespace std;
const int N=1e7+10,mod=1e9+7;
struct matrix
{
	int n,m,g[5][5];
	inline void init(int _n,int _m)
	{
		memset(g,0,sizeof(g));
		n=_n;
		m=_m;
	}
	inline matrix operator *(const matrix &b)
	{
		matrix res;
		res.init(n,b.m);
		for(re int i=1; i<=res.n; i++)
			for(re int j=1; j<=res.m; j++)
				for(re int k=1; k<=b.n; k++)
					res.g[i][j]=(res.g[i][j]+1ll*g[i][k]*b.g[k][j])%mod;
		return res;
	}
	template<typename T>
	inline matrix operator ^ (T b)
	{
		matrix res,a=(*this);
		res.init(4,4);
		for(re int i=1; i<=4; i++)
			res.g[i][i]=1;
		while(b)
		{
			if(b&1)
				res=res*a;
			a=a*a;
			b>>=1;
		}
		return res;
	}
	inline void print()
	{
		for(re int i=1; i<=n; i++)
		{
			for(re int j=1; j<=m; j++)
				printf("%d ",g[i][j]);
			putchar('
');
		}
		return;
	}
} a,b,f;
int n;
char s[5010];
int main()
{
	a.init(1,4);
	a.g[1][4]=1;
	b.init(4,4);
	b.g[1][1]=2;
	b.g[2][1]=b.g[2][2]=1;
	b.g[3][2]=2,b.g[3][3]=1;
	b.g[4][3]=b.g[4][4]=1;
	f.init(4,4);
	for(re int i=1; i<=4; i++)
		f.g[i][i]=1;
	scanf("%s",s+1);
	for(re int i=1; s[i]; i++)
		f=(f^10)*(b^((int)s[i]-'0'));
	printf("%d
",(a*f).g[1][1]);
	return 0;
}