【BZOJ3328】PYXFIB（單位根反演，矩陣快速冪）

【BZOJ3328】PYXFIB（單位根反演，矩陣快速冪）

題面

BZOJ

題解

首先要求的式子是：\(\displaystyle \sum_{i=0}^n [k|i]{n\choose i}f_i\)。
斐波那契數列若是要快速算顯然就只能對應着一個矩陣，因此咱們就直接默認\(f_i\)是一個矩陣的形式。
若是沒有\([k|i]\)這個東西這個玩意看着就很像一個二項式定義的展開。
實際上二項式定理對於矩陣而言顯然是對的，不妨設斐波那契數列的轉移矩陣是\(A\)，
那麼\(\displaystyle (A+I)^n=\sum_{i=0}^n {n\choose i}A^i\)，其中\(I\)是單位矩陣。
由於有前面那個東西，因此用單位根反演直接拆：
\[\begin{aligned} \sum_{i=0}^n[k|i]{n\choose i}f_i&=\sum_{i=0}^n{n\choose i}f_i\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{ji}\\ &=\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^n{n\choose i}f_i(\omega^j)^i \end{aligned}\]
彷佛把\(\omega\)一塊兒丟進矩陣裏面就能夠直接計算了。
（也就是後面那個東西是\((A\omega^j+I)^n\)）

#include<iostream>
#include<cstdio>
using namespace std;
#define ll long long
inline ll read()
{
    ll x=0;bool t=false;char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=true,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return t?-x:x;
}
ll n;int K,MOD,g,w,ans;
int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
struct Matrix
{
    int s[2][2];
    void clear(){s[0][0]=s[0][1]=s[1][0]=s[1][1]=0;}
    void init(){s[0][0]=s[1][1]=1;s[0][1]=s[1][0]=0;}
    int*operator[](int x){return s[x];}
}A,I;
Matrix operator*(Matrix a,Matrix b)
{
    Matrix c;c.clear();
    for(int i=0;i<2;++i)
        for(int j=0;j<2;++j)
            for(int k=0;k<2;++k)
                c[i][j]=(c[i][j]+1ll*a[i][k]*b[k][j])%MOD;
    return c;
}
Matrix operator+(Matrix a,Matrix b)
{
    for(int i=0;i<2;++i)
        for(int j=0;j<2;++j)
            a[i][j]=(a[i][j]+b[i][j])%MOD;
    return a;
}
Matrix operator*(Matrix a,int b)
{
    for(int i=0;i<2;++i)
        for(int j=0;j<2;++j)
            a[i][j]=1ll*a[i][j]*b%MOD;
    return a;
}
Matrix fpow(Matrix a,ll b)
{
    Matrix s;s.init();
    while(b){if(b&1)s=s*a;a=a*a;b>>=1;}
    return s;
}
int fac[100],tot;
int Getg()
{
    int x=MOD-1;tot=0;
    for(int i=2;i*i<=x;++i)
        if(x%i==0)
        {
            fac[++tot]=i;
            while(x%i==0)x/=i;
        }
    if(x>1)fac[++tot]=x;
    for(int i=2;;++i)
    {
        bool fl=true;
        for(int j=1;j<=tot;++j)
            if(fpow(i,(MOD-1)/fac[j])==1){fl=false;break;}
        if(fl)return i;
    }
}
int main()
{
    int T=read();I.init();
    while(T--)
    {
        n=read();K=read();MOD=read();
        g=Getg();w=fpow(g,(MOD-1)/K);ans=0;
        for(int i=0,tw=1;i<K;++i,tw=1ll*tw*w%MOD)
        {
            A[0][0]=1;A[0][1]=1;A[1][0]=1;A[1][1]=0;
            A=fpow(A*tw+I,n);
            ans=(0ll+ans+A[0][1]+A[1][1])%MOD;
        }
        ans=1ll*ans*fpow(K,MOD-2)%MOD;
        printf("%d\n",ans);
    }
    return 0;
}