T446080 百战名 讲解
zrl123456
·
·
题解
题目
递推题。
题目简述：给你 a,b,c,p,n，使 x_1,x_2 为 ax^2+bx+c=0 的两根，求 (x_1^n+x_2^n)\bmod p。
众所周知，一元二次方程的公式是 \displaystyle x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}，但是最后的小数取模就很难弄，我们应该换个思路。
我们往递推的角度想，设 f_i=(x_1^i+x_2^i)\bmod p，则：
\begin{aligned}f_i&=(x_1^i+x_2^i)\bmod p\\&=(x_1^i+x_2^i+x_1^{i-1}x_2+x_1x_2^{i-1}-x_1^{i-1}x_2-x_1x_2^{i-1})\bmod p\\&=((x_1+x_2)(x_1^{i-1}+x_2^{i-1})-x_1x_2(x_1^{i-2}+x_2^{i-2}))\bmod p\\&=((x_1+x_2)f(i-1)-x_1x_2f(i-2))\bmod p\end{aligned}
这样我们就推出了状态转移方程，那么（\Delta=b^2-4ac）：
\begin{aligned}x_1+x_2&=\frac{-b+\sqrt\Delta}{2a}+\frac{-b-\sqrt\Delta}{2a}\\&=\frac{-b+\sqrt\Delta-b-\sqrt\Delta}{2a}\\&=\frac{-2b}{2a}\\&=\frac{-b}{a}\end{aligned}
我们用 exgcd（费马小定理也可以）求逆，exgcd 详见 P1495 【模板】中国剩余定理（CRT）/ 曹冲养猪 讲解。
\begin{aligned}x_1x_2&=(\frac{-b+\sqrt\Delta}{2a})(\frac{-b-\sqrt\Delta}{2a})\\&=\frac{(-b+\sqrt\Delta)(-b-\sqrt\Delta)}{4a^2}\\&=\frac{b^2-b\sqrt\Delta+b\sqrt\Delta-\Delta}{4a^2}\\&=\frac{b^2-\Delta}{4a^2}\\&=\frac{b^2-(b^2-4ac)}{4a^2}\\&=\frac{4ac}{4a^2}\\&=\frac{c}{a}\end{aligned}
同上用 exgcd 求 \displaystyle\frac{c}{a}。
这样我们成功得到了矩阵 res：
\begin{bmatrix}\frac{-b}{a}&-\frac{c}{a}\\1&0\end{bmatrix}
然后我们得到 res^{n-2}，再来一个矩阵 ans：
\begin{bmatrix}f_2\\f_1\end{bmatrix}
$$\begin{aligned}f_2&=x_1^2+x_2^2\\&=x_1^2+x_2^2+2x_1x_2-2x_1x_2\\&=(x_1+x_2)^2-2x_1x_2\\&=(\frac{-b}{a})^2-\frac{2c}{a}\end{aligned}$$
这样 $res^{n-2}\times ans$ 的 $a_{1,1}$ 就是答案。
代码：
```cpp
#include<bits/stdc++.h>
#define int long long
#define inl inline
#define INF LLONG_MAX
#define rep(i,x,y) for(int i=x;i<=y;++i)
#define rer(i,x,y,cmp) for(int i=x;i<=y&&cmp;++i)
#define per(i,x,y) for(int i=x;i>=y;--i)
#define FAST ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define endl '\n'
using namespace std;
int a,b,c,p,n,x,y,sum,ts,numm;
struct matrix{
int a[5][5];
inl matrix init(){
matrix res;
memset(res.a,0,sizeof(res.a));
return res;
}
inl matrix init2(){
matrix res=init();
res.a[1][1]=res.a[2][2]=1;
return res;
}
inl friend matrix operator*(matrix a,matrix b){
matrix res=res.init();
rep(i,1,2)
rep(j,1,2)
rep(k,1,2) (res.a[i][j]+=a.a[i][k]*b.a[k][j])%=p;
return res;
}
inl friend matrix operator^(matrix a,int b){
matrix num=a,ans=ans.init2();
while(b){
if(b&1) ans=ans*num;
num=num*num;
b>>=1;
}
return ans;
}
}ans,num;
inl void exgcd(int a,int b,int p){
if(p==0){
x=a;
y=0;
return;
}
exgcd(a,p,b%p);
int tx=x;
x=y;
y=tx-b/p*y;
}
signed main(){
FAST;
cin>>a>>b>>c>>p>>n;
exgcd(-b,a,p);
sum=((x%p)+p)%p;
exgcd(c,a,p);
ts=((x%p)+p)%p;
numm=(sum*sum%p-(ts<<1)%p+p)%p;
ans=num=ans.init();
num.a[1][1]=numm;
num.a[2][1]=sum;
ans.a[1][1]=sum;
ans.a[1][2]=((-ts)%p+p)%p;
ans.a[2][1]=1;
if(n==1){
cout<<sum;
return 0;
}
if(n==2){
cout<<numm;
return 0;
}
ans=ans^(n-2);
ans=ans*num;
cout<<ans.a[1][1];
return 0;
}
```