大致题意: 求(sum_{i=0}^{lfloorfrac nk floor}C_n^{i imes k} imes F_{i imes k})。
前言
斐波那契?(nle 10^{18})?矩乘啊!
然后就不会了。。。
单位根反演
考虑枚举(i imes k),其实可以看作是在([0,n])范围内枚举(i)并要求(k|i),即:
[sum_{i=0}^n[k|i] imes C_n^i imes F_i
]
看到([k|i]),我们就可以套上单位根反演的公式([k|n]=frac 1ksum_{i=0}^{k-1}omega_k^{ni})得到:
[sum_{i=0}^n(frac 1ksum_{j=0}^{k-1}omega_k^{ij}) imes C_n^{i} imes F_i
]
我们调整枚举顺序,改为先枚举(j),就得到:
[frac 1ksum_{j=0}^{k-1}sum_{i=0}^nC_n^i imes F_i imes(omega_k^j)^i
]
一个众所周知的事实,(F_i)可以表示为矩阵快速幂,即:(([1,1])表示右下角的格子)
[F_i=egin{bmatrix}0 & 1\1 & 1end{bmatrix}^i[1,1]
]
那么代入原式就有:
[(frac 1ksum_{j=0}^{k-1}sum_{i=0}^nC_n^i imes (egin{bmatrix}0 & 1\1 & 1end{bmatrix} imesomega_k^j)^i)[1,1]
]
观察这个式子,我们发现,它可以直接套用二项式定理得到:(其中(E)为单位矩阵)
[frac 1ksum_{j=0}^{k-1}(omega_k^jegin{bmatrix}0&1\1&1end{bmatrix}+E)^i[1,1]
]
由于(k)很小,直接枚举(k)然后做矩乘就可以了。
代码
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define K 20000
#define LL long long
using namespace std;
LL n;int k,g,X;
struct M
{
int a[2][2];I M() {a[0][0]=a[0][1]=a[1][0]=a[1][1]=0;}
I int* operator [] (CI x) {return a[x];}
I M operator + (Con M& o) Con//矩阵加法
{
M t;t[0][0]=(a[0][0]+o.a[0][0])%X,t[0][1]=(a[0][1]+o.a[0][1])%X,
t[1][0]=(a[1][0]+o.a[1][0])%X,t[1][1]=(a[1][1]+o.a[1][1])%X;return t;
}
I M operator * (CI x) Con//矩阵乘数
{
M t;t[0][0]=1LL*a[0][0]*x%X,t[0][1]=1LL*a[0][1]*x%X,
t[1][0]=1LL*a[1][0]*x%X,t[1][1]=1LL*a[1][1]*x%X;return t;
}
I M operator * (Con M& o) Con//矩阵乘法
{
M t;t[0][0]=(1LL*a[0][0]*o.a[0][0]+1LL*a[0][1]*o.a[1][0])%X,
t[0][1]=(1LL*a[0][0]*o.a[0][1]+1LL*a[0][1]*o.a[1][1])%X,
t[1][0]=(1LL*a[1][0]*o.a[0][0]+1LL*a[1][1]*o.a[1][0])%X,
t[1][1]=(1LL*a[1][0]*o.a[0][1]+1LL*a[1][1]*o.a[1][1])%X;return t;
}
I M operator ^ (LL y) Con//矩阵快速幂
{
M x=*this,t;t[0][0]=t[1][1]=1;W(y) y&1&&(t=t*x,0),x=x*x,y>>=1;return t;
}
}E,U,S;
I int QP(RI x,RI y,CI X) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
int p[K+5];I void GetGR(CI x)//求原根
{
RI i,j,t=0,s=x-1;for(i=2;i*i<=s;++i) if(!(s%i)) {p[++t]=i;W(!(s%i)) s/=i;}//求出质因数
for(s^1&&(p[++t]=s),i=2;i<=x;++i)//枚举数
{
for(j=1;j<=t;++j) if(QP(i,(x-1)/p[j],x)==1) break;if(j>t) return (void)(g=i);//判断是否为原根
}
}
int main()
{
E[0][0]=E[1][1]=U[0][1]=U[1][0]=U[1][1]=1;//初始化单位矩阵以及斐波那契矩阵
RI Tt,i,w,t;scanf("%d",&Tt);W(Tt--)
{
scanf("%lld%d%d",&n,&k,&X),S=M();
for(GetGR(X),w=QP(g,(X-1)/k,X),t=1,i=0;i^k;t=1LL*t*w%X,++i) S=S+(((U*t)+E)^n);//枚举k,累加矩阵,w为单位根
printf("%d
",1LL*QP(k,X-2,X)*S[1][1]%X);//输出答案,注意乘1/k
}return 0;
}