【BZOJ3328】PYXFIB（矩乘+单位根反演）

点此看题面

大致题意： 求\(\sum_{i=0}^{\lfloor\frac nk\rfloor}C_n^{i\times k}\times F_{i\times k}\)。

前言

斐波那契？\(n\le 10^{18}\)？矩乘啊！

然后就不会了。。。

单位根反演

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

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

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

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

我们调整枚举顺序，改为先枚举\(j\)，就得到：

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

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

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

那么代入原式就有：

\[(\frac 1k\sum_{j=0}^{k-1}\sum_{i=0}^nC_n^i\times (\begin{bmatrix}0 & 1\\1 & 1\end{bmatrix}\times\omega_k^j)^i)[1,1] \]

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

\[\frac 1k\sum_{j=0}^{k-1}(\omega_k^j\begin{bmatrix}0&1\\1&1\end{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\n",1LL*QP(k,X-2,X)*S[1][1]%X);//输出答案，注意乘1/k
	}return 0;
}