8

My task is to modify Sergiu Dotenco's Well Equidistributed Long-period Linear (WELL) algorithm code to not use boost (not saying boost is bad, but due to some company's policy i have to remove it).

now, Sergiu's WELL is using boost's mpl library, there are quite some logic behind it. So one way is to read up all those, then naturally i would be able to finish the task. The other way is, replacing bit by bit with some best guess.

I'm on the 2nd way to hope this try-and-error approach would be faster. So far I've successfully replaced boost::mpl::if_ and if_c with std::conditional, but hit error when try to update IsPowerOfTwo and Power2Modulo etc, that's why i'm seeking help there.

Below is the code, how to rewrite it without boost, but only c++17?

/**
    * Conditional expression of type (r & (r - 1)) == 0 which allows to check
    * whether a number @f$r@f$ is of type @f$2^n@f$.
    */
    typedef boost::mpl::equal_to<
        boost::mpl::bitand_<
        boost::mpl::_,
        boost::mpl::minus<boost::mpl::_, boost::mpl::int_<1>
        >
        >,
        boost::mpl::int_<0>
    > IsPowerOfTwo;

    template<class UIntType, UIntType r>
    struct Power2Modulo
    {
        typedef typename boost::mpl::apply<
            IsPowerOfTwo,
            boost::mpl::integral_c<UIntType, r>
        >::type type;

        BOOST_STATIC_ASSERT(type::value);

        template<class T>
        static T calc(T value)
        {
            return value & (r - 1);
        }
    };

If possible, pls give a short example on how to call it? I tried to instantiate IsPowerOfTwo or Power2Modulo in main with

Detail::IsPowerOfTwo p0;     

or

Detail::Power2Modulo<int, 3> p1;

but got compilation error.

I asked a relevant question before and got some suggestion. However, not familiar to metaprogramming and boost, I don't quite get it.

KARTHIKEYAN.A
  • 18,210
  • 6
  • 124
  • 133
athos
  • 6,120
  • 5
  • 51
  • 95
  • 1
    code review is for ppl to review code written by me right? this question is to seek help rewrite a piece of code from a 3rd party. @JakeFreeman – athos Dec 13 '17 at 02:10
  • 1
    Simplest way to check if the number is power of two is `((n & (n - 1)) == 0)`, while modulo of power of two would be `(n & (n - 1))` – miradham Dec 13 '17 at 02:19
  • @miradham this is a question on compilation time metaprogramming, not runtime coding. – athos Dec 13 '17 at 02:20
  • @KenWhite my task is to modify Sergiu Dotenco's WELL code to not use boost (not saying boost is bad, but due to some company's policy i have to remove it). now, Sergiu's WELL is using boost's mpl library, there are quite some logic behind it. So one way is to read up all those, then naturally i would be able to finish the task. The other way is, replacing bit by bit with some best guess. So far I've successfully replaced `boost::mpl::if_` with `std::conditional`, but hit error when try to update `IsPowerOfTwo` and `Power2Modulo` etc, that's why i'm seeking help there. – athos Dec 13 '17 at 02:29
  • This is classical X/Y problem. If you want to know how to test for power-of-2 do **not** ask about misusing some implementation that you don't want to be using. Ask about _what you want to achieve_, perhaps with your current attempt. (Also note that `int` is not unsigned, so it probably doesn't fit UIntType requirements?) – sehe Dec 13 '17 at 08:36
  • @sehe i actually asked "how to **rewrite without using boost** a piece of code that does A", the question is about *rewrite*, the target is to mimic that and rewrite a whole library. Unfortunately, people keep on telling me how to *do A* in a different way or the question is duplicate as there are questions already on how to *do A*. funny... – athos Dec 13 '17 at 08:50
  • If you want to reimplement the template on any c++11 + capable compiler: https://stackoverflow.com/a/19399478/85371. Mmm. I've understood it now. Al you need is the power-of-2 check, and you have it. – sehe Dec 13 '17 at 08:53
  • @sehe Sergiu's implementation of WELL depends on boost MPL, all these topics have much logic inside. So, I fancy if the target is not to fully understand Sergiu's code, WELL algorithm and boost MPL, but just remove the dependency on MPL, an easier way is find an example on rewrite a small piece, then i carry on for the rest. – athos Dec 13 '17 at 08:57
  • @sehe no jobs is easy. that's why i'm trying to minimum the effort. – athos Dec 13 '17 at 09:07
  • 1
    @sehe by minimum effort i mean, there might be thousands of lines of code in MPL, but Sergiu only make use of around 10 functions of it, in multiple places of course. I already successfully changed 3 of them using STL (if_c, if_, and int_), if the codes above can be also replaced using STL, the job can be done. – athos Dec 13 '17 at 09:14
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/161092/discussion-between-sehe-and-athos). – sehe Dec 13 '17 at 09:15

