【题解】bzoj3328 PYXFIB

PYXFIB

( ext{Solution:})

[sum_{i=0}^{n/k}inom{n}{i imes k}F_{i imes k} \ =sum_{i=0}^n inom{n}{i}F_i[k|i] \ =sum_{i=0}^ninom{n}{i}F_{i}frac{1}{k}sum_{j=0}^{k-1}omega_k^{ij} \ =frac{1}{k}sum_{j=0}^{k-1}(1+omega_k^j)^nF_i ]

考虑将 (F) 转化掉，设初始矩阵为 (B=[1,1],) 转移矩阵为 (A,)​

[A=egin{bmatrix} 1& 1\ 1&0end{bmatrix} ]

[=B(frac{1}{k}sum_{j=0}^{k-1}(1+omega_k^j)^n A^i) \ =B(frac{1}{k}sum_{j=0}^{k-1}(I+omega_k^jA)^n) ]

直接矩阵快速幂即可。

注意一些细节：

• 原根的求解

根据结论最小原根不超过 (p^{0.25},) 暴力枚举原根并且暴力枚举 (varphi(p)=p-1) 的素因数，判断是否存在一个素因数 (x, o g^{frac{p-1}{x}}equiv 1(mod p).) 存在则不合法。

• 关于题目中说到 (pmod k=1)

这个意思就是说其 (k) 次单位根存在。不需要扩域。

• 第一步的转化式子

观察到只有整除的地方有贡献，所以考虑把整除去掉，直接枚举，判断哪些数有贡献，从而得到单位根反演的形式。

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e5+10;
int n,p,k,T;
inline int Add(int x,int y){return (x+y)%p;}
inline int Mul(int x,int y){return 1ll*x*y%p;}
inline int qpow(int x,int y){
    int res=1;
    while(y){
        if(y&1)res=Mul(res,x);
        x=Mul(x,x);y>>=1;
    }
    return res;
}
bool vis[N];
int cnt,pr[N];
vector<int>vec;
void pre(){
    for(int i=2;i<=100000;++i){
        if(!vis[i])pr[++cnt]=i;
        for(int j=1;j<=cnt&&i*pr[j]<=100000;++j){
            vis[i*pr[j]]=1;
            if(i%pr[j]==0)break;
        }
    }
}
struct Matrix{
    int a[3][3];
    Matrix(){memset(a,0,sizeof a);}
    Matrix operator*(const Matrix &B)const{
        Matrix res;
        for(int i=0;i<3;++i)
            for(int j=0;j<3;++j)
                for(int k=0;k<3;++k)
                    res.a[i][j]=Add(res.a[i][j],Mul(a[i][k],B.a[k][j]));
        return res;
    }
    Matrix operator+(const Matrix &B)const{
        Matrix res;
        for(int i=0;i<3;++i)
            for(int j=0;j<3;++j)
                res.a[i][j]=Add(a[i][j],B.a[i][j]);
        return res;
    }
};
Matrix qpow(Matrix A,int b){
    Matrix res;
    for(int i=0;i<3;++i)res.a[i][i]=1;
    while(b){
        if(b&1)res=res*A;
        A=A*A;b>>=1;
    }
    return res;
}
bool check(int x){
    for(auto i:vec){
        int c=qpow(x,(p-1)/i);
        if(c==1)return false;
    }
    return true;
}
inline int getinv(int x){return qpow(x,p-2);}
void print(Matrix x){
    for(int i=0;i<3;++i)
        for(int j=0;j<3;++j)
            printf("%lld%c",x.a[i][j],j==2?'
':' ');
    puts("");
}
signed main(){
    pre();
    scanf("%lld",&T);
    while(T--){
        vec.clear();
        scanf("%lld%lld%lld",&n,&k,&p);
        for(int i=1;i<=cnt;++i)
            if((p-1)%pr[i]==0)
                vec.push_back(pr[i]);
        Matrix A,B;
        B.a[1][1]=1;B.a[1][2]=1;
        // print(B);
        int g;
        for(g=2;g<p;++g)if(check(g))break;
        // cout<<g<<endl;
        int Wk=qpow(g,(p-1)/k);
        Matrix Ans;
        for(int i=0;i<k;++i){
            for(int j=0;j<3;++j)
                for(int jj=0;jj<3;++jj)
                    A.a[j][jj]=0;
            A.a[1][1]=A.a[1][2]=A.a[2][1]=qpow(Wk,i);
            for(int j=0;j<3;++j)A.a[j][j]++;
            // print(A);
            Ans=Ans+qpow(A,n);
        }
        int kn=getinv(k);
        for(int i=0;i<3;++i)
            for(int j=0;j<3;++j)
                Ans.a[i][j]=Mul(Ans.a[i][j],kn);
        // print(Ans);
        Ans=B*Ans;
        // print(Ans);
        printf("%lld
",Ans.a[1][2]);
    }
    return 0;
}