【bzoj3328】PYXFIB 题解

题目大意

求 ∑i=0⌊nk⌋CiknFik (mod p)

其中 F 表示斐波那契数列，F[0]=F[1]=1

      n<=1e18， k<=20000， p<=1e9 且保证 p 为质数、p mod k=1

这个题。。。

      这个式子涉及了很多东西，需要选手能灵活运用各种数论知识。。。
      其次还要有强大的构造能力。。。

      所以考场上碰到就打暴力吧，现在就当是学些技巧了。

题解

      首先斐波那契数列可以用矩阵来表示，设为 A=| 11 10 |，则 F(i)=Ai 的左上角。
      于是原式等于

∑i=0⌊nk⌋CiknAik (mod p)

      上标为 ik 很难看，于是换种写法：

∑i=0n [i mod k=0]CinAi (mod p)

      现在考虑怎么替换 [i mod k=0] 这个条件。
      令 p 的原根为 g，设 w=gp−1k，我们来观察一下 ∑k−1j=0wij的性质。
      1、i mod k=0：则 wij 恒等于 1，因此原式 =k；
      2、i mod k≠0：则由等比数列求和得：原式 =wik−1wi−1=1−1wi−1=0；
      综上，[i mod k=0] 可以替换成 1k∑k−1j=0wij
      （超厉害。。。）
      于是现在原式变成：

∑i=0nCinAi1k∑j=0k−1wij (mod p)

      发现 Cin 乘上一堆上标为 i 的东西很像二项式展开，于是交换一下顺序变成：

1k∑j=0k−1∑i=0nCinAiwij (mod p)

=1k∑j=0k−1(Awj+I)n (mod p)，其中 I 表示单位矩阵

      然后，枚举 j 就做完了。

代码

#include<cmath>#include<cstdio>#include<cstring>#define fo(i,a,b) for(int i=a;i<=b;i++)using namespace std;typedef long long LL;const int maxsqrtp=5e4+5;struct Arr{    LL n[2][2];};LL n,p;int k;LL mi(LL x,LL y){    LL re=1;    for(; y; y>>=1, (x*=x)%=p) if (y&1) (re*=x)%=p;    return re;}int g,p0,pr[maxsqrtp];void find_g(){    int sqrtp=sqrt(p), pp=p-1;    p0=0;    fo(i,2,sqrtp) if (pp%i==0)    {        pr[++p0]=i;        while (pp%i==0) pp/=i;    }    if (pp>1) pr[++p0]=pp;    for(g=2; ; g++)    {        bool pd=1;        fo(i,1,p0) if (mi(g,(p-1)/pr[i])==1) {pd=0; break;}        if (pd) return;    }}Arr F,I;Arr mul(Arr a,Arr b){    Arr re;    re.n[0][0]=(a.n[0][0]*b.n[0][0]+a.n[0][1]*b.n[1][0])%p;    re.n[0][1]=(a.n[0][0]*b.n[0][1]+a.n[0][1]*b.n[1][1])%p;    re.n[1][0]=(a.n[1][0]*b.n[0][0]+a.n[1][1]*b.n[1][0])%p;    re.n[1][1]=(a.n[1][0]*b.n[0][1]+a.n[1][1]*b.n[1][1])%p;    return re;}Arr mulnum(Arr a,LL b){    (a.n[0][0]*=b)%=p;    (a.n[0][1]*=b)%=p;    (a.n[1][0]*=b)%=p;    (a.n[1][1]*=b)%=p;    return a;}Arr plusI(Arr a){    (a.n[0][0]+=1)%=p;    (a.n[1][1]+=1)%=p;    return a;}Arr Mi(Arr x,LL y){    Arr re=I;    for(; y; y>>=1, x=mul(x,x)) if (y&1) re=mul(re,x);    return re;}int T;int main(){    F.n[0][0]=F.n[0][1]=F.n[1][0]=1;    F.n[1][1]=0;    I.n[0][0]=I.n[1][1]=1;    scanf("%d",&T);    while (T--)    {        scanf("%lld %d %lld",&n,&k,&p);        find_g();        LL w0=mi(g,(p-1)/k);        LL ans=0, w=1;        fo(j,0,k-1)        {            Arr t=Mi(plusI(mulnum(F,w)),n);            (ans+=t.n[0][0])%=p;            (w*=w0)%=p;        }        printf("%lld\n",ans*mi(k,p-2)%p);    }}