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

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

题面

BZOJ

题解

首先要求的式子是：(displaystyle sum_{i=0}^n [k|i]{nchoose i}f_i)。
斐波那契数列如果要快速算显然就只能对应着一个矩阵，所以我们就直接默认(f_i)是一个矩阵的形式。
如果没有([k|i])这个东西这个玩意看着就很像一个二项式定义的展开。
实际上二项式定理对于矩阵而言显然是对的，不妨设斐波那契数列的转移矩阵是(A)，
那么(displaystyle (A+I)^n=sum_{i=0}^n {nchoose i}A^i)，其中(I)是单位矩阵。
因为有前面那个东西，所以用单位根反演直接拆：

[egin{aligned} sum_{i=0}^n[k|i]{nchoose i}f_i&=sum_{i=0}^n{nchoose i}f_ifrac{1}{k}sum_{j=0}^{k-1}omega_k^{ji}\ &=frac{1}{k}sum_{j=0}^{k-1}sum_{i=0}^n{nchoose i}f_i(omega^j)^i end{aligned}]

似乎把(omega)一起丢进矩阵里面就可以直接计算了。
（也就是后面那个东西是((Aomega^j+I)^n)）

#include<iostream>
#include<cstdio>
using namespace std;
#define ll long long
inline ll read()
{
	ll x=0;bool t=false;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
	if(ch=='-')t=true,ch=getchar();
	while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
	return t?-x:x;
}
ll n;int K,MOD,g,w,ans;
int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
struct Matrix
{
	int s[2][2];
	void clear(){s[0][0]=s[0][1]=s[1][0]=s[1][1]=0;}
	void init(){s[0][0]=s[1][1]=1;s[0][1]=s[1][0]=0;}
	int*operator[](int x){return s[x];}
}A,I;
Matrix operator*(Matrix a,Matrix b)
{
	Matrix c;c.clear();
	for(int i=0;i<2;++i)
		for(int j=0;j<2;++j)
			for(int k=0;k<2;++k)
				c[i][j]=(c[i][j]+1ll*a[i][k]*b[k][j])%MOD;
	return c;
}
Matrix operator+(Matrix a,Matrix b)
{
	for(int i=0;i<2;++i)
		for(int j=0;j<2;++j)
			a[i][j]=(a[i][j]+b[i][j])%MOD;
	return a;
}
Matrix operator*(Matrix a,int b)
{
	for(int i=0;i<2;++i)
		for(int j=0;j<2;++j)
			a[i][j]=1ll*a[i][j]*b%MOD;
	return a;
}
Matrix fpow(Matrix a,ll b)
{
	Matrix s;s.init();
	while(b){if(b&1)s=s*a;a=a*a;b>>=1;}
	return s;
}
int fac[100],tot;
int Getg()
{
	int x=MOD-1;tot=0;
	for(int i=2;i*i<=x;++i)
		if(x%i==0)
		{
			fac[++tot]=i;
			while(x%i==0)x/=i;
		}
	if(x>1)fac[++tot]=x;
	for(int i=2;;++i)
	{
		bool fl=true;
		for(int j=1;j<=tot;++j)
			if(fpow(i,(MOD-1)/fac[j])==1){fl=false;break;}
		if(fl)return i;
	}
}
int main()
{
	int T=read();I.init();
	while(T--)
	{
		n=read();K=read();MOD=read();
		g=Getg();w=fpow(g,(MOD-1)/K);ans=0;
		for(int i=0,tw=1;i<K;++i,tw=1ll*tw*w%MOD)
		{
			A[0][0]=1;A[0][1]=1;A[1][0]=1;A[1][1]=0;
			A=fpow(A*tw+I,n);
			ans=(0ll+ans+A[0][1]+A[1][1])%MOD;
		}
		ans=1ll*ans*fpow(K,MOD-2)%MOD;
		printf("%d
",ans);
	}
	return 0;
}