理论部分
二次剩余
在数论中，整数 $X$ 对整数 $p$ 的二次剩余是指 $X^2$ 除以 $p$ 的余数。
当存在某个 $X$，使得式子 $X^2 equiv d(mod p)$ 成立时，称“ $d$ 是模 $p$ 的二次剩余”
当对任意 $X$，$X^2 equiv d(mod p)$ 都不成立时，称“ $d$ 是模 $p$ 的二次非剩余”
矩阵的相似对角化
相似矩阵：对于矩阵 $A$ 和 $B$，若存在可逆矩阵 $P$，使得 $P^{-1}AP=B$，则称 $A$ 相似于 $B$，记作 $Asim B$。
矩阵可对角化是指与对角矩阵相似，即 $P^{-1}AP = diag(lambda _1, lambda _2,...,lambda _n)$
定理：$n$ 阶矩阵 $A$ 与对角矩阵相似的充要条件是 $A$ 有 $n$ 个线性无关的特征向量。简单的判断就是 $|A| eq 0$.
分析
来源于 acdreamers 的广义Fibonacci数列找循环节，这里只是梳理一下。
将递推式表示成转移矩阵的形式，相当于只需找到转移矩阵的周期。
即求最小的 $n$，使得
$${egin{bmatrix} a & b\ 1 & 0 end{bmatrix}}^n(mod p) = egin{bmatrix} 1 & 0\ 0 & 1 end{bmatrix}$$
我们可以找到一个满足条件的 $n$，然后枚举其因子。
设 $A$ 相似于 $D$，$A$ 的特征根分别为 $lambda _1,lambda _2$.
则D为
$D = egin{bmatrix} lambda _1 & 0\ 0 & lambda _2 end{bmatrix}$
因为 $A^n = T^{-1}D^nT$，所以 $D^n equiv I(mod p)$ 时，$D^n equiv I(mod p)$
由于$$D^n = {egin{bmatrix} lambda _1 & 0\ 0 & lambda _2 end{bmatrix}}^n = egin{bmatrix} lambda _1^n & 0\ 0 & lambda _2^n end{bmatrix}(mod p) = I$$
于是 $lambda _1^nequiv 1(mod p), lambda _2^nequiv 1(mod p)$.
如何找到 $n$ 呢？
直接上结论吧！（看不懂了555）
设 $Delta = a^2+4b$
- 若 $Delta$ 是 $p$ 的二次剩余，$n=p-1$
- 若 $Delta$ 是 $p$ 的二次非剩余，$n=(p+1)(p-1)$
- 若 $Delta=p$，则无法相似化。幸好题目不会出现这种情况
上述结论还是针对 $p$ 为素数，$p$ 为合数呢？
首先，质因数分解 $p=p_1^{a_1}p_2^{a_2}...p_k^{a_k}$,
设 $g(x)$ 为模 $x$ 时循环节长度，有结论 $g(p_i^{a_i}) = p_i^{a_1-1} g(p_i)$
所以
$$
egin{aligned}
g(x) & =g(p_1^{a_1}p_2^{a_2}...p_k^{a_k})\
&= p_1^{a_1-1} g(p_1) cdot p_2^{a_2-1} g(p_2) ... p_k^{a_k-1} g(p_k)\
&= (p / p_1p_2... p_k) g(p_1)g(p_2) ... g(p_k)
end{aligned}$$
实现
二次剩余我又不会判断，直接当作 $(p+1)(p-1)$，大不了多个倍数。
按道理最小的循环节不是 $n$，而是 $n$ 的因数，感觉分解再判断挺麻烦的，不管了，大不了多个倍数。
#include<bits/stdc++.h> using namespace std; typedef long long ll; const int N=1000000+10; ll x0,x1,a,b,n,mod; char s[N]; const int maxn = 100000 + 10; //p最为2e9，不会有两个超过1e5的质因数 int prime[maxn], pcnt; //prime[i]表示第i个素数 bool is_prime[maxn + 1]; //is_prime[i]为true表示i是素数 int sieve(int n) { int cnt = 0; for (int i = 0; i <= n; i++) is_prime[i] = true; is_prime[0] = is_prime[1] = false; for (ll i = 2; i <= n; i++) { if (is_prime[i]) { prime[cnt++] = i; for (ll j = i * i; j <= n; j += i) is_prime[j] = false; //i * i可能爆int } } return cnt; } ll solve(ll x){ ll ans1=1,ans2=1,xx=x; for(int i=0;i<pcnt;i++){ if(1ll*prime[i]*prime[i]>x) break; if(x%prime[i]==0){ ans1*=(prime[i]-1)*(prime[i]+1); ans2*=prime[i]; while(x%prime[i]==0) x/=prime[i]; } } if(x>1){ ans1*=(x-1)*(x+1); ans2*=x; } return xx/ans2*ans1; } ll qmul(ll x,ll y,ll p){ //快速乘 x%=p; y%=p; ll ans=0; while(y){ if(y&1){ ans+=x; if(ans>=p) ans-=p; //这样写不能有负数 } x<<=1; if(x>=p) x-=p; y>>=1; } return ans; } struct Mat{ int r,c; ll m[3][3]; Mat(){ memset(m,0,sizeof(m)); } }; Mat mmul(Mat x,Mat y,ll p){ Mat ans; ans.r=x.r; ans.c=y.c; for(int i=0;i<x.r;i++) for(int k=0;k<x.c;k++) for(int j=0;j<y.c;j++){ ans.m[i][j]+=qmul(x.m[i][k],y.m[k][j],p); if(ans.m[i][j]>=p) ans.m[i][j]-=p; } return ans; } Mat mpow(Mat x,ll y,ll p){ Mat ans; ans.r=x.r; ans.c=x.c; for(int i=0;i<ans.c;i++) ans.m[i][i]=1; while(y){ if(y&1) ans=mmul(ans,x,p); x=mmul(x,x,p); y>>=1; } return ans; } int main(){ pcnt = sieve(100000); while(scanf("%lld%lld%lld%lld",&x0,&x1,&a,&b) == 4){ scanf("%s%lld",s,&mod); ll lop=solve(mod); //循环节长度 n=0; int lens=strlen(s); for(int i=0;i<lens;i++){ n=qmul(n,10,lop)+s[i]-'0'; if(n>=lop) n-=lop; } Mat A,T; A.r=2; A.c=1; A.m[0][0]=x1; A.m[1][0]=x0; T.r=2; T.c=2; T.m[0][0]=a; T.m[0][1]=b; T.m[1][0]=1; if(n>1){ T=mpow(T,n-1,mod); A=mmul(T,A,mod); } printf("%lld ",A.m[0][0]); } return 0; }
参考链接：
1. 分析： https://blog.csdn.net/acdreamers/article/details/25616461