glest-source/source/shared_lib/sources/streflop/Random.cpp

/*
    streflop: STandalone REproducible FLOating-Point
    Nicolas Brodu, 2006
    Code released according to the GNU Lesser General Public License

    Heavily relies on GNU Libm, itself depending on netlib fplibm, GNU MP, and IBM MP lib.
    Uses SoftFloat too.

    Please read the history and copyright information in the documentation provided with the source code
*/

// Include time(0) function to get a seed based on system time
#include <time.h>
#include <iostream>
using namespace std;
#include "streflop.h"

// Include endian-specific code
#undef __BYTE_ORDER
#undef __FLOAT_WORD_ORDER
#include "System.h"

namespace streflop {

//////////////////////////////////////////////////////////////////////
// Code stolen and adapted from mt19937ar.c
//////////////////////////////////////////////////////////////////////

/*
   A C-program for MT19937, with initialization improved 2002/1/26.
   Coded by Takuji Nishimura and Makoto Matsumoto.

   Before using, initialize the state by using init_genrand(seed)  
   or init_by_array(init_key, key_length).

   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
   All rights reserved.                          
   Copyright (C) 2005, Mutsuo Saito,
   All rights reserved.                          

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:

     1. Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.

     2. Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.

     3. The names of its contributors may not be used to endorse or promote 
        products derived from this software without specific prior written 
        permission.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


   Any feedback is very welcome.
   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
*/

#if STREFLOP_RANDOM_GEN_SIZE == 32

/* Period parameters */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
#define LOWER_MASK 0x7fffffffUL /* least significant r bits */

/* initializes mt[N] with a seed */
inline void init_genrand(SizedUnsignedInteger<32>::Type s, RandomState& state)
{
    state.seed = s;
    state.mt[0]= s; // & 0xffffffffUL; // NB060508: unnecessary with the use of sized types
    for (state.mti=1; state.mti<N; state.mti++) {
        state.mt[state.mti] = 
        (1812433253UL * (state.mt[state.mti-1] ^ (state.mt[state.mti-1] >> 30)) + state.mti);
        /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
        /* In the previous versions, MSBs of the seed affect   */
        /* only MSBs of the array mt[].                        */
        /* 2002/01/09 modified by Makoto Matsumoto             */
        //mt[mti] &= 0xffffffffUL;  // NB060508: unnecessary with the use of sized types
        /* for >32 bit machines */
    }
}

/* generates a random number on [0,0xffffffff]-interval */
inline SizedUnsignedInteger<32>::Type genrand_int(RandomState& state)
{
    SizedUnsignedInteger<32>::Type y;
    static SizedUnsignedInteger<32>::Type mag01[2]={0x0UL, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */

    if (state.mti >= N) { /* generate N words at one time */
        int kk;

        //if (state.mti == N+1)   /* if init_genrand() has not been called, */
            //init_genrand(5489UL, state); /* a default initial seed is used */

        for (kk=0;kk<N-M;kk++) {
            y = (state.mt[kk]&UPPER_MASK)|(state.mt[kk+1]&LOWER_MASK);
            state.mt[kk] = state.mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        for (;kk<N-1;kk++) {
            y = (state.mt[kk]&UPPER_MASK)|(state.mt[kk+1]&LOWER_MASK);
            state.mt[kk] = state.mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        y = (state.mt[N-1]&UPPER_MASK)|(state.mt[0]&LOWER_MASK);
        state.mt[N-1] = state.mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];

        state.mti = 0;
    }
  
    y = state.mt[state.mti++];

    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return y;
}

#else

//////////////////////////////////////////////////////////////////////
// End of code adapted from mt19937ar.c
// Now adapt code from the 64-bit version in mt19937-64.c
//////////////////////////////////////////////////////////////////////

/* 
   A C-program for MT19937-64 (2004/9/29 version).
   Coded by Takuji Nishimura and Makoto Matsumoto.

   This is a 64-bit version of Mersenne Twister pseudorandom number
   generator.

   Before using, initialize the state by using init_genrand64(seed)  
   or init_by_array64(init_key, key_length).

   Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,
   All rights reserved.                          

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:

     1. Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.

     2. Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.

     3. The names of its contributors may not be used to endorse or promote 
        products derived from this software without specific prior written 
        permission.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

   References:
   T. Nishimura, ``Tables of 64-bit Mersenne Twisters''
     ACM Transactions on Modeling and 
     Computer Simulation 10. (2000) 348--357.
   M. Matsumoto and T. Nishimura,
     ``Mersenne Twister: a 623-dimensionally equidistributed
       uniform pseudorandom number generator''
     ACM Transactions on Modeling and 
     Computer Simulation 8. (Jan. 1998) 3--30.

   Any feedback is very welcome.
   http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html
   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)
*/

#define NN 312
#define MM 156
#define MATRIX_A 0xB5026F5AA96619E9ULL
#define UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */
#define LM 0x7FFFFFFFULL /* Least significant 31 bits */

/* initializes mt[NN] with a seed */
inline void init_genrand(SizedUnsignedInteger<64>::Type seed, RandomState& state)
{
    state.seed = seed;
    state.mt[0] = seed;
    for (state.mti=1; state.mti<NN; state.mti++)
        state.mt[state.mti] =  (SizedUnsignedInteger<64>::Type(6364136223846793005ULL) * (state.mt[state.mti-1] ^ (state.mt[state.mti-1] >> 62)) + state.mti);
}

/* generates a random number on [0, 2^64-1]-interval */
inline SizedUnsignedInteger<64>::Type genrand_int(RandomState& state)
{

    SizedUnsignedInteger<64>::Type x;
    static SizedUnsignedInteger<64>::Type mag01[2]={0ULL, MATRIX_A};

    if (state.mti >= NN) { /* generate NN words at one time */

        /* if init_genrand64() has not been called, */
        /* a default initial seed is used     */
        //if (state.mti == NN+1)
            //init_genrand64(5489ULL, state);
    	int i = 0;
        for (i=0;i<NN-MM;i++) {
            x = (state.mt[i]&UM)|(state.mt[i+1]&LM);
            state.mt[i] = state.mt[i+MM] ^ (x>>1) ^ mag01[(int)(x&1ULL)];
        }
        for (;i<NN-1;i++) {
            x = (state.mt[i]&UM)|(state.mt[i+1]&LM);
            state.mt[i] = state.mt[i+(MM-NN)] ^ (x>>1) ^ mag01[(int)(x&1ULL)];
        }
        x = (state.mt[NN-1]&UM)|(state.mt[0]&LM);
        state.mt[NN-1] = state.mt[MM-1] ^ (x>>1) ^ mag01[(int)(x&1ULL)];

        state.mti = 0;
    }
  
    x = state.mt[state.mti++];

    x ^= (x >> 29) & 0x5555555555555555ULL;
    x ^= (x << 17) & 0x71D67FFFEDA60000ULL;
    x ^= (x << 37) & 0xFFF7EEE000000000ULL;
    x ^= (x >> 43);

    return x;
}
#endif

//////////////////////////////////////////////////////////////////////
// End of code adapted from mt19937-64.c
//////////////////////////////////////////////////////////////////////

// Bit getter utilities
template<int nbits> struct Accessor {
    typedef typename SizedUnsignedInteger<nbits>::Type Type;
    static inline Type getRandomInt(RandomState& state) {
        return static_cast<Type>(genrand_int(state));
    }
};

// Specialize for 32 bits generator case
#if STREFLOP_RANDOM_GEN_SIZE == 32
template<> struct Accessor<64> {
    typedef SizedUnsignedInteger<64>::Type Type;
    static inline Type getRandomInt(RandomState& state) {
        return static_cast<Type>(genrand_int(state)) | (static_cast<Type>(genrand_int(state)) << 32);
    }
};
#endif


// This code inspired from a trick found in Richard J. Wagner's Mersene class and the optimization
// by Magnus Jonsson, also found at http://aggregate.org.
// The trick consists in:
// - draw random numbers in as close as possible as the target
// - reject these which are out of range
// The goal is to avoid operator %

template<int bits_size> struct RandomIntRestrictor {
};

// for 8-bits
template<> struct RandomIntRestrictor<8> {
    typedef SizedUnsignedInteger<8>::Type Type;

    static inline Type getRestrictedRandomInt(Type n, RandomState& state)
    {
        // First propagate leading 1 to all other bits
        Type mask = n;
        mask |= mask >> 1;
        mask |= mask >> 2;
        mask |= mask >> 4;

        // Draw only that number of bits, until a number is in the desired range [0,n]
        // Worse case is number of loops proba decreasing in 1/2^nloops
        Type ret;
        do {
            ret = Accessor<8>::getRandomInt(state) & mask;
        } while( ret > n );

        return ret;
    }
};

// for 16-bits
template<> struct RandomIntRestrictor<16> {
    typedef SizedUnsignedInteger<16>::Type Type;

    static inline Type getRestrictedRandomInt(Type n, RandomState& state)
    {
        // First propagate leading 1 to all other bits
        Type mask = n;
        mask |= mask >> 1;
        mask |= mask >> 2;
        mask |= mask >> 4;
        mask |= mask >> 8;

        // Draw only that number of bits, until a number is in the desired range [0,n]
        // Worse case is number of loops proba decreasing in 1/2^nloops
        Type ret;
        do {
            ret = Accessor<16>::getRandomInt(state) & mask;
        } while( ret > n );

        return ret;
    }
};

// for 32-bits
template<> struct RandomIntRestrictor<32> {
    typedef SizedUnsignedInteger<32>::Type Type;

    static inline Type getRestrictedRandomInt(Type n, RandomState& state)
    {
        // First propagate leading 1 to all other bits
        Type mask = n;
        mask |= mask >> 1;
        mask |= mask >> 2;
        mask |= mask >> 4;
        mask |= mask >> 8;
        mask |= mask >> 16;

        // Draw only that number of bits, until a number is in the desired range [0,n]
        // Worse case is number of loops proba decreasing in 1/2^nloops
        Type ret;
        do {
            ret = Accessor<32>::getRandomInt(state) & mask;
        } while( ret > n );

        return ret;
    }
};

// for 64-bits
template<> struct RandomIntRestrictor<64> {
    typedef SizedUnsignedInteger<64>::Type Type;

    static inline Type getRestrictedRandomInt(Type n, RandomState& state)
    {
        // First propagate leading 1 to all other bits
        Type mask = n;
        mask |= mask >> 1;
        mask |= mask >> 2;
        mask |= mask >> 4;
        mask |= mask >> 8;
        mask |= mask >> 16;
        mask |= mask >> 32;

        // Draw only that number of bits, until a number is in the desired range [0,n]
        // Worse case is number of loops proba descreasing in 1/2^nloops
        Type ret;
        do {
            ret = Accessor<64>::getRandomInt(state) & mask;
        } while( ret > n );

        return ret;
    }
};


// Now implement the Random functions

/* range checker.
   But why should clean caller code be slowed down by checks?
   => caller code should provide good arguments, or be fixed

/// Common function that will be used for all random integer types
template<typename an_int_type, int num_bounds, int min_excluded> struct RandomSwitcher {
    // random selection
    static inline an_int_type getMinMaxRandom(an_int_type min, an_int_type max, RandomState& state) {
        // Works in both signed and unsigned arithmetic
        if (max<=min) {
            if (max==min) return min; // easy
            return RandomSwitcher<an_int_type,num_bounds,(num_bounds==1)?(!min_excluded):min_excluded>::getMinMaxRandom(max, min, state);
        }

        // now max > min
        an_int_type range = max - min + 1;

        // overflow, asking the whole range of integers
        if (range==0) return getRandomBits<a_type>(sizeof(an_int_type)*STREFLOP_INTEGER_TYPES_CHAR_BITS, state);

        // range >= 2, num_bounds <= 2, OK to subtract 2 - num_bounds
        range -= 2-num_bounds;
        // ask the impossible, like RandomEE(3,4)
        if (range == 0) return min;

        // This gives the number of possible integers to choose from (at least 1)
        // Now generate the number
        an_int_type ret = min + min_excluded + static_cast<an_int_type>(RandomIntRestrictor< sizeof(an_int_type)*STREFLOP_INTEGER_TYPES_CHAR_BITS >::getRestrictedRandomInt(range,state));
        return ret;
    }
};
*/

#define SPECIALIZE_RANDOM_FOR_TYPE(a_type,use_signed) \
template<> a_type Random<a_type>(RandomState& state) { \
    return Accessor<sizeof(a_type)*STREFLOP_INTEGER_TYPES_CHAR_BITS>::getRandomInt(state); \
} \
template<> a_type Random<true, true, a_type>(a_type min, a_type max, RandomState& state) { \
    return static_cast<a_type>(RandomIntRestrictor< sizeof(a_type)*STREFLOP_INTEGER_TYPES_CHAR_BITS >::getRestrictedRandomInt(max-min,state)) + min; \
} \
template<> a_type Random<true, false, a_type>(a_type min, a_type max, RandomState& state) { \
    return static_cast<a_type>(RandomIntRestrictor< sizeof(a_type)*STREFLOP_INTEGER_TYPES_CHAR_BITS >::getRestrictedRandomInt(max-min-1,state)) + min; \
} \
template<> a_type Random<false, true, a_type>(a_type min, a_type max, RandomState& state) { \
    return max - static_cast<a_type>(RandomIntRestrictor< sizeof(a_type)*STREFLOP_INTEGER_TYPES_CHAR_BITS >::getRestrictedRandomInt(max-min-1,state)); \
} \
template<> a_type Random<false, false, a_type>(a_type min, a_type max, RandomState& state) { \
    return static_cast<a_type>(RandomIntRestrictor< sizeof(a_type)*STREFLOP_INTEGER_TYPES_CHAR_BITS >::getRestrictedRandomInt(max-min-2,state)) + min + 1; \
}
SPECIALIZE_RANDOM_FOR_TYPE(char,true)
SPECIALIZE_RANDOM_FOR_TYPE(unsigned char,false)
SPECIALIZE_RANDOM_FOR_TYPE(short,true)
SPECIALIZE_RANDOM_FOR_TYPE(unsigned short,false)
SPECIALIZE_RANDOM_FOR_TYPE(int,true)
SPECIALIZE_RANDOM_FOR_TYPE(unsigned int,false)
SPECIALIZE_RANDOM_FOR_TYPE(long,true)
SPECIALIZE_RANDOM_FOR_TYPE(unsigned long,false)
SPECIALIZE_RANDOM_FOR_TYPE(long long,true)
SPECIALIZE_RANDOM_FOR_TYPE(unsigned long long,false)


// Random for float types is even more dependent on size


/*

EXPERIMENTAL:

    The ULP_Random functions use a uniform distribution in the ULP space.
    This means, each possibly representable float is given exactly the same weight for the
    random choice. This is an exponential scale where there are as many numbers between
    successive powers of 2 (i.e as many numbers between 1 and 2 as there are between 1024 and 2048).
    This effectively corresponds to the maximum machine precision, but it is unfortunately
    not what one means by "uniform" in the traditional sense of the term (it is uniform in terms
    of bit patterns, so for the computer!)

    Algorithm: treat floats as binary pattern, thanks to IEEE754 ordered property
    => this gives a "range" like max-min for the integers, where the unit is ulp
    => this gives the maximum precision for floats, because a random number is chosen between exactly how many
    numbers can be represented with that float format in that range


#define SPECIALIZE_RANDOM_FOR_SIMPLE_ULP(X,Y) \
Simple ULP_Random ## X ## Y(Simple min, Simple max) { \
 \
    // Convert to signed integer for quick test of bit sign \
    SizedInteger<32>::Type imin = *reinterpret_cast<SizedUnsignedInteger<32>::Type*>(&min); \
    SizedInteger<32>::Type imax = *reinterpret_cast<SizedUnsignedInteger<32>::Type*>(&max); \
 \
    // Infinity is a perfectly fine number, with a bit pattern contiguous to the min and max floats \
    // This makes sense if excluding bounds: RandomEE(-inf,+inf) returns a possible float at random \
 \
    // Rule out NaNs \
    if (imin&0x7fffffff > 0x7f800000) return SimpleNaN; \
    if (imax&0x7fffffff > 0x7f800000) return SimpleNaN; \
 \
    // Convert to 2-complement representation \
    if (imin<0) imin = 0x80000000 - imin; \
    if (imax<0) imax = 0x80000000 - imax; \
 \
    // Now magically get an integer at random in that range \
    // This gives EXACTLY one choice per representable float \
    // It is non-uniform in the float space, but uniform in the bit-pattern space \
    SizedInteger<32>::Type iret = Random ## X ## Y(imin,imax); \
 \
    // convert back to 2-complement to IEEE754 \
    if (iret<0) iret = 0x80000000 - iret; \
 \
    // cast to float \
    return *reinterpret_cast<Simple*>(&iret); \
}
SPECIALIZE_RANDOM_FOR_SIMPLE_ULP(E,I)
SPECIALIZE_RANDOM_FOR_SIMPLE_ULP(E,E)
SPECIALIZE_RANDOM_FOR_SIMPLE_ULP(I,E)
SPECIALIZE_RANDOM_FOR_SIMPLE_ULP(I,I)
*/


// Return a random float
template<> Simple Random<Simple>(RandomState& state) {
    // Generate bits
    SizedUnsignedInteger<32>::Type ret = Accessor<32>::getRandomInt(state);

    // Discard NaNs and Inf, ignore sign
    while ((ret & 0x7fffffff) >= 0x7f800000) ret = Accessor<32>::getRandomInt(state);

    // cast to float
    return *reinterpret_cast<Simple*>(&ret);
}

// Random in 1..2 - ideal IE case
template<> Simple Random12<true,false,Simple>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Generate bits
    SizedUnsignedInteger<32>::Type r12 = Accessor<32>::getRandomInt(state);

    // Simple precision keeps only 23 bits
    r12 &= 0x007FFFFF;

    // Insert exponent so it's in the [1.0-2.0) range
    r12 |= 0x3F800000;

    return *reinterpret_cast<Simple*>(&r12);
}

// Random in 1..2 - near ideal EI case
template<> Simple Random12<false,true,Simple>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Generate bits
    SizedUnsignedInteger<32>::Type r12 = Accessor<32>::getRandomInt(state);

    // Simple precision keeps only 23 bits
    r12 &= 0x007FFFFF;

    // Insert exponent so it's in the [1.0-2.0) range
    r12 |= 0x3F800000;

    // Bitwise add 1 so it's in the (1.0-2.0] range
    r12 += 1;

    return *reinterpret_cast<Simple*>(&r12);
}

// Random in 1..2 - need to include both bounds
template<> Simple Random12<true,true,Simple>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Generate bits
    SizedUnsignedInteger<32>::Type r12 = Accessor<32>::getRandomInt(state);

