使用 C 语言实现的 Monte Carlo 方法求解圆周率

Question

我写了一个函数，输入一个长长整型值n 作为迭代次数。该函数应该给出一个良好的 pi估计值. 然而，当使用大数值 n时，所有值都趋近于 3.000 而不是 3.1415，所以我不确定出了什么问题。

我做错了什么吗？

以下是我的代码：

double estimate_pi(long long n){
    double randomx, randomy, equation, pi;
    long long i, incircle = 0;

    for(i = 0; i < n; i++){
        randomx = (double)(rand() % (1+1-0) + 0);
        randomy = (double)(rand() % (1+1-0) + 0);

        equation = randomx * randomx + randomy * randomy;

        if(equation <= 1){
            incircle++;
        }
    }

    pi = (long double)4 * (long double)incircle / (long double)n;

    return pi;
}

在主函数中，打印出10个pi的值：

int main(void){

    long long N;
    double pi_approx;
    int i;

    printf("Input a value of N: ");
    if(scanf("%ld", &N) != 1){
        printf("Error, input must be an integer!\n");
        exit(EXIT_SUCCESS);
    } 
    if(N < 1){
        printf("Error, the integer must be positive!\n");
        exit(EXIT_SUCCESS); 
    }

    srand(time(NULL));
    for(i = 0; i < 10; i++){
        pi_approx = estimate_pi(N);
        printf("%.10f\n", pi_approx);
    }
    return 0;
}

- Pengibaby

Answer 1

它按照应有的方式运作。问题在于实现。

C语言中rand()函数返回0到RAND_MAX之间的整数。关键字是整数。

然后你计算该整数模2的结果，可能是0或1。这给您留下了4个可能的点：(0,0)，(0,1)，(1,0)，(1,1)。

在这4个点中，只有1个在半径为1的圆外部：(1,1)。也就是说，在4个可能的点中，有3个在圆内。

您应该替换代码使用浮点值，而不是整数，这样可以计算出圆内外点的比例。

- Guille

---

非常感谢您的解释，这帮助我更好地理解了！^^ - Pengibaby

Answer 2

您需要使用浮点数随机化，或者使用一个非常大半径的圆形来实现。

因此，不是用：

    randomx = (double)(rand() % (1+1-0) + 0);
    randomy = (double)(rand() % (1+1-0) + 0);

你使用

    randomx = rand();
    randomy = rand();

您需要判断它是否落在半径为RAND_MAX的圆内。

   #define RMAX ((double)RAND_MAX*(double)RAND_MAX)
   equation <= RMAX;

你需要处理细节。阅读 man 3 rand 以了解 rand() 函数返回整数。

- alinsoar

---

非常感谢！这是我第一次编写C程序，所以还不太熟悉。但这确实帮了我很多！^^ - Pengibaby

Answer 3

你的randomx和randomy变量被限制为整数值，因为rand()函数返回一个整数。请点击此处查看实现效果。因此，你的两个变量将是 0 或 1 中的一个, 因此你的点将随机落在 (0,0)，(1,0)，(0,1) 或 (1,1) 中的一个，其中有 3/4 的概率在圆内。这就是结果为 3 的原因。

如果想要生成 0 到 1 之间的随机数，可以参考如何在 C 中生成随机浮点数。

- Louen

---

谢谢你提供的链接和解释！^^ - Pengibaby

Answer 4

为了让你的方法奏效，你需要生成从区间[0,1](或者近似这个区间)上均匀分布的double值。但是你现在所生成的是来自集合{0, 1}的随机整数，并将它们转换为double类型。这并不能得到任何类似于你目标分布的结果。最初的回答：

- John Bollinger

---

谢谢您的解释！^^ - Pengibaby

Answer 5

其他人已经指出了你的错误。我提供了一个优化版本基于另一个SO线程

   int main(void)
    {
        long points = 1000000000; // Some input
        long m = 0;
        unsigned long HAUSNUMERO = 1;
        double DIV1byMAXbyMAX = 1. / RAND_MAX / RAND_MAX;
    
        unsigned int aThreadSpecificSEED_x = HAUSNUMERO + 1
        unsigned int aThreadSpecificSEED_y = HAUSNUMERO - 1
        for(long i = 0; i < points; i++)
        {
            double x = rand_r( &aThreadSpecificSEED_x );
            double y = rand_r( &aThreadSpecificSEED_y );
            m += (1  >= ( x * x + y * y ) * DIV1byMAXbyMAX);
        }
        printf("Pi is roughly %lf\n", (double) 4*m / (double) points);
    }

- dreamcrash