使用球体体积进行 Pi 估计答案

【问题标题】：Pi estimation using sphere volume使用球体体积进行 Pi 估计
【发布时间】：2021-06-12 12:21:13
【问题描述】：

我的任务是计算 pi 的近似值，精度至少为 10^-6。 Monte Carlo 算法不提供所需的精度。我只需要通过球体的体积来使用计算。你有什么建议？我很高兴看到 CUDA 或纯 C++ 中的代码示例。谢谢。

【问题讨论】：

原则上，如果您选择了正确的数据类型 RNG 并让它运行足够长的时间，则通过 Monte Carlo 进行的估计应该能够提供任何准确性。那么它对您来说还不够高效吗？
@PaulG。亲爱的，我明白你的意思。我的意思是蒙特卡洛算法在这种情况下效率很低。此外，由于 x、y、z 值的随机性，相同迭代次数的结果差异很大。
是否有充分的理由在 3D 中而不是在 2D 中使用圆形区域？这只是一个挑战吗？
@PaulG。只是一个挑战:)

【解决方案1】：

从字面上看，您可以使用Darboux integral 来测量球体一个八分圆的体积。此代码示例测量半径为 2 的圆的一个象限的面积。

#include <iomanip>
#include <iostream>
#include <queue>
#include <vector>

enum class Relationship { kContained, kDisjoint, kIntersecting };

static const double kRadius = 2;

class Rectangle {
public:
  Rectangle(const double x0, const double y0, const double x1, const double y1)
      : x0_(x0), y0_(y0), x1_(x1), y1_(y1) {}

  double Area() const { return (x1_ - x0_) * (y1_ - y0_); }

  Relationship RelationshipToCircle() const {
    if (x1_ * x1_ + y1_ * y1_ < kRadius * kRadius) {
      return Relationship::kContained;
    }
    if (x0_ * x0_ + y0_ * y0_ > kRadius * kRadius) {
      return Relationship::kDisjoint;
    }
    return Relationship::kIntersecting;
  }

  std::vector<Rectangle> Subdivide() const {
    const double xm = 0.5 * (x0_ + x1_);
    const double ym = 0.5 * (y0_ + y1_);
    return {
        {x0_, y0_, xm, ym},
        {x0_, ym, xm, y1_},
        {xm, y0_, x1_, ym},
        {xm, ym, x1_, y1_},
    };
  }

private:
  double x0_;
  double y0_;
  double x1_;
  double y1_;
};

int main() {
  const Rectangle unit_rectangle{0, 0, kRadius, kRadius};
  double lower_bound = 0;
  double upper_bound = unit_rectangle.Area();
  std::queue<Rectangle> rectangles;
  rectangles.push(unit_rectangle);
  while (upper_bound - lower_bound > 1e-6) {
    const Rectangle r = rectangles.front();
    rectangles.pop();
    switch (r.RelationshipToCircle()) {
    case Relationship::kContained:
      lower_bound += r.Area();
      break;
    case Relationship::kDisjoint:
      upper_bound -= r.Area();
      break;
    case Relationship::kIntersecting:
      for (const Rectangle s : r.Subdivide()) {
        rectangles.push(s);
      }
      break;
    }
  }
  std::cout << std::setprecision(10) << 0.5 * (lower_bound + upper_bound)
            << '\n';
}

【讨论】：

是的，我草拟了代码并检查了它——它运行得非常快速和准确。我想这就是我需要的。谢谢！

【解决方案2】：

泰勒级数可用于计算精确到小数点后 5 位的 pi 值。

#include<iostream>
using namespace std;
double pi(int n) {
    double sum = 0.0;
    int sign = 1;
    for (int i = 0; i < n; ++i) {           
        sum += sign/(2.0*i+1.0);
        sign *= -1;
    }
    return 4.0*sum;
}
int main(){
   cout << pi(10000000) << endl;
}

【讨论】：

我在这个答案中看不到任何“球体的体积”。
@trincot 对，我对使用球体体积的解决方案感兴趣。否则，我本可以使用出色的 Chudnovsky 算法。