/*
 * Random numbers generators.
 * 2021-09-25 Udo Krebelder <udo@krebelder.at>
 *
 * Compile with 64 bit WinDDK 7.10 compiler:
 *   cl -nologo -I../../MiniCRT/current -W3 -O2 -Ob1 -GS- -c rand.c
 */

#include "winutil.h"

#ifdef _MSC_VER
  #define rotl(x, r) _rotl64(x, r)
#else
  #define rotl(x, r) (((x) << r) | ((x) >> (64 - r)))
#endif

/*
 * This is Java 8 SplittableRandom generator without the split feature.
 * http://dx.doi.org/10.1145/2714064.2660195
 * http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
 * https://www.pcg-random.org/posts/bugs-in-splitmix.html
 * https://www.pcg-random.org/posts/birthday-test.html
 *
 * It is a very fast generator passing TestU01 BigCrush,
 * and it can be useful if you absolutely want 64 bits of state.
 * It fails the Birthday test (but values are not unique).
 */
STATIC INLINE UINT64 SplitmixExInline(SPLITMIX_STATE *state)
{
    UINT64 x;

    state->data += 0x9e3779b97f4a7c15;
    x = state->data;

    x ^= x >> 30;
    x *= 0xbf58476d1ce4e5b9;
    x ^= x >> 27;
    x *= 0x94d049bb133111eb;
    x ^= x >> 31;

    return x;
}


UINT64 SplitmixRandEx(SPLITMIX_STATE *state)
{
    return SplitmixExInline(state);
}


/* Default splitmix state */
#ifndef SPLITMIX_ZERO_SEED
#define SPLITMIX_ZERO_SEED 0x0992ccaf6a6fca05  /* 31 ones */
#endif


VOID SplitmixSeedEx(SPLITMIX_STATE *state, UINT64 seed)
{
    state->data = seed ^ SPLITMIX_ZERO_SEED;
}


/*
 * State for splitmix.  The state can be seeded with any value.
 */
SPLITMIX_STATE g_SplitmixState = { SPLITMIX_ZERO_SEED };


UINT64 SplitmixRand(VOID)
{
    return SplitmixExInline(&g_SplitmixState);
}


VOID SplitmixSeed(UINT64 seed)
{
    SplitmixSeedEx(&g_SplitmixState, seed);
}


/*
 * xoroshiro128++
 * https://prng.di.unimi.it/
 *
 * This is xoshiro128++ 1.0, the default random number generator.
 * It has excellent speed, and a state size of 128 bits.
 *
 * The state must be seeded so that it is not everywhere zero.
 * If you have a 64-bit seed, we suggest to seed a splitmix generator
 * and use its output to fill state.
 */
STATIC INLINE UINT64 RandExInline(RAND_STATE *state)
{
    UINT64 s0, s1, t, result;

    s0 = state->data[0];
    s1 = state->data[1];
    t = s0 + s1;
    result = rotl(t, 17) + s0;

    s1 ^= s0;
    state->data[0] = rotl(s0, 49) ^ s1 ^ (s1 << 21);
    state->data[1] = rotl(s1, 28);

    return result;
}


UINT64 RandEx(RAND_STATE *state)
{
    return RandExInline(state);
}


VOID RandSeedEx(RAND_STATE *state, UINT64 seed)
{
    SPLITMIX_STATE s;

    SplitmixSeedEx(&s, seed);
    state->data[0] = SplitmixRandEx(&s);
    state->data[1] = SplitmixRandEx(&s);
}


/*
 * This is the jump function for the generator.
 * It is equivalent to 2^64 calls to RandEx(state).
 * It can be used to generate 2^64 non-overlapping
 * subsequences for parallel computations.
 */
VOID RandJumpEx(RAND_STATE *state)
{
    STATIC CONST UINT64 JUMP[] = {
        0x2bd7a6a6e99c2ddc, 0x0992ccaf6a6fca05
    };

    UINT   i, b;
    UINT64 s0 = 0;
    UINT64 s1 = 0;

    for (i = 0; i < COUNTOF(JUMP); i++) {
        for (b = 0; b < 64; b++) {
            if (JUMP[i] & (UINT64)1 << b) {
                s0 ^= state->data[0];
                s1 ^= state->data[1];
            }
            RandEx(state);
        }
    }

    state->data[0] = s0;
    state->data[1] = s1;
}


/*
 * State for Rand. The state should be seeded with
 * equal number of zeros and ones.
 */
RAND_STATE g_RandState = {
    0x2bd7a6a6e99c2ddc, 0x0992ccaf6a6fca05
};


UINT64 Rand(VOID)
{
    return RandExInline(&g_RandState);
}


VOID RandSeed(UINT64 seed)
{
    RandSeedEx(&g_RandState, seed);
}


VOID RandJump(VOID)
{
    RandJumpEx(&g_RandState);
}


/*
 * Generate a random number between min and max-1 from the
 * uniform distributed 64 bit random number argument.
 */
