Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "MooseRandomPerturbation.h"
      11             : 
      12         138 : MooseRandomPerturbation::MooseRandomPerturbation(uint64_t seed, unsigned int n, unsigned int rounds)
      13         138 :   : _k0(static_cast<uint32_t>(seed)),
      14         138 :     _k1(static_cast<uint32_t>(seed >> 32)),
      15         138 :     _n(n),
      16         138 :     _half_bits((ceilLog2(n) + 1) / 2),
      17         138 :     _half_mask((uint32_t(1) << _half_bits) - 1),
      18         138 :     _rounds(rounds)
      19             : {
      20             :   mooseAssert(_n > 0, "n must be > 0");
      21             :   mooseAssert(_rounds > 0, "rounds must be greater than 0.");
      22         138 : }
      23             : 
      24             : uint32_t
      25       52462 : MooseRandomPerturbation::permute(uint32_t x) const
      26             : {
      27             :   mooseAssert(x < _n, "x must be < n");
      28             : 
      29       52462 :   uint32_t y = x;
      30             :   do
      31             :   {
      32       76732 :     y = permutePadded(y);
      33       76732 :   } while (y >= _n);
      34             : 
      35       52462 :   return y;
      36             : }
      37             : 
      38             : uint32_t
      39       10162 : MooseRandomPerturbation::invert(uint32_t y) const
      40             : {
      41             :   mooseAssert(y < _n, "y must be < n");
      42             : 
      43       10162 :   uint32_t x = y;
      44             :   do
      45             :   {
      46       18134 :     x = invertPadded(x);
      47       18134 :   } while (x >= _n);
      48             : 
      49       10162 :   return x;
      50             : }
      51             : 
      52             : uint32_t
      53       76732 : MooseRandomPerturbation::permutePadded(uint32_t x) const
      54             : {
      55             :   // Split x into the upper (L) and lower (R) half_bits-wide halves.
      56       76732 :   uint32_t L = x >> _half_bits;
      57       76732 :   uint32_t R = x & _half_mask;
      58             : 
      59      690588 :   for (unsigned int r = 0; r < _rounds; ++r)
      60             :   {
      61             :     // Standard balanced Feistel step: new_L = R, new_R = L XOR F(R).
      62      613856 :     const uint32_t F = roundFunction(R, r);
      63      613856 :     const uint32_t newL = R;
      64      613856 :     const uint32_t newR = (L ^ F) & _half_mask;
      65      613856 :     L = newL;
      66      613856 :     R = newR;
      67             :   }
      68             : 
      69       76732 :   return ((L & _half_mask) << _half_bits) | (R & _half_mask);
      70             : }
      71             : 
      72             : uint32_t
      73       18134 : MooseRandomPerturbation::invertPadded(uint32_t y) const
      74             : {
      75       18134 :   uint32_t L = (y >> _half_bits) & _half_mask;
      76       18134 :   uint32_t R = y & _half_mask;
      77             : 
      78             :   // Run rounds in reverse; recover prev_R = L then prev_L = R XOR F(prev_R).
      79      163206 :   for (int r = static_cast<int>(_rounds) - 1; r >= 0; --r)
      80             :   {
      81      145072 :     const uint32_t prevR = L;
      82      145072 :     const uint32_t F = roundFunction(prevR, static_cast<unsigned int>(r));
      83      145072 :     const uint32_t prevL = (R ^ F) & _half_mask;
      84      145072 :     L = prevL;
      85      145072 :     R = prevR;
      86             :   }
      87             : 
      88       18134 :   return ((L & _half_mask) << _half_bits) | (R & _half_mask);
      89             : }
      90             : 
      91             : uint32_t
      92      758928 : MooseRandomPerturbation::roundFunction(uint32_t half, unsigned int round) const
      93             : {
      94      758928 :   uint32_t x = half;
      95      758928 :   x ^= _k0;
      96      758928 :   x += 0x9e3779b9U * static_cast<uint32_t>(round + 1); // round-dependent constant
      97      758928 :   x ^= _k1;
      98      758928 :   x = mix32(x);
      99      758928 :   return x & _half_mask;
     100             : }
     101             : 
     102             : unsigned int
     103         138 : MooseRandomPerturbation::ceilLog2(uint32_t n)
     104             : {
     105             :   mooseAssert(n > 0, "n must be > 0");
     106             : 
     107         138 :   unsigned int bits = 0;
     108         138 :   uint32_t v = n - 1;
     109        1758 :   while (v > 0)
     110             :   {
     111        1620 :     ++bits;
     112        1620 :     v >>= 1;
     113             :   }
     114         138 :   return bits;
     115             : }
     116             : 
     117             : uint32_t
     118      758928 : MooseRandomPerturbation::mix32(uint32_t x)
     119             : {
     120      758928 :   x ^= x >> 16;
     121      758928 :   x *= 0x7feb352dU;
     122      758928 :   x ^= x >> 15;
     123      758928 :   x *= 0x846ca68bU;
     124      758928 :   x ^= x >> 16;
     125      758928 :   return x;
     126             : }