BZOJ3328：PYXFIB（数论+矩阵快速幂）-CSDN博客

本文链接：https://blog.csdn.net/q582116859/article/details/79103984

我觉得这是一道非常细的数学题。

题面
题意：模P下，求

\sum k | i, i < n C i n * F i

$\sum_{k|i,i<n} C_n^i*F_i$
F为斐波那契数，

n≤1e18，p≤1e9，k≤2e4，k|(p−1) $n≤1e18，p≤1e9，k≤2e4，k|(p-1)$ 。

其实数学题也是有套路的（根据本题现学）
看到菲波那契就应该是矩阵
看到整除P-1就应该是原根有关
看到组合数就应该是二项式

先搞个算斐波那契的矩阵 $A$ ={1,1}{1,0}，和单位矩阵 $I$ 。

然后把二项式定理套到矩阵里去
$(I+A)^n=\sum_{i=1}^nC_n^i*A^i$

而原式为 $\sum_{i=1}^nC_n^i*A^i*[i％k==0]$

先搞个原根g，满足g的0~p-2次幂对应了1~p-1。
再搞个单位根 $w=g^{\frac{p-1}{k}}$
和ntt的某步推导差不多， $[i％k==0]=\frac{1}{k}\sum_{j=0}^{k-1}w^{ij}$

原式变为 $\frac{1}{k}\sum_{i=1}^nC_n^i*A^i*\sum_{j=0}^{k-1}w^{ij}$

换下枚举顺序

1 k * \sum j = 0 k - 1 w j \sum i = 1 n C i n * A i w i

$\frac{1}{k}*\sum_{j=0}^{k-1}w^j\sum_{i=1}^nC_n^i*A^iw^i$
显而易见的二项式定理形式

= 1 k * \sum j = 0 k - 1 w j * w n \sum i = 1 n C i n * A i (1 w) n - i

$=\frac{1}{k}*\sum_{j=0}^{k-1}w^j*w^n\sum_{i=1}^nC_n^i*A^i(\frac{1}{w})^{n-i}$
然后后半部分上矩阵快速幂就行了。

#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>

using namespace std;
#define mmst(a, b) memset(a, b, sizeof(a))
#define mmcp(a, b) memcpy(a, b, sizeof(b))

typedef long long LL;

const int N=20020;

LL T,n,k,p,ans;
LL w,wn;

struct yy
{
    LL a[2][2];
}f,g,I;

yy operator *(yy x,yy y)
{
    yy re=f;
    for(int i=0;i<=1;i++)
    for(int j=0;j<=1;j++)
    for(int k=0;k<=1;k++)
    re.a[i][j]=(re.a[i][j]+x.a[i][k]*y.a[k][j])%p;
    return re;
}

LL cheng(LL a,LL b)
{
    LL res=1;
    for(;b;b>>=1,a=a*a%p)
    if(b&1)
    res=res*a%p;
    return res;
}

bool ok(LL x)
{
    for(LL i=2;i*i<=p;i++)
    if((p-1)%i==0&&cheng(x,(p-1)/i)==1)
    return 0;
    return 1;
}

LL find(LL p)
{
    if(p==2)
    return 1;
    LL res=2;
    for(;!ok(res);)
    res++;
    return cheng(res,(p-1)/k);
}

LL F(LL x)
{
    g.a[0][0]=x+1;
    g.a[1][1]=x;
    g.a[0][1]=g.a[1][0]=1;
    LL tu=n;
    yy hy=I;
    for(;tu;tu>>=1,g=g*g)
    if(tu&1)
    hy=hy*g;

    x=cheng(cheng(x,n),p-2);
    return x*hy.a[0][0]%p;
}

int main()
{
    I.a[0][0]=I.a[1][1]=1;

    cin>>T;
    while(T--)
    {
        ans=0;
        cin>>n>>k>>p;

        wn=find(p);
        w=cheng(cheng(wn,k-1),p-2);
        for(int i=0;i<k;i++)
        {
            ans=(ans+F(w))%p;
            w=w*wn%p;
        }

        ans=ans*cheng(k,p-2)%p;
        cout<<ans<<endl;
    }

    return 0;
}