STATIC INLINE UINT64 RandRangeEx_Inline(UINT64 rnd_num, UINT64 min, UINT64 max)
{
    UINT64 hi, lo;

    lo = _umul128(rnd_num, max, &hi);
    UNUSED(lo);
    return hi + min;
}


UINT64 RandRangeEx(UINT64 rnd_num, UINT64 min, UINT64 max)
{
    return RandRangeEx_Inline(rnd_num, min, max);
}


UINT64 RandRange(UINT64 min, UINT64 max)
{
    UINT64 rnd_num = RandExInline(&g_RandState);
    return RandRangeEx_Inline(rnd_num, min, max);
}


/* For compatibiliy with random */
INT random(VOID)
{
    return (UINT)RandRange(0, 0x7fffffff);
}


void srandom(UINT seed)
{
    RandSeed(seed);
}


/*
 * Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs".
 * The state word must be initialized to non-zero.
 * Note that the generator fail the Birthday test because
 * values do not repeat, which is sometimes useful.
 */
INLINE UINT XorShift32ExInline(XORSHIFT32_STATE *state)
{
    UINT x = state->data;

    x ^= x << 13;
    x ^= x >> 17;
    x ^= x << 5;
    state->data = x;

    return x;
}


UINT XorShift32RandEx(XORSHIFT32_STATE *state)
{
    return XorShift32ExInline(state);
}


#ifndef XORSHIFT32_ZERO_SEED
#define XORSHIFT32_ZERO_SEED 0x6a6fca05
#endif


VOID XorShift32SeedEx(XORSHIFT32_STATE *state, UINT seed)
{
    state->data = seed ^ XORSHIFT32_ZERO_SEED;
}


/*
 * State for XorShift32.
 * The state should be seeded something not zero.
 */
XORSHIFT32_STATE g_XorShiftState = { XORSHIFT32_ZERO_SEED };


VOID XorShift32Seed(UINT seed)
{
    XorShift32SeedEx(&g_XorShiftState, seed);
}


UINT XorShift32Rand(VOID)
{
    return XorShift32ExInline(&g_XorShiftState);
}


#if 0
/*
 * 64 bit XorShiftStar generator
 * https://de.wikipedia.org/wiki/Xorshift
 * https://en.wikipedia.org/wiki/xorshift
 * https://stackoverflow.com/questions/53886131/how-does-xorshift32-works
 * https://www.pcg-random.org/posts/birthday-test.html
 */

typedef struct {
    UINT64      data;
} XORSHIFTSTAR_STATE;


UINT64 XorShiftStarEx(XORSHIFTSTAR_STATE *state)
{
    UINT64 x = state->data;

    x ^= x >> 12;
    x ^= x << 25;
    x ^= x >> 27;
    state->data = x;

    return x * 0x2545f4914f6cdd1d;
}
#endif


#if 0
/*
 * Implementation of a 32-bit KISS generator which uses no multiply instructions
 * https://eprint.iacr.org/2011/007.pdf
 */ 
UINT JKISS32(VOID)
{
    STATIC UINT x = 123456789, y = 234567891, z = 345678912, w = 456789123, c = 0; 

    UINT t;

    y ^= y << 5;
    y ^= y >> 7;
    y ^= y << 22; 

    t = z + w + c;
    z = w;
    c = t < 0;
    w = t & 2147483647; 

    x += 1411392427; 

    return x + y + w; 
}
#endif


/*
 * KISS - George Marsaglia
 * https://de.wikipedia.org/wiki/KISS_(Zufallszahlengenerator)
 * https://www.thecodingforums.com/threads/64-bit-kiss-rngs.673657/
 * 
 */

#if 0
#include <stdio.h>

static UINT64 x=1234567890987654321ULL,
              c=123456123456123456ULL,
              y=362436362436362436ULL,
              z=1066149217761810ULL, t;

#define MWC (t=(x<<58)+c, c=(x>>6), x+=t, c+=(x<t), x)
#define XSH ( y^=(y<<13), y^=(y>>17), y^=(y<<43) )
#define CNG ( z=6906969069LL*z+1234567 )
#define KISS (MWC+XSH+CNG)

int main(void)
{
    int i;

    for(i=0; i < 100000000; i++)
        t = KISS;
    (t == 1666297717051644203ULL) ?
        printf("100 million uses of KISS OK"): printf("Fail");
}
#endif


#if 0
typedef struct {
    UINT64 x, y, z, c;
} KISS_STATE;


KISS_STATE g_KissState = {
    1066149217761810,
    362436362436362436,
    1234567890987654321,
    123456123456123456,
};


UINT64 Kiss64(VOID)
{
    KISS_STATE *s = &g_KissState;
    UINT64 t;

    /* Linearer Kongruenzgenerator */
    s->x = 6906969069ULL * s->x + 1234567;

    /* Xorshift */
    s->y ^= s->y << 13;
    s->y ^= s->y >> 17;
    s->y ^= s->y << 43;

    /* Multiply-with-carry */
    t = (s->z << 58) + s->c;
    s->c = s->z >> 6;
    s->z += t;
    s->c += s->z < t;

    return s->x + s->y + s->z;
}
#endif
