Pi estimation using sphere volume

My task is to calculate the approximate value of pi with an accuracy of at least 10^-6. The Monte Carlo algorithm does not provide the required accuracy. I need to use the calculation only through the volume of the sphere. What do you advise? I would be glad to see examples of code in CUDA or pure C++. Thank you.

To be completely literal about it, you could use a Darboux integral to measure the volume of one octant of the sphere. This code sample measures the area of one quadrant of a circle of radius 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';
}