\(\text{Solution:}\)
\[\sum_{i=0}^{n/k}\binom{n}{i\times k}F_{i\times k} \\ =\sum_{i=0}^n \binom{n}{i}F_i[k|i] \\ =\sum_{i=0}^n\binom{n}{i}F_{i}\frac{1}{k}\sum_{j=0}^{k-1}\omega_k^{ij} \\ =\frac{1}{k}\sum_{j=0}^{k-1}(1+\omega_k^j)^nF_i \]考虑将 \(F\) 转化掉，设初始矩阵为 \(B=[1,1],\) 转移矩阵为 \(A,\)
\[A=\begin{bmatrix} 1& 1\\ 1&0\end{bmatrix} \] \[=B(\frac{1}{k}\sum_{j=0}^{k-1}(1+\omega_k^j)^n A^i) \\ =B(\frac{1}{k}\sum_{j=0}^{k-1}(I+\omega_k^jA)^n) \]直接矩阵快速幂即可。
注意一些细节：
- 原根的求解
根据结论最小原根不超过 \(p^{0.25},\) 暴力枚举原根并且暴力枚举 \(\varphi(p)=p-1\) 的素因数，判断是否存在一个素因数 \(x,\to g^{\frac{p-1}{x}}\equiv 1(\bmod p).\) 存在则不合法。
- 关于题目中说到 \(p\bmod k=1\)
这个意思就是说其 \(k\) 次单位根存在。不需要扩域。
- 第一步的转化式子
观察到只有整除的地方有贡献，所以考虑把整除去掉，直接枚举，判断哪些数有贡献，从而得到单位根反演的形式。
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e5+10;
int n,p,k,T;
inline int Add(int x,int y){return (x+y)%p;}
inline int Mul(int x,int y){return 1ll*x*y%p;}
inline int qpow(int x,int y){
int res=1;
while(y){
if(y&1)res=Mul(res,x);
x=Mul(x,x);y>>=1;
}
return res;
}
bool vis[N];
int cnt,pr[N];
vector<int>vec;
void pre(){
for(int i=2;i<=100000;++i){
if(!vis[i])pr[++cnt]=i;
for(int j=1;j<=cnt&&i*pr[j]<=100000;++j){
vis[i*pr[j]]=1;
if(i%pr[j]==0)break;
}
}
}
struct Matrix{
int a[3][3];
Matrix(){memset(a,0,sizeof a);}
Matrix operator*(const Matrix &B)const{
Matrix res;
for(int i=0;i<3;++i)
for(int j=0;j<3;++j)
for(int k=0;k<3;++k)
res.a[i][j]=Add(res.a[i][j],Mul(a[i][k],B.a[k][j]));
return res;
}
Matrix operator+(const Matrix &B)const{
Matrix res;
for(int i=0;i<3;++i)
for(int j=0;j<3;++j)
res.a[i][j]=Add(a[i][j],B.a[i][j]);
return res;
}
};
Matrix qpow(Matrix A,int b){
Matrix res;
for(int i=0;i<3;++i)res.a[i][i]=1;
while(b){
if(b&1)res=res*A;
A=A*A;b>>=1;
}
return res;
}
bool check(int x){
for(auto i:vec){
int c=qpow(x,(p-1)/i);
if(c==1)return false;
}
return true;
}
inline int getinv(int x){return qpow(x,p-2);}
void print(Matrix x){
for(int i=0;i<3;++i)
for(int j=0;j<3;++j)
printf("%lld%c",x.a[i][j],j==2?'\n':' ');
puts("");
}
signed main(){
pre();
scanf("%lld",&T);
while(T--){
vec.clear();
scanf("%lld%lld%lld",&n,&k,&p);
for(int i=1;i<=cnt;++i)
if((p-1)%pr[i]==0)
vec.push_back(pr[i]);
Matrix A,B;
B.a[1][1]=1;B.a[1][2]=1;
// print(B);
int g;
for(g=2;g<p;++g)if(check(g))break;
// cout<<g<<endl;
int Wk=qpow(g,(p-1)/k);
Matrix Ans;
for(int i=0;i<k;++i){
for(int j=0;j<3;++j)
for(int jj=0;jj<3;++jj)
A.a[j][jj]=0;
A.a[1][1]=A.a[1][2]=A.a[2][1]=qpow(Wk,i);
for(int j=0;j<3;++j)A.a[j][j]++;
// print(A);
Ans=Ans+qpow(A,n);
}
int kn=getinv(k);
for(int i=0;i<3;++i)
for(int j=0;j<3;++j)
Ans.a[i][j]=Mul(Ans.a[i][j],kn);
// print(Ans);
Ans=B*Ans;
// print(Ans);
printf("%lld\n",Ans.a[1][2]);
}
return 0;
}