円周率を正確に求めるプログラム

simple-pi.cpp

/*
 * 円周率を正確に求めるプログラム
 *
 * - 半径が10^nの円の第一象限を縦横に10^n分割してその面積の下限と上限を求めます
 */

#include <cstdint>
#include <cmath>
#include <iostream>

using INT = std::int32_t;
constexpr INT N = 4;
constexpr INT M = std::pow(10, N);
constexpr INT MM = M * M;

INT calc_lower() {
    // 下限を求める。各分割の左上が円に含まれるならば積分する。
    INT sum = 0;
    for (INT i = 1; i <= M; ++i)
        for (INT j = 1; j <= M; ++j)
            if (i * i + j * j <= MM)
                sum += 1;
    return sum * 4;
}

INT calc_upper() {
    // 上限を求める。各分割の右下が円に含まれないならば積分する。
    INT sum = 0;
    for (INT i = 0; i < M; ++i)
        for (INT j = 0; j < M; ++j)
            if (i * i + j * j > MM)
                sum += 1;
    return (MM - sum) * 4;
}

int main() {
    std::cout << "n = " << N << std::endl;
    std::cout << "lower = " << calc_lower() << " / " << MM << std::endl;
    std::cout << "upper = " << calc_upper() << " / " << MM << std::endl;
    return 0;
}