自动搬运

查看原文

来自洛谷，原作者为

	VinstaG173 Si Dieu n'existait pas, il faudrait l'inventer.

搬运于2025-08-24 23:12:22，当前版本为作者最后更新于2025-04-05 19:33:36，作者可能在搬运后再次修改，您可在原文处查看最新版

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

以下是正文

---

首先我们有一个贪心策略：每当剩余未确定结果的情况多于 $n$ 个时，将其中的 $cn$ 个确定结果（每种菜各分配 $c$ 种可能），剩下的情况继续掷骰子。

则我们的过程会是这样的：掷了 $v$ 次骰子时会有 $r_v=k^v\bmod{n}$ 种情况未确定结果，故第 $v+1$ 次会有 $n\lfloor\frac{kr_v}n\rfloor$ 种情况被确定结果。那么期望次数就是

$$\sum_{v\ge1}\dfrac{n\lfloor\frac{kr_v}n\rfloor}{k^{v+1}}(v+1). $$

由 Euler 定理，$r_v=k^v\bmod{n}$ 会进入循环，设出现循环的最小位置为 $k^b\equiv k^{b+d}\pmod{n}$，则对于 $v\ge b$ 有：

$$\begin{aligned} E&=\sum_{v\ge b}\dfrac{n\lfloor\frac{kr_v}n\rfloor}{k^{v+1}}(v+1)\\ &=\sum_{v=b+1}^{b+d}n\left\lfloor\frac{kr_{v-1}}n\right\rfloor\sum_{i\ge0}\dfrac{v+id}{k^{v+id}}\\ &=\sum_{v=b+1}^{b+d}\dfrac{vn\lfloor\frac{kr_{v-1}}n\rfloor}{k^v}+\dfrac1{k^d}\left(E+d\sum_{v\ge b}\dfrac{n\lfloor\frac{kr_v}n\rfloor}{k^{v+1}}\right)\\ &=\sum_{v=b+1}^{b+d}\dfrac{vn\lfloor\frac{kr_{v-1}}n\rfloor}{k^v}+\dfrac1{k^d}\left(E+\dfrac{dr_b}{k^b}\right), \end{aligned}$$

最后一步是由于括号里的和式本质上是 $\ge b+1$ 次才确定结果的概率，根据贪心策略其值即为第 $b$ 次没有确定结果的概率 $\frac{r_b}{k^b}$。

从而我们可以列出方程求解：

$$\begin{aligned} \left(1-\dfrac1{k^d}\right)E&=\dfrac{dr_b}{k^{b+d}}+\sum_{v=b+1}^{b+d}\dfrac{vn\lfloor\frac{kr_{v-1}}n\rfloor}{k^v},\\ E&=\left(1-\dfrac1{k^d}\right)^{-1}\left(\dfrac{dr_b}{k^{b+d}}+\sum_{v=b+1}^{b+d}\dfrac{vn\lfloor\frac{kr_{v-1}}n\rfloor}{k^v}\right). \end{aligned}$$

因此我们可以用桶记录 $r_v$ 求出循环，预处理 $k^{-v}$ 后可 $O(d)$ 求出 $E$，由于 $E$ 是 $v\ge b$ 时的答案，因此对于 $v<b$ 直接 $O(b)$ 计算即可得到最终结果。由 Euler 定理 $b,d=O(\varphi(n))$，故总时间复杂度 $O(\varphi(n))$，空间复杂度 $O(n)$。注意判定 $n\mid k^b$ 的情况。

Code:

int n,k,m,b,c,d;
inline ll qpw(ll x,int v){
	ll r=1;while(v){
		(v&1)&&(r=r*x%m);
		x=x*x%m,v>>=1;
	}return r;
}

ll x,y,tmp,ans;
ll inv[3000007];
int vis[3000007];
inline void solve(){
	cin>>n>>k>>m;x=inv[0]=1;
	inv[1]=qpw(k,m-2);
	for(c=1;!vis[x];++c,x=x*k%n){
		vis[x]=c,inv[c+1]=inv[c]*inv[1]%m;
	}d=c-vis[x];b=vis[x]-1;
	if(x){
		for(int i=1;i<=d;++i,x=x*k%n){
			tmp=(tmp+(x*k/n*n)%m*inv[i+b]%m*(i+b))%m;
		}tmp=(tmp+inv[d]*d%m*x%m*inv[b])%m;
		tmp=tmp*qpw(m+1-inv[d],m-2)%m;
	}y=1;
	for(int i=1;i<=b;++i,y=y*k%n){
		ans=(ans+(y*k/n*n)%m*inv[i]%m*i)%m;
	}ans=(ans+tmp)%m;
	cout<<ans<<"\n";
}