DRMを利用した円周率計算 (Zの計算フラグの導入)

chudnovsky.cc

#include "computer.h"

void Chudnovsky::setXYZ(int64_t k, mpz_class& x, mpz_class& y, mpz_class& z) {
  static constexpr int64_t A = 13591409;
  static constexpr int64_t B = 545140134;
  static constexpr int64_t C = 640320;

  if (k == 0) {
    x = 1;
  } else {
    x = k * C;
    x *= k * C;
    x *= k * C / 24;
  }
  y = A + B * k;
  z = -(6 * k + 5);
  z *= 6 * k + 1;
  z *= 2 * k + 1;
}

void Chudnovsky::postProcess(mpz_class& x, mpz_class& y) {
  x *= 426880;
  mpf_sqrt_ui(pi_.get_mpf_t(), 10005);
  pi_ *= x;
  pi_ /= y;
}

computer.cc

#include "computer.h"

#include <cmath>
#include <cstdio>
#include <fstream>
#include <iomanip>
#include <iostream>

void Computer::drm(const int64_t n0,
                   const int64_t n1,
                   mpz_class& x0,
                   mpz_class& y0,
                   mpz_class& z0,
                   bool need_z) {
  if (n0 + 1 == n1) {
    setXYZ(n0, x0, y0, z0);
    return;
  }

  mpz_class x1, y1, z1;
  int64_t m = (n0 + n1) / 2;
  drm(n0, m, x0, y0, z0, true);
  drm(m, n1, x1, y1, z1, need_z);

  // y0 = x1 * y0 + y1 * z0;
  y0 *= x1;
  y1 *= z0;
  y0 += y1;

  x0 *= x1;
  if (need_z)
    z0 *= z1;
}

void Computer::compute() {
  const int64_t precision = digits_ * std::log2(10) + 10;
  pi_.set_prec(precision);

  const int64_t n = terms(digits_);
  mpz_class x, y, z;
  drm(0, n, x, y, z, false);
  postProcess(x, y);
}

void Computer::output() {
  char filename[50];
  std::snprintf(filename, 45, "pi_peria_%s.out", name());
  std::ofstream ofs(filename);
  if (!ofs.is_open()) {
    std::cerr << "Could not open the file\n";
    return;
  }
  ofs << std::setprecision(digits_ + 2) << pi_;
}

void Computer::check(const char* answer_filename) {
  if (!answer_filename)
    return;

  const int64_t precision = digits_ * std::log2(10) + 10;
  std::ifstream ifs(answer_filename);
  if (!ifs.is_open())
    return;
  mpf_class answer(0, precision);
  ifs >> answer;
  answer -= pi_;

  static constexpr int kBase = 10;
  static constexpr int kDigits = 10;
  mpf_out_str(stderr, kBase, kDigits, answer.get_mpf_t());
}

computer.h

#pragma once

#include <cstdint>

#include <gmpxx.h>

class Computer {
 public:
  Computer(int64_t digits) : digits_(digits) {}
  virtual ~Computer() = default;

  void compute();
  void output();
  void check(const char*);

 protected:
  mpf_class pi_;

 private:
  void drm(const int64_t n0,
           const int64_t n1,
           mpz_class& x0,
           mpz_class& y0,
           mpz_class& z0,
           bool need_z);
  virtual void setXYZ(int64_t k, mpz_class& x, mpz_class& y, mpz_class& z) = 0;
  virtual void postProcess(mpz_class& x, mpz_class& y) = 0;
  virtual int64_t terms(int64_t digits) const = 0;
  virtual const char* name() const = 0;

  const int64_t digits_;
};

class Chudnovsky : public Computer {
 public:
  Chudnovsky(int64_t digits) : Computer(digits) {}
  ~Chudnovsky() override = default;

 private:
  void setXYZ(int64_t k, mpz_class& x, mpz_class& y, mpz_class& z) override;
  void postProcess(mpz_class& x, mpz_class& y) override;
  int64_t terms(int64_t digits) const override { return digits / 14; }
  const char* name() const override { return "chudnovsky"; }
};

class Ramanujan : public Computer {
 public:
  Ramanujan(int64_t digits) : Computer(digits) {}
  ~Ramanujan() override = default;

 private:
  void setXYZ(int64_t k, mpz_class& x, mpz_class& y, mpz_class& z) override;
  void postProcess(mpz_class& x, mpz_class& y) override;
  int64_t terms(int64_t digits) const override { return digits / 7.98; }
  const char* name() const override { return "ramanujan"; }
};

main.cc

#include <chrono>
#include <cstdio>
#include <memory>
#include <string>

#include "computer.h"

constexpr int kDefaultDigits = 10000;

class Timer {
  using Clock = std::chrono::system_clock;
  using Ms = std::chrono::milliseconds;

 public:
  Timer(const std::string& name) : name_(name), start_(Clock::now()) {}
  ~Timer() {
    auto duration = Clock::now() - start_;
    std::printf("%7s: %.3f\n", name_.c_str(),
                std::chrono::duration_cast<Ms>(duration).count() * 1e-3);
  }

 private:
  std::string name_;
  Clock::time_point start_;
};

int main(int argc, char* argv[]) {
  int64_t digits = (argc > 1) ? std::atoi(argv[1]) : kDefaultDigits;
  if (digits <= 0)
    digits = kDefaultDigits;
  char* answer_file = (argc > 2) ? argv[2] : nullptr;

  std::unique_ptr<Computer> computer(new Chudnovsky(digits));
  // std::unique_ptr<Computer> computer(new Ramanujan(digits));
  {
    Timer all("All");
    {
      Timer timer("Compute");
      computer->compute();
    }
    computer->output();
  }
  computer->check(answer_file);

  return 0;
}

ramanujan.cc

#include "computer.h"

void Ramanujan::setXYZ(int64_t k, mpz_class& x, mpz_class& y, mpz_class& z) {
  static constexpr int64_t A = 1103;
  static constexpr int64_t B = 26390;
  static constexpr int64_t C = 396;

  if (k == 0) {
    x = 1;
  } else {
    x = k * C / 2;
    x *= k * C / 2;
    x *= k * C / 2;
    x *= k * C;
  }
  y = A + B * k;
  z = 4 * k + 1;
  z *= 2 * k + 1;
  z *= 4 * k + 3;
  z *= k + 1;
}

void Ramanujan::postProcess(mpz_class& x, mpz_class& y) {
  x *= 9801;
  y *= 4;
  mpf_sqrt_ui(pi_.get_mpf_t(), 2);
  pi_ *= x;
  pi_ /= y;
}