【HDU5950】Recursive sequence（矩阵快速幂）

BUPT2017 wintertraining(15) #6F

题意

\(f(1)=a,f(2)=b,f(i)=2*(f(i-2)+f(i-1)+i^4)\)
给定n，a，b ，\(N,a,b < 2^{31}\)，求f(n)% 2147493647。

题解

\[f[i]=(f[i-1]+2*f[i-2]+i^4)*2\\ i^4=(i-1)^4+4*(i-1)^3+6*(i-1)^2+4*(i-1)+1 \]

我们可以构造出矩阵乘法

\[\left[ \begin{matrix} f_{i}\\ f_{i-1}\\ i^4\\ i^3\\ i^2\\ i\\ 1\\ \end{matrix} \right] = \left[ \begin{matrix} 1&2&1&4&6&4&1\\ 1&0&0&0&0&0&0\\ 0&0&1&4&6&4&1\\ 0&0&0&1&3&3&1\\ 0&0&0&0&1&2&1\\ 0&0&0&0&0&1&1\\ 0&0&0&0&0&0&1\\ \end{matrix} \right] * \left[ \begin{matrix} f_{i-1}\\ f_{i-2}\\ (i-1)^4\\ (i-1)^3\\ (i-1)^2\\ i-1\\ 1\\ \end{matrix} \right] \]

B为\([f_2,f_1,2^4,2^3,2^2,2,1]^T\)于是\(f(n)=A^{n-2}*B\)的第一项。
有了递推关系，再用矩阵快速幂解决就好了。

代码

#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
#include <iostream>
using namespace std;
const ll mod=2147493647;
struct Mat{
	int r,c;
	ll a[10][10];
	Mat(int _r,int _c){
		r=_r;c=_c;
		memset(a,0,sizeof a);
	}
	Mat operator *(const Mat &b)const{
		Mat c(r,b.c);
		for(int i=0;i<r;i++)
		for(int j=0;j<b.c;j++)
		for(int k=0;k<b.r;k++){
			c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j]%mod)%mod;
		}
		return c;
	}
}A(7,7),B(7,1);

Mat qpow(Mat a,int b){
	Mat c(a.r,a.c);
	for(int i=0;i<a.r;i++)c.a[i][i]=1;
	while(b){
		if(b&1)c=c*a;
		b>>=1;
		a=a*a;
	}
	return c;
}
int main() {
	int at[10][10]={{1,2,1,4,6,4,1},
				  {1,0,0,0,0,0,0},
				  {0,0,1,4,6,4,1},
				  {0,0,0,1,3,3,1},
				  {0,0,0,0,1,2,1},
				  {0,0,0,0,0,1,1},
				  {0,0,0,0,0,0,1}};
	for(int i=0;i<7;i++)for(int j=0;j<7;j++)A.a[i][j]=at[i][j];
	int t,n,a,b;
	cin>>t;
	while(t--){
		scanf("%d%d%d",&n,&a,&b);
		B.a[0][0]=b;B.a[1][0]=a;
		B.a[6][0]=1;
		for(int i=5;i>1;i--)B.a[i][0]=B.a[i+1][0]*2;
		if(n==1){
			printf("%d\n",a);
		}else if(n==2){
			printf("%d\n",b);
		}else{
			Mat C=qpow(A,n-2)*B;
			printf("%lld\n",C.a[0][0]);
		}
	}
	return 0;
}