2019牛客暑期多校训练营（第五场）B.generator 1(广义斐波那契数列找循环节)

题意

已知$x_0,x_1,a,b(1\le 10^9)$，求
$$x_n=a{x}_{n-1}+b{x}_{n-2}(modm)$$
$1\le n\le 10^{10^6},10^9<m\le 2\times 10^9$

思路

由于题目中给的$n$特别大，所以直接做矩阵快速幂不行，当然可以用10进制快速幂做。这里介绍一下找循环节的办法。
显然的在模意义下广义斐波那契必有循环节，那么循环节和什么有关呢，对于广义斐波那契数列

$$\{\begin{array}{c}f[1]=c\\ f[2]=d\\ f[n]=af[n-1]+bf[n-2],n>2\end{array}$$
我们可以尝试固定$a,b$改变$c,d$，会发现数列的循环节其实没有变化的，要说证明的话我们可以构造一下矩阵
$$\left[\begin{array}{c}{f}_n\\ f_{n-1}\end{array}\right]=\left[\begin{array}{cc}a& b\\ 1& 0\end{array}\right]\left[\begin{array}{c}{f}_{n-1}\\ f_{n-2}\end{array}\right](modm)$$
而循环节就是在模意义下${\left[\begin{array}{cc}a& b\\ 1& 0\end{array}\right]}^k=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right](modm)$,这里的$k$就是循环节，所以实际上只会和$a,b$有关和$c,d$是无关的

接下来有这么一个定理若模数为$m=p_1^{a_1}{p}_2^{a_2}\cdots p_k^{a_k}$，$g(m)$表示模$m$下循环节，那么就有
$$g(m)=lcm(g(p_1^{a_1}),g(p_2^{a_2}),\cdots ,g(p_k^{a_k}))$$
$lcm$为最小公倍数
并且对于$g(p^k)$，$p$为质数,有$g(p^k)=g(p)\ast p^{k-1}$
那么我们现在就可以将给的模数$m$唯一分解了然后对于其每一个质因数求解循环节，现在考虑对于质因数$p$怎么求解循环节

我们可以对其的通项下手进行分析，对于这种2阶线性常系数齐次递推关系，我们可以用特征方程的方式来求解其的通项即
$$x^2-ax-b=0$$
考虑用求根公式来求解$\Delta =a^2+4b$，考虑$\Delta$在$(modp)$的情况$p>2$

1. $\Delta \equiv 0(modp)$，即有重根的情况，就有$f[n]=A^n(Bn+C)$
   根据求根公式$$A=\frac{a+\sqrt{\Delta }}{2}=\frac{a}{2}$$
   有$A=a\ast inv2(modp)$那我们就可以用费马小定理来得到循环节了
   $$A^{p-1}\equiv A^0\equiv 1(modp)$$
   又有$p\equiv 0(modp)$
   所以循环节必是$lcm(p-1,p)=p(p-1)$的倍数
2. $\Delta ̸\equiv 0(modp)$，即有两个不同的根的情况有$f[n]=A{B}^n+C{D}^n$
   根据求根公式我们求的是$\sqrt{\Delta }$的所以我们得考虑$\Delta$是否是模$p$意义下的二次剩余
   若$\Delta$为二次剩余即$x^2\equiv \Delta (modp)$，$x$有解，那么表示在模$p$意义下$\sqrt{\Delta }$可以表示，那么
   $$B=\frac{a+\sqrt{\Delta }}{2},D=\frac{a+\sqrt{\Delta }}{2}$$在模意义下都是整数，根据费马小定理或者是扩展欧拉降幂，$B^n,C^n$在模$p$的情况下循环节就是$p-1$的倍数
   若$\Delta$不为二次剩余即$x^2\equiv \Delta (modp)$，$x$无解，那么我们可以考虑将$B^{p+1}$展开
   
   $$B^{p+1}=(\frac{a+\sqrt{\Delta }}{2}{)}^{p+1}=\frac{1}{2^{p+1}}\sum _{i=0}^{p+1}{C}_{p+1}^i{a}^{p+1-i}{\sqrt{\Delta }}^i$$
   $$\equiv \frac{1}{4}({\sqrt{\Delta }}^{p+1}+a{\sqrt{\Delta }}^p+0+\cdots +0+a\sqrt{\Delta }+a^2)(modp)$$
   
   上式根据组合数$C_{p+1}^i=\frac{(p+1)(p)(p-1)\cdots }{1\ast 2\ast 3\ast \cdots }$，所以中间部分都为$0$，$\frac{1}{2^{p+1}}$根据逆元可以变为$\frac{1}{4}$
   再根据欧拉准则
   $${\Delta }^{\frac{p-1}{2}}\equiv -1(modp)$$
   带入上式有
   $$B^{p+1}\equiv \frac{1}{4}(-\Delta -a\sqrt{\Delta }+a\sqrt{\Delta }+a^2)(modp)$$
   有$\Delta =a^2+4b$，所以有
   $$B^{p+1}\equiv -b(modp)$$
   再根据费马小定理
   $$(B^{p+1}{)}^{p-1}\equiv (-b{)}^{p-1}\equiv 1(modp)$$
   所以
   $$B^{(p+1)(p-1)}\equiv 1(modp)$$
   同理
   $$D^{(p+1)(p-1)}\equiv 1(modp)$$
   也就是说循环节必是$(p+1)(p-1)$的倍数

