[bzoj3328]PYXFIB

题目大意

求
$\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor}C_n^{i*k}*F_{i*k}$
F表示斐波那契数列

一点思路

这个形式很眼熟？
$\sum_{i=1}^nC_n^i*x^i=(1+x)^n$
那么同样，我们知道$F_i=A^i[0][0]$
A是一个2*2的矩阵，满足A[0][0]=A[0][1]=A[1][0]=1，而A[1][1]=0。
$\sum_{i=1}^nC_n^i*F_i=(I+A)^n[0][0]$
I就是单位矩阵。

继续本题

然而k的倍数才能算？？
注意到本题k是p-1的因数。
p是个质数，那么我们可以找到一个原根g。
找原根的方法是暴力枚举并检验，检验只需要检验p-1的每一个因数的次幂。
设$w=g^{(p-1)/k}$
那么$w^k=1$
有什么用呢？
注意以下等式成立
$\sum_{j=0}^{k-1}w^{ij}=[i\%k=0]*k$
为啥啊？
假如$i\%k=0$，显然正确。
否则，式子等于等比数列求和，为$\frac{1-w^{ik}}{1-w^i}$
分子为0！
因此证明了这个结论。
那么这条等式怎么用呢？
$\frac{1}{k}\sum_{i=0}^nC_n^i*A^i[0][0]*(\sum_{j=0}^{k-1}w^{ij})$就是答案了！
我们可以设$F(x)=x^{-n}*(xI+A)^n=\sum_{i=0}^nC_n^ix^{-i}A^i[0][0]$
那么答案就是$\frac{1}{k}[\sum_{j=0}^{k-1}F(w^{-j})]$

#include<cstdio>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
int o[2][2],one[2][2],a[2][2],b[2][2],c[2][2];
int ys[100000],sta[80];
int i,j,k,l,t,m,p,w,ans,ca,top;
ll n;
int qsm(int x,ll y){
    if (!y) return 1;
    int t=qsm(x,y/2);
    t=(ll)t*t%p;
    if (y%2) t=(ll)t*x%p;
    return t;
}
void work(int n){
    top=0;
    int i,t=floor(sqrt(n));
    fo(i,2,t)
        if (n%i==0){
            ys[++top]=i;
            if (n/i>i) ys[++top]=n/i;
        }
}
int getg(int n){
    work(n-1);
    int i,j;
    bool czy;
    fo(i,2,n){
        czy=1;
        fo(j,1,top)
            if (qsm(i,ys[j])==1){
                czy=0;
                break;
            }
        if (czy) break;
    }
    return i;
}
int getny(int x){
    return qsm(x,p-2);
}
int F(int x){
    a[0][0]=a[1][1]=x;
    a[0][1]=a[1][0]=0;
    int i,j,k;
    fo(i,0,1)
        fo(j,0,1)
            a[i][j]+=b[i][j];
    fo(i,0,1)
        fo(j,0,1)
            c[i][j]=one[i][j];
    top=0;
    ll t=n;
    while (t){
        sta[++top]=t%2;
        t/=2;
    }
    while (top){
        fo(i,0,1)
            fo(j,0,1)
                o[i][j]=0;
        fo(k,0,1)
            fo(i,0,1)
                fo(j,0,1)
                    o[i][j]=(o[i][j]+(ll)c[i][k]*c[k][j]%p)%p;
        fo(i,0,1)
            fo(j,0,1)
                c[i][j]=o[i][j];
        if (sta[top]){
            fo(i,0,1)
                fo(j,0,1)
                    o[i][j]=0;
            fo(k,0,1)
                fo(i,0,1)
                    fo(j,0,1)
                        o[i][j]=(o[i][j]+(ll)c[i][k]*a[k][j]%p)%p;
            fo(i,0,1)
                fo(j,0,1)
                    c[i][j]=o[i][j];
        }
        top--;
    }
    return (ll)getny(qsm(x,n))*c[0][0]%p;
}
int main(){
    one[0][0]=one[1][1]=1;
    b[0][0]=b[0][1]=b[1][0]=1;
    scanf("%d",&ca);
    while (ca--){
        scanf("%lld%d%d",&n,&k,&p);
        w=getg(p);
        w=qsm(w,(p-1)/k);
        ans=0;
        fo(i,0,k-1)
            (ans+=F(getny(qsm(w,i))))%=p;
        ans=(ll)ans*getny(k)%p;
        (ans+=p)%=p;
        printf("%d\n",ans);
    }
}