使用球体积的 Pi 估计

Pi estimation using sphere volume

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

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

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

从字面上看,您可以使用 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';
}