解决完$p$为奇素数的情况，再考虑$p=2$的情况，我们可以通过打表打法发现对于$p=2$时任意的$a,b$可以循环节可以为$3$的倍数，但是$p=2^k$时，循环节并不遵从$g(2^k)=g(2)\ast 2^{k-1}$，其实本质在于求根公式中有一个除以$2$的操作，所以做不了同余式
对于$p=2^k$时，例如$p=4,a=2,b=1$时循环节为$4$，所以循环节并不符合$g(2^k)=g(2)\ast 2^{k-1}=6$，而是4，那么我们其实可以取这两者的最小公倍数，那么其实可以发现对于$p=2^k$的情况循环节满足是$g(2^k)=3\ast 2^k$.这样就可以找到循环节了。

顺便提一下如何判定一个数是否是模$p$的二次剩余
$a^{\frac{p-1}{2}}\equiv 1(modp)$时$a$就是模$p$下的二次剩余

以上题解可以用一个例子加深理解
$$f[n]=f[n-1]+f[n-2]\%1000000009$$
$$f_n=\frac{\sqrt{5}}{5}[(\frac{1+\sqrt{5}}{2}{)}^n-(\frac{1-\sqrt{5}}{2}{)}^n]$$
$$383008016^2\equiv 5(mod1000000009)$$
$$383008016\equiv \sqrt{5}(mod1000000009)$$
$$\frac{\sqrt{5}}{5}\equiv 276601605(mod1000000009)$$
$$\frac{1+\sqrt{5}}{2}\equiv 691504013(mod1000000009)$$
$$\frac{1-\sqrt{5}}{2}\equiv 308495997(mod1000000009)$$
$$f_n\equiv 276601605(691504013^n-308495997^n)(mod1000000009)$$
本题由于求得的循环节可能很大所以在对输入的数进行乘10取模的时候可能会爆LL，所以用快速乘代替

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 2;
LL mod = 1000000007;
int cnt,ct;
LL pri[100005];
LL num[100005];
struct Matrix
{
    LL m[N][N];
} ;
Matrix A;
Matrix I = {1, 0, 0, 1};
char str[1000005];
map<int,int>y;
int flag=0;
Matrix multi(Matrix a,Matrix b)
{
    Matrix c;
    for(int i=0; i<N; i++)
    {
        for(int j=0; j<N; j++)
        {
            c.m[i][j]  =0;
            for(int k=0; k<N; k++)
            {
                c.m[i][j] =(c.m[i][j]+a.m[i][k] * b.m[k][j]%mod)%mod;
            }
        }
    }
    return c;
}
Matrix power(Matrix A,LL n)
{
    Matrix ans = I, p = A;
    while(n>0)
    {
        if(n & 1)
        {
            ans = multi(ans,p);
            n--;
        }
        n >>= 1;
        p = multi(p,p);
    }
    return ans;
}
LL quick_mod(LL a,LL b,LL MOD)
{
    LL ans = 1;
    a %= MOD;
    while(b>0)
    {
        if(b & 1)
        {
            ans = ans * a % MOD;
            b--;
        }
        b >>= 1;
        a = a * a % MOD;
    }
    return ans;
}
LL quick_mul(LL a,LL b,LL MOD)
{
    LL ans=0;
    while(b)
    {
        if(b&1)
        {
            ans=ans+a;
            if(ans>MOD)
            {
                ans%=MOD;
                flag=1;
            }
        }
        a=a+a;
        if(a>MOD)
        {
            a%=MOD;
            flag=1;
        }
        b=b/2;
    }
    return ans;
}
LL Legendre(LL a,LL p,LL MOD)
{
    LL t = quick_mod(a,(p-1)>>1,MOD);
    if(t == 1) return 1;
    return -1;
}
void yinzi(LL x)
{
    for(int i=2; i*i<=x; i++)
    {
        if(x%i==0)
        {
            while(x%i==0)
            {
                y[i]++;
                x=x/i;
            }
        }
    }
    if(x>1)
    {
        y[x]++;
    }
}
LL lcm(LL a,LL b)
{
    return a/__gcd(a,b)*b;
}
LL get_len(LL a,LL b,LL mod)
{
    yinzi(mod);
    for(map<int,int>::iterator i=y.begin(); i!=y.end(); i++)
    {
        pri[cnt]=i->first;
        num[cnt]=i->second;
        cnt++;
    }
    LL ans=1;
    LL c=a*a+4*b;
    for(int i=0; i<cnt; i++)
    {
        LL p=1;
        if(pri[i]==2)
        {
            p*=3;
            for(int j=0; j<=num[i]; j++)
                p*=pri[i];
        }
        else if(c%pri[i]==0)
        {
            p*=pri[i]*(pri[i]-1);
            for(int j=1; j<num[i]; j++)
                p*=pri[i];
        }
        else if(Legendre(c,pri[i],pri[i])==1)
        {
            p*=(pri[i]-1);
            for(int j=1; j<num[i]; j++)
                p*=pri[i];
        }
        else
        {
            p*=(pri[i]-1)*(pri[i]+1);
            for(int j=1; j<num[i]; j++)
                p*=pri[i];
        }
        ans=lcm(ans,p);
    }
    return ans;
}
int main()
{
    LL a,b,c,d;
    cin>>c>>d>>a>>b;
    cin>>str>>mod;
    LL p=get_len(a,b,mod);
    LL n=0;
    int len=strlen(str);
    for(int i=0; i<len; i++)
    {
        n=quick_mul(n,10,p)+str[i]-'0';
    }
    if(flag)
        n+=p;
    A.m[0][0]=a;
    A.m[0][1]=b;
    A.m[1][0]=1;
    A.m[1][1]=0;
    Matrix T=power(A,n-1);
    long long ans=(T.m[0][0]*d%mod+T.m[0][1]*c%mod)%mod;
    printf("%lld\n",ans);
    return 0;
}