include/opencascade/MathUtils_Random.hxx

0001 // Copyright (c) 2025 OPEN CASCADE SAS
0002 //
0003 // This file is part of Open CASCADE Technology software library.
0004 //
0005 // This library is free software; you can redistribute it and/or modify it under
0006 // the terms of the GNU Lesser General Public License version 2.1 as published
0007 // by the Free Software Foundation, with special exception defined in the file
0008 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0009 // distribution for complete text of the license and disclaimer of any warranty.
0010 //
0011 // Alternatively, this file may be used under the terms of Open CASCADE
0012 // commercial license or contractual agreement.
0013
0014 #ifndef _MathUtils_Random_HeaderFile
0015 #define _MathUtils_Random_HeaderFile
0016
0017 #include <cstdint>
0018
0019 namespace MathUtils
0020 {
0021
0022 //! High-quality pseudo-random number generator based on xoshiro256**.
0023 //!
0024 //! xoshiro256** is a general-purpose PRNG designed by David Blackman
0025 //! and Sebastiano Vigna. It has:
0026 //! - 256-bit state (period 2^256 - 1)
0027 //! - Passes all BigCrush statistical tests
0028 //! - Very fast on 64-bit hardware
0029 //! - Equidistributed to 4 dimensions
0030 //!
0031 //! Suitable for Monte Carlo methods, stochastic optimization,
0032 //! and any application requiring high-quality randomness.
0033 class RandomGenerator
0034 {
0035 public:
0036   //! Initialize with a seed value.
0037   //! Uses SplitMix64 to expand a single seed into the full 256-bit state,
0038   //! ensuring good initialization even from poor seeds.
0039   //! @param theSeed seed value (default 1)
0040   explicit RandomGenerator(const uint64_t theSeed = 1) { SetSeed(theSeed); }
0041
0042   //! Re-seed the generator.
0043   //! @param theSeed seed value
0044   void SetSeed(const uint64_t theSeed)
0045   {
0046     // Use SplitMix64 to initialize the state from a single seed.
0047     // This ensures good state initialization even from sequential seeds.
0048     uint64_t aSeedState = theSeed;
0049     myState[0]          = splitmix64(aSeedState);
0050     myState[1]          = splitmix64(aSeedState);
0051     myState[2]          = splitmix64(aSeedState);
0052     myState[3]          = splitmix64(aSeedState);
0053   }
0054
0055   //! Generate next 64-bit unsigned integer.
0056   //! @return pseudo-random value in [0, 2^64)
0057   uint64_t NextInt()
0058   {
0059     const uint64_t aResult = rotl(myState[1] * 5, 7) * 9;
0060     const uint64_t aT      = myState[1] << 17;
0061
0062     myState[2] ^= myState[0];
0063     myState[3] ^= myState[1];
0064     myState[1] ^= myState[2];
0065     myState[0] ^= myState[3];
0066
0067     myState[2] ^= aT;
0068     myState[3] = rotl(myState[3], 45);
0069
0070     return aResult;
0071   }
0072
0073   //! Generate next double in [0, 1).
0074   //! Uses 53 bits of randomness for full double precision.
0075   //! @return pseudo-random value in [0, 1)
0076   double NextReal() { return static_cast<double>(NextInt() >> 11) * (1.0 / 9007199254740992.0); }
0077
0078 private:
0079   //! 64-bit left rotation.
0080   static constexpr uint64_t rotl(const uint64_t theX, const int theK)
0081   {
0082     return (theX << theK) | (theX >> (64 - theK));
0083   }
0084
0085   //! SplitMix64 step for seed expansion.
0086   //! Takes a reference to the seed state and advances it.
0087   static uint64_t splitmix64(uint64_t& theState)
0088   {
0089     theState += 0x9e3779b97f4a7c15ULL;
0090     uint64_t aZ = theState;
0091     aZ          = (aZ ^ (aZ >> 30)) * 0xbf58476d1ce4e5b9ULL;
0092     aZ          = (aZ ^ (aZ >> 27)) * 0x94d049bb133111ebULL;
0093     return aZ ^ (aZ >> 31);
0094   }
0095
0096   uint64_t myState[4]; //!< 256-bit generator state
0097 };
0098
0099 } // namespace MathUtils
0100
0101 #endif // _MathUtils_Random_HeaderFile