基于随机数互质概率求π的程序未收敛到π，原因何在？

嘿，我来帮你揪出问题所在！你提到用两个随机数互质的概率（6/π²）来近似π，但结果收敛到3.912而不是我们熟悉的≈3.1416，而且换了rand()和Mersenne Twister两种随机数生成器结果都一样——这说明问题大概率出在逻辑细节上，而不是随机数生成器本身。我来列几个最可能的坑，你对应排查下：

1. 你大概率统计了「不互质」的数对比例，而非互质的 #

咱们反推一下：根据公式π≈√(6/概率)，如果得到3.912，那代入计算的概率是6/(3.912²)≈0.392，这正好是不互质数对的理论比例（1 - 6/π²≈0.3921）。所以很可能你写代码时搞反了判断条件：

• 错误写法：if (gcd(a,b) != 1) count++，然后用count/total作为概率代入公式
• 正确写法：应该统计gcd(a,b) == 1的数对数量，用这个比例来计算

2. 随机数包含0，拉低了互质比例 #

0和任何正整数的最大公约数都是那个正整数本身（除非那个数是1），所以如果你的随机数生成范围包含0，会出现大量不互质的数对，导致统计出的互质比例远低于理论值0.6079，最终让π的近似值偏大。比如：

• 错误的分布：std::uniform_int_distribution<>(0, 100000)
• 正确的分布：应该从1开始生成正整数，比如std::uniform_int_distribution<>(1, 1000000)

3. GCD（最大公约数）函数实现错误 #

如果你的GCD逻辑写错了，会导致互质判断完全错误。比如：

• 终止条件搞反：把while (b != 0)写成while (a != 0)
• 取模运算出错：比如写成a = b % a而不是b = a % b
• 没有处理大数情况（不过一般小范围测试也能体现问题）

这里给你一个修正后的完整示例代码，你可以对比下自己的代码：

#include <iostream>
#include <random>
#include <cmath>

// 正确的迭代式GCD实现
long long gcd(long long a, long long b) {
    while (b != 0) {
        long long temp = b;
        b = a % b;
        a = temp;
    }
    return a;
}

int main() {
    const int total_trials = 1000000; // 尽量增大试验次数，收敛更准确
    int coprime_count = 0;

    // 使用Mersenne Twister生成器，从1开始生成正整数
    std::mt19937 rng(std::random_device{}());
    std::uniform_int_distribution<long long> dist(1, 1000000);

    for (int i = 0; i < total_trials; ++i) {
        long long num1 = dist(rng);
        long long num2 = dist(rng);
        if (gcd(num1, num2) == 1) {
            coprime_count++;
        }
    }

    double coprime_prob = static_cast<double>(coprime_count) / total_trials;
    double pi_estimate = sqrt(6.0 / coprime_prob);

    std::cout << "近似π值: " << pi_estimate << std::endl;
    return 0;
}

最后总结 #

你得到3.912的核心原因几乎肯定是前两点：要么搞反了互质/不互质的统计逻辑，要么随机数包含0。先检查这两个地方，应该就能让结果收敛到π附近了！

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