std::poisson_distribution在λ>100时精度崩溃，因double累乘PMF导致下溢/上溢；应改用指数间隔法（数值稳定）或λ>1000时用正态近似。
为什么标准库的 std::poisson_distribution 在高 λ 下精度会崩？
当 λ（均值）超过约 100，std::poisson_distribution 的默认实现（Knuth 方法或 inversion-based）开始出现显著偏差：生成的分布尾部过轻、均值漂移、甚至卡死在极低概率分支。根本原因是其内部使用 double 累乘计算概率质量函数（PMF），而 exp(-λ) * λ^k / k! 在 k ≈ λ 附近时，λ^k 和 k! 都极大，double 无法精确表示中间结果，导致下溢/上溢和舍入误差雪崩。
用 std::gamma_distribution + 指数间隔法绕过 PMF 计算
泊松过程的本质是事件间隔服从独立同分布的指数分布。生成 n 个泊松事件等价于：从时间 0 开始，不断累加 std::exponential_distribution<double>(1.0/λ)</double> 的采样值，直到总和超过观测窗口长度（通常设为 1）。事件数即为所求——这完全避开阶乘和幂运算，只依赖对数和加法，数值稳定。
- 实际只需一个
std::gamma_distribution<double>(k, 1.0/λ)</double>（k 是整数）就能直接生成第 k 个事件发生时间，但更常用的是循环累加指数变量 - λ 很大时，用中心极限定理近似：当 λ > 1000，可直接用
std::normal_distribution<double>(λ, sqrt(λ))</double>，误差通常 std::round() 后再 cast 到 int - 若必须严格服从泊松且 λ 极大（如 1e6），建议切换到任意精度库（如 Boost.Multiprecision）配合对数伽马函数重写 PMF inversion，但性能下降 10x+，一般场景不必要
手写 rejection sampling 保精度但控制成本
对中等 λ（50–500），可基于正态近似构造 proposal 分布，再用 rejection sampling 修正偏差。关键不是“通用”，而是让拒绝率可控：
- proposal 用
std::normal_distribution<double>(λ, std::sqrt(λ))</double>，但采样后强制截断到 [0, 2*λ] 区间，避免 tail 太重 - acceptance ratio 计算必须用 log-PDF：比较
log(poisson_pmf(k, λ)) - log(normal_pdf(k, λ, sqrt(λ)))，所有项都用 lgamma、log 等函数，避免直接算阶乘 - 预计算
lgamma(k+1)时，k 最大不超过 2*λ，可用 vector 缓存，避免重复调用 lgamma；若 λ 动态变化，缓存策略要重置
实操建议：按 λ 规模选方案
别硬套一个方法。λ 决定数值瓶颈在哪，也决定你愿意为精度牺牲多少速度：
立即学习“C++免费学习笔记（深入）”；
- λ ≤ 20：直接用
std::poisson_distribution<int>(λ)</int>，标准库足够准 - 20 std::exponential_distribution），简单、稳定、无额外依赖
- λ > 1000：优先尝试正态近似 +
std::round()，验证前 10000 个样本的均值与方差是否满足 |mean−λ| - 所有方案都必须用高质量随机引擎，如
std::mt19937_64，别用std::minstd_rand——后者低阶位周期短，在 rejection sampling 中易暴露相关性
最常被忽略的是：即使算法数学正确，若没校验生成序列的实际统计量（均值、方差、零计数比例），就等于没验证。尤其在 λ 接近 1000 时，不同编译器对 std::lgamma 的实现略有差异，会导致跨平台结果偏移。