    // Keep 2^23 + 1 possibilities, discard others
    // Note: %= operator is nicely converted into reciprocal multiply and shift by compiler
    r12 %= 0x00800001;
    // Choose to avoid % operator by having about 1/2 chance of rejection. Not faster.
//    while ((r12 &= 0x00FFFFFF) > 0x00800000) r12 = Accessor<32>::getRandomInt(state);

/* could also use multiply by reciprocal and then find remainder. For div by 0x00800001:
; dividend: register other than EAX or memory location
MOV    EAX, 0FFFFFE01h
MUL    dividend
SHR    EDX, 23
; quotient now in EDX
*/

    // bitwise add exponent so it's in the [1.0-2.0] range
    r12 += 0x3F800000;

    return *reinterpret_cast<Simple*>(&r12);
}

// Random in 1..2 - need to exclude both bounds
template<> Simple Random12<false,false,Simple>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Generate bits
    SizedUnsignedInteger<32>::Type r12 = Accessor<32>::getRandomInt(state);

    // Keep 2^23 - 1 possibilities
    //r12 %= 0x007FFFFF;
    // Choose to avoid % operator by having very small chance of rejection
    // Could we find a branchless version for % by 2^N-1 ?
    while ((r12 &= 0x007FFFFF) == 0x007FFFFF) r12 = Accessor<32>::getRandomInt(state);

    // bitwise add exponent so it's in the (1.0-2.0) range
    r12 += 0x3F800001;

    return *reinterpret_cast<Simple*>(&r12);
}


