bzoj3328 PYXFIB 矩阵快速幂+单位根反演

Description

---

给定$n,k,p$，保证$k\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$
求$$\sum _{i=0}^{\lfloor \frac{n}{k}\rfloor }\left(\genfrac{}{}{0ex}{}{n}{ki}\right)F(ki)(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
其中F(i)是斐波那契数列第i项
$n\le 10^{18},k\le 20000,p\le 10^9$保证p为质数

Solution

---

很妙的姿势，感觉这个不太好出题啊，做法貌似很套路

我们直接枚举ki，那么柿子就是$$\sum _{i=0}^n[k|i](\genfrac{}{}{0ex}{}{n}{i})F(i)(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
然后单位根反演有个结论，就是$$[n|k]=\frac{1}{n}\sum _{i=0}^{n-1}{w_n}^{ki}$$
这是一个公比鬼畜的等比数列求和，证明的话我们分情况算一下就行了

往上面代就有$$\sum _{i=0}^n\frac{1}{k}\sum _{j=0}^{k-1}{w_k}^{ij}\left(\genfrac{}{}{0ex}{}{n}{i}\right)F(i)$$
交换求和符号然后交换字母，那么就是$$\frac{1}{k}\sum _{i=0}^{k-1}\sum _{j=0}^n{w_k}^{ij}F(i)\left(\genfrac{}{}{0ex}{}{n}{j}\right)$$
根据斐波那契数列的矩阵推法，我们设$$A=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$$
那么$F(i)=A^i[0,0]$。这里矩阵可能长成别的样子，同样也是对的
于是那个就可以二项式弄回去了，也就是$$\frac{1}{k}\sum _{i=0}^{k-1}{\left({w_k}^{i}A+I\right)}^n$$
这里$I$是单位矩阵，$A\ast I=A$

注意到这里的答案是膜意义下的，那么$$w_k=g^{\frac{p-1}{k}}$$这里g是p的原根，这也是为什么题目保证了p膜k余1

于是我们求一个原根，做快速幂和矩阵快速幂就可以得到答案了

Code

---

#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <vector>

#define rep(i,st,ed) for (int i=st;i<=ed;++i)
#define fill(x,t) memset(x,t,sizeof(x))

typedef long long LL;

int v[233333],MOD;

void upd(LL &x,LL v) {
	x+=v; (x>=MOD)?(x-=MOD):0;
}

struct Mat {
	LL r[2][2];

	Mat() {fill(r,0);}

	LL* operator [](int x) {
		return r[x];
	}

	Mat operator *(LL x) {
		Mat C=*this;
		C[0][0]=C[0][0]*x%MOD;
		C[0][1]=C[0][1]*x%MOD;
		C[1][0]=C[1][0]*x%MOD;
		C[1][1]=C[1][1]*x%MOD;
		return C;
	}

	Mat operator +(Mat B) {
		Mat C=*this;
		upd(C[0][0],B[0][0]);
		upd(C[0][1],B[0][1]);
		upd(C[1][0],B[1][0]);
		upd(C[1][1],B[1][1]);
		return C;
	}

	Mat operator *(Mat B) {
		Mat A=*this,C;
		C[0][0]=(A[0][0]*B[0][0]+A[0][1]*B[1][0])%MOD;
		C[0][1]=(A[0][0]*B[0][1]+A[0][1]*B[1][1])%MOD;
		C[1][0]=(A[1][0]*B[0][0]+A[1][1]*B[1][0])%MOD;
		C[1][1]=(A[1][0]*B[0][1]+A[1][1]*B[1][1])%MOD;
		return C;
	}
} A,I,B,S;

LL read() {
	LL x=0,v=1; char ch=getchar();
	for (;ch<'0'||ch>'9';v=(ch=='-')?(-1):v,ch=getchar());
	for (;ch<='9'&&ch>='0';x=x*10LL+ch-'0',ch=getchar());
	return x*v;
}

LL ksm(LL x,LL dep) {
	LL res=1;
	for (;dep;dep>>=1,x=x*x%MOD) {
		(dep&1)?(res=res*x%MOD):0;
	}
	return res;
}

Mat ksm(Mat A,LL dep) {
	Mat C=I;
	for (;dep;dep>>=1,A=A*A) {
		if (dep&1) C=C*A;
	}
	return C;
}

int get_g(int MOD) {
	int tmp=MOD-1; v[0]=0;
	for (int i=2;i*i<=tmp;++i) {
		if (tmp%i==0) {
			while (tmp%i==0) tmp/=i;
			v[++v[0]]=i;
		}
	}
	if (tmp!=1) v[++v[0]]=tmp;
	rep(g,2,MOD) {
		bool flag=false;
		rep(i,1,v[0]) {
			if (ksm(g,(MOD-1)/v[i])==1) {
				flag=true; break;
			}
		}
		if (!flag) return g;
	}
}

int main(void) {
	A[0][0]=A[0][1]=A[1][0]=I[0][0]=I[1][1]=1;
	int T; scanf("%d",&T);
	for (;T--;) {
		LL n,k; scanf("%lld%lld%d",&n,&k,&MOD);
		LL g=get_g(MOD),wn=ksm(g,(MOD-1)/k),w=1;
		LL ans=0;
		rep(i,0,k-1) {
			S=ksm(A*w+I,n);
			upd(ans,S[0][0]);
			w=w*wn%MOD;
		}
		ans=ans*ksm(k,MOD-2)%MOD;
		printf("%lld\n", ans);
	}
	return 0;
}