洛谷 P5175 数列 做题记录

纯纯数学题。

看到 \(n\le 10^{18}\) 不难想到矩乘，但是 \(\log_2 10^{18} \approx 60\)，再加上 \(T=30000\) 的多测，运算量已经来到了 \(1.8 \times 10^6\)，所以我们最多有一个边长为 \(\sqrt[3]{\frac{1.5\times 10^8}{6 \times 10^6}} \approx 4\) 的矩阵。

\[\because a_i=xa_{i-1}+ya_{i-2} \]

\[\therefore {a_i}^2=(xa_{i-1}+ya_{i-2})^2=x^2{a_{i-1}}^2+y^2{x_{i-2}}^2+2xya_{i-1}a_{i-2} \]

\[\therefore a_{i}a_{i-1}=a_{i-1}(xa_{i-1}+ya_{i-2})=x{a_{i-1}}^2+ya_{i-1}a_{i-2} \]

至此，我们想要得到 \({a_i}^2\) 的所有东西都有了，直接用矩乘加速就好了。

矩阵长这样：

\[\begin{bmatrix} {a_{i-1}}^2 & {a_{i-2}}^2 & a_{i-1}a_{i-2} & sum_i-1 \end{bmatrix} \times \begin{bmatrix} x^2 & 1 & x & x^2\\ y^2 & 0 & 0 & y^2\\ 2xy & 0 & y & 2xy \\ 0 & 0 & 0 & 1 \end{bmatrix} \]

点击查看代码

struct node {
	int x[5][5];
}qwq,awa;
const int mod=1e9+7;
node mul(node a,node b) {
	node res; m0(res.x);
	For(i,1,4) For(j,1,4) For(k,1,4) (res.x[i][j]+=a.x[i][k]*b.x[k][j]%mod)%=mod;
	return res;
}
node qpo(node a,int b) {
	node res;m0(res.x);
	For(i,1,4) res.x[i][i]=1;
	for(;b;b>>=1,a=mul(a,a)) if(b&1) res=mul(res,a);
	return res;
}
void work() {
	int n,a1,a2;
	in3(n,a1,a2);
	int x,y;
	in2(x,y);
	if(n==1) {printf("%lld\n",a1*a1%mod); return ;}
	if(n==2) {printf("%lld\n",(a1*a1%mod+a2*a2%mod)%mod); return ;}
	m0(qwq.x);m0(awa.x);
	qwq.x[1][1]=x*x%mod,qwq.x[2][1]=y*y%mod,qwq.x[3][1]=2*x*y%mod;
	qwq.x[1][2]=1,qwq.x[1][3]=x,qwq.x[3][3]=y;
	qwq.x[1][4]=x*x%mod,qwq.x[2][4]=y*y%mod,qwq.x[3][4]=2*x*y%mod,qwq.x[4][4]=1;
	awa.x[1][1]=a2*a2%mod,awa.x[1][2]=a1*a1%mod,awa.x[1][3]=a1*a2%mod,awa.x[1][4]=(a1*a1%mod+a2*a2%mod)%mod;
	awa=mul(awa,qpo(qwq,n-2));
	printf("%lld\n",awa.x[1][4]);
}