2 Answers2

17

So, I looked at that library, and created a no-boost fork adapting the WELL pseudo-random-number-generator to pure c++11.

See here on my github: https://github.com/sehe/well-random (the default branch is no-boost).

What is well-random?

well-random is a c++11 fork from random, a collection of various pseudo-random number generators and distributions that were intended to accompany the Boost Random Number Library.

This fork currently only adopted the WELL generator and its tests.

Getting started

The no-boost branch no longer requires any boost library. Instead it requires c++11. To compile the tests make sure first CMake 2.8 is installed, then enter :

$ cmake . -DCMAKE_BUILD_TYPE=Release

in your terminal or command prompt on Windows inside project's directory to generate the appropriate configuration that can be used to compile the tests using make/nmake or inside an IDE.

What Was Refactored

  1. BOOST_STATIC_ASSERT to STATIC_ASSERT (this becomes obsolete with c++17: http://en.cppreference.com/w/cpp/language/static_assert)
  2. BOOST_STATIC_CONSTANT to static constexpr
  3. BOOST_PREVENT_MACRO_SUBSTITUTION -> PREVENT_MACRO_SUBSTITUTION (trivial macro)
  4. BOOST_THROW_EXCEPTION dropped. NOTE This implies the code cannot be compiled with exception support disabled.
  5. All things related to Boost Test

    • BOOST_CHECK -> CHECK

      #define MESSAGE_PREAMBLE() (std::cerr << __FILE__ << ":" << __LINE__ << " ")
      
      
      #define CHECK(test) do { if (!(test)) MESSAGE_PREAMBLE() << #test << "\n"; } while (0)
      
    • BOOST_CHECK_EQUAL -> CHECK_EQUAL

      #define CHECK_EQUAL(expected,actual) do { \
          auto&& _e = expected; \
          auto&& _a = actual; \
          if (_e != _a) \
              MESSAGE_PREAMBLE() << "expected:" << #expected << " = " << _e << "\n" \
                        << "\tactual:" << #actual << " = " << _a << "\n"; \
      } while (0)
      
    • BOOST_AUTO_TEST_CASE - dropped. The test driver is main now:

      int main() {
          //CHECK_EQUAL(16, Detail::shift<2>(64));
          //CHECK_EQUAL(64, Detail::shift<-2>(16));
          //CHECK_EQUAL(32, Detail::shift<0>(32));
          //CHECK(Detail::is_powerof2(512u));
          //CHECK(not Detail::is_powerof2(0u));
      
          WellTestCase<Well512a,   0x2b3fe99e>::run();
          WellTestCase<Well521a,   0xc9878363>::run();
          WellTestCase<Well521b,   0xb75867f6>::run();
          WellTestCase<Well607a,   0x7b5043ea>::run();
          WellTestCase<Well607b,   0xaedee7da>::run();
          WellTestCase<Well800a,   0x2bfe686f>::run();
          WellTestCase<Well800b,   0xf009e1bd>::run();
          WellTestCase<Well1024a,  0xd07f528c>::run();
          WellTestCase<Well1024b,  0x867f7993>::run();
          WellTestCase<Well19937a, 0xb33a2cd5>::run();
          WellTestCase<Well19937b, 0x191de86a>::run();
          WellTestCase<Well19937c, 0x243eaed5>::run();
          WellTestCase<Well21701a, 0x7365a269>::run();
          WellTestCase<Well23209a, 0x807dacb >::run();
          WellTestCase<Well23209b, 0xf1a77751>::run();
          WellTestCase<Well44497a, 0xfdd7c07b>::run();
          WellTestCase<Well44497b, 0x9406547b>::run();
      }
      
  6. boost::ref -> std::ref (from <functional>)

  7. Boost Range helpers replaced by standard c++ (boost::size, boost::end for arrays)

  8. using ulong_long_type = unsigned long long;

  9. Conditional operators shift and mod have been re-implemented with straight-up SFINAE based on std::enable_if instead of using MPL meta-programming:

    template<class UIntType, unsigned N>
    struct Left
    {
        static UIntType shift(UIntType a)
        {
            return a << N;
        }
    };
    
    template<class UIntType, unsigned N>
    struct Right
    {
        static UIntType shift(UIntType a)
        {
            return a >> N;
        }
    };
    
    template<int N, class UIntType>
    inline UIntType shift(UIntType a)
    {
        return boost::mpl::if_c<(N < 0),
                    Left<UIntType, -N>,
                    Right<UIntType, N>
                >::type::shift(a);
    }
    

    became:

    template <typename UIntType, signed N, typename Enable = void> struct Shift;
    
    template <typename UIntType, signed N>
        struct Shift<UIntType, N, typename std::enable_if<(N>=0)>::type> {
            static UIntType apply(UIntType a) { return a >> N; }
        };
    
    template <typename UIntType, signed N>
        struct Shift<UIntType, N, typename std::enable_if<(N<0)>::type> {
            static UIntType apply(UIntType a) { return a << -N; }
        };
    
    template<int N, class UIntType>
    inline UIntType shift(UIntType a) { return Shift<UIntType, N>::apply(a); }
    
  10. Likewise, the Modulo switch (Power2Modulo and GenericModulo) that looked like this:

    /**
     * Conditional expression of type (r & (r - 1)) == 0 which allows to check
     * whether a number @f$r@f$ is of type @f$2^n@f$.
     */
    typedef boost::mpl::equal_to<
                boost::mpl::bitand_<
                    boost::mpl::_,
                    boost::mpl::minus<boost::mpl::_, boost::mpl::int_<1>
                >
            >,
            boost::mpl::int_<0>
        > IsPowerOfTwo;
    
    template<class UIntType, UIntType r>
    struct Power2Modulo
    {
        typedef typename boost::mpl::apply<
                IsPowerOfTwo,
                boost::mpl::integral_c<UIntType, r>
            >::type type;
    
        BOOST_STATIC_ASSERT(type::value);
    
        template<class T>
        static T calc(T value)
        {
            return value & (r - 1);
        }
    };
    
    template<class UIntType, UIntType r>
    struct GenericModulo
    {
        /**
         * @brief Determines @a value modulo @a r.
         *
         * @pre value >= 0 and value < 2 * r
         * @post value >= 0 and value < r
         */
        template<class T>
        static T calc(T value)
        {
            BOOST_STATIC_ASSERT(!std::numeric_limits<UIntType>::is_signed);
            assert(value < 2 * r);
    
            if (value >= r)
                value -= r;
    
            return value;
        }
    };
    
    template<class UIntType, UIntType r>
    struct Modulo
    {
        typedef typename boost::mpl::apply<
                IsPowerOfTwo,
                boost::mpl::integral_c<UIntType, r>
            >::type rIsPowerOfTwo;
    
        static UIntType calc(UIntType value)
        {
            // Use the bitwise AND for power 2 modulo arithmetic, or subtraction
            // otherwise. Subtraction is about two times faster than direct modulo
            // calculation.
            return boost::mpl::if_<
                        rIsPowerOfTwo,
                            Power2Modulo<UIntType, r>,
                            GenericModulo<UIntType, r>
                    >::type::calc(value);
        }
    };
    

    became much simpler with a little bit of c++11 (constexpr!) goodness:

    template <typename T, typename = typename std::enable_if<!std::is_signed<T>()>::type>
    constexpr static bool is_powerof2(T v) { return v && ((v & (v - 1)) == 0); }
    
    template<class UIntType, UIntType r>
    struct Modulo {
        template<class T> static T calc(T value) { return calc(value, std::integral_constant<bool, is_powerof2(r)>{}); }
        /**
         * @brief Determines @a value modulo @a r.
         *
         * @pre value >= 0 and value < 2 * r
         * @post value >= 0 and value < r
         */
        template<class T> static T calc(T value, std::true_type) { return value & (r - 1); }
        template<class T> static T calc(T value, std::false_type) {
            STATIC_ASSERT(!std::numeric_limits<UIntType>::is_signed);
            assert(value < 2 * r);
    
            if (value >= r)
                value -= r;
    
            return value;
        }
    };
    
  11. <boost/cstdint.hpp> -> <cstdint> (replacing ::boost by ::std for uint_least32_t and uint32_t)

  12. Well_quoted type function replaced by an alias template (template<...> using T = ... see http://en.cppreference.com/w/cpp/language/type_alias ad 2)

  13. typedefs rewritten as type aliases.

Full Listing

Live On Coliru

// Copyright (c) Sergiu Dotenco 2010, 2011, 2012
// Copyright (c) Seth Heeren - made independent of BOOST using C++11 - 2017
//
// Distributed under the Boost Software License, Version 1.0. (See accompanying
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

/**
 * @brief Implementation of the Well Equidistributed Long-period Linear (WELL)
 *        pseudo-random number generator.
 * @file well.hpp
 */

#ifndef WELL_HPP
#define WELL_HPP

#include <algorithm>
#include <cassert>
#include <cstddef>
#include <iomanip>
#include <istream>
#include <limits>
#include <ostream>
#include <functional>
#include <stdexcept>

#define STATIC_ASSERT(x) static_assert(x, #x)
#define PREVENT_MACRO_SUBSTITUTION

//! @cond hide_private

namespace Detail {
    using ulong_long_type = unsigned long long;

    template <typename UIntType, signed N, typename Enable = void> struct Shift;

    template <typename UIntType, signed N>
        struct Shift<UIntType, N, typename std::enable_if<(N>=0)>::type> {
            static UIntType apply(UIntType a) { return a >> N; }
        };

    template <typename UIntType, signed N>
        struct Shift<UIntType, N, typename std::enable_if<(N<0)>::type> {
            static UIntType apply(UIntType a) { return a << -N; }
        };

    template<int N, class UIntType>
    inline UIntType shift(UIntType a) {
        return Shift<UIntType, N>::apply(a);
    }

/**
 * @name Transformation matrices @f$M0,\dotsc,M6@f$ from Table I
 * @{
 */

struct M0
{
    template<class T>
    static T transform(T)
    {
        return T(0);
    }
};

struct M1
{
    template<class T>
    static T transform(T x)
    {
        return x;
    }
};

template<int N>
struct M2
{
    template<class T>
    static T transform(T x)
    {
        return shift<N>(x);
    }
};

template<int N>
struct M3
{
    template<class T>
    static T transform(T x)
    {
        return x ^ shift<N>(x);
    }
};

template<std::uint_least32_t a>
struct M4
{
    template<class T>
    static T transform(T x)
    {
        T result = x >> 1;

        if ((x & 1) == 1)
            result ^= a;

        return result;
    }
};

template<int N, std::uint_least32_t b>
struct M5
{
    template<class T>
    static T transform(T x)
    {
        return x ^ (shift<N>(x) & b);
    }
};

template
<
    std::size_t w,
    std::uint_least32_t q,
    std::uint_least32_t a,
    std::uint_least32_t ds,
    std::uint_least32_t dt
>
struct M6
{
    template<class T>
    static T transform(T x)
    {
        T result = ((x << q) ^ (x >> (w - q))) & ds;

        if ((x & dt) != 0)
            result ^= a;

        return result;
    }
};

//! @}

template <typename T, typename = typename std::enable_if<!std::is_signed<T>()>::type>
constexpr static bool is_powerof2(T v) { return v && ((v & (v - 1)) == 0); }

template<class UIntType, UIntType r>
struct Modulo {
    template<class T> static T calc(T value) { return calc(value, std::integral_constant<bool, is_powerof2(r)>{}); }
    /**
     * @brief Determines @a value modulo @a r.
     *
     * @pre value >= 0 and value < 2 * r
     * @post value >= 0 and value < r
     */
    template<class T> static T calc(T value, std::true_type) { return value & (r - 1); }
    template<class T> static T calc(T value, std::false_type) {
        STATIC_ASSERT(!std::numeric_limits<UIntType>::is_signed);
        assert(value < 2 * r);

        if (value >= r)
            value -= r;

        return value;
    }
};

template<std::uint_least32_t b, std::uint_least32_t c>
struct MatsumotoKuritaTempering
{
    template<std::size_t r, class UIntType, std::size_t N>
    static UIntType apply(UIntType x, UIntType (&)[N], std::size_t)
    {
        x ^= (x << 7) & b;
        x ^= (x << 15) & c;

        return x;
    }
};

template<std::uint_least32_t mask>
struct HaraseTempering
{
    template<std::size_t r, class UIntType, std::size_t N>
    static UIntType apply(UIntType x, UIntType (&s)[N], std::size_t m2)
    {
        return x ^ (s[Modulo<UIntType, r>::calc(m2 + 1)] & mask);
    }
};

struct NoTempering
{
    template<std::size_t r, class UIntType, std::size_t N>
    static UIntType apply(UIntType x, UIntType (&)[N], std::size_t)
    {
        return x;
    }
};

} // namespace Detail

//! @endcond

/**
 * @brief Well Equidistributed Long-period Linear (WELL) pseudo-random number
 *        generator.
 *
 * The implementation is based on the "Improved Long-Period Generators Based on
 * Linear Recurrences Modulo 2" paper by Francois Panneton, Pierre L'Ecuyer and
 * Makoto Matsumoto from ACM Transactions on Mathematical Software, 32 (1,
 * March) 2006, pp. 1-16.
 *
 * @tparam UIntType The unsigned integer type.
 * @tparam w Word size.
 * @tparam r State size.
 */
template
<
    class UIntType,
    std::size_t w,
    std::size_t r,
    std::size_t p,
    std::size_t m1,
    std::size_t m2,
    std::size_t m3,
    class T0,
    class T1,
    class T2,
    class T3,
    class T4,
    class T5,
    class T6,
    class T7,
    class Tempering // mpl pluggable
>
class Well
{
    STATIC_ASSERT(!std::numeric_limits<UIntType>::is_signed);
    STATIC_ASSERT(w <= static_cast<std::size_t>(std::numeric_limits<UIntType>::digits));
    STATIC_ASSERT(r > 0 && p < w);
    STATIC_ASSERT(m1 > 0 && m1 < r);
    STATIC_ASSERT(m2 > 0 && m2 < r);
    STATIC_ASSERT(m3 > 0 && m3 < r);

public:
    //! The unsigned integer type.
    typedef UIntType result_type;

    //! Word size.
    static constexpr std::size_t word_size = w;
    //! State size.
    static constexpr std::size_t state_size = r;
    //! Number of mask bits.
    static constexpr std::size_t mask_bits = p;
    //! Default seed value.
    static constexpr UIntType default_seed = 5489U;

    /**
     * @brief Initializes the class using the specified seed @a value.
     *
     * @param value The seed value to be used for state initialization.
     */
    explicit Well(result_type value = default_seed)
    {
        seed(value);
    }

    template<class InputIterator>
    Well(InputIterator& first, InputIterator last)
    {
        seed(first, last);
    }

    template<class Generator>
    explicit Well(Generator& g)
    {
        seed(g);
    }

    template<class Generator>
    void seed(Generator& g)
    {
        // Ensure std::generate_n doesn't copy the generator g by using
        // std::reference_wrapper
        std::generate_n(state_, state_size, std::ref(g));
    }

    void seed(result_type value = default_seed)
    {
        if (value == 0U)
            value = default_seed;

        state_[0] = value;

        std::size_t i = 1;
        UIntType *const s = state_;

        // Same generator used to seed Mersenne twister
        for ( ; i != state_size; ++i)
            s[i] = (1812433253U * (s[i - 1] ^ (s[i - 1] >> (w - 2))) + i);

        index_ = i;
    }

    template<class InputIterator>
    void seed(InputIterator& first, InputIterator last)
    {
        index_ = 0;
        std::size_t i = 0;

        for ( ; i != state_size && first != last; ++i, ++first)
            state_[i] = *first;

        if (first == last && i != state_size)
            throw std::invalid_argument("Seed sequence too short");
    }

    /**
     * @brief Generates a random number.
     */
    result_type operator()()
    {
        const UIntType upper_mask = ~0U << p;
        const UIntType lower_mask = ~upper_mask;

        // v[i,j] = state[(r-i+j) mod r]
        std::size_t i = index_;
        // Equivalent to r-i but allows to avoid negative values in the
        // following two expressions
        std::size_t j = i + r;
        std::size_t k = mod(j - 1); // [i,r-1]
        std::size_t l = mod(j - 2); // [i,r-2]

        std::size_t im1 = i + m1;
        std::size_t im2 = i + m2;
        std::size_t im3 = i + m3;

        UIntType z0, z1, z2, z3, z4;

        z0 = (state_[k] & upper_mask) | (state_[l] & lower_mask);
        z1 = T0::transform(state_[i]) ^
             T1::transform(state(im1));
        z2 = T2::transform(state(im2)) ^
             T3::transform(state(im3));
        z3 = z1 ^ z2;
        z4 = T4::transform(z0) ^ T5::transform(z1) ^
             T6::transform(z2) ^ T7::transform(z3);

        state_[i] = z3; // v[i+1,1]
        state_[k] = z4; // v[i+1,0]

        index_ = k;

        return Tempering::template apply<r>(z4, state_, im2);
    }

    result_type min PREVENT_MACRO_SUBSTITUTION () const
    {
        return 0U;
    }

    result_type max PREVENT_MACRO_SUBSTITUTION () const
    {
        return ~0U >> (std::numeric_limits<UIntType>::digits - w);
    }

    void discard(Detail::ulong_long_type z)
    {
        while (z-- > 0) {
            operator()();
        }
    }

    /**
     * @brief Compares the state of two generators for equality.
     */
    friend bool operator==(const Well& lhs, const Well& rhs)
    {
        for (std::size_t i = 0; i != state_size; ++i)
            if (lhs.compute(i) != rhs.compute(i))
                return false;

        return true;
    }

    /**
     * @brief Compares the state of two generators for inequality.
     */
    friend bool operator!=(const Well& lhs, const Well& rhs)
    {
        return !(lhs == rhs);
    }

    /**
     * @brief Writes the state to the specified stream.
     */
    template<class E, class T>
    friend std::basic_ostream<E, T>&
        operator<<(std::basic_ostream<E, T>& out, const Well& well)
    {
        E space = out.widen(' ');

        for (std::size_t i = 0; i != state_size; ++i)
            out << well.compute(i) << space;

        return out;
    }

    /**
     * @brief Reads the generator state from the specified input stream.
     */
    template<class E, class T>
    friend std::basic_istream<E, T>&
        operator>>(std::basic_istream<E, T>& in, Well& well)
    {
        for (std::size_t i = 0; i != state_size; ++i)
            in >> well.state_[i] >> std::ws;

        well.index_ = state_size;

        return in;
    }

private:
    template<class T>
    static T mod(T value)
    {
        return Detail::Modulo<T, r>::calc(value);
    }

    UIntType state(std::size_t index) const
    {
        return state_[mod(index)];
    }

    UIntType compute(std::size_t index) const
    {
        return state_[(index_ + index + r) % r];
    }

    UIntType state_[r];
    std::size_t index_;
};

namespace Detail {
    /**
     * @name Base definitions with pluggable tempering method
     * @{
     */

    template <typename Tempering>
    using Well512a_base = Well<
        std::uint32_t, 32, 16, 0, 13, 9, 5, M3<-16>, M3<-15>, M3<11>, M0, M3<-2>, M3<-18>, M2<-28>,
        M5<-5, 0xda442d24>, Tempering>;

    template <typename Tempering>
    using Well521a_base = Well<
        std::uint32_t, 32, 17, 23, 13, 11, 10, M3<-13>, M3<-15>, M1, M2<-21>,
        M3<-13>, M2<1>, M0, M3<11>, Tempering>;

    template <typename Tempering>
    using Well521b_base = Well<
        std::uint32_t, 32, 17, 23, 11, 10, 7, M3<-21>, M3<6>, M0, M3<-13>, M3<13>,
        M2<-10>, M2<-5>, M3<13>, Tempering>;

    template <typename Tempering>
    using Well607a_base = Well<
        std::uint32_t, 32, 19, 1, 16, 15, 14, M3<19>, M3<11>, M3<-14>, M1, M3<18>,
        M1, M0, M3<-5>, Tempering>;

    template <typename Tempering>
    using Well607b_base = Well<
        std::uint32_t, 32, 19, 1, 16, 18, 13, M3<-18>, M3<-14>, M0, M3<18>,
        M3<-24>, M3<5>, M3<-1>, M0, Tempering>;

    template <typename Tempering>
    using Well800a_base = Well<
        std::uint32_t, 32, 25, 0, 14, 18, 17, M1, M3<-15>, M3<10>, M3<-11>, M3<16>,
        M2<20>, M1, M3<-28>, Tempering>;

    template <typename Tempering>
    using Well800b_base = Well<
        std::uint32_t, 32, 25, 0, 9, 4, 22, M3<-29>, M2<-14>, M1, M2<19>, M1,
        M3<10>, M4<0xd3e43ffd>, M3<-25>, Tempering>;

    template <typename Tempering>
    using Well1024a_base = Well<
        std::uint32_t, 32, 32, 0, 3, 24, 10, M1, M3<8>, M3<-19>, M3<-14>, M3<-11>,
        M3<-7>, M3<-13>, M0, Tempering>;

    template <typename Tempering>
    using Well1024b_base = Well<
        std::uint32_t, 32, 32, 0, 22, 25, 26, M3<-21>, M3<17>, M4<0x8bdcb91e>,
        M3<15>, M3<-14>, M3<-21>, M1, M0, Tempering>;

    template <typename Tempering>
    using Well19937a_base = Well<
        std::uint32_t, 32, 624, 31, 70, 179, 449, M3<-25>, M3<27>, M2<9>, M3<1>,
        M1, M3<-9>, M3<-21>, M3<21>, Tempering>;

    template <typename Tempering>
    using Well19937b_base = Well<
        std::uint32_t, 32, 624, 31, 203, 613, 123, M3<7>, M1, M3<12>, M3<-10>,
        M3<-19>, M2<-11>, M3<4>, M3<-10>, Tempering>;

    template <typename Tempering>
    using Well21701a_base = Well<
        std::uint32_t, 32, 679, 27, 151, 327, 84, M1, M3<-26>, M3<19>, M0, M3<27>,
        M3<-11>, M6<32, 15, 0x86a9d87e, 0xffffffef, 0x00200000>, M3<-16>,
        Tempering>;

    template <typename Tempering>
    using Well23209a_base = Well<
        std::uint32_t, 32, 726, 23, 667, 43, 462, M3<28>, M1, M3<18>, M3<3>,
        M3<21>, M3<-17>, M3<-28>, M3<-1>, Tempering>;

    template <typename Tempering>
    using Well23209b_base = Well<
        std::uint32_t, 32, 726, 23, 610, 175, 662, M4<0xa8c296d1>, M1, M6<32, 15,
        0x5d6b45cc, 0xfffeffff, 0x00000002>, M3<-24>, M3<-26>, M1, M0, M3<16>,
        Tempering>;

    template <typename Tempering>
    using Well44497a_base = Well<
        std::uint32_t, 32, 1391, 15, 23, 481, 229, M3<-24>, M3<30>, M3<-10>,
        M2<-26>, M1, M3<20>, M6<32, 9, 0xb729fcec, 0xfbffffff, 0x00020000>, M1, Tempering>;
    //! @}

} // namespace Detail

using Well512a   = Detail::Well512a_base<Detail::NoTempering>;
using Well521a   = Detail::Well521a_base<Detail::NoTempering>;
using Well521b   = Detail::Well521b_base<Detail::NoTempering>;
using Well607a   = Detail::Well607a_base<Detail::NoTempering>;
using Well607b   = Detail::Well607b_base<Detail::NoTempering>;
using Well800a   = Detail::Well800a_base<Detail::NoTempering>;
using Well800b   = Detail::Well800b_base<Detail::NoTempering>;
using Well1024a  = Detail::Well1024a_base<Detail::NoTempering>;
using Well1024b  = Detail::Well1024b_base<Detail::NoTempering>;
using Well19937a = Detail::Well19937a_base<Detail::NoTempering>;
using Well19937b = Detail::Well19937b_base<Detail::NoTempering>;
using Well19937c = Detail::Well19937a_base<Detail::MatsumotoKuritaTempering<0xe46e1700, 0x9b868000>>;
using Well21701a = Detail::Well21701a_base<Detail::NoTempering>;
using Well23209a = Detail::Well23209a_base<Detail::NoTempering>;
using Well23209b = Detail::Well23209b_base<Detail::NoTempering>;
using Well44497a = Detail::Well44497a_base<Detail::NoTempering>;
using Well44497b = Detail::Well44497a_base<Detail::MatsumotoKuritaTempering<0x93dd1400, 0xfa118000>>;

/**
 * @name Maximally equidistributed versions using Harase's tempering method
 * @{
 */
using Well800a_ME   = Detail::Well800a_base<Detail::HaraseTempering<0x4880>>;
using Well800b_ME   = Detail::Well800b_base<Detail::HaraseTempering<0x17030806>>;
using Well19937a_ME = Detail::Well19937a_base<Detail::HaraseTempering<0x4118000>>;
using Well19937b_ME = Detail::Well19937b_base<Detail::HaraseTempering<0x30200010>>;
using Well21701a_ME = Detail::Well21701a_base<Detail::HaraseTempering<0x1002>>;
using Well23209a_ME = Detail::Well23209a_base<Detail::HaraseTempering<0x5100000>>;
using Well23209b_ME = Detail::Well23209b_base<Detail::HaraseTempering<0x34000300>>;
using Well44497a_ME = Detail::Well44497a_base<Detail::HaraseTempering<0x48000000>>;
//! @}

#endif // WELL_HPP

// Copyright (c) Sergiu Dotenco 2010
// Copyright (c) Seth Heeren - made independent of BOOST using C++11 - 2017
//
// Distributed under the Boost Software License, Version 1.0. (See accompanying
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

/**
 * @brief WELL PRNG implementation unit test.
 * @file welltest.cpp
 */

#include <algorithm>
#include <memory>
#include <iostream>

// #include "well.hpp"

#define MESSAGE_PREAMBLE() (std::cerr << __FILE__ << ":" << __LINE__ << " ")

#define CHECK_EQUAL(expected,actual) do { \
    auto&& _e = expected; \
    auto&& _a = actual; \
    if (_e != _a) \
        MESSAGE_PREAMBLE() << "expected:" << #expected << " = " << _e << "\n" \
                  << "\tactual:" << #actual << " = " << _a << "\n"; \
} while (0)

