Twitterのモンテカルロ法の速度比較のコード（C++版）

mc.cpp

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

#define USE_MOD 0

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

int main()
{
    auto seed = 20231226;
    auto const 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++;
        }
    }

    std::cout.setf(std::ios::fixed, std::ios::floatfield);
    std::cout << "pi = "
              << std::setprecision(16)
              << 4.0 * static_cast<double>(num_inside) / static_cast<double>(num_points)
              << 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));
}
}