一本通-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) \]

则
我们可以推出初始矩阵

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

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

\[\begin{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)^{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('\n');
		}
		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\n",(a*f).g[1][1]);
	return 0;
}