#define CHECK(test) do { if (!(test)) MESSAGE_PREAMBLE() << #test << "\n"; } while (0)

/**
 * @brief Generic WELL test case.
 *
 * The test case performs the following checks:
 * -# The last generated value is equal to the value generate by the reference
 *    implementation after @f$10^9@f$ iterations. The generator is seeded using
 *    an array filled with 1s.
 * -# The @c min and @c max methods of the @ref Well generator return 0 and
 *    @f$2^{32}-1@f$ respectively.
 *
 * @tparam RandomNumberGenerator WELL PRNG implementation type.
 * @tparam Expected The expected result after @f$10^9@f$ iterations.
 */
template
<
    class RandomNumberGenerator,
    typename RandomNumberGenerator::result_type Expected
>
class WellTestCase
{
    RandomNumberGenerator rng;

    typedef typename RandomNumberGenerator::result_type result_type;

    result_type generate()
    {
        unsigned state[RandomNumberGenerator::state_size];
        std::uninitialized_fill_n(state, RandomNumberGenerator::state_size, 1);

        unsigned* p = state;
        rng.seed(p, p + RandomNumberGenerator::state_size);

        result_type x = 0;

        int iterations = 1000000000;

        while (iterations-- > 0)
            x = rng();

        return x;
    }

public:
    static void run()
    {
        WellTestCase c;

        CHECK_EQUAL(c.generate(), Expected);
        CHECK_EQUAL(c.rng.min(), 0U);
        CHECK_EQUAL(c.rng.max(), ~0U);
        CHECK_EQUAL(c.rng, c.rng);
        CHECK(c.rng == c.rng);
    }
};

