Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to rewrite this code without using boost?

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.

like image 289
athos Avatar asked Dec 13 '17 01:12

athos


2 Answers

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
}
like image 116
sehe Avatar answered Oct 08 '22 13:10

sehe


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/).

like image 28
Rumburak Avatar answered Oct 08 '22 14:10

Rumburak