蒙特卡洛模拟求圆周率

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

/* 该函数模拟了在一个半径为r的圆内投掷n个点的情况，并返回样本中落在圆内的点的数量 */
int monte_carlo_pi(int n, float r)
{
    int i;
    float x, y, distance;
    int p = 0; /* 落在圆内的点的数量 */

    for (i = 0; i < n; i++) {
        /* 对每个点进行均匀随机分布，即随机生成一个二维点 (x, y)，且每个点被生成的概率相等 */
        x = (float) rand() / RAND_MAX * 2 * r - r;  /* 在[-r, r]范围内随机生成x坐标 */
        y = (float) rand() / RAND_MAX * 2 * r - r;  /* 在[-r, r]范围内随机生成y坐标 */

        /* 判断该点是否在圆内 */
        distance = sqrt(x*x + y*y);  /* 计算点到原点的距离 */
        if (distance <= r) {   /* 如果距离小于等于r，则表示该点在圆内*/
            p++;  /* 投掷的点数加1 */
        }
    }

    return p;
}

/* 主函数，对输入的参数进行处理，调用monte_carlo_pi函数模拟多次投掷情况，并输出结果 */
int main(void)
{
    int i;  /* 循环计数器 */
    int n = 1000000;  /* 每次模拟中投掷的点的数量 */
    float r = 1.0;  /* 圆的半径 */
    double mu = 0, s = 0, p_mean, p_stddev, err;  /* p的均值、p的方差、样本均值、样本标准差、误差程度 */

    srand(time(NULL));  /* 设置随机数种子 */

    /* 利用不同的随机种子进行10次模拟 */
    for (i = 0; i < 10; i++) {
        /* 进行一次投掷，并累计落在圆内的点的数量 */
        int p = monte_carlo_pi(n, r);  /* 落在圆内的点的数量 */
        mu += p;  /* 各次投掷所得到的p的和，后面再除以10即为均值 */
        s += p * p;  /* 各次投掷所得到的p的平方的和，后面再除以10即为方差 */
    }

    /* 计算得到的p的平均值和标准差 */
    mu /= (n * 10); /* 根据大数定理，采样数量越多，样本均值越接近总体均值*/

    // 计算样本标准差
    s = s / (n * 10 - 1); // 样本方差：S^2=(∑(Xi-X)²)/(N-1)
    p_stddev = sqrt(s); // 样本标准差： S = √S^2

    err = 1.96 * p_stddev / sqrt(n * 10); /* 根据95%置信区间计算误差程度，1.96为正态分布的值 */

    /* 输出结果 */
    printf("样本数量：%d\n", n * 10);  /* 总共采样了多少个样本 */
    printf("圆周率的估计值：%lf\n", 4 * mu);  /* 各次投掷所得到的p的均值的均值 */
    printf("误差程度：+-%lf\n", 4 * err); /* 置信区间为平均值加减误差程度 */

    return 0;  /* 程序正常结束*/
}

样本数量：10000000
圆周率的估计值：3.141101
误差程度：+-0.045014