Middle Square Weyl Sequence PRNG in C++14

msws.hpp

// PRNG: Middle Square Weyl Sequence
//
// This is a rewrite of mswsrngv3 from https://mswsrng.wixsite.com/rand in
// C++14. This file is licensed under MIT license.
//
// Use msws::Msws as a UniformRandomBitGenerator and RandomNumberEngine.
// #define MSWS_NOIOS to remove formatted I/O support with C++ input/output
// streams.
//
// MIT License
//
// Copyright (c) 2020 EAirPeter
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

#pragma once

#ifndef MSWS_HPP_
#define MSWS_HPP_

#include <cstdint>

#ifndef MSWS_NOIOS
#include <iosfwd>
#endif // ifndef MSWS_NOIOS

namespace msws {
using ::std::uint32_t;
using ::std::uint64_t;

namespace detail {
  constexpr uint64_t sconst[]{
      0x37e1c9b5e1a2b843, 0x56e9d7a3d6234c87, 0xc361be549a24e8c7,
      0xd25b9768a1582d7b, 0x18b2547d3de29b67, 0xc1752836875c29ad,
      0x4e85ba61e814cd25, 0x17489dc6729386c1, 0x7c1563ad89c2a65d,
      0xcdb798e4ed82c675, 0xd98b72e4b4e682c1, 0xdacb7524e4b3927d,
      0x53a8e9d7d1b5c827, 0xe28459db142e98a7, 0x72c1b3461e4569db,
      0x1864e2d745e3b169, 0x6a2c143bdec97213, 0xb5e1d923d741a985,
      0xb4875e967bc63d19, 0x92b64d5a82db4697, 0x7cae812d896eb1a5,
      0xb53827d41769542d, 0x6d89b42c68a31db5, 0x75e26d434e2986d5,
      0x7c82643d293cb865, 0x64c3bd82e8637a95, 0x2895c34d9dc83e61,
      0xa7d58c34dea35721, 0x3dbc5e687c8e61d5, 0xb468a235e6d2b193,
  };
} // namespace detail

// Satisfies:
//  * RandomNumberEngine
//  * UniformRandomBitGenerator
class Msws {
public:
  using result_type = uint32_t;

  static constexpr uint32_t min() noexcept { return 0; }

  static constexpr uint32_t max() noexcept { return UINT32_MAX; }

private:
  struct InternalInit {};

  explicit constexpr Msws(uint64_t x, uint64_t w, uint64_t s) noexcept
      : x(x), w(w), s(s) {}

  explicit constexpr Msws(InternalInit, uint64_t init) noexcept
      : x(init), w(init), s(init) {}

  static constexpr Msws InitLocalMsws(uint64_t n) noexcept {
    uint64_t r = n / 100000000;
    uint64_t t = n % 100000000;
    uint64_t s = detail::sconst[r % 30];
    r /= 30;
    uint64_t xw = t * s + r * s * 100000000;
    return Msws(xw, xw, s);
  }

  static constexpr uint64_t InitInternalState(uint64_t n) noexcept {
    Msws msws = InitLocalMsws(n);

    uint64_t u = (msws() & 0x7) << 1 | 1;
    uint32_t v = uint32_t(1) << u;

    uint32_t c = 0;
    for (int m = 60; m > 0;) {
      uint32_t j = msws();
      for (int i = 0; i < 32; i += 4) {
        uint32_t k = (j >> i) & 0xf;
        if (k && !(c & (uint32_t(1) << k))) {
          c |= uint32_t(1) << k;
          u |= uint64_t(k) << m;
          m -= 4;
          if (m == 24 || m == 28)
            c = uint32_t(1) << k | v;
          if (!m)
            break;
        }
      }
    }
    return u;
  }

  template<class SeedSeq> static uint64_t GetSeedFromSeq(SeedSeq& seq) {
    typename SeedSeq::result_type seed[1];
    seq.generate(seed, seed + 1);
    return uint64_t(seed[0]);
  }

public:
  constexpr Msws() noexcept = default;
  constexpr Msws(const Msws&) noexcept = default;
  constexpr Msws(Msws&&) noexcept = default;

  constexpr Msws& operator=(const Msws&) noexcept = default;
  constexpr Msws& operator=(Msws&&) noexcept = default;

  ~Msws() noexcept = default;

  explicit constexpr Msws(uint64_t seed) noexcept
      : Msws(InternalInit{}, InitInternalState(seed)) {}

  template<class SeedSeq>
  explicit Msws(SeedSeq& seq) : Msws(GetSeedFromSeq(seq)) {}

  constexpr void seed() noexcept { *this = Msws(); }
  constexpr void seed(uint64_t seed) noexcept { *this = Msws(seed); }

  template<class SeedSeq> constexpr void seed(SeedSeq& seq) {
    seed(GetSeedFromSeq(seq));
  }

  constexpr uint32_t operator()() noexcept {
    x *= x;
    x += w += s;
    return uint32_t(x = (x >> 32) | (x << 32));
  }

  constexpr void discard(unsigned long long nToSkip) noexcept {
    while (nToSkip--)
      (void)(*this)();
  }

  constexpr bool operator==(const Msws& rhs) const noexcept {
    return x == rhs.x && w == rhs.w && s == rhs.s;
  }

  constexpr bool operator!=(const Msws& rhs) const noexcept {
    return !(*this == rhs);
  }

#ifndef MSWS_NOIOS
  template<class CharT, class Traits>
  friend ::std::basic_istream<CharT, Traits>&
  operator>>(::std::basic_istream<CharT, Traits>& is, Msws& m) {
    return is >> m.x >> m.w >> m.s;
  }

  template<class CharT, class Traits>
  friend ::std::basic_ostream<CharT, Traits>&
  operator<<(::std::basic_ostream<CharT, Traits>& os, const Msws& m) {
    return os << m.x << ' ' << m.w << ' ' << m.s;
  }
#endif // ifndef MSWS_NOIOS

private:
  uint64_t x = 0, w = 0, s = 0xb5ad4eceda1ce2a9;
};

} // namespace msws

#endif // ifndef MSWS_HPP_