Monte Carlo 用C求圆周率的方法

Question

我编写了一个函数，它接受一个 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函数中，打印pi的10个值：

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

Answer 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() returns 整数。

Answer 2

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

现场观看here。

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

想要0到1之间的随机数可以查How to generate random float number in C

Answer 3

它正常工作。问题是实施。

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

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

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

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

Answer 4

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

Answer 5

其他人已经指出了你的错误。我提供的是优化版

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

Monte Carlo 用C求圆周率的方法

Monte Carlo Method of finding pi using C

c

pi

montecarlo