BZOJ3328: PYXFIB

我们要求这个东西 $\sum_{i=0}^nC_n^iFi[i\ mod\ k=0]$
先不考虑$k|i$这个条件，即只计算$\sum_{i=0}^nC_n^iFi$
设出Fib数列的转移矩阵$A$，那么$Fi$就是$A^i$的左上角
再考虑二项式定理$(a+b)^n=\sum_{i=0}^na^ib^{n-i}$，在矩阵运算中同样适用（满足分配律，结合律）
那么我们可以得到

$$\sum_{i=0}^nC_n^iFi=\sum_{i=0}^n(_i^n)A^iI^{n-i}=(A+I)^n$$ $I$是单位矩阵

---

接下来考虑构造一个$[i\ mod\ k=0]$的函数
注意到 $k\equiv1\ Mod\ p$，即 $k|p-1,phi(p)=p-1$
因此我们求出 $p$的原根$g$，令 $w=g^{\frac{p-1}{k}}$，那么 $w^i\equiv1\mod p$当且仅当 $k|i$
并且有这么个柿子 $\sum_{j=0}^{k-1}w^{ij}=k[i\ mod\ k=0]$
证明：当 $k|i$时， $w^i=1$，显然和为 $k$，否则是个以1为首项，$w^i$为公比的等比数列， $\dfrac{1-w^{ik}}{1-w^i}$，分子为0，和为0
那么我们就构造出了我们想要的函数：
$\frac{1}{k}\sum_{j=0}^{k-1}w^{ij}=[i\ mod\ k=0]$

---

那么最后的答案就是
$\sum_{i=0}^n(^n_i)Fi[i\ mod\ k=0]$
=$\sum_{i=0}^n(^n_i)A^i\frac{1}{k}\sum_{j=0}^{k-1}w^{ij}$
=$\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^n(^n_i)A^iI^{n-i}w^{ji}$
=$\frac{1}{k}\sum_{j=0}^{k-1}w^{nj}\sum_{i=0}^n(^n_i)A^i(Iw^{-j})^{n-i}$
=$\frac{1}{k}\sum_{j=0}^{k-1}w^{nj}(A+Iw^{-j})^n$
注意这里$w^{ij}$只能和$I$结合到一起不能$Aw$这样子，因为一般矩阵$A$是不满足交换律的，但是单位矩阵$I$特殊，他满足交换律，所以可以弄成$Iw^{-j}$这样子

code：

#include<set>
#include<map>
#include<deque>
#include<queue>
#include<stack>
#include<cmath>
#include<ctime>
#include<bitset>
#include<string>
#include<vector>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<climits>
#include<complex>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;

const int maxn = 110;

int pwmod(int x,ll k,int mod)
{
    int re=1;
    for(;k;k>>=1,x=(ll)x*x%mod) if(k&1ll)
        re=(ll)re*x%mod;
    return re;
}
int Inv(int x,int mod){return pwmod(x,mod-2,mod);}
int t[maxn],tp;
int get_g(int P)
{
    int p=P-1;
    tp=0;
    for(int i=2;i*i<=p;i++) if(p%i==0)
    {
        t[++tp]=i; while(p%i==0) p/=i;
    }
    if(p>1) t[++tp]=p;
    p=P-1;
    for(int i=2;;i++)
    {
        int ok=1;
        for(int j=1;j<=tp;j++)
            if(pwmod(i,p/t[j],P)==1) { ok=0;break; }
        if(ok) return i;
    }
}

ll n;
int K,mod,ans;
inline void add(int &a,const int &b){a+=b;if(a>=mod)a-=mod;}
struct mat
{
    int a[2][2];
    mat(){memset(a,0,sizeof a);}
    friend inline mat operator *(const mat &x,const mat &y)
    {
        mat re;
        for(int i=0;i<2;i++) for(int k=0;k<2;k++)
            for(int j=0;j<2;j++) add(re.a[i][j],(ll)x.a[i][k]*y.a[k][j]%mod);
        return re;
    }
    friend inline mat operator ^(const mat &x,const int &y)
    {
        mat re;
        for(int i=0;i<2;i++) for(int j=0;j<2;j++) re.a[i][j]=(ll)x.a[i][j]*y%mod;
        return re;
    }
    friend inline mat operator +(const mat &x,const mat &y)
    {
        mat re;
        for(int i=0;i<2;i++) for(int j=0;j<2;j++) re.a[i][j]=(x.a[i][j]+y.a[i][j])%mod;
        return re;
    }
}A,one,W;
mat pw(mat x,ll k)
{
    mat re=one;
    for(;k;k>>=1,x=x*x) if(k&1ll)
        re=re*x;
    return re;
}

int main()
{
    //freopen("tmp.in","r",stdin);
    //freopen("tmp.out","w",stdout);

    A.a[0][0]=A.a[0][1]=A.a[1][0]=1;
    one.a[0][0]=one.a[1][1]=1;

    int T;scanf("%d",&T);
    while(T--)
    {
        ans=0;

        scanf("%lld%d%d",&n,&K,&mod);
        int g=get_g(mod);
        int w=pwmod(g,(mod-1)/K,mod),inv=Inv(w,mod),wj=1;
        W=one;
        for(int j=0;j<K;j++)
        {
            mat temp=A+W,tn=pw(temp,n);
            add(ans,(ll)pwmod(wj,n,mod)*tn.a[0][0]%mod);
            if(ans<0) ans+=mod;
            W=W^inv; wj=(ll)wj*w%mod;
        }

        printf("%lld\n",(ll)ans*Inv(K,mod)%mod);
    }

    return 0;
}