2019牛客多校第五场 generator 1——广义斐波那契循环节&&矩阵快速幂

理论部分

二次剩余

在数论中,整数 $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;
}
View Code

参考链接:

1. 分析: https://blog.csdn.net/acdreamers/article/details/25616461

2. 代码:https://www.cnblogs.com/LMCC1108/p/11286388.html