// Ranlux24 PRNG. (1.06)
// Chaotic with long period but slow due to discard.

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

// Utility.

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

#define UNLIKELY(cond) \
    __builtin_expect(!!(cond), 0)

#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;
}

// Engine Base.

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

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

void engineBase_discard(EngineBase *self, long n)
{
    return self->discard(self, n);
}

void engineBase_delete(EngineBase *self)
{
    return self->delete(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;
}

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

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

    if (UNLIKELY(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];
}

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

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

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

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.next = subtractWithCarryEngine_next;
    self->base.discard = subtractWithCarryEngine_discard;
    self->base.delete = subtractWithCarryEngine_delete;
    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] = nrand48(state) & self->m;
    return (EngineBase *) self;
}

// Discard Block.

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

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

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

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

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

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

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

    DiscardBlockEngine *self = NEW_T(DiscardBlockEngine, 1);

    self->base.next = discardBlockEngine_next;
    self->base.discard = discardBlockEngine_discard;
    self->base.delete = discardBlockEngine_delete;
    self->owned = *unique; *unique = 0; // Transfer ownership.
    self->p = p;
    self->r = r;
    self->i = r;
    return (EngineBase *) self;
}

// Ranlux24.

EngineBase *ranlux24_new(unsigned long seed)
{
    EngineBase *ranlux24_base = subtractWithCarryEngine_new(24, 10, 24, seed);
    return discardBlockEngine_new(&ranlux24_base, 223, 23);
}

// 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;
}

int main(void)
{
    unsigned long seed = time(0);
    EngineBase *engine = ranlux24_new(seed);

    for (int i = 0; i < 24; i++)
    {
        long r = engineBase_next(engine);
        printf("%ld\n", r);
    }

    long n = 10000000;
    double t = clock_now();
    engineBase_discard(engine, n);
    printf("Elapsed: %.9fs\n", clock_now() - t);

    engineBase_delete(engine);
    return 0;
}