// Ranlux24 PRNG. (2.01)
// 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;
}

// Ranlux24.

typedef struct {
    int *x, c, i, j;
} Ranlux24;

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

static int subtract_with_carry_next(Ranlux24 *self)
{
    self->i += 1;
    int *x = self->x, i = self->i;

    if (UNLIKELY(i >= 24))
    {
        int c = self->c;
        for (i = 0; i < 10; i++)
            x[i] = next(x, i + 14, i, &c);
        for (i = 10; i < 24; i++)
            x[i] = next(x, i - 10, i, &c);
        self->c = c;
        self->i = i = 0;
    }
    return x[i];
}

static void subtract_with_carry_discard(Ranlux24 *self, int n)
{
    for (int i = 0; i < n; i++)
        subtract_with_carry_next(self);
}

int ranlux24_next(Ranlux24 *self)
{
    if (self->j == 0)
    {
        subtract_with_carry_discard(self, 223 - 23);
        self->j = 23;
    }
    self->j -= 1;
    return subtract_with_carry_next(self);
}

void ranlux24_discard(Ranlux24 *self, int n)
{
    for (int i = 0; i < n; i++)
        ranlux24_next(self);
}

void ranlux24_delete(Ranlux24 *self)
{
    if (self != 0)
    {
        free(self->x);
        free(self);
    }
}

Ranlux24 *ranlux24_new(unsigned int seed)
{
    Ranlux24 *self = NEW_T(Ranlux24, 1);
    self->x = NEW_T(int, 24);
    self->c = 0;
    self->i = 24-1;
    self->j = 23;

    unsigned short rs[] = {0x330E, seed, seed>>16};
    for (int i = 0; i < 24; i++)
        self->x[i] = nrand48(rs) & 0xFFFFFF;
    return self;
}

// 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 int seed = time(0);
    Ranlux24 *engine = ranlux24_new(seed);

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

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

    ranlux24_delete(engine);
    return 0;
}