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

mc_xorshift.cpp

#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 __AVX2_AVAILABLE__
#include "SIMDInstructionSet.h"
#include "Xoshiro256Plus.h"

namespace {
inline double mcpi();
}

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

namespace {
double mcpi()
{
    using Xoshiro256PlusAVX2 = SEFUtility::RNG::Xoshiro256Plus<SIMDInstructionSet::AVX2>;

    auto constexpr seed = 20231226;
    auto constexpr num_points = 1000000000;
    auto num_inside = 0;
    Xoshiro256PlusAVX2 avx_rng(seed);

    auto const loopnum = num_points / 2;
    for (auto i = 0; i < loopnum; i++) {
        auto const next4_avx = avx_rng.dnext4();
        auto const x = next4_avx[0];
        auto const y = next4_avx[1];
        
        auto const r2 = x * x + y * y;
        if (r2 < 1.0) {
            num_inside++;
        }

        auto const x_2 = next4_avx[2];
        auto const y_2 = next4_avx[3];
        
        auto const r2_2 = x_2 * x_2 + y_2 * y_2;
        if (r2_2 < 1.0) {
            num_inside++;
        }
    }

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