BZOJ 3328: PYXFIB

前置知识单位根反演自己去浅谈单位根反演

考虑先给原来的式子变个形,\(\sum_{i=0}^{\lfloor \frac{n}{k}\rfloor} C_{n}^{ik}\cdot F_{ik}=\sum_{i=0}^n [k|i] C_n^i\cdot F_i\)

然后先把\(F_i\)做出来,我们令\(A\)斐波那契的转移矩阵,\(I\)单位矩阵

\(F_n=A^n\)的左上角,为了简化以下的式子,我们直接记\(F_n=A^n\)

则用一下单位根反演

\[\sum_{i=0}^n [k|i] C_n^i\cdot F_i \]

\[=\sum_{i=0}^n \frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{ij}\cdot C_n^i\cdot F_i \]

\[=\frac{1}{k}\sum_{j=0}^{k-1} \sum_{i=0}^nC_n^i(\omega_k^j\cdot A)^i \]

用二项式定理合起来:

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

那么就可以\(O(T\cdot k\log^3 n)\)计算了

最后提一下求一个任意质数\(p\)的原根的方法,对\(p-1\)分解质因数,然后枚举原根\(g\),若存在\(g^{\frac{p-1}{p_i}}=1\)则说明这个\(g\)不合法

#include<cstdio>
#include<cstring>
#define RI register int
#define CI const int&
using namespace std;
typedef long long LL;
const int N=2;
LL n; int t,k,mod,g,ans,pri[100],cnt,ret;
inline int quick_pow(int x,int p=mod-2,int mul=1)
{
    for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
inline void inc(int& x,CI y)
{
    if ((x+=y)>=mod) x-=mod;
}
inline int getroot(CI p)
{
    RI i,j; int x=p-1; for (i=2;i*i<=x;++i) if (x%i==0)
    for (pri[++cnt]=i;x%i==0;x/=i); for (i=2;;++i)
    {
        bool flag=1; for (j=1;j<=cnt;++j)
        if (quick_pow(i,(p-1)/pri[j])==1) { flag=0; break; }
        if (flag) return quick_pow(i,(p-1)/k);
    }
}
struct Matrix
{
    int n,mat[N][N];
    inline Matrix(void)
    {
        n=2; memset(mat,0,sizeof(mat));
    }
    inline int* operator [] (CI x) { return mat[x]; }
    friend inline Matrix operator * (Matrix A,Matrix B)
    {
        Matrix C; for (RI i=0,j,k;i<A.n;++i)
        for (j=0;j<B.n;++j) for (k=0;k<A.n;++k)
        inc(C[i][j],1LL*A[i][k]*B[k][j]%mod); return C;
    }
    friend inline Matrix operator ^ (Matrix A,LL p)
    {
        Matrix T; for (RI i=0;i<T.n;++i) T[i][i]=1;
        for (;p;p>>=1LL,A=A*A) if (p&1) T=T*A; return T;
    }
};
inline int F(CI bs)
{
    Matrix P; P[0][0]=P[0][1]=P[1][0]=bs; P[1][1]=1;
    inc(P[0][0],1); P=P^n; return P[0][0];
}
int main()
{
    for (scanf("%d",&t);t;--t)
    {
        scanf("%lld%d%d",&n,&k,&mod); g=getroot(mod);
        ans=0; ret=1; for (RI i=0;i<k;++i) inc(ans,F(ret)),ret=1LL*ret*g%mod;
        printf("%d\n",1LL*ans*quick_pow(k)%mod);
    }
    return 0;
}
posted @ 2019-10-23 16:03  空気力学の詩  阅读(290)  评论(1编辑  收藏  举报