Twitterのモンテカルロ法のC++版の速度比較コード（線形合同法自前実装）

mc.cpp

#include <cstdint> // for std::int32_t
#include <iomanip> // for std::setprecision
#include <ios> // for std::ios::fixed, std::ios::floatfield
#include <iostream> // for std::cout, std::endl
#include <utility> // for std::make_pair, std::pair

#define USE_MOD 1

namespace {
inline std::pair<std::int32_t, double> LCGs(std::int32_t seed);
inline double mcpi();
}

int main()
{
    std::cout.setf(std::ios::fixed, std::ios::floatfield);
    std::cout << "pi = "
              << std::setprecision(16)
              << mcpi()
              << std::endl;
}

namespace {
// Linear congruential generators (LCGs)
// Parameters are provided by Park and Miller
// See https://c-faq.com/lib/rand.html
std::pair<std::int32_t, double> LCGs(std::int32_t seed)
{
    auto const a = 48271;
    auto const m = 2147483647;
    auto const q = m / a;
#if USE_MOD
    auto const r = m % a;
#else
    auto const r = m - q * a;
#endif

    auto const hi = seed / q;
#if USE_MOD
    auto const lo = seed % q;
#else
    auto const lo = seed - hi * q;
#endif
    auto const test = a * lo - r * hi;

    if (test > 0) {
        seed = test;
    } else {
        seed = test + m;
    }

    return std::make_pair(seed, static_cast<double>(seed) / static_cast<double>(m));
}

double mcpi()
{
    auto seed = 20231226;
    auto constexpr num_points = 1000000000;
    auto num_inside = 0;

    for (auto i = 0; i < num_points; i++) {
        auto const [seedtmp, x] = LCGs(seed);
        auto const [seedtmp2, y] = LCGs(seedtmp);

        seed = seedtmp2;

        auto const r2 = x * x + y * y;
        if (r2 < 1.0) {
            num_inside++;
        }
    }

    return 4.0 * static_cast<double>(num_inside) / static_cast<double>(num_points);
}
}