random/engine/RanluxppRngEngine.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // SPDX-FileCopyrightText: 2020 CERN
0003 // SPDX-FileCopyrightText: 2025 Celeritas contributors
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 /*!
0007  * \file corecel/random/engine/RanluxppRngEngine.hh
0008  *
0009  * Original source:
0010  * https://github.com/apt-sim/AdePT/blob/master/include/AdePT/copcore/Ranluxpp.h
0011  */
0012 //---------------------------------------------------------------------------//
0013 #pragma once
0014
0015 #include "corecel/Assert.hh"
0016 #include "corecel/random/data/RanluxppRngData.hh"
0017 #include "corecel/random/distribution/GenerateCanonical.hh"
0018 #include "corecel/random/distribution/detail/GenerateCanonical32.hh"
0019 #include "corecel/sys/ThreadId.hh"
0020
0021 #include "detail/RanluxppImpl.hh"
0022
0023 namespace celeritas
0024 {
0025 //---------------------------------------------------------------------------//
0026 /*!
0027  * Implements the RANLUX++ random number generator engine with modifications.
0028  *
0029  * This implementation of RANLUX++
0030  * \citep{sibidanov-revisionsubtractwithborrow-2017,
0031  * https://doi.org/10.1016/j.cpc.2017.09.005} is based directly on the
0032  * implementation of
0033  * \citet{hahnfeld-ranlux-2021, https://doi.org/10.1051/epjconf/202125103008}.
0034  * The RANLUX++ algorithm
0035  * is an optimization of the RANLUX generator \citep{james-ranluxfortran-1994,
0036  * https://doi.org/10.1016/0010-4655(94)90233-X}, based on work by Luscher's
0037  * modification \citep{luscher-portablehighquality-1994,
0038  * http://arxiv.org/abs/hep-lat/9309020} of Marsaglia and Zaman's RCARRY
0039  * \citep{james-reviewpseudorandom-1990,
0040  * https://doi.org/10.1016/0010-4655(90)90032-V}. As discussed in the RANLUX
0041  * theory paper, the algorithm is essentially a linear congruential generator
0042  * (LCG) with a huge state.
0043  *
0044  * The underlying RCARRY algorithm used an array of 24 24-bit integer words,
0045  * which with today's large integer sizes can be written as 9 64-bit integers.
0046  * A given state is used to extract 12 samples, and the lower 32
0047  * bits of each is used as entropy.
0048  *
0049  * \todo The decision to discard the high bits rather than the low bits from
0050  * each word is likely undesirable at least for the last number in the state,
0051  * since it is known that LCG integers have highly correlated sequential LSBs.
0052  * Also instead, the class could be adapted to provide 18 32-bit samples
0053  * instead of discarding 16 out of every 48 bits.
0054  */
0055 class RanluxppRngEngine
0056 {
0057   public:
0058     //@{
0059     //! Public types.
0060     using result_type = std::uint32_t;
0061     using Initializer_t = RanluxppInitializer;
0062     using ParamsRef = NativeCRef<RanluxppRngParamsData>;
0063     using StateRef = NativeRef<RanluxppRngStateData>;
0064     using RngStateInitializer_t = RanluxppRngStateInitializer;
0065     //@}
0066
0067   public:
0068     //! Instantiate with optional default seed.
0069     CELER_FUNCTION RanluxppRngEngine(
0070         ParamsRef const& params, StateRef const& state, TrackSlotId tid)
0071         : params_(params)
0072     {
0073         CELER_EXPECT(tid < state.state.size());
0074         state_ = &state.state[tid];
0075     }
0076
0077     //! Lowest value potentially generated.
0078     static CELER_CONSTEXPR_FUNCTION result_type min() { return 0u; }
0079
0080     //! Highest value potentially generated.
0081     static CELER_CONSTEXPR_FUNCTION result_type max()
0082     {
0083         return celeritas::numeric_limits<result_type>::max();
0084     }
0085
0086     // Initialize state with the given seed initializer
0087     inline CELER_FUNCTION RanluxppRngEngine& operator=(
0088         Initializer_t const& init);
0089
0090     // Initialize state with the given state initializer
0091     inline CELER_FUNCTION RanluxppRngEngine& operator=(
0092         RngStateInitializer_t const& state_init);
0093
0094     // Generate a 32-bit random integer.
0095     inline CELER_FUNCTION result_type operator()();
0096
0097     // Advance the state \c count times.
0098     inline CELER_FUNCTION void discard(RanluxppUInt count);
0099
0100     // Initialize a state for a new spawned RNG.
0101     inline CELER_FUNCTION RngStateInitializer_t branch();
0102
0103     // Force the state to produce the next block of random bits
0104     inline CELER_FUNCTION void advance() { this->advance(*state_); }
0105
0106   private:
0107     /// IMPLEMENTATION ///
0108
0109     // Produce the next block of random bits for the given state.
0110     inline CELER_FUNCTION void advance(RanluxppRngState& state);
0111
0112     /// DATA ///
0113     static constexpr int offset_ = 48;
0114     ParamsRef const& params_;
0115     RanluxppRngState* state_;
0116 };
0117
0118 //---------------------------------------------------------------------------//
0119 /*!
0120  * Specialization of GenerateCanonical for RanluxppRngEngine.
0121  */
0122 template<class RealType>
0123 struct GenerateCanonical<RanluxppRngEngine, RealType>
0124 {
0125     //!@{
0126     //! \name Type aliases
0127     using real_type = RealType;
0128     using result_type = RealType;
0129     //!@}
0130
0131     //! Declare that we use the 32-bit canonical generator
0132     static constexpr auto policy = GenerateCanonicalPolicy::builtin32;
0133
0134     //! Sample a random number on [0, 1)
0135     CELER_FORCEINLINE_FUNCTION result_type operator()(RanluxppRngEngine& rng)
0136     {
0137         return detail::GenerateCanonical32<RealType>()(rng);
0138     }
0139 };
0140
0141 //---------------------------------------------------------------------------//
0142 // INLINE FUNCTIONS
0143 //---------------------------------------------------------------------------//
0144 /*!
0145  * Initialize state for the given seed and subsequence.
0146  */
0147 inline CELER_FUNCTION RanluxppRngEngine& RanluxppRngEngine::operator=(
0148     Initializer_t const& init)
0149 {
0150     // Skip forward (2^96) * (init.seed + init.subsequence) states
0151     RanluxppArray9 new_a_seed = celeritas::detail::compute_power_modulus(
0152         params_.advance_sequence, init.seed + init.subsequence);
0153
0154     // Convert to Ranluxpp number and save state
0155     state_->value = celeritas::detail::to_ranlux(new_a_seed);
0156     state_->position = 0;
0157
0158     // Skip forward another init.offset samples
0159     if (init.offset > 0)
0160     {
0161         this->discard(init.offset);
0162     }
0163     return *this;
0164 }
0165
0166 //---------------------------------------------------------------------------//
0167 /*!
0168  * Initialize state for the given state initializer.
0169  */
0170 inline CELER_FUNCTION RanluxppRngEngine& RanluxppRngEngine::operator=(
0171     RngStateInitializer_t const& state_init)
0172 {
0173     (*state_).value = state_init.value;
0174     return *this;
0175 }
0176
0177 //---------------------------------------------------------------------------//
0178 /*!
0179  * Skip `n` random numbers without generating them.
0180  */
0181 CELER_FUNCTION void RanluxppRngEngine::discard(RanluxppUInt n)
0182 {
0183     CELER_EXPECT(n > 0);
0184     CELER_ASSERT(params_.max_position > 0);
0185
0186     int left = (params_.max_position - state_->position) / offset_;
0187     CELER_ASSERT(left >= 0);
0188     if (n < static_cast<RanluxppUInt>(left))
0189     {
0190         // Just skip the next few entries in the currently
0191         // available bits.
0192         state_->position += n * offset_;
0193         CELER_ASSERT(state_->position <= params_.max_position);
0194         return;
0195     }
0196
0197     n -= left;
0198     // Need to advance and possibly skip multiple blocks (each block is 576
0199     // random bits, or 12 48-bit samples)
0200     constexpr int n_per_state = ParamsRef::max_position / offset_;
0201     int skip = n / n_per_state;
0202
0203     RanluxppArray9 a_skip = celeritas::detail::compute_power_modulus(
0204         params_.advance_state, skip + 1);
0205
0206     RanluxppArray9 lcg = celeritas::detail::to_lcg(state_->value);
0207     lcg = celeritas::detail::compute_mod_multiply(a_skip, lcg);
0208     state_->value = celeritas::detail::to_ranlux(lcg);
0209
0210     // Potentially skip numbers in the freshly generated block.
0211     int remaining = n - skip * n_per_state;
0212     CELER_ASSERT(remaining >= 0);
0213     state_->position = remaining * offset_;
0214     CELER_ENSURE(state_->position <= params_.max_position);
0215 }
0216
0217 //---------------------------------------------------------------------------//
0218 /*!
0219  * Return the next random bits, generating a new block if necessary.
0220  */
0221 CELER_FUNCTION auto RanluxppRngEngine::operator()() -> result_type
0222 {
0223     if (state_->position + offset_ > params_.max_position)
0224     {
0225         this->advance(*state_);
0226     }
0227
0228     // Extract the Nth word from the state
0229     int idx = state_->position / 64;
0230     int offset = state_->position % 64;
0231     int num_bits = 64 - offset;
0232
0233     RanluxppUInt bits = state_->value.number[idx] >> offset;
0234     if (num_bits < offset_)
0235     {
0236         bits |= state_->value.number[idx + 1] << num_bits;
0237     }
0238     bits &= ((RanluxppUInt(1) << offset_) - 1);
0239
0240     state_->position += offset_;
0241
0242     CELER_ENSURE(state_->position <= params_.max_position);
0243     return bits & 0xffffffffu;
0244 }
0245
0246 //---------------------------------------------------------------------------//
0247 /*!
0248  * Initialize a state for a new spawned RNG.
0249  *
0250  * \par Branching
0251  * Branching is performed in two steps.  First, the state of the new RNG
0252  * (\f$x^{\prime}\f$) is initialized as
0253  * \f[
0254  *      x^{i,\prime}_j = x^i_j ^ x^{i+1}_j \, .
0255  * \f]
0256  * Second, to decorrelate the new RNG from this RNG, the new RNG is
0257  * advanced forward to the next block
0258  */
0259 CELER_FUNCTION RanluxppRngEngine::RngStateInitializer_t
0260 RanluxppRngEngine::branch()
0261 {
0262     // Create a new state
0263     RanluxppRngState new_state = *state_;
0264
0265     // Advance the RNG
0266     this->advance(*state_);
0267
0268     // XOR the new state with the updated and advanced state
0269     for (auto i : celeritas::range(9))
0270     {
0271         new_state.value.number[i] ^= (*state_).value.number[i];
0272     }
0273
0274     // Advance the new rng to decorrelate it
0275     this->advance(new_state);
0276
0277     // Create and return the state initializer
0278     RngStateInitializer_t initializer{new_state.value};
0279
0280     return initializer;
0281 }
0282
0283 //---------------------------------------------------------------------------//
0284 /*!
0285  * Advance to the next state (block of random bits).
0286  */
0287 CELER_FUNCTION void RanluxppRngEngine::advance(RanluxppRngState& state)
0288 {
0289     RanluxppArray9 lcg = celeritas::detail::to_lcg(state.value);
0290     lcg = celeritas::detail::compute_mod_multiply(params_.advance_state, lcg);
0291     state.value = celeritas::detail::to_ranlux(lcg);
0292     state.position = 0;
0293 }
0294
0295 //---------------------------------------------------------------------------//
0296 }  // end namespace celeritas