基于球体积法实现1e-6精度π值估算的方案及代码咨询

核心原理 #

基于三维单位球体积公式 $V=\frac{4}{3}\pi R^3$（$R=1$ 时 $V=\frac{4}{3}\pi$），通过数值积分计算单位球体积后反推 $\pi=\frac{3V}{4}$，该方法收敛速度远高于蒙特卡洛算法，可轻松达到1e-6及更高精度。
我们仅计算第一卦限的单位球体积（$x,y,z \in [0,1]$ 且满足 $x^2+y2+z^2 \leq 1$），乘以8即可得到整个单位球的体积，无需处理负数逻辑。

推荐实现方案：复合辛普森积分法 #

复合辛普森积分收敛阶为 $O(h^4)$，仅需每维200个采样点即可将误差控制在1e-9级别，完全满足精度要求：

• 内层积分：固定x值，计算 $f(y)=\sqrt{1-x^2-y2}$ 从 $y=0$ 到 $y=\sqrt{1-x^2}$ 的积分
• 外层积分：对x从0到1做辛普森积分，得到第一卦限体积
• 最终计算：总球体积=第一卦限体积8，$\pi=\frac{3总球体积}{4}$

纯C++实现代码 #

#define _USE_MATH_DEFINES
#include <iostream>
#include <cmath>
#include <iomanip>

// 计算固定x时y方向的辛普森积分
double simpson_y(double x, int sample_cnt) {
    double y_max = sqrt(1.0 - x * x);
    double step = y_max / sample_cnt;
    double sum = sqrt(1.0 - x * x); // y=0时的函数值，y=y_max时函数值为0无需加
    for (int i = 1; i < sample_cnt; ++i) {
        double y = i * step;
        double f_val = sqrt(1.0 - x * x - y * y);
        sum += (i % 2 == 1) ? 4 * f_val : 2 * f_val;
    }
    return sum * step / 3.0;
}

double calc_pi(int sample_cnt = 200) {
    double step = 1.0 / sample_cnt;
    double sum = simpson_y(0.0, sample_cnt) + simpson_y(1.0, sample_cnt);
    for (int i = 1; i < sample_cnt; ++i) {
        double x = i * step;
        double inner_val = simpson_y(x, sample_cnt);
        sum += (i % 2 == 1) ? 4 * inner_val : 2 * inner_val;
    }
    double total_volume = 8 * sum * step / 3.0;
    return 3 * total_volume / 4.0;
}

int main() {
    double pi = calc_pi(200);
    std::cout << std::fixed << std::setprecision(10) << "计算得到的π值：" << pi << std::endl;
    std::cout << "与标准π的误差：" << std::abs(pi - M_PI) << std::endl;
    return 0;
}

CUDA并行实现代码 #

当采样点数量较大时，可通过CUDA并行计算每个x对应的内层积分，大幅提升计算速度：

#define _USE_MATH_DEFINES
#include <iostream>
#include <cmath>
#include <iomanip>
#include <cuda_runtime.h>

__device__ double dev_sqrt(double x) {
    return __dsqrt_rn(x);
}

// 每个线程计算一个x对应的内层y方向积分
__global__ void calc_inner_integral(double* d_out, int sample_cnt) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= sample_cnt) return;
    double x = double(idx) / sample_cnt;
    double y_max = dev_sqrt(1.0 - x * x);
    double step = y_max / sample_cnt;
    double sum = dev_sqrt(1.0 - x * x);
    for (int j = 1; j < sample_cnt; ++j) {
        double y = j * step;
        double f_val = dev_sqrt(1.0 - x * x - y * y);
        sum += (j % 2 == 1) ? 4 * f_val : 2 * f_val;
    }
    d_out[idx] = sum * step / 3.0;
}

int main() {
    const int sample_cnt = 2000;
    double* h_vals = new double[sample_cnt];
    double* d_vals;
    cudaMalloc(&d_vals, sample_cnt * sizeof(double));

    dim3 block_size(256);
    dim3 grid_size((sample_cnt + block_size.x - 1) / block_size.x);
    calc_inner_integral<<<grid_size, block_size>>>(d_vals, sample_cnt);
    cudaMemcpy(h_vals, d_vals, sample_cnt * sizeof(double), cudaMemcpyDeviceToHost);

    // 外层辛普森积分在CPU端计算
    double step = 1.0 / sample_cnt;
    double sum = h_vals[0] + h_vals[sample_cnt - 1];
    for (int i = 1; i < sample_cnt - 1; ++i) {
        sum += (i % 2 == 1) ? 4 * h_vals[i] : 2 * h_vals[i];
    }
    double total_volume = 8 * sum * step / 3.0;
    double pi = 3 * total_volume / 4.0;

    std::cout << std::fixed << std::setprecision(10) << "CUDA计算得到的π值：" << pi << std::endl;
    std::cout << "与标准π的误差：" << std::abs(pi - M_PI) << std::endl;

    delete[] h_vals;
    cudaFree(d_vals);
    return 0;
}

可选低复杂度方案 #

若不需要太高计算效率，也可采用均匀网格采样法：

• 将x/y/z三个轴各拆分为1000个均匀采样点，总采样量为1e9
• 统计满足 $x^2+y2+z^2 \leq 1$ 的点数量count
• 总球体积 $V=8*count/(1000^3)$，计算得到的π误差约为5e-7，刚好满足1e-6的精度要求

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