使用 C 求 pi 的蒙特卡洛方法

Posted 2023-02-16

【中文标题】使用 C 求 pi 的蒙特卡洛方法

【英文标题】：Monte Carlo Method of finding pi using C

【发布时间】：2019-08-04 12:54:48

【问题描述】：

我编写了一个函数，它接收一个 long long 值 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;

在main函数中，打印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;

【参考方案1】：

您需要使用浮点随机化，或者使用半径很大的圆。

所以不是

    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() 返回整数。

【讨论】：

非常感谢！这是我第一次写C程序，所以对它还很陌生。但这真的很有帮助！ ^^

【参考方案2】：

您的randomx 和randomy 变量被限制为一个整数 值，因为rand() 函数返回一个整数。

现场观看here。

因此，您的两个变量中的每一个都是 1 或 0，因此您的点将是 (0,0)、(1,0)、(0,1)、(1,1)、有 3:4 的机会进入圈子。因此你的结果是 3。

如果你想要一个介于 0 和 1 之间的随机数，可以查找 How to generate random float number in C。

【讨论】：

感谢您的链接和您的解释！ ^^

【参考方案3】：

它可以正常工作。问题在于实施。

C rand() 函数 returns an integer 在 0 到 RAND_MAX 的范围内。那里的关键字是整数。

然后您计算该整数模 2 的结果，它可以是 0 或 1。这会留下 4 个可能的点：(0,0)、(0,1)、(1,0)、(1, 1).

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

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

【讨论】：

非常感谢您的解释，这帮助我更了解了这一点！ ^^

【参考方案4】：

为了您的工作方法，您需要从区间 [0,1]（或大约如此）上的均匀分布中生成 double 值。您正在生成从二元素集合 0, 1 中抽取的随机整数，并将它们转换为类型 double。这不会像您为您的目的所需的分发那样远程产生任何东西。

【讨论】：

感谢您的解释！ ^^

【参考方案5】：

其他人已经指出了你的错误。我提供优化版based on another SO Thead

   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);