P5175 题解

题意简述

给出数列 $a_n(1\le n\le {10}^{18})$ 的两项 $a_1,a_2$ 与递推公式 $a_n=x{a}_{n-1}+y{a}_{n-2}$，求：

$$S_n=\sum _{k=1}^n{a}_k^2\mathrm{mod}({10}^9+7)$$

题目分析

一看见 $1\le n\le {10}^{18}$，就大概能知道要用 $O(\log n)$ 级别的算法。再一看递推，就知道要用矩阵快速幂了。

题目所求的答案是前 $n$ 项的 $a_k^2$ 之和，如果我们能快速求出每个 $a_k^2$，那么再用矩阵快速幂递推求和是十分简单的。

注意到 $a_{k+1}=x{a}_k+y{a}_{k-1}$，那么就有 $a_{k+1}^2=x^2{a}_k^2+y^2{a}_{k-1}^2+2xy{a}_k{a}_{k-1}$，如果我们维护了 $a_k^2$ 与 $a_{k-1}^2$，就可以递推地求出 $a_{k+1}^2,a_{k+2}^2$ 等等。

还没完，注意到 $a_{k+1}$ 展开后还有一项 $2xy{a}_k{a}_{k-1}$，去掉系数就还剩下 $a_k{a}_{k-1}$。维护它应该怎么处理呢？将 $a_k$ 拆开得到 $a_k{a}_{k-1}=(x{a}_{k-1}+y{a}_{k-2})a_{k-1}=x{a}_{k-1}^2+y{a}_{k-1}{a}_{k-2}$，其中 $a_{k-1}^2$ 已经维护了，所剩的只有 $y{a}_{k-1}{a}_{k-2}$。这就意味着我们维护每一个 $a_k{a}_{k-1}$ 就可以求出 $a_{k+1}{a}_k$ 和 $a_{k+2}{a}_{k+1}$ 等等。

综上，我们使用答案矩阵
$\left[\begin{array}{cccc}{S}_k& a_k^2& a_{k-1}^2& a_k{a}_{k-1}\end{array}\right]$ 进行维护。初始值为
$\left[\begin{array}{cccc}{S}_2=a_1^2+a_2^2& a_2^2& a_1^2& a_2{a}_1\end{array}\right]$。

为了将其递推到 $\left[\begin{array}{cccc}{S}_{k+1}& a_{k+1}^2& a_k^2& a_{k+1}{a}_k\end{array}\right]$，需要有一个转移矩阵：

$$\left[\begin{array}{cccc}1& 0& 0& 0\\ x^2& x^2& 1& x\\ y^2& y^2& 0& 0\\ 2xy& 2xy& 0& y\end{array}\right]$$

我们一列一列地进行解释。

对于第 $2$ 列，有 $a_{k+1}^2=0S_k+x^2{a}_k^2+y^2{a}_{k-1}^2+2xy{a}_k{a}_{k-1}$。

那么对于第 $1$ 列，有 $S_{k+1}=S_k+a_{k+1}^2=1S_k+x^2{a}_k^2+y^2{a}_{k-1}^2+2xy{a}_k{a}_{k-1}$。

对于第 $3$ 列，有 $a_k^2=0S_k+1a_k^2+0a_{k-1}^2+0a_k{a}_{k-1}$，没什么可解释的。

对于第 $4$ 列，有 $a_{k+1}{a}_k=0S_k+x{a}_k^2+0a_{k-1}^2+y{a}_k{a}_{k-1}$。

现在我们只需要计算出转移矩阵的 $n-2$ 次幂，再乘到答案矩阵上就好了。

需要注意的是要特判 $n=1$，否则会挂掉。然后就是该开 long long 的地方一定要开，不然也会挂。

分析一下时间复杂度。一次矩阵乘法是 $O(4^3)=O(64)$，因此矩阵快速幂的时间复杂度为 $O(64\log n)$，总复杂度为 $O(64T\log n)$，在 $T\le 3\times {10}^4,1\le n\le {10}^{18}$，时限 $1.50s$ 的情况下用个快读快写就可以过了。

代码实现

#include<bits/stdc++.h>
using namespace std;
const long long mod=1000000007;
int t;
long long n,x,y,a1,a2;
void rd(long long &x)
{
	x=0ll;
	char c=getchar();
	for(;c>'9'||c<'0';c=getchar());
	for(;c<='9'&&c>='0';c=getchar())
		x=(x<<3ll)+(x<<1ll)+c-'0';
}//long long快读 
void rd(int &x)
{
	x=0;
	char c=getchar();
	for(;c>'9'||c<'0';c=getchar());
	for(;c<='9'&&c>='0';c=getchar())
		x=(x<<3)+(x<<1)+c-'0';
}//int快读 
void pt(long long x)
{
	if(x>=10ll)
		pt(x/10ll);
	putchar(x%10ll+'0');
}//快写 
struct matrix
{
	int n,m;
	long long a[5][5];
	void init(long long k)
	{
		for(int i=1;i<=n;i++)
			for(int j=1;j<=m;j++)
				if(i!=j)
					a[i][j]=0ll;
				else
					a[i][j]=k;
	}//初始化 k=0零矩阵，k=1单位矩阵 
	friend matrix operator *(matrix a,matrix b)
	{
		matrix c;
		c.n=a.n,c.m=a.m;
		c.init(0ll);
		for(int i=1;i<=a.n;i++)
			for(int k=1;k<=a.m;k++)
				for(int j=1;j<=b.m;j++)
					c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%mod)%mod;
		return c;
	}//矩阵乘法 
	friend matrix operator ^(matrix a,long long b)
	{
		matrix res;
		res.n=res.m=a.n;
		res.init(1ll);
		for(;b;b>>=1ll)
		{
			if(b&1ll)
				res=res*a;
			a=a*a;
		}
		return res;
	}//矩阵快速幂 
}ans,trans;
int main()
{
	rd(t);
	while(t--)
	{
		rd(n),rd(a1),rd(a2),rd(x),rd(y);
		if(n==1ll)//特判 n=1，不然会 T 得很惨
		{
			pt(a1*a1%mod);
			putchar('\n');
			continue;
		}
		ans.n=1,ans.m=4;
		ans.a[1][1]=(a2*a2%mod+a1*a1%mod)%mod,ans.a[1][2]=a2*a2%mod,ans.a[1][3]=a1*a1%mod,ans.a[1][4]=a1*a2%mod;//初始化答案矩阵 
		trans.n=trans.m=4;
		trans.init(0ll);
		trans.a[2][1]=trans.a[2][2]=x*x%mod,trans.a[3][1]=trans.a[3][2]=y*y%mod,trans.a[4][1]=trans.a[4][2]=2ll*x*y%mod;
		trans.a[1][1]=trans.a[2][3]=1ll;
		trans.a[2][4]=x,trans.a[4][4]=y;//初始化转移矩阵 
		trans=trans^(n-2ll);//转移矩阵快速幂 
		ans=ans*trans;//答案矩阵乘上转移矩阵 
		pt(ans.a[1][1]);
		putchar('\n');
	}
	return 0;
}