/**
 * @brief Defines the actual test case.
 *
 * @param name The name of the test case.
 * @param type WELL pseudo-random generator type.
 * @param expected The expected result after @f$10^9@f$ iterations.
 *
 * @hideinitializer
 */
int main() {
    CHECK_EQUAL(16, Detail::shift<2>(64));
    CHECK_EQUAL(64, Detail::shift<-2>(16));
    CHECK_EQUAL(32, Detail::shift<0>(32));
    CHECK(Detail::is_powerof2(512u));
    CHECK(not Detail::is_powerof2(0u));

    WellTestCase<Well512a,   0x2b3fe99e>::run();

#ifndef COLIRU // stay in execution time limits
    WellTestCase<Well521a,   0xc9878363>::run();
    WellTestCase<Well521b,   0xb75867f6>::run();
    WellTestCase<Well607a,   0x7b5043ea>::run();
    WellTestCase<Well607b,   0xaedee7da>::run();
    WellTestCase<Well800a,   0x2bfe686f>::run();
    WellTestCase<Well800b,   0xf009e1bd>::run();
    WellTestCase<Well1024a,  0xd07f528c>::run();
    WellTestCase<Well1024b,  0x867f7993>::run();
    WellTestCase<Well19937a, 0xb33a2cd5>::run();
    WellTestCase<Well19937b, 0x191de86a>::run();
    WellTestCase<Well19937c, 0x243eaed5>::run();
    WellTestCase<Well21701a, 0x7365a269>::run();
    WellTestCase<Well23209a, 0x807dacb >::run();
    WellTestCase<Well23209b, 0xf1a77751>::run();
    WellTestCase<Well44497a, 0xfdd7c07b>::run();
    WellTestCase<Well44497b, 0x9406547b>::run();
#endif
}
sehe
  • 374,641
  • 47
  • 450
  • 633
  • `template ()>::type>` has a compilation error: "E0065 expected a ';'". Seems shall change to `template ()>>`. – athos Dec 15 '17 at 02:40
  • if I add `Well44497a_ME well44497a_me; auto r=well44497a_me();` in `main()`, there'll be a compilation error : C2027 use of undefined type 'Detail::Modulo' – athos Dec 15 '17 at 04:06
  • after replacing your `Modulo` part with Caleth's solution in https://stackoverflow.com/a/47711587/826203 , it can compile. Still, I wonder why it doesn't pass, as your `Modulo` seems is using the same idea of `is_powerof2`, i.e. using `enable_if` for SFINAE. – athos Dec 15 '17 at 06:43
  • 2
    My latest [commit](https://github.com/sehe/well-random/commit/a4ad83d9f3c2f5715e565ea364fbf89a4fbd9354) fixes the compatibility for MSVC++. See it **[Live On MSVC](http://rextester.com/ROYB40593)** as well as **[Live On GCC](http://coliru.stacked-crooked.com/a/8fefbef5ae28410b)** and **[Clang](http://coliru.stacked-crooked.com/a/fc860f46a6ea6b18)**. They were MSVC conformance issues (see the commit message). – sehe Dec 15 '17 at 08:14
8

Using C++17, this code becomes way simpler and error messages are friendlier on the eye.

This is a sample implementation of Power2Modulo:

#include <type_traits>

template<class UIntType, UIntType r>
struct Power2Modulo
{
  static_assert(std::is_unsigned_v<UIntType>);
  static_assert((r & (r - 1)) == 0,
     "The second parameter of this struct is required to be a power of 2");

  template<class T>
  [[nodiscard]] static constexpr T calc(T value)
  {
    return value & (r - 1);
  }
};

You can use it like this:

int main()
{ 
  /* This code fails to compile with friendly error message
  Power2Modulo<unsigned, 12> x;
  */

  // Using the static function
  using Mod16 = Power2Modulo<unsigned, 16>;
  static_assert(Mod16::calc(15) == 15);
  static_assert(Mod16::calc(16) == 0);
  static_assert(Mod16::calc(17) == 1);

  // Using it like a member function
  Power2Modulo<unsigned, 4> mod4;
  static_assert(mod4.calc(15) == 3);
  static_assert(mod4.calc(16) == 0);
  static_assert(mod4.calc(17) == 1);
}

Tested with clang-6 and gcc-8 and VisualC++ (via http://webcompiler.cloudapp.net/).

Rumburak
  • 3,416
  • 16
  • 27