本文章由 WyOJ Shojo 从洛谷专栏拉取，原发布时间为 2024-04-26 16:39:03

题目

递推题。
题目简述：给你 $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}\
&0\end{bmatrix}$$ 然后我们得到 $res^{n-2}$，再来一个矩阵 $ans$： $$\begin{bmatrix}f_2\f_1\end{bmatrix}$$ $\displaystyle f_1=\frac{-b}{a}$，我们再求 $f_2$： $$\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}$ 就是答案。

代码：

#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;
}