///////// Double versions  ///////////

// Return a random float
template<> Double Random<Double>(RandomState& state) {
    // Generate bits
    SizedUnsignedInteger<64>::Type ret = Accessor<64>::getRandomInt(state);

    // Discard NaNs and Inf, ignore sign
    while ((ret & 0x7fffffffffffffffULL) >= 0x7ff0000000000000ULL) ret = Accessor<64>::getRandomInt(state);

    // cast to Double
    return *reinterpret_cast<Double*>(&ret);
}


// Random in a 1..2 - ideal IE case
template<> Double Random12<true,false,Double>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Generate bits
    SizedUnsignedInteger<64>::Type r12 = Accessor<64>::getRandomInt(state);

    // Double precision keeps only 52 bits
    r12 &= 0x000FFFFFFFFFFFFFULL;

    // Insert exponent so it's in the [1.0-2.0) range
    r12 |= 0x3FF0000000000000ULL;

    // scale from 1-2 interval to the desired interval
    return *reinterpret_cast<Double*>(&r12);
}

// Random in a 1..2 - near ideal EI case
template<> Double Random12<false,true,Double>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Generate bits
    SizedUnsignedInteger<64>::Type r12 = Accessor<64>::getRandomInt(state);

    // Double precision keeps only 52 bits
    r12 &= 0x000FFFFFFFFFFFFFULL;

    // Insert exponent so it's in the [1.0-2.0) range
    r12 |= 0x3FF0000000000000ULL;

    // Bitwise add 1 so it's in the (1.0-2.0] range
    r12 += 1;

    // scale from 1-2 interval to the desired interval
    return *reinterpret_cast<Double*>(&r12);
}

// Random in a 1..2 - need to include both bounds
template<> Double Random12<true,true,Double>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Generate bits
    SizedUnsignedInteger<64>::Type r12 = Accessor<64>::getRandomInt(state);

    // Keep 2^52 + 1 possibilities
    // See comment in Simple version
    // allow %= only for 64-bit register machines
#if STREFLOP_RANDOM_GEN_SIZE == 64
    r12 %= 0x0010000000000001ULL;
#else
    // Choose to avoid % operator by having about 1/2 chance of rejection. Is this faster?
    while ((r12 &= 0x001FFFFFFFFFFFFFULL) > 0x0010000000000000ULL) r12 = Accessor<64>::getRandomInt(state);
#endif

    // bitwise add exponent so it's in the [1.0-2.0] range
    r12 += 0x3FF0000000000000ULL;

    return *reinterpret_cast<Double*>(&r12);
}

// Random in a 1..2 - need to exclude both bounds
template<> Double Random12<false,false,Double>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Generate bits
    SizedUnsignedInteger<64>::Type r12 = Accessor<64>::getRandomInt(state);

    // Keep 2^52 - 1 possibilities
    // See comment in Simple version
    // r12 %= 0x000FFFFFFFFFFFFFULL;
    // Choose to avoid % operator by having very small chance of rejection
    while ((r12 &= 0x000FFFFFFFFFFFFFULL) == 0x000FFFFFFFFFFFFFULL) r12 = Accessor<64>::getRandomInt(state);

    // bitwise add exponent so it's in the (1.0-2.0) range
    r12 += 0x3FF0000000000001ULL;

    return *reinterpret_cast<Double*>(&r12);
}


#ifdef Extended

// little endian
#if __FLOAT_WORD_ORDER == 1234

/// Little endian is fine, always the same offsets whatever the memory model
template<int Nbytes> struct ExtendedConverter {
    // Sign and exponent
    static inline SizedUnsignedInteger<16>::Type* sexpPtr(Extended* e) {
        return reinterpret_cast<SizedUnsignedInteger<16>::Type*>(reinterpret_cast<char*>(e)+8);
    }
    // Mantissa
    static inline SizedUnsignedInteger<64>::Type* mPtr(Extended* e) {
        return reinterpret_cast<SizedUnsignedInteger<64>::Type*>(e);
    }
};

// Big endian
#elif __FLOAT_WORD_ORDER == 4321
template<int Nbytes> struct ExtendedConverter {
}

/// Extended is softfloat
template<> struct ExtendedConverter<10> {
    // Sign and exponent
    static inline SizedUnsignedInteger<16>::Type* sexpPtr(Extended* e) {
        return reinterpret_cast<SizedUnsignedInteger<16>::Type*>(e);
    }
    // Mantissa
    static inline SizedUnsignedInteger<64>::Type* mPtr(Extended* e) {
        return reinterpret_cast<SizedUnsignedInteger<64>::Type*>(reinterpret_cast<char*>(e)+2);
    }
};

/// Default gcc model for x87 - 32 bits (or with -m96bit-long-double)
template<> struct ExtendedConverter<12> {
    // Sign and exponent
    static inline SizedUnsignedInteger<16>::Type* sexpPtr(Extended* e) {
        return reinterpret_cast<SizedUnsignedInteger<16>::Type*>(reinterpret_cast<char*>(e)+2);
    }
    // Mantissa
    static inline SizedUnsignedInteger<64>::Type* mPtr(Extended* e) {
        return reinterpret_cast<SizedUnsignedInteger<64>::Type*>(reinterpret_cast<char*>(e)+4);
    }
};

/// Default gcc model for x87 - 64 bits (or with -m128bit-long-double)
template<> struct ExtendedConverter<16> {
    // Sign and exponent
    static inline SizedUnsignedInteger<16>::Type* sexpPtr(Extended* e) {
        return reinterpret_cast<SizedUnsignedInteger<16>::Type*>(reinterpret_cast<char*>(e)+6);
    }
    // Mantissa
    static inline SizedUnsignedInteger<64>::Type* mPtr(Extended* e) {
        return reinterpret_cast<SizedUnsignedInteger<64>::Type*>(reinterpret_cast<char*>(e)+8);
    }
};
#else
#error unknown byte order
#endif


// Return a random float
template<> Extended Random<Extended>(RandomState& state) {
    // Work directly on Extended bits
    Extended ret;

    // Generate 63 bits for the mantissa
    *ExtendedConverter<sizeof(Extended)>::mPtr(&ret) = Accessor<64>::getRandomInt(state);

    // Generate 16 bits for the exponent
    *ExtendedConverter<sizeof(Extended)>::sexpPtr(&ret) = Accessor<16>::getRandomInt(state);

    // Discard NaNs and Inf, ignore sign
    while ((*ExtendedConverter<sizeof(Extended)>::sexpPtr(&ret) & 0x7fff) == 0x7fff)
        *ExtendedConverter<sizeof(Extended)>::sexpPtr(&ret) = Accessor<16>::getRandomInt(state);

    // x87 extended format oddity: the first mantissa bit (always 1 for normal numbers) is visible, not hidden!
    if ((*ExtendedConverter<sizeof(Extended)>::sexpPtr(&ret) & 0x7fff) != 0)
        *ExtendedConverter<sizeof(Extended)>::mPtr(&ret) |= 0x8000000000000000ULL;

    return ret;
}


// Random in 1..2 - ideal IE case
template<> Extended Random12<true,false,Extended>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Work directly on Extended bits
    Extended r12;

    // Generate 63 bits for the mantissa
    *ExtendedConverter<sizeof(Extended)>::mPtr(&r12) = Accessor<64>::getRandomInt(state);

    // x87 extended format oddity: the first mantissa bit (always 1 for normal numbers) is visible, not hidden!
    // Since in this case this is a number between 1 and 2, this bit has to be set explicitly
    *ExtendedConverter<sizeof(Extended)>::mPtr(&r12) |= 0x8000000000000000ULL;

    // Set exponent so it's in the [1.0-2.0) range
    *ExtendedConverter<sizeof(Extended)>::sexpPtr(&r12) = 0x3FFF;

    // scale from 1-2 interval to the desired interval
    return r12;
}

// Random in 1..2 - near ideal EI case
template<> Extended Random12<false,true,Extended>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Work directly on Extended bits
    Extended r12;

    // Generate 63 bits for the mantissa
    *ExtendedConverter<sizeof(Extended)>::mPtr(&r12) = Accessor<64>::getRandomInt(state);

    // x87 extended format oddity: the first mantissa bit (always 1 for normal numbers) is visible, not hidden!
    // Since in this case this is a number between 1 and 2, this bit has to be set explicitly
    *ExtendedConverter<sizeof(Extended)>::mPtr(&r12) |= 0x8000000000000000LL;

    // Set exponent so it's in the (1.0-2.0] range
    // Replace 1.0 by 2.0
    if (*ExtendedConverter<sizeof(Extended)>::mPtr(&r12) == 0x8000000000000000LL)
        *ExtendedConverter<sizeof(Extended)>::sexpPtr(&r12) = 0x4000;
    else
        *ExtendedConverter<sizeof(Extended)>::sexpPtr(&r12) = 0x3FFF;

    return r12;
}

// Random in 1..2 - need to include both bounds
template<> Extended Random12<true,true,Extended>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Work directly on Extended bits
    Extended r12;

    // Generate 64 bits for the mantissa
    do {
        *ExtendedConverter<sizeof(Extended)>::mPtr(&r12) = Accessor<64>::getRandomInt(state);
    }
    // Include both bounds : repeat the random get till it's in the desired range
    // Keep the possibility for 2.0 in addition to all [1.0-2.0)
    while (*ExtendedConverter<sizeof(Extended)>::mPtr(&r12) > 0x8000000000000000LL);

    // Set exponent so it's in the [1.0-2.0] range
    // Check if 2.0 was selected
    if (*ExtendedConverter<sizeof(Extended)>::mPtr(&r12) == 0x8000000000000000LL)
        *ExtendedConverter<sizeof(Extended)>::sexpPtr(&r12) = 0x4000;
    // otherwise, set the correct mantissa bit and the correct exponent
    else {
        // x87 extended format oddity: the first mantissa bit (always 1 for normal numbers) is visible, not hidden!
        // Since in this case this is a number between 1 and 2, this bit has to be set explicitly
        *ExtendedConverter<sizeof(Extended)>::mPtr(&r12) |= 0x8000000000000000LL;
        *ExtendedConverter<sizeof(Extended)>::sexpPtr(&r12) = 0x3FFF;
    }

    return r12;
}

// Random in 1..2 - need to exclude both bounds
template<> Extended Random12<false,false,Extended>(RandomState& state) {

    // Get uniform number between 1 and 2 at max precision

    // Work directly on Extended bits
    Extended r12;

    // Generate 63 bits for the mantissa
    do {
        *ExtendedConverter<sizeof(Extended)>::mPtr(&r12) = Accessor<64>::getRandomInt(state);
        // x87 extended format oddity: the first mantissa bit (always 1 for normal numbers) is visible, not hidden!
        // Since in this case this is a number between 1 and 2, this bit has to be set explicitly
        *ExtendedConverter<sizeof(Extended)>::mPtr(&r12) |= 0x8000000000000000LL;
    }
    // Exclude both bounds : repeat the random get till it's in the desired range
    // Remove the possibility for 1.0 compared to all [1.0-2.0)
    while (*ExtendedConverter<sizeof(Extended)>::mPtr(&r12) == 0x8000000000000000LL);

    // Set exponent so it's in the (1.0-2.0) range
    *ExtendedConverter<sizeof(Extended)>::sexpPtr(&r12) = 0x3FFF;

    // scale from 1-2 interval to the desired interval
    return r12;
}

#endif

// This is a way to hide the implementation from the header
// And also to ensure there is only one template instanciation, instead of duplicating
// the code in all object files
template<typename FloatType> static inline FloatType NRandom_Generic(FloatType *secondary, RandomState& state) {

    FloatType x, y, d;
    // Generate a point strictly inside the unit circle
    do {
        // Use IE as this is the most convenient for the generation,
        // any will do since the bounds are rejected anyway by unit circle strict comparison
        x = RandomIE(FloatType(-1.0), FloatType(1.0), state);
        y = RandomIE(FloatType(-1.0), FloatType(1.0), state);
        d = x*x + y*y;
    } while (d>=FloatType(1.0));
    // Convert from x/y to uniform
    FloatType conv = sqrt( FloatType(-2.0) * log(d) / d );
    // Use one result as secondary independent variable, if user wishes
    if (secondary) *secondary = x * conv;
    // return the other
    return y * conv;
}

// Specialize for the Float types declared in the header
template<> Simple NRandom(Simple *secondary, RandomState& state) {
    return NRandom_Generic<Simple>(secondary,state);
}

/*
template<> Double NRandom(Double *secondary, RandomState& state) {
    return NRandom_Generic<Double>(secondary,state);
}
#if defined(Extended)
template<> Extended NRandom(Extended *secondary, RandomState& state) {
    return NRandom_Generic<Extended>(secondary,state);
}
#endif
*/

// May save one mul
template<typename FloatType> static inline FloatType NRandom_Generic(FloatType mean, FloatType std_dev, FloatType *secondary, RandomState& state) {

    FloatType x, y, d;
    // Generate a point strictly inside the unit circle
    do {
        // Use IE as this is the most convenient for the generation,
        // any will do since the bounds are rejected anyway by unit circle strict comparison
        x = RandomIE(FloatType(-1.0), FloatType(1.0), state);
        y = RandomIE(FloatType(-1.0), FloatType(1.0), state);
        d = x*x + y*y;
    } while (d>=FloatType(1.0));
    // Convert from x/y to uniform
    FloatType conv = sqrt( FloatType(-2.0) * log(d) / d ) * std_dev;
    // Use one result as secondary independent variable, if user wishes
    if (secondary) *secondary = x * conv + mean;
    // return the other
    return y * conv + mean;
}

// Specialize for the Float types declared in the header
template<> Simple NRandom(Simple mean, Simple std_dev, Simple *secondary, RandomState& state) {
    return NRandom_Generic<Simple>(mean, std_dev, secondary,state);
}

/*
template<> Double NRandom(Double mean, Double std_dev, Double *secondary, RandomState& state) {
    return NRandom_Generic<Double>(mean, std_dev, secondary,state);
}
#if defined(Extended)
template<> Extended NRandom(Extended mean, Extended std_dev, Extended *secondary, RandomState& state) {
    return NRandom_Generic<Extended>(mean, std_dev, secondary,state);
}
#endif
*/

SizedUnsignedInteger<32>::Type RandomInit(RandomState& state) {
    return RandomInit(SizedUnsignedInteger<32>::Type(time(0)));
}

SizedUnsignedInteger<32>::Type RandomInit(SizedUnsignedInteger<32>::Type seed, RandomState& state) {
    init_genrand(seed,state);
    return state.seed;
}

SizedUnsignedInteger<32>::Type RandomSeed(RandomState& state) {
    return state.seed;
}

// Default state holder, so single threaded applications don't bother setting up an object
RandomState DefaultRandomState;

} // end streflop namespace