BZOJ 3328: PYXFIB 单位根反演+矩阵乘法+二项式定理

如果写过 LJJ 学二项式那道题的话这道题就不难了.

#include <bits/stdc++.h>

#define ll long long

#define setIO(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)

using namespace std;

int K,bu[10000],G;

ll Mod,N;

struct M

{

    ll a[2][2];

    M () { memset(a,0,sizeof(a));}

    ll * operator [] (const int &x) { return a[x]; }

    M operator * (const M &b) const

    {

        M c;

        int i,j,k;

        for(i=0;i<2;++i)  for(j=0;j<2;++j)  for(k=0;k<2;++k)  c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j]%Mod)%Mod;

        return c;

    }

}A,W;

inline ll pow_N(ll x,ll y)

{

    ll tmp=1ll;

    for(;y;y>>=1,x=x*x%Mod)     if(y&1)    tmp=tmp*x%Mod;

    return tmp;

}

inline void pow_M(ll y)

{

    while(y)

    {

        if(y&1)   A=A*W;

        y>>=1;

        W=W*W;

    }

}

inline int get_G()

{

    ll tmp=Mod-1;

    int i,j,cnt=0;

    for(i=2;i*i<=tmp;++i)

    {

        if(tmp%i==0)

        {

            bu[++cnt]=i;

            while(tmp%i==0)    tmp/=i;

        }

    }

    if(tmp>1)    bu[++cnt]=tmp;

    for(G=2;;++G)

    {

        int flag=1;

        for(j=1;j<=cnt;++j)    if(pow_N(G,(Mod-1)/bu[j])==1)   { flag=0; break; }

        if(flag==1)   break;

    }

}

void solve()

{

    int i,j;

    scanf("%lld%d%lld",&N,&K,&Mod);

    get_G();

    ll wn=pow_N(G,(Mod-1)/K),t,ans=0;

    for(i=0;i<=K-1;++i)

    {

        A[0][0]=1,A[1][1]=A[0][1]=A[1][0]=0;

        t=pow_N(wn,Mod-1-i);

        W[0][0]=W[1][0]=W[0][1]=1,W[0][0]+=t,W[1][1]=t;

        pow_M(N);

        t=A[0][0];

        ans=(ans+t*pow_N(wn,N%(Mod-1)*i))%Mod;

    }

    printf("%lld\n",ans*pow_N(K,Mod-2)%Mod);

}

int main()

{

    //  setIO("input");

    int i,j,T;

    scanf("%d",&T);

    while(T--)   solve();

    return 0;

}