【BZOJ3328】PYXFIB （单位根反演）（矩阵快速幂）

传送门

---

解析：

显然我们要求的是：
$$Ans=\sum _{i=0}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)F_i[k\mid i]$$

我们考虑将$F$表示成矩阵的形式（以下用$A$表示该矩阵）。即我们求得答案矩阵后还原数据就行了。

那么现在求：$$Ans=\sum _{i=0}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)A^i[k\mid i]$$

由于满足$mod\%k=1$，所以在$\%mod$意义下存在$k$次单位根。

直接单位根反演得到：
$$\begin{array}{rlrl}Ans& =& & \frac{1}{k}\sum _{i=0}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)A^i\sum _{j=0}^{k-1}{\omega }_k^{ij}\\ & =& & \frac{1}{k}\sum _{j=0}^{k-1}\sum _{i=0}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)(A\cdot {\omega }_k^j{)}^i\\ & =& & \frac{1}{k}\sum _{j=0}^{k-1}(A\ast {\omega }_k^j+I{)}^i\end{array}$$

直接上矩阵快速幂算就行了。

---

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const

using std::cerr;

int mod;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline int dec(int a,int b){return (a-=b)<0?a+mod:a;}
inline void Inc(int &a,int b){(a+=b)>=mod&&(a-=mod);}
inline void Dec(int &a,int b){(a-=b)<0&&(a+=mod);}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline int power(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));
	return res;
}

struct matrix{
	int a[2][2];
	
	int *operator[](int pos){return a[pos];}
	cs int *operator[](int pos)cs{return a[pos];}
	
	friend matrix operator+(cs matrix &A,cs matrix &B){
		matrix C;
		C[0][0]=add(A[0][0],B[0][0]),C[0][1]=add(A[0][1],B[0][1]);
		C[1][0]=add(A[1][0],B[1][0]),C[1][1]=add(A[1][1],B[1][1]);
		return C;
	}
	
	friend matrix operator*(cs matrix &A,cs matrix &B){
		matrix C;
		C[0][0]=add(mul(A[0][0],B[0][0]),mul(A[0][1],B[1][0]));
		C[0][1]=add(mul(A[0][0],B[0][1]),mul(A[0][1],B[1][1]));
		C[1][0]=add(mul(A[1][0],B[0][0]),mul(A[1][1],B[1][0]));
		C[1][1]=add(mul(A[1][0],B[0][1]),mul(A[1][1],B[1][1]));
		return C;
	}
}A;

cs matrix I=(matrix){1,0,0,1};

inline matrix power(matrix A,ll b){
	matrix res=I;
	for(;b;b>>=1,A=A*A)if(b&1)res=res*A;
	return res;
}

int factor[100],cnt;
inline int find_g(){
	int x=mod-1;cnt=0;
	for(int re i=2;i*i<=x;++i)
	if(x%i==0){
		factor[++cnt]=i;
		while(x%i==0)x/=i;
	}
	if(x>1)factor[++cnt]=x;
	for(int re i=2;;++i){
		bool flag=true;
		for(int re j=1;j<=cnt;++j)if(power(i,(mod-1)/factor[j])==1){flag=false;break;}
		if(flag)return i;
	}
}

int T;
ll n;
int k,ans;
signed main(){
//	freopen("pyxfib.in","r",stdin);
	scanf("%d",&T);
	while(T--){
		ans=0;
		scanf("%lld%d%d",&n,&k,&mod);
		int g=find_g(),w=power(g,(mod-1)/k);
		for(int re i=0,now=1;i<k;++i,now=mul(now,w)){
			A[0][0]=A[0][1]=A[1][0]=now;A[1][1]=0;
			A=power(A+I,n);
			Inc(ans,add(A[0][1],A[1][1]));
		}
		ans=mul(ans,power(k,mod-2));
		std::cout<<ans<<"\n";
	}
	return 0;
}