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

mc_dsfmt.cpp

#include <array>    // for std::array
#include <iomanip>  // for std::setprecision
#include <ios>      // for std::ios::fixed, std::ios::floatfield
#include <iostream> // for std::cout, std::endl
#define HAVE_SSE2
#define HAVE_AVX2
#include "dSFMT.h"
#include "dSFMT-2203-avx512.h"

namespace {
static auto constexpr ARRAYSIZE = 64;
static auto constexpr RANDSIZE = ARRAYSIZE / 2;
inline double mcpi();
}

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

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

    dsfmt_t dsfmt;
    dsfmt_init_gen_rand(&dsfmt, seed);

    alignas(64) std::array<double, ARRAYSIZE> randarray;
    auto const loopnum = num_points / (RANDSIZE / 2);

    for (auto i = 0; i < loopnum; i++) {
        dsfmt_fill_array_close_open(&dsfmt, randarray.data(), RANDSIZE);

        for (auto j = 0; j < RANDSIZE; j += 2) {
            auto const x = randarray[j];
            auto const y = randarray[j + 1];
        
            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);
}
}