【bzoj3328】PYXFIB【单位根反演】【矩阵快速幂】_【省选基础组合数与反演单位根反演,矩阵乘法】pyxfib-CSDN博客

本文链接：https://blog.csdn.net/qq_42835815/article/details/102768010

传送门

入门小板子~

令F(i)表示斐波那契第i项，则所求即为

$\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_{i=0}^n\sum_{j=0}^{k-1}\binom{n}{i}F(i)(\omega_k^j)^i$

都知道斐波那契可以矩阵快速幂求，而这道题最开始那个组合数和斐波那契加起来正好是一个二项式展开形式。

所以我们大可以用矩阵来表示这个东东。

令 $A=\begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}$ 即初始矩阵， $I=\begin{bmatrix} 1 &0 \\ 0& 1 \end{bmatrix}$ 即单位矩阵，则这个二项式就是 (A+I)^n 的展开形式。

原式如果这样表示：

$=\frac{1}{k}\sum_{i=0}^n\sum_{j=0}^{k-1}\binom{n}{i}A^i(\omega_k^j)^i$

我们发现可以合并后两项。

$=\frac{1}{k}\sum_{i=0}^n\sum_{j=0}^{k-1}\binom{n}{i}(A\omega_k^j)^i$

所以所求即为

$\frac{1}{k}\sum_{j=0}^{k-1}(A\omega_k^j+I)^n$

咕。

#include<bits/stdc++.h>
using namespace std;
#define in read()
#define int long long
int in{
	int cnt=0,f=1;char ch=0;
	while(!isdigit(ch)){
		ch=getchar();if(ch=='-')f=-1;
	}
	while(isdigit(ch)){
		cnt=cnt*10+ch-48;
		ch=getchar();
	}return cnt*f;
}
int mod,n,k,g,w;
int fac[233],cnt;
int ksm(int a,int b){
	int sum=1;
	while(b){
		if(b&1)sum=sum*a%mod;a=a*a%mod;b>>=1;
	}return sum;
}
void getg(){
	cnt=0;int x=mod-1;
	for(int i=2;i*i<=mod;i++){
		if(x%i==0){
			fac[++cnt]=i;
			while(x%i==0)x/=i;
		}
	}if(x>1)fac[++cnt]=x;
	for(int i=2;;i++){
		int flag=1;
		for(int j=1;j<=cnt;j++){
			if(ksm(i,(mod-1)/fac[j])==1){
			
				flag=0;break;
			}
		}
		if(flag){
			g=i;return;
		}
	}
}
struct node{
	int a[3][3];
	node operator *(const node &b){
		node ans;memset(ans.a,0,sizeof(ans));
		for(int i=1;i<=2;i++){
			for(int j=1;j<=2;j++){
				for(int k=1;k<=2;k++){
					ans.a[i][j]=(ans.a[i][j]+a[i][k]*b.a[k][j]%mod)%mod;
				}
			}
		}return ans;
	}
	node operator +(const node &b){
		node c;memset(c.a,0,sizeof(c.a));
		for(int i=1;i<=2;i++){
			for(int j=1;j<=2;j++)c.a[i][j]=(a[i][j]+b.a[i][j])%mod;
		}
		return c;
	}
};
node I;
node operator *(const node &a,const int &b){
	node c;memset(c.a,0,sizeof(c.a));
	for(int i=1;i<=2;i++)for(int j=1;j<=2;j++)c.a[i][j]=a.a[i][j]*b%mod;
	return c;
}

node jzksm(node a,int b){
	node sum=I;
	while(b){
		if(b&1)sum=sum*a;
		a=a*a;b>>=1;
	}return sum;
}int ans;
signed main(){
	int t=in;I.a[2][2]=I.a[1][1]=1;
	while(t--){
		n=in;k=in;mod=in;getg();//cout<<g<<endl;
		ans=0;
		w=ksm(g,(mod-1)/k);
		for(int i=0,t=1;i<k;i++,t=t*w%mod){
			node A;A.a[1][1]=A.a[2][1]=A.a[1][2]=1;A.a[2][2]=0;
			A=jzksm(A*t+I,n);
			ans=(ans+A.a[1][1])%mod;
		}
		ans=ans*ksm(k,mod-2)%mod;cout<<ans<<'\n';
	}
	return 0;
}