hdu-4686 Jzzhu and Sequences 【矩阵快速幂】

其中a₀ = A0
a_i = a_i-1*AX+AY
b₀ = B0
b_i = b_i-1*BX+BY

最后的结果mod 1,000,000,007

n<=10^18.

ai*bI=Ax*Bx*ai-1*bi-1+Ax*By*ai-1+AyBxbi-1+AyBy

递推式可以构造由矩阵快速幂快速求取。

我理解的递推式构造矩阵的形式可以先这样想（说错了板砖轻拍）:

A*[每个递推式的递推项 ]=[每个递推项前进一项]

【f(n-1),a(n-1),b(n-1),1,s(n-2)】*A=【f(n),a(n),b(n),1,s(n-1)】

A那么由此构造A:

   | axbx    0   0   0   1   |

   | axby   ax   0   0   0  |

   | aybx    0   bx  0   0  |

   | ayay   ay  by   1   0 |

   | 0         0    0    0   1 |

然后 【f(1), a(1) ,b(1), 1, s(0)】*A^(n-1) = 【f(n), a(n), b(n), 1, s(n-1)】

矩阵快速幂求得A^(n-1) 再处理 s(n)=s(n-1)+f(n);

我的代码的矩阵式A*[]=[]形式的。呵呵。A矩阵是我写的矩阵的逆置矩阵。

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <vector>
#define L 5
#define MOD 1000000007
using namespace std;
typedef long long int llint;
typedef vector<llint> vec;
typedef vector<vec> mat;
mat mul(mat a,mat b)
{//矩阵乘法
	mat c(L,vec(L));
	for(int i=0;i<L;++i)
		for(int j=0;j<L;++j)
		{
			c[i][j]=0;
			for(int k=0;k<L;++k)
			{
				c[i][j]+=a[i][k]*b[k][j];
				c[i][j]%=MOD;
			}
		}
	return c;
}
mat quickpow(mat p,llint n)
{
	mat res(L,vec(L));
	//单位矩阵
	for(int i=0;i<L;++i)
		for(int j=0;j<L;++j)
			if(i==j)
				res[i][j]=1;
			else res[i][j]=0;
	while(n>0)
	{
		if(n&1)
			res=mul(p,res);
		n=n>>1;
		p=mul(p,p);
		
	}
	return res;
}
int main()
{
	llint n;
	llint a0,ax,ay;
	llint b0,bx,by;
	while(~scanf("%lld",&n))
	{
		scanf("%lld%lld%lld",&a0,&ax,&ay);
		scanf("%lld%lld%lld",&b0,&bx,&by);
		if(n==0)
			printf("%lld\n",0);
		else if(n==1)
			printf("%lld\n", (a0*b0)%MOD);
		else {
			
			mat p(L,vec(L));
			for(int i=0;i<L;++i)
				for(int j=0;j<L;++j)
					p[i][j]=0;
			//初始A数组
			p[0][0]=(ax*bx)%MOD,p[0][1]=(ax*by)%MOD,
			p[0][2]=(ay*bx)%MOD,p[0][3]=(ay*by)%MOD;
			p[1][1]=ax%MOD,p[1][3]=ay%MOD;
			p[2][2]=bx%MOD,p[2][3]=by%MOD;
			p[3][3]=1;
			p[4][0]=1,p[4][4]=1;
			mat f=quickpow(p,n-1);
			//f(n)
			llint f0=(a0*b0)%MOD;
			llint fn=0,fs=0;
			fn=(fn+f[0][0]*f0)%MOD;
			fn=(fn+f[0][1]*a0)%MOD;
			fn=(fn+f[0][2]*b0)%MOD;
			fn=(fn+f[0][3])%MOD;
			//s(n-1)
			fs=(fs+f[4][0]*f0)%MOD;
			fs=(fs+f[4][1]*a0)%MOD;
			fs=(fs+f[4][2]*b0)%MOD;
			fs=(fs+f[4][3])%MOD;

			printf("%lld\n",(fs+fn)%MOD);
		}
	}
	return 0;
}