自动搬运

查看原文

来自洛谷，原作者为

	Purslane AFO

搬运于2025-08-24 23:00:04，当前版本为作者最后更新于2024-06-29 23:27:57，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

Solution

单位根反演！

---

式子有点过于不优美了，转化为：

$$\sum_{i=0}^n [k \mid i] \dbinom{n}{i} F_i$$

对 $[k \mid i]$ 使用单位根反演，设 $\omega$ 是一个 $k$ 次本原单位根，则

$$\begin{aligned} & \sum_{i=0}^n [k \mid i] \dbinom{n}{i} F_i \\ =& \sum_{i=0}^n (\dfrac{1}{k}\sum_{j=0}^{k-1} \omega^{ij}) \dbinom{n}{i} F_i \\ =& \dfrac{1}{k}\sum_{j=0}^{k-1} \sum_{i=0}^n (\omega^j)^i \dbinom{n}{i}F_i \end{aligned} $$

这个形式很像二项式定理。

设 $M = \left[ \begin{matrix}1 & 1\\ 1 & 0 \end{matrix} \right]$，那么 $F_i = \left[ \begin{matrix} 1&0\end{matrix} \right] (M^i \left[ \begin{matrix}1 \\ 0 \end{matrix} \right])$，所以我们本质上对每个 $j$ 求出

$$\sum_{i=0}^n \dbinom{n}{i}(\omega^jM)^i$$

即可。矩阵和 $I$ 相乘是有交换律的，因此可以使用二项式定理，得到 $(I + \omega^j M)^n$，容易使用矩阵快速幂优化。

---

所以 $\omega$ 是多少？现在是模意义，我们要求 $x^k \equiv 1 \pmod p$ 的一个根，其中 $x \neq 1$。

设 $g$ 是 $p$ 的一个原根，有 $(g^{\frac{p-1}{k}})^k \equiv 1 \pmod p$，令 $\omega \equiv g^{\frac{p-1}{k}} \pmod p$ 即可。

由于题目保证 $k \mid p-1$，这个高次同余方程有解。

Fun Fact：我不会枚举一个数得所有质因数，看来实际 CCF 评级只有 $3$ 级。

#include<bits/stdc++.h>
#define ll long long
#define ffor(i,a,b) for(int i=(a);i<=(b);i++)
#define roff(i,a,b) for(int i=(a);i>=(b);i--)
using namespace std;
int T,k,p;
ll n;
vector<int> frac;
inline int qpow(int base,int d) {
	int ans=1;
	while(d) {
		if(d&1) ans=1ll*ans*base%p;
		base=1ll*base*base%p,d>>=1;
	}
	return ans;
}
inline int check_gen(const int v) {
	ffor(i,0,(int)frac.size()-1) if(qpow(v,(p-1)/frac[i])==1) return 0;
	return 1;
}
int calc_gen(void) {
	ffor(i,2,p) if(check_gen(i)) return i;	
}
struct Matrix {int v[2][2];};
Matrix operator *(Matrix A,Matrix B) {
	Matrix C; memset(C.v,0,sizeof(C.v));
	ffor(i,0,1) ffor(j,0,1) ffor(k,0,1) C.v[i][k]=(C.v[i][k]+1ll*A.v[i][j]*B.v[j][k])%p;
	return C;	
}
Matrix operator ^(Matrix A,ll n) {
	Matrix C; C.v[0][0]=C.v[1][1]=1,C.v[0][1]=C.v[1][0]=0;
	while(n) {
		if(n&1) C=C*A;
		A=A*A,n>>=1;
	}
	return C;
}
signed main() {
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	cin>>T;
	while(T--) {
		cin>>n>>k>>p,frac.clear();
		int v=p-1;
		ffor(i,2,v/i) if(v%i==0) {
			frac.push_back(i);
			while(v%i==0) v/=i;	
		}
		if(v>1) frac.push_back(v);
		int ans=0,g=calc_gen(),w=qpow(g,(p-1)/k);
		ffor(j,0,k-1) {
			Matrix mt;
			mt.v[0][0]=mt.v[0][1]=mt.v[1][0]=qpow(w,j),mt.v[1][1]=0;
			mt.v[0][0]++,mt.v[1][1]++;
			mt=mt^n;
			ans=(ans+mt.v[0][0])%p;
		}
		ans=1ll*ans*qpow(k,p-2)%p;
		cout<<ans<<'\n';
	}
	return 0;
}