// Pseudorandom number generation. (1.10)

#include <stdlib.h>
#include <limits.h>
#include <assert.h>
#include <stdio.h>
#include <time.h>

//
// Utility.
//

#define BITS(x) \
    ((int) sizeof(x) * CHAR_BIT)

#define NEW_T(T, n) \
    ((T *) allocate(n * sizeof(T)))

void *allocate(size_t n)
{
    void *r = malloc(n);
    assert(r != 0 || n == 0);
    return r;
}

unsigned long uniform_ulong(unsigned short *state)
{
    unsigned long r = 0;
    for (int i = 0; i < BITS(r); i += 32)
        r += (unsigned long) jrand48(state) << i;
    return r;
}

//
// Engine Base.
//

typedef struct {
    struct EngineBaseMethod *method;
} EngineBase;

struct EngineBaseMethod {
    long (*next)    (EngineBase *);
    void (*discard) (EngineBase *, long);
    void (*delete)  (EngineBase *);
};

long engine_next(EngineBase *self)
{
    return self->method->next(self);
}

void engine_discard(EngineBase *self, long n)
{
    self->method->discard(self, n);
}

void engine_delete(EngineBase *self)
{
    self->method->delete(self);
}

//
// Linear Congruential.
//

typedef struct {
    EngineBase base;
    long x, a, c, m; int d;
} LinearCongruentialEngine;

static long linearCongruentialEngine_next(EngineBase *base)
{
    LinearCongruentialEngine *self = (LinearCongruentialEngine *) base;

    self->x = ((unsigned long long) self->x * self->a + self->c) % self->m;
    return self->x >> self->d;
}

static void linearCongruentialEngine_discard(EngineBase *base, long n)
{
    for (long i = 0; i < n; i++)
        linearCongruentialEngine_next(base);
}

static void linearCongruentialEngine_delete(EngineBase *base)
{
    free(base);
}

static struct EngineBaseMethod linearCongruentialEngine_method = {
    linearCongruentialEngine_next,
    linearCongruentialEngine_discard,
    linearCongruentialEngine_delete
};

EngineBase *linearCongruentialEngine_new(long a, long c, long m, int d,
        unsigned long seed)
{
    assert(m != 0);
    assert(0 <= d && d < BITS(long));

    LinearCongruentialEngine *self = NEW_T(LinearCongruentialEngine, 1);

    self->base.method = &linearCongruentialEngine_method;
    self->x = seed % m;
    self->a = a;
    self->c = c;
    self->m = m;
    self->d = d;
    return (EngineBase *) self;
}

//
// Subtract With Carry.
//

typedef struct {
    EngineBase base;
    long *x, m; int c, i, r, s;
} SubtractWithCarryEngine;

static inline long next(long *x, long mask, int i_s, int i_r, int *carry)
{
    long y = x[i_s] - x[i_r] - *carry;
    *carry = -(y >> (BITS(y) - 1));
    return y & mask;
}

static long subtractWithCarryEngine_next(EngineBase *base)
{
    SubtractWithCarryEngine *self = (SubtractWithCarryEngine *) base;

    self->i += 1;
    long *x = self->x;
    int   i = self->i, r = self->r;

    if (i == r)
    {
        long m = self->m;
        int  c = self->c, s = self->s, t = r - s;

        for (i = 0; i < s; i++)
            x[i] = next(x, m, i + t, i, &c);
        for (i = s; i < r; i++)
            x[i] = next(x, m, i - s, i, &c);
        self->c = c;
        self->i = i = 0;
    }
    return x[i];
}

static void subtractWithCarryEngine_discard(EngineBase *base, long n)
{
    for (long i = 0; i < n; i++)
        subtractWithCarryEngine_next(base);
}

static void subtractWithCarryEngine_delete(EngineBase *base)
{
    SubtractWithCarryEngine *self = (SubtractWithCarryEngine *) base;

    if (self != 0)
    {
        free(self->x);
        free(self);
    }
}

static struct EngineBaseMethod subtractWithCarryEngine_method = {
    subtractWithCarryEngine_next,
    subtractWithCarryEngine_discard,
    subtractWithCarryEngine_delete
};

EngineBase *subtractWithCarryEngine_new(int w, int s, int r, unsigned long seed)
{
    assert(0 < w && w < BITS(long));
    assert(0 < s && s < r);

    SubtractWithCarryEngine *self = NEW_T(SubtractWithCarryEngine, 1);

    self->base.method = &subtractWithCarryEngine_method;
    self->x = NEW_T(long, r);
    self->m = -1UL >> (BITS(long) - w);
    self->c = 0;
    self->i = r-1;
    self->r = r;
    self->s = s;

    unsigned short state[3] = {0x330E, seed, seed>>16};

    for (int i = 0; i < r; i++)
        self->x[i] = uniform_ulong(state) & self->m;
    return (EngineBase *) self;
}

//
// Discard Block.
//

typedef struct {
    EngineBase base, *owned;
    int p, r, i;
} DiscardBlockEngine;

static long discardBlockEngine_next(EngineBase *base)
{
    DiscardBlockEngine *self = (DiscardBlockEngine *) base;

    if (self->i == 0)
    {
        engine_discard(self->owned, self->p - self->r);
        self->i = self->r;
    }
    self->i -= 1;
    return engine_next(self->owned);
}

static void discardBlockEngine_discard(EngineBase *base, long n)
{
    for (long i = 0; i < n; i++)
        discardBlockEngine_next(base);
}

static void discardBlockEngine_delete(EngineBase *base)
{
    DiscardBlockEngine *self = (DiscardBlockEngine *) base;

    if (self != 0)
    {
        engine_delete(self->owned);
        free(self);
    }
}

static struct EngineBaseMethod discardBlockEngine_method = {
    discardBlockEngine_next,
    discardBlockEngine_discard,
    discardBlockEngine_delete,
};

EngineBase *discardBlockEngine_new(EngineBase **unique, int p, int r)
{
    assert(0 < r && r <= p);

    DiscardBlockEngine *self = NEW_T(DiscardBlockEngine, 1);

    self->base.method = &discardBlockEngine_method;
    self->owned = *unique; *unique = 0; // Transfer ownership.
    self->p = p;
    self->r = r;
    self->i = r;
    return (EngineBase *) self;
}

//
// Predefined Engines.
//

EngineBase *default_rand_new(unsigned long seed)
{
    assert(BITS(long) >= 50);

    // 32-bit with a period of 2^48.
    seed = seed * 0x10000L + 0x330E;
    return linearCongruentialEngine_new(25214903917L, 11, 1L << 48, 16, seed);
}

EngineBase *minstd_rand_new(unsigned long seed)
{
    // 31-bit on the closed interval [1, 2^31-2].
    if (seed % 2147483647L == 0)
        seed = 1;
    return linearCongruentialEngine_new(48271L, 0, 2147483647L, 0, seed);
}

EngineBase *ranlux24_new(unsigned long seed)
{
    // 24-bit chaotic with long period.
    EngineBase *ranlux24_base = subtractWithCarryEngine_new(24, 10, 24, seed);
    return discardBlockEngine_new(&ranlux24_base, 223, 23);
}

EngineBase *ranlux48_new(unsigned long seed)
{
    // 48-bit chaotic with long period.
    EngineBase *ranlux48_base = subtractWithCarryEngine_new(48, 5, 12, seed);
    return discardBlockEngine_new(&ranlux48_base, 389, 11);
}

//
// Main.
//

double clock_now(void)
{
    struct timespec now;
    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &now);
    return now.tv_sec + now.tv_nsec / 1.0E+09;
}

void show(int n, EngineBase *(*new_)(unsigned long), unsigned long seed)
{
    EngineBase *engine = new_(seed);

    for (int i = 0; i < n; i++)
        printf("%ld\n", engine_next(engine));

    double t = clock_now();
    engine_discard(engine, 1000000L);
    printf("Elapsed: %.9fs\n", clock_now() - t);

    engine_delete(engine);
}

int main(void)
{
    unsigned long seed = time(0);
    show(8, default_rand_new, seed);
    show(8, minstd_rand_new, seed);
    show(8, ranlux24_new, seed);
    show(8, ranlux48_new, seed);
    return 0;
}