【bzoj3328】PYXFIB 题解

题目大意

  求：

$$\sum_{i=0}^{\lfloor \frac{n}{k} \rfloor} C_{n}^{ik}F_{ik} \pmod p$$

  其中 $F$ 表示斐波那契数列，$F[0]=F[1]=1$。

  $n \leq 10^{18},\ k \leq 20000,\ p \leq 10^9$ 且保证 $p$ 为质数、$p\bmod k=1$

题解

  首先斐波那契数列可以用矩阵来表示，设为 $A=\begin{bmatrix} 1&1 \\ 1&0\end{bmatrix}$，则 $F(i)=A^i$ 的左上角。
  于是原式等于

$$\sum_{i=0}^{\lfloor \frac{n}{k} \rfloor} C_n^{ik}A^{ik} \pmod p$$

  上标为 $ik$ 很难看，于是换种写法：

$$\sum_{i=0}^n ~[i \bmod k=0]C_n^i A^i \pmod p$$

  现在考虑怎么替换 $[i \bmod k=0]$ 这个条件。（据说这个东西叫单位根反演？）
  令 $p$ 的原根为 $g$，设 $w=g^{\frac{p-1}{k}}$，我们来观察一下 $\sum_{j=0}^{k-1}w^{ij}$的性质。
  1、$i \bmod k=0$：则 $w^{ij}$ 恒等于 $1$，因此原式 $=k$；
  2、$i \bmod k\not=0$：则由等比数列求和得：原式 $=\frac{w^{ik}-1}{w^i-1}=\frac{1-1}{w^i-1}=0$；
  综上，$[i \bmod k=0]$ 可以替换成 $\frac{1}{k}\sum_{j=0}^{k-1}w^{ij}$
  （超厉害。。。）
  于是现在原式变成：

$$\sum_{i=0}^n C_n^i A^i \frac{1}{k}\sum_{j=0}^{k-1}w^{ij} \pmod p$$

  发现 $C_n^i$ 乘上一堆上标为 $i$ 的东西很像二项式展开，于是交换一下顺序变成：

$$\begin{aligned} &\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^n C_n^i A^i w^{ij} \pmod p \\ =&\frac{1}{k}\sum_{j=0}^{k-1}(Aw^j+I)^n \pmod p \end{aligned}$$

  其中 $I$ 表示单位矩阵。

  然后，枚举 $j$ 就做完了。

代码

#include<cmath>
#include<cstdio>
#include<cstring>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;

typedef long long LL;

const int maxsqrtp=5e4+5;

struct Arr{
	LL n[2][2];
};

LL n,p;
int k;

LL mi(LL x,LL y)
{
	LL re=1;
	for(; y; y>>=1, (x*=x)%=p) if (y&1) (re*=x)%=p;
	return re;
}

int g,p0,pr[maxsqrtp];
void find_g()
{
	int sqrtp=sqrt(p), pp=p-1;
	p0=0;
	fo(i,2,sqrtp) if (pp%i==0)
	{
		pr[++p0]=i;
		while (pp%i==0) pp/=i;
	}
	if (pp>1) pr[++p0]=pp;
	for(g=2; ; g++)
	{
		bool pd=1;
		fo(i,1,p0) if (mi(g,(p-1)/pr[i])==1) {pd=0; break;}
		if (pd) return;
	}
}

Arr F,I;
Arr mul(Arr a,Arr b)
{
	Arr re;
	re.n[0][0]=(a.n[0][0]*b.n[0][0]+a.n[0][1]*b.n[1][0])%p;
	re.n[0][1]=(a.n[0][0]*b.n[0][1]+a.n[0][1]*b.n[1][1])%p;
	re.n[1][0]=(a.n[1][0]*b.n[0][0]+a.n[1][1]*b.n[1][0])%p;
	re.n[1][1]=(a.n[1][0]*b.n[0][1]+a.n[1][1]*b.n[1][1])%p;
	return re;
}
Arr mulnum(Arr a,LL b)
{
	(a.n[0][0]*=b)%=p;
	(a.n[0][1]*=b)%=p;
	(a.n[1][0]*=b)%=p;
	(a.n[1][1]*=b)%=p;
	return a;
}
Arr plusI(Arr a)
{
	(a.n[0][0]+=1)%=p;
	(a.n[1][1]+=1)%=p;
	return a;
}
Arr Mi(Arr x,LL y)
{
	Arr re=I;
	for(; y; y>>=1, x=mul(x,x)) if (y&1) re=mul(re,x);
	return re;
}

int T;
int main()
{
	F.n[0][0]=F.n[0][1]=F.n[1][0]=1;
	F.n[1][1]=0;
	I.n[0][0]=I.n[1][1]=1;
	
	scanf("%d",&T);
	while (T--)
	{
		scanf("%lld %d %lld",&n,&k,&p);
		
		find_g();
		LL w0=mi(g,(p-1)/k);
		
		LL ans=0, w=1;
		fo(j,0,k-1)
		{
			Arr t=Mi(plusI(mulnum(F,w)),n);
			(ans+=t.n[0][0])%=p;
			(w*=w0)%=p;
		}
		
		printf("%lld\n",ans*mi(k,p-2)%p);
	}
}