题面
题解
我们设\(A=\begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix}\)，那么\(A^n\)的左上角就是\(F\)的第\(n\)项
所以我们现在转化为求
\[\sum_{i=0}^n[k|i]{n\choose i}A^i
\]
\]
把\([k|i]\)单位根反演一下
\[\begin{aligned}
ans
&=\sum_{i=0}^n[k|i]{n\choose i}A^i\\
&={1\over k}\sum_{i=0}^n{n\choose i}A^i\sum_{j=0}^{k-1}\omega^{ij}_k\\
&={1\over k}\sum_{j=0}^{k-1}\sum_{i=0}^n{n\choose i}A^i\omega^{ij}_k\\
&={1\over k}\sum_{j=0}^{k-1}\left(A\omega_k^j+E\right)^n\\
\end{aligned}
\]
ans
&=\sum_{i=0}^n[k|i]{n\choose i}A^i\\
&={1\over k}\sum_{i=0}^n{n\choose i}A^i\sum_{j=0}^{k-1}\omega^{ij}_k\\
&={1\over k}\sum_{j=0}^{k-1}\sum_{i=0}^n{n\choose i}A^i\omega^{ij}_k\\
&={1\over k}\sum_{j=0}^{k-1}\left(A\omega_k^j+E\right)^n\\
\end{aligned}
\]
其中\(E\)表示单位矩阵
注意，这里能用二项式定理是因为\(A,E\)之间满足交换律，如果是普通的矩阵是不能这么做的
然后求个原根就行了（话说我今天才知道该怎么求原根）
//minamoto
#include<bits/stdc++.h>
#define R register
#define ll long long
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
const int N=5e5+5;
int p[N],P,g,tot,k,w,wn;ll n;
inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}
inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
int ksm(R int x,R int y){
R int res=1;
for(;y;y>>=1,x=mul(x,x))if(y&1)res=mul(res,x);
return res;
}
struct Matrix{
int a[2][2];
inline Matrix(){a[0][0]=a[0][1]=a[1][0]=a[1][1]=0;}
inline int* operator [](const int &x){return a[x];}
Matrix operator *(Matrix b){
Matrix res;
res[0][0]=(1ll*a[0][0]*b[0][0]+1ll*a[0][1]*b[1][0])%P;
res[0][1]=(1ll*a[0][0]*b[0][1]+1ll*a[0][1]*b[1][1])%P;
res[1][0]=(1ll*a[1][0]*b[0][0]+1ll*a[1][1]*b[1][0])%P;
res[1][1]=(1ll*a[1][0]*b[0][1]+1ll*a[1][1]*b[1][1])%P;
return res;
}
Matrix operator +(Matrix b){
Matrix res;
res[0][0]=add(a[0][0],b[0][0]);
res[0][1]=add(a[0][1],b[0][1]);
res[1][0]=add(a[1][0],b[1][0]);
res[1][1]=add(a[1][1],b[1][1]);
return res;
}
Matrix operator *(const int &x){
Matrix res;
res[0][0]=mul(a[0][0],x);
res[0][1]=mul(a[0][1],x);
res[1][0]=mul(a[1][0],x);
res[1][1]=mul(a[1][1],x);
return res;
}
}E,A,ans;
Matrix ksm(Matrix x,ll y){
Matrix res;res[0][0]=res[1][1]=1;
for(;y;y>>=1,x=x*x)if(y&1)res=res*x;
return res;
}
void getrt(){
int phi=P-1,x=phi;tot=0;
for(R int i=2;1ll*i*i<=x;++i)if(x%i==0){
p[++tot]=i;
while(x%i==0)x/=i;
}
if(x>1)p[++tot]=x;
fp(i,2,phi){
x=0;
fp(j,1,tot)if(ksm(i,phi/p[j])==1){x=1;break;}
if(!x)return g=i,void();
}
}
int main(){
// freopen("testdata.in","r",stdin);
E[0][0]=E[1][1]=A[0][0]=A[0][1]=A[1][0]=1;
int T;scanf("%d",&T);
while(T--){
scanf("%lld%d%d",&n,&k,&P);
ans=Matrix(),getrt(),w=ksm(g,(P-1)/k),wn=1;
fp(i,0,k-1)ans=ans+ksm(A*wn+E,n),wn=mul(wn,w);
printf("%d\n",mul(ans[0][0],ksm(k,P-2)));
}
return 0;
}