如何在常数时间内完成N次非规则骰子投掷模拟？

首先明确：不存在完全与n无关的严格常数时间算法——因为最终要输出每个面的出现次数，总共有k个（k是骰子面数）结果，至少需要O(k)时间来生成这些数值。但我们可以做到时间复杂度与n无关（仅和骰子面数k相关），这在n极大时，比原O(n)的方法快几个数量级。

核心思路：利用多项分布直接生成计数 #

n次独立的非规则骰子投掷，每个面的出现次数服从多项分布（Multinomial Distribution）。我们不需要模拟每一次投掷，而是直接基于这个分布生成各个面的总计数，时间复杂度为O(k)。

实现方式1：用numpy的multinomial函数（最简便高效） #

numpy内置了多项分布的抽样函数，底层实现经过高度优化，适合大n场景：

import numpy as np
from collections import Counter

def fast_roll(dist, n):
    # 直接生成每个面的出现次数
    counts = np.random.multinomial(n, dist)
    # 转换为Counter格式，和原函数输出一致
    return Counter({side: cnt for side, cnt in zip(range(1, len(dist)+1), counts)})

# 测试示例
print(fast_roll([0.1, 0.3, 0.4, 0.2], 10000000))

运行这个函数，n=1e7时的速度会比原random.choices快很多，因为它不需要循环1e7次，只需要处理k个面的抽样逻辑。

实现方式2：手动基于二项分布递推（不用numpy） #

多项分布可以分解为一系列二项分布：

1. 先从n次试验中，抽样第一个面出现的次数X₁~Binomial(n, p₁)
2. 剩下的n-X₁次试验中，抽样第二个面出现的次数X₂~Binomial(n-X₁, p₂/(1-p₁))
3. 以此类推，直到最后一个面，次数为剩余的试验次数

代码示例：

import random
from collections import Counter

def multinomial_roll(dist, n):
    remaining = n
    counts = []
    # 归一化概率，处理浮点误差
    total_p = sum(dist)
    norm_dist = [p / total_p for p in dist]
    
    for p in norm_dist[:-1]:
        if remaining <= 0:
            counts.append(0)
            continue
        # 用Python内置的二项分布抽样（3.10+支持）
        cnt = random.binomial(remaining, p)
        counts.append(cnt)
        remaining -= cnt
    # 最后一个面的次数为剩余值
    counts.append(remaining)
    
    return Counter({i+1: cnt for i, cnt in enumerate(counts)})

# 测试示例
print(multinomial_roll([0.1, 0.3, 0.4, 0.2], 10000000))

这个方法同样不需要循环n次，仅需遍历k个面，时间复杂度为O(k)。

为什么原方法慢？ #

原random.choices需要对n次投掷逐一模拟，每次都要做一次加权随机选择，时间复杂度是O(n)。当n达到1e7量级时，循环次数带来的开销会非常大；而基于多项分布的方法直接生成最终计数，完全避开了n次循环，效率提升极其明显。

内容的提问来源于stack exchange，提问作者Max