题解：P10664 BZOJ3328 PYXFIB

· · 题解

Solution

单位根反演！

式子有点过于不优美了，转化为：

\sum_{i=0}^n [k \mid i] \dbinom{n}{i} F_i

对 [k \mid i] 使用单位根反演，设 \omega 是一个 k 次本原单位根，则

\begin{aligned} & \sum_{i=0}^n [k \mid i] \dbinom{n}{i} F_i \\ =& \sum_{i=0}^n (\dfrac{1}{k}\sum_{j=0}^{k-1} \omega^{ij}) \dbinom{n}{i} F_i \\ =& \dfrac{1}{k}\sum_{j=0}^{k-1} \sum_{i=0}^n (\omega^j)^i \dbinom{n}{i}F_i \end{aligned}

这个形式很像二项式定理。

设 M = \left[ \begin{matrix}1 & 1\\ 1 & 0 \end{matrix} \right]，那么 F_i = \left[ \begin{matrix} 1&0\end{matrix} \right] (M^i \left[ \begin{matrix}1 \\ 0 \end{matrix} \right])，所以我们本质上对每个 j 求出

\sum_{i=0}^n \dbinom{n}{i}(\omega^jM)^i

即可。矩阵和 I 相乘是有交换律的，因此可以使用二项式定理，得到 (I + \omega^j M)^n，容易使用矩阵快速幂优化。

所以 \omega 是多少？现在是模意义，我们要求 x^k \equiv 1 \pmod p 的一个根，其中 x \neq 1。

设 g 是 p 的一个原根，有 (g^{\frac{p-1}{k}})^k \equiv 1 \pmod p，令 \omega \equiv g^{\frac{p-1}{k}} \pmod p 即可。

由于题目保证 k \mid p-1，这个高次同余方程有解。

Fun Fact：我不会枚举一个数得所有质因数，看来实际 CCF 评级只有 3 级。

#include<bits/stdc++.h>
#define ll long long
#define ffor(i,a,b) for(int i=(a);i<=(b);i++)
#define roff(i,a,b) for(int i=(a);i>=(b);i--)
using namespace std;
int T,k,p;
ll n;
vector<int> frac;
inline int qpow(int base,int d) {
    int ans=1;
    while(d) {
        if(d&1) ans=1ll*ans*base%p;
        base=1ll*base*base%p,d>>=1;
    }
    return ans;
}
inline int check_gen(const int v) {
    ffor(i,0,(int)frac.size()-1) if(qpow(v,(p-1)/frac[i])==1) return 0;
    return 1;
}
int calc_gen(void) {
    ffor(i,2,p) if(check_gen(i)) return i;  
}
struct Matrix {int v[2][2];};
Matrix operator *(Matrix A,Matrix B) {
    Matrix C; memset(C.v,0,sizeof(C.v));
    ffor(i,0,1) ffor(j,0,1) ffor(k,0,1) C.v[i][k]=(C.v[i][k]+1ll*A.v[i][j]*B.v[j][k])%p;
    return C;   
}
Matrix operator ^(Matrix A,ll n) {
    Matrix C; C.v[0][0]=C.v[1][1]=1,C.v[0][1]=C.v[1][0]=0;
    while(n) {
        if(n&1) C=C*A;
        A=A*A,n>>=1;
    }
    return C;
}
signed main() {
    ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
    cin>>T;
    while(T--) {
        cin>>n>>k>>p,frac.clear();
        int v=p-1;
        ffor(i,2,v/i) if(v%i==0) {
            frac.push_back(i);
            while(v%i==0) v/=i; 
        }
        if(v>1) frac.push_back(v);
        int ans=0,g=calc_gen(),w=qpow(g,(p-1)/k);
        ffor(j,0,k-1) {
            Matrix mt;
            mt.v[0][0]=mt.v[0][1]=mt.v[1][0]=qpow(w,j),mt.v[1][1]=0;
            mt.v[0][0]++,mt.v[1][1]++;
            mt=mt^n;
            ans=(ans+mt.v[0][0])%p;
        }
        ans=1ll*ans*qpow(k,p-2)%p;
        cout<<ans<<'\n';
    }
    return 0;
}