【BZOJ3328】PYXFIB(矩乘+单位根反演)

点此看题面

大致题意:(sum_{i=0}^{lfloorfrac nk floor}C_n^{i imes k} imes F_{i imes k})

前言

斐波那契?(nle 10^{18})?矩乘啊!

然后就不会了。。。

单位根反演

考虑枚举(i imes k),其实可以看作是在([0,n])范围内枚举(i)并要求(k|i),即:

[sum_{i=0}^n[k|i] imes C_n^i imes F_i ]

看到([k|i]),我们就可以套上单位根反演的公式([k|n]=frac 1ksum_{i=0}^{k-1}omega_k^{ni})得到:

[sum_{i=0}^n(frac 1ksum_{j=0}^{k-1}omega_k^{ij}) imes C_n^{i} imes F_i ]

我们调整枚举顺序,改为先枚举(j),就得到:

[frac 1ksum_{j=0}^{k-1}sum_{i=0}^nC_n^i imes F_i imes(omega_k^j)^i ]

一个众所周知的事实,(F_i)可以表示为矩阵快速幂,即:(([1,1])表示右下角的格子)

[F_i=egin{bmatrix}0 & 1\1 & 1end{bmatrix}^i[1,1] ]

那么代入原式就有:

[(frac 1ksum_{j=0}^{k-1}sum_{i=0}^nC_n^i imes (egin{bmatrix}0 & 1\1 & 1end{bmatrix} imesomega_k^j)^i)[1,1] ]

观察这个式子,我们发现,它可以直接套用二项式定理得到:(其中(E)为单位矩阵)

[frac 1ksum_{j=0}^{k-1}(omega_k^jegin{bmatrix}0&1\1&1end{bmatrix}+E)^i[1,1] ]

由于(k)很小,直接枚举(k)然后做矩乘就可以了。

代码

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define K 20000
#define LL long long
using namespace std;
LL n;int k,g,X;
struct M
{
	int a[2][2];I M() {a[0][0]=a[0][1]=a[1][0]=a[1][1]=0;}
	I int* operator [] (CI x) {return a[x];}
	I M operator + (Con M& o) Con//矩阵加法
	{
		M t;t[0][0]=(a[0][0]+o.a[0][0])%X,t[0][1]=(a[0][1]+o.a[0][1])%X,
		t[1][0]=(a[1][0]+o.a[1][0])%X,t[1][1]=(a[1][1]+o.a[1][1])%X;return t;
	}
	I M operator * (CI x) Con//矩阵乘数
	{
		M t;t[0][0]=1LL*a[0][0]*x%X,t[0][1]=1LL*a[0][1]*x%X,
		t[1][0]=1LL*a[1][0]*x%X,t[1][1]=1LL*a[1][1]*x%X;return t;
	}
	I M operator * (Con M& o) Con//矩阵乘法
	{
		M t;t[0][0]=(1LL*a[0][0]*o.a[0][0]+1LL*a[0][1]*o.a[1][0])%X,
		t[0][1]=(1LL*a[0][0]*o.a[0][1]+1LL*a[0][1]*o.a[1][1])%X,
		t[1][0]=(1LL*a[1][0]*o.a[0][0]+1LL*a[1][1]*o.a[1][0])%X,
		t[1][1]=(1LL*a[1][0]*o.a[0][1]+1LL*a[1][1]*o.a[1][1])%X;return t;
	}
	I M operator ^ (LL y) Con//矩阵快速幂
	{
		M x=*this,t;t[0][0]=t[1][1]=1;W(y) y&1&&(t=t*x,0),x=x*x,y>>=1;return t;
	}
}E,U,S;
I int QP(RI x,RI y,CI X) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
int p[K+5];I void GetGR(CI x)//求原根
{
	RI i,j,t=0,s=x-1;for(i=2;i*i<=s;++i) if(!(s%i)) {p[++t]=i;W(!(s%i)) s/=i;}//求出质因数
	for(s^1&&(p[++t]=s),i=2;i<=x;++i)//枚举数
	{
		for(j=1;j<=t;++j) if(QP(i,(x-1)/p[j],x)==1) break;if(j>t) return (void)(g=i);//判断是否为原根
	}
}
int main()
{
	E[0][0]=E[1][1]=U[0][1]=U[1][0]=U[1][1]=1;//初始化单位矩阵以及斐波那契矩阵
	RI Tt,i,w,t;scanf("%d",&Tt);W(Tt--)
	{
		scanf("%lld%d%d",&n,&k,&X),S=M();
		for(GetGR(X),w=QP(g,(X-1)/k,X),t=1,i=0;i^k;t=1LL*t*w%X,++i) S=S+(((U*t)+E)^n);//枚举k,累加矩阵,w为单位根
		printf("%d
",1LL*QP(k,X-2,X)*S[1][1]%X);//输出答案,注意乘1/k
	}return 0;
}
原文地址:https://www.cnblogs.com/chenxiaoran666/p/BZOJ3328.html