src/math/GaloisLFSR.hpp

   1 #ifndef BLOBS_MATH_GALOISLFSR_HPP_
   2 #define BLOBS_MATH_GALOISLFSR_HPP_
   3
   4 #include <cassert>
   5 #include <cstdint>
   6 #include <limits>
   7
   8
   9 namespace blobs {
  10 namespace math {
  11
  12 class GaloisLFSR {
  13
  14 public:
  15         // seed should be non-zero
  16         explicit GaloisLFSR(std::uint64_t seed) noexcept
  17         : state(seed) {
  18                 if (state == 0) {
  19                         state = 1;
  20                 }
  21         }
  22
  23         // get the next bit
  24         bool operator ()() noexcept {
  25                 bool result = state & 1;
  26                 state >>= 1;
  27                 if (result) {
  28                         state |= 0x8000000000000000;
  29                         state ^= mask;
  30                 } else {
  31                         state &= 0x7FFFFFFFFFFFFFFF;
  32                 }
  33                 return result;
  34         }
  35
  36         template<class T>
  37         T operator ()(T &out) noexcept {
  38                 constexpr int num_bits =
  39                         std::numeric_limits<T>::digits +
  40                         std::numeric_limits<T>::is_signed;
  41                 for (int i = 0; i < num_bits; ++i) {
  42                         operator ()();
  43                 }
  44                 return out = static_cast<T>(state);
  45         }
  46
  47         /// special case for randrom(boolean), since static_cast<bool>(0b10) == true
  48         bool operator ()(bool &out) noexcept {
  49                 return out = operator ()();
  50         }
  51
  52         template<class T>
  53         T Next() noexcept {
  54                 T next;
  55                 return (*this)(next);
  56         }
  57
  58         double SNorm() noexcept {
  59                 return (2.0 * UNorm()) - 1.0;
  60         }
  61
  62         double UNorm() noexcept {
  63                 return double(Next<std::uint64_t>()) / double(std::numeric_limits<std::uint64_t>::max());
  64         }
  65
  66         template<class Container>
  67         typename Container::reference From(Container &c) {
  68                 assert(c.size() > 0);
  69                 return c[Next<typename Container::size_type>() % c.size()];
  70         }
  71         template<class Container>
  72         typename Container::const_reference From(const Container &c) {
  73                 assert(c.size() > 0);
  74                 return c[Next<typename Container::size_type>() % c.size()];
  75         }
  76
  77 private:
  78         std::uint64_t state;
  79         // bits 64, 63, 61, and 60 set to 1 (counting from 1 lo to hi)
  80         static constexpr std::uint64_t mask = 0xD800000000000000;
  81
  82 };
  83
  84 }
  85 }
  86
  87 #endif