random: adds new random generator

Author: quesnel

lib/include/irritator/random.hpp

@@ -0,0 +1,149 @@
+// Copyright (c) 2025 INRAE Distributed under the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef ORG_VLEPROJECT_IRRITATOR_RANDOM_HPP
+#define ORG_VLEPROJECT_IRRITATOR_RANDOM_HPP
+
+#include <array>
+#include <bit>
+#include <concepts>
+#include <cstdint>
+#include <limits>
+
+#ifdef _MSC_VER
+#include <intrin.h>
+#endif
+
+/**
+ * A Pseudo Random Generator based on Salmon et al. (Random123).
+ *
+ * The interface of @c philox_64 in fully compatible with the standard random
+ * header.
+ */
+class philox_64
+{
+public:
+    using result_type = uint64_t;
+
+    static constexpr uint64_t PHILOX_M0 = 0xD2B74407B1CADAC9ULL;
+    static constexpr uint64_t PHILOX_W0 = 0x9E3779B97F4A7C15ULL;
+    static constexpr int      ROUNDS    = 10;
+
+    /**
+     * @brief Ctor
+     * @param seed Global simulation @c seed.
+     * @param index Unique identifiant for model
+     * @param step First step
+     */
+    constexpr Philox64(const uint64_t seed,
+                       const uint64_t index,
+                       const uint64_t step = 0) noexcept
+      : m_key({ seed })
+      , m_counter({ index, step })
+      , m_buffer_pos(2) // To force buffer computation
+    {}
+
+    /** Build the next random number. Defined by the STL requirements for random
+     * number generators. */
+    constexpr result_type operator()() noexcept
+    {
+        if (m_buffer_pos >= 2) {
+            refill_buffer();
+        }
+        returnm_buffer[m_buffer_pos++];
+    }
+
+    /** Discard @c z random numbers. Useful to jump ahead in the random
+     * sequence. */
+    constexpr void discard(const uint64_t z) noexcept
+    {
+        m_counter[1] += z;
+        m_buffer_pos = 2; // Invalid the buffer.
+    }
+
+    /** A state setter to give full random access for example to replay
+     * simulation for example @c index is the @c model_id and @c step is the
+     * number of the random.  */
+    constexpr void set_state(const uint64_t index, const uint64_t step) noexcept
+    {
+        m_counter[0] = index;
+        m_counter[1] = step;
+        m_buffer_pos = 2; // Invalid the buffer.
+    }
+
+    /** Give the minimum integer of the random generator. Defined by the STL
+     * requirements for random number generators. */
+    static constexpr result_type min() noexcept
+    {
+        return std::numeric_limits<result_type>::min();
+    }
+
+    /** Give the maximun integer of the random generator. Defined by the STL
+     * requirements for random number generators. */
+    static constexpr result_type max() noexcept
+    {
+        return std::numeric_limits<result_type>::max();
+    }
+
+private:
+    std::array<uint64_t, 1> m_key; /**< Global seed. */
+    std::array<uint64_t, 2>
+      m_counter; /**< m_counter[0] the model_id and m_counter[1] the step. */
+    std::array<uint64_t, 2>
+        m_buffer; /**< Buffer to store result and next result. */
+    int m_buffer_pos;
+
+    // A 128-bit multiplication (high and low parts)
+    static constexpr void mulhilo(uint64_t  a,
+                                  uint64_t  b,
+                                  uint64_t& lo,
+                                  uint64_t& hi) noexcept
+    {
+#if defined(__SIZEOF_INT128__) && !defined(__wasm__)
+        __uint128_t product =
+          static_cast<__uint128_t>(a) * static_cast<__uint128_t>(b);
+        lo = static_cast<uint64_t>(product);
+        hi = static_cast<uint64_t>(product >> 64);
+#elif defined(_MSC_VER) && defined(_WIN64)
+        lo = _umul128(a, b, &hi);
+#else
+        /* First calculate all of the cross products. */
+        uint64_t lo_lo = (lhs & 0xffffffff) * (rhs & 0xffffffff);
+        uint64_t hi_lo = (lhs >> 32) * (rhs & 0xffffffff);
+        uint64_t lo_hi = (lhs & 0xffffffff) * (rhs >> 32);
+        uint64_t hi_hi = (lhs >> 32) * (rhs >> 32);
+
+        /* Now add the products together. These will never overflow. */
+        uint64_t cross = (lo_lo >> 32) + (hi_lo & 0xffffffff) + lo_hi;
+        uint64_t upper = (hi_lo >> 32) + (cross >> 32) + hi_hi;
+
+        hi = upper;
+        lo = (cross << 32) | (lo_lo & 0xffffffff);
+#endif
+    }
+
+    constexpr void refill_buffer() noexcept
+    {
+        uint64_t ctr0 = m_counter[0];
+        uint64_t ctr1 = m_counter[1];
+        uint64_t key0 = m_key[0];
+
+        for (int i = 0; i < ROUNDS; ++i) {
+            uint64_t hi, lo;
+            mulhilo(PHILOX_M0, ctr0, lo, hi);
+
+            ctr0 = hi ^ key0 ^ ctr1;
+            ctr1 = lo;
+
+            key0 += PHILOX_W0;
+        }
+
+        m_buffer[0] = ctr0;
+        m_buffer[1] = ctr1;
+
+        // Finaly increase the counter and reset the buffer position.
+        m_counter[1]++;
+        m_buffer_pos = 0;
+    }
+};