pi.cpp

// PI

#include <algorithm>
#include <chrono>
#include <format>
#include <limits>
#include <print>
#include <ranges>
#include <span>
#include <string>

using number = unsigned long long;

constexpr int FIG     = 8;
constexpr number BASE = 0xffffffffULL;
constexpr number NN   = 100000000ULL;

// C = 1.0 / (4 * log10(2)) / FIG;
template<int X> class PI {
private:
    static constexpr double C_COEFF = 0.830482023722 / FIG;
    static constexpr int M = X / FIG + 2;
    static constexpr int N = static_cast<int>(X * C_COEFF + 2);

    number ans[N];

public:
    void add(number* a, number* b);
    void sub(number* a, number* b);
    void div(number* a, number d);
    void div_into(number* dst, const number* src, number d);
    void mul(number* a, number d);
    void atan(number* a, number m);
    bool has_bit(const number* a);
    void decimal(number* a, number* w);
    void print();
    void calc();
};

// a <- arctan(1 / m)
template<int X>
void PI<X>::atan(number* a, number m)
{
    number b[N] = { 1 };
    number c[N];
    div(b, m);
    std::ranges::copy(std::span(b, N), a);
    const number m2 = m * m;
    // div_into eliminates the O(N) copy (copy + div -> single pass)
    for (int i = 1; ; i++) {
        div(b, m2);
        div_into(c, b, static_cast<number>(i * 2 + 1));
        if (!has_bit(c)) break;
        if (i & 1) sub(a, c);
        else add(a, c);
    }
}

// a <- a + b
template<int X>
void PI<X>::add(number* a, number* b)
{
    for (int i = 0; i < N; i++) {
        number x = a[i] + b[i];
        if (x & ~BASE) {
            a[i] = x & BASE;
            int j;
            for (j = i - 1; a[j] == BASE; j--) a[j] = 0;
            a[j]++;
        } else a[i] = x;
    }
}

// a <- a - b  (a > b)
template<int X>
void PI<X>::sub(number* a, number* b)
{
    for (int i = 0; i < N; i++) {
        if (a[i] < b[i]) {
            a[i] = (BASE + 1) + a[i] - b[i];
            int j;
            for (j = i - 1; a[j] == 0; j--) a[j] = BASE;
            a[j]--;
        } else a[i] -= b[i];
    }
}

// a <- a / d  (d <= BASE)
template<int X>
void PI<X>::div(number* a, number d)
{
    number res = 0;
    for (int i = 0; i < N; i++) {
        res <<= (FIG * 4);
        number x = a[i] + res;
        number q = x / d;
        a[i] = q;
        res = x - q * d;
    }
}

// dst <- src / d  (d <= BASE)
// Combines copy + div into a single pass, avoiding O(N) intermediate copy
template<int X>
void PI<X>::div_into(number* dst, const number* src, number d)
{
    number res = 0;
    for (int i = 0; i < N; i++) {
        res <<= (FIG * 4);
        number x = src[i] + res;
        number q = x / d;
        dst[i] = q;
        res = x - q * d;
    }
}

// a <- a * d  (d <= BASE)
template<int X>
void PI<X>::mul(number* a, number d)
{
    number q = 0;
    for (int i = N - 1; i >= 0; i--) {
        number x = a[i] * d + q;
        a[i] = x & BASE;
        q = x >> (FIG * 4);
    }
}

template<int X>
bool PI<X>::has_bit(const number* a)
{
    // Scan from back: small values have nonzero bits at high indices (low significance),
    // so back-to-front finds them faster as the series converges
    for (int i = N - 1; i >= 0; i--) if (a[i]) return true;
    return false;
}

// hex to decimal number
template<int X>
void PI<X>::decimal(number* a, number* w)
{
    number b[N];
    std::ranges::copy(std::span(a, N), b);
    w[0] = b[0];
    b[0] = 0;
    for (int i = 1; i < M; i++) {
        mul(b, NN);
        w[i] = b[0];
        b[0] = 0;
    }
}

// print answer
template<int X>
void PI<X>::print()
{
    number w[M];
    decimal(ans, w);
    // Buffer all output to avoid per-element std::print overhead
    std::string output;
    output.reserve(M * (FIG + 1) + 20);
    output += std::format("{:>4}.", w[0]);
    for (int i = 1; i < M; i++) {
        output += std::format("{:08}", w[i]);
        if (i % 6 == 0) output += "\n     ";
        else output += ' ';
    }
    output += '\n';
    std::print("{}", output);
}

// calculate PI
// Machin's formula: 4 * arctan(1/5) - arctan(1/239) = PI / 4
template<int X>
void PI<X>::calc()
{
    number x[N];
    atan(ans, 5);
    mul(ans, 4);
    atan(x, 239);
    sub(ans, x);
    mul(ans, 4);
}

int main()
{
    PI<10000> pi;

    const auto start = std::chrono::steady_clock::now();
    pi.calc();
    const auto end = std::chrono::steady_clock::now();
    pi.print();

    const std::chrono::duration<double> elapsed = end - start;
    std::print(stderr, "Time: {}\n", elapsed.count());
}