[BZOJ3328]PYXFIB/[JZOJ5150]却不悔付此华年

题目大意

给定$n,K,P$，求

$$\sum_{i=0}^{\left\lfloor\frac n k\right\rfloor}{n\choose iK}{fib}_{iK}$$
对 $P$取模的结果。
其中$fib_i$是斐波那契数列的第 $i$项。
一个测试点$T$组数据。

其中$1\le n\le{10}^{18},1\le k\le2\times10^4,1\le p\le10^9,1\le T\le20$
保证$p\equiv1\pmod K$

---

题目分析

好一道脑洞题。
题目相当于求

$$\sum_{i=0}^n{n\choose i}{fib}_i[i\ \mathrm{mod}\ K=0]$$
把斐波那契数列写成矩阵形式，令
$$A=\begin{pmatrix} 1&1\\ 1&0 \end{pmatrix}$$
则 ${fib}_i=A^i_{11}$，题目变成求
$$\sum_{i=0}^n{n\choose i}A^i[i\ \mathrm{mod}\ K=0]$$
先不考虑这个 $[i\ \mathrm{mod}\ K=0]$，然后发现这个和二项式定理很像
$$\sum_{i=0}^n{n\choose i}x^i=(x+1)^n$$
因为单位矩阵 $I$满足可交换性，于是我们有
$$\sum_{i=0}^n{n\choose i}A^i=\left(A+I\right)^n$$
现在考虑怎么解决 $[i\ \mathrm{mod}\ K=0]$，考虑构造一个函数 $F(x)$乘在 $A^i$后面，使其指数为 $i$且可以等效代替$[i\ \mathrm{mod}\ K=0]$。
怎么构造呢？考虑使用单位根的性质来构造。
令 $g$为模$p$意义下的原根，设 $\omega=g^{\frac{P-1}k}$根据单位根的求和引理，我们有：
$$\sum_{j=0}^{K-1}\omega^{ij}= \begin{cases} 0, &i\not\equiv0\pmod K\\[2ex] K, &i\equiv0\pmod K \end{cases}$$
这个怎么证明？你等比数列求和一下，把 $\omega$用 $g$表示一下就好了。
于是答案等于
$$\frac1K\sum_{i=0}^n{n\choose i}A^i\sum_{j=0}^{K-1}\omega^{ij}$$
考虑把 $j$提出去
$$\frac1K\sum_{j=0}^{K-1}\sum_{i=0}^n{n\choose i}A^i\omega^{ij}=\frac1K\sum_{j=0}^{K-1}\left(\omega^{ij}A+I\right)^n$$
枚举 $j$，然后后面的部分可以直接矩阵快速幂，然后就做完了ヽ(≧∀≦)ﾉ
原根暴力找（一般比较小），判断的话只用检查$P-1$除以 $P-1$的质因子的次幂。
时间复杂度 $O\left(\log^2 P+K\log n\right)$。

---

代码实现

#include <iostream>
#include <cstring>
#include <cstdio>

using namespace std;

typedef long long LL;

const int L=31622;

int pri[L+5],p[L+5],a[L+5];
bool mark[L+5];
int T,K,P,g,cnt;
LL n;

struct matrix
{
    int num[2][2];
    int r,c;

    inline matrix operator*(matrix const mat)const
    {
        matrix ret;ret.r=r,ret.c=mat.c,memset(ret.num,0,sizeof ret.num);
        for (int i=0;i<ret.r;++i)
            for (int k=0;k<c;++k)
                if (num[i][k])
                    for (int j=0;j<ret.c;++j)
                        (ret.num[i][j]+=1ll*num[i][k]*mat.num[k][j]%P)%=P;
        return ret;
    }

    inline matrix operator+(matrix const mat)const
    {
        matrix ret;ret.r=r,ret.c=c;
        for (int i=0;i<ret.r;++i)
            for (int j=0;j<ret.c;++j)
                ret.num[i][j]=(num[i][j]+mat.num[i][j])%P;
        return ret;
    }

    inline matrix operator*(int const k)const
    {
        matrix ret;ret.r=r,ret.c=c;
        for (int i=0;i<ret.r;++i)
            for (int j=0;j<ret.c;++j)
                ret.num[i][j]=1ll*num[i][j]*k%P;
        return ret;
    }
}A,I,res;

inline matrix operator*=(matrix &x,matrix y){return x=x*y;}
inline matrix operator+=(matrix &x,matrix y){return x=x+y;}
inline matrix operator*=(matrix &x,int y){return x=x*y;}

inline matrix operator^(matrix x,LL y)
{
    matrix ret=I;
    for (;y;y>>=1,x*=x) if (y&1) ret*=x;
    return ret;
}

int quick_power(int x,int y)
{
    int ret=1;
    for (;y;y>>=1,x=1ll*x*x%P) if (y&1) ret=1ll*ret*x%P;
    return ret;
}

void pre()
{
    mark[1]=1;
    for (int i=2;i<=L;++i)
    {
        if (!mark[i]) pri[++pri[0]]=i;
        for (int j=1;j<=pri[0];++j)
        {
            if (1ll*i*pri[j]>L) break;
            mark[i*pri[j]]=1;
            if (!(i%pri[j])) break;
        }
    }
    A.r=A.c=2,A.num[0][0]=A.num[0][1]=A.num[1][0]=1;
    I.r=I.c=2,I.num[0][0]=I.num[1][1]=1;
}

void calc()
{
    cnt=0;
    int x=P-1;
    for (int i=1;i<=pri[0]&&1ll*pri[i]*pri[i]<=P-1&&x>1;++i)
        if (!(x%pri[i]))
        {
            p[++cnt]=pri[i],a[cnt]=0;
            for (;!(x%pri[i]);x/=pri[i]) ++a[cnt];
        }
    if (x>1) p[++cnt]=x,a[cnt]=1;
    g=0;
    for (int i=2;i<P&&!g;++i)
    {
        bool judge=1;
        for (int j=1;j<=cnt&&judge;++j) judge&=quick_power(i,(P-1)/p[j])!=1;
        if (judge) g=i;
    }
    int w=quick_power(g,(P-1)/K);
    res.r=res.c=2,memset(res.num,0,sizeof res.num);
    for (int j=0;j<K;++j) res+=((A*quick_power(w,j)+I)^n);
    res*=quick_power(K,P-2);
}

int main()
{
    freopen("you.in","r",stdin),freopen("you.out","w",stdout);
    for (pre(),scanf("%d",&T);T--;) scanf("%lld%d%d",&n,&K,&P),calc(),printf("%d\n",res.num[0][0]);
    fclose(stdin),fclose(stdout);
    return 0;
}