source: libcfa/src/bits/random.hfa @ dd46fd3

// 
// Cforall Version 1.0.0 Copyright (C) 2022 University of Waterloo
// 
// The contents of this file are covered under the licence agreement in the
// file "LICENCE" distributed with Cforall.
//
// random.hfa -- 
// 
// Author           : Peter A. Buhr
// Created On       : Fri Jan 14 07:18:11 2022
// Last Modified By : Peter A. Buhr
// Last Modified On : Wed Nov 30 18:32:25 2022
// Update Count     : 111
// 

#pragma once

#include <stdint.h>

#define GLUE2( x, y ) x##y
#define GLUE( x, y ) GLUE2( x, y )

// Set default PRNG for architecture size.
#ifdef __x86_64__                                                                               // 64-bit architecture
#define LEHMER64
#define XORSHIFT_6_21_7
//#define XOSHIRO256PP
//#define XOSHIRO128PP
#else                                                                                                   // 32-bit architecture
#define LEHMER64
#define XORSHIFT_6_21_7
#endif // __x86_64__

// C/CFA PRNG name and random-state.

#ifdef LEHMER64
#define PRNG_NAME_64 lehmer64
#define PRNG_STATE_64_T __uint128_t
#endif // LEHMER64

#ifdef XORSHIFT_6_21_7
#define PRNG_NAME_32 xorshift_6_21_7
#define PRNG_STATE_32_T uint32_t
#endif // XORSHIFT_6_21_7

#ifdef XOSHIRO256PP
#define PRNG_NAME_64 xoshiro256pp
#define PRNG_STATE_64_T struct GLUE(PRNG_NAME_64,_t)
PRNG_STATE_64_T { uint64_t s[4]; };
#endif // XOSHIRO256PP

#ifdef XOSHIRO128PP
#define PRNG_NAME_32 xoshiro128pp
#define PRNG_STATE_32_T struct GLUE(PRNG_NAME_32,_t)
PRNG_STATE_32_T { uint32_t s[4]; };
#endif // XOSHIRO128PP

#define PRNG_SET_SEED_64 GLUE(PRNG_NAME_64,_set_seed)
#define PRNG_SET_SEED_32 GLUE(PRNG_NAME_32,_set_seed)


// Default PRNG used by runtime.
#ifdef __x86_64__                                                                               // 64-bit architecture
#define PRNG_NAME PRNG_NAME_64
#define PRNG_STATE_T PRNG_STATE_64_T
#else                                                                                                   // 32-bit architecture
#define PRNG_NAME PRNG_NAME_32
#define PRNG_STATE_T PRNG_STATE_32_T
#endif // __x86_64__

#define PRNG_SET_SEED GLUE(PRNG_NAME,_set_seed)


#ifdef __cforall                                                                                // don't include in C code (invoke.h)

// https://prng.di.unimi.it/xoshiro128plusplus.c
//
// This is xoshiro128++ 1.0, one of our 32-bit all-purpose, rock-solid generators. It has excellent speed, a state size
// (128 bits) that is large enough for mild parallelism, and it passes all tests we are aware of.
//
// For generating just single-precision (i.e., 32-bit) floating-point numbers, xoshiro128+ is even faster.
//
// The state must be seeded so that it is not everywhere zero.

#ifndef XOSHIRO128PP
struct xoshiro128pp_t { uint32_t s[4]; };
#endif // ! XOSHIRO128PP

static inline uint32_t xoshiro128pp( xoshiro128pp_t & rs ) with(rs) {
        inline uint32_t rotl( const uint32_t x, int k ) {
                return (x << k) | (x >> (32 - k));
        }

        const uint32_t result = rotl( s[0] + s[3], 7 ) + s[0];
        const uint32_t t = s[1] << 9;

        s[2] ^= s[0];
        s[3] ^= s[1];
        s[1] ^= s[2];
        s[0] ^= s[3];
        s[2] ^= t;
        s[3] = rotl( s[3], 11 );
        return result;
}

static inline void xoshiro128pp_set_seed( xoshiro128pp_t & state, uint32_t seed ) {
        state = (xoshiro128pp_t){ {seed, seed, seed, seed} };
} // xoshiro128pp_set_seed

// This is xoshiro256++ 1.0, one of our all-purpose, rock-solid generators.  It has excellent (sub-ns) speed, a state
// (256 bits) that is large enough for any parallel application, and it passes all tests we are aware of.
//
// For generating just floating-point numbers, xoshiro256+ is even faster.
//
// The state must be seeded so that it is not everywhere zero. If you have a 64-bit seed, we suggest to seed a
// splitmix64 generator and use its output to fill s.

#ifndef XOSHIRO256PP
struct xoshiro256pp_t { uint64_t s[4]; };
#endif // ! XOSHIRO256PP

static inline uint64_t xoshiro256pp( xoshiro256pp_t & rs ) with(rs) {
        inline uint64_t rotl(const uint64_t x, int k) {
                return (x << k) | (x >> (64 - k));
        }

        const uint64_t result = rotl( s[0] + s[3], 23 ) + s[0];
        const uint64_t t = s[1] << 17;

        s[2] ^= s[0];
        s[3] ^= s[1];
        s[1] ^= s[2];
        s[0] ^= s[3];
        s[2] ^= t;
        s[3] = rotl( s[3], 45 );
        return result;
}

static inline void xoshiro256pp_set_seed( xoshiro256pp_t & state,  uint64_t seed ) {
        state = (xoshiro256pp_t){ {seed, seed, seed, seed} };
} // xoshiro256pp_set_seed

#ifdef __SIZEOF_INT128__
        // Pipelined to allow out-of-order overlap with reduced dependencies. Critically, the current random state is
        // returned (copied), and then compute and store the next random value.
        //--------------------------------------------------
        static inline uint64_t lehmer64( __uint128_t & state ) {
                __uint128_t ret = state;
                state *= 0xda942042e4dd58b5;
                return ret >> 64;
        } // lehmer64

        static inline void lehmer64_set_seed( __uint128_t & state, uint64_t seed ) {
                state = seed;
        } // lehmer64_set_seed

        //--------------------------------------------------
        static inline uint64_t wyhash64( uint64_t & state ) {
                state += 0x60bee2bee120fc15;
                __uint128_t tmp;
                tmp = (__uint128_t) state * 0xa3b195354a39b70d;
                uint64_t m1 = (tmp >> 64) ^ tmp;
                tmp = (__uint128_t)m1 * 0x1b03738712fad5c9;
                uint64_t m2 = (tmp >> 64) ^ tmp;
                return m2;
        }

        static inline void wyhash64_set_seed( __uint128_t & state, uint64_t seed ) {
                state = seed;
        } // lehmer64_set_seed
#endif // __SIZEOF_INT128__

//--------------------------------------------------
static inline uint64_t xorshift_13_7_17( uint64_t & state ) {
        uint64_t ret = state;
        state ^= state << 13;
        state ^= state >> 7;
        state ^= state << 17;
        return ret;
}

static inline void xorshift_13_7_17_set_seed( uint64_t & state, uint32_t seed ) {
        state = seed;
}

//--------------------------------------------------
static inline uint32_t xorshift_6_21_7( uint32_t & state ) {
        uint32_t ret = state;
        state ^= state << 6;
        state ^= state >> 21;
        state ^= state << 7;
        return ret;
} // xorshift_6_21_7

static inline void xorshift_6_21_7_set_seed( uint32_t & state, uint32_t seed ) {
        state = seed;
}

//--------------------------------------------------
typedef struct {
        uint32_t a, b, c, d;
        uint32_t counter;
} xorwow__state_t;

// The state array must be initialized to not be all zero in the first four words.
static inline uint32_t xorwow( xorwow__state_t & state ) {
        // Algorithm "xorwow" from p. 5 of Marsaglia, "Xorshift RNGs".
        uint32_t ret = state.a + state.counter;
        uint32_t t = state.d;

        uint32_t const s = state.a;
        state.d = state.c;
        state.c = state.b;
        state.b = s;

        t ^= t >> 2;
        t ^= t << 1;
        t ^= s ^ (s << 4);
        state.a = t;

        state.counter += 362437;
        return ret;
}

// Used in __tls_rand_fwd
//--------------------------------------------------
#define M  (1_l64u << 48_l64u)
#define A  (25214903917_l64u)
#define AI (18446708753438544741_l64u)
#define C  (11_l64u)
#define D  (16_l64u)

// Bi-directional LCG random-number generator
static inline uint32_t LCGBI_fwd( uint64_t & state ) {
        state = (A * state + C) & (M - 1);
        return state >> D;
}

static inline uint32_t LCGBI_bck( uint64_t & state ) {
        unsigned int r = state >> D;
        state = AI * (state - C) & (M - 1);
        return r;
}

#undef M
#undef A
#undef AI
#undef C
#undef D

#endif // __cforall