// Ranlux24 PRNG. (1.03)
// 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)
{
return r;
}
// Engine Base.
typedef struct EngineBase {
int (*next) (struct EngineBase*);
void (*discard) (struct EngineBase*, int);
void (*delete) (struct EngineBase*);
} EngineBase;
int EngineBase_next(EngineBase *self)
{
return self->next(self);
}
void EngineBase_discard(EngineBase *self, int n)
{
return self->discard(self, n);
}
void EngineBase_delete(EngineBase *self)
{
return self->delete(self);
}
// Subtract With Carry.
typedef struct {
EngineBase base;
int *x, c, i, m, r, s;
} SubtractWithCarryEngine;
static inline int next(int *x, int mask, int i_s, int i_r, int *carry)
{
int y = x[i_s] - x[i_r] - *carry;
*carry = -(y >> (BITS(int) - 1));
return y & mask;
}
int SubtractWithCarryEngine_next(EngineBase *base)
{
SubtractWithCarryEngine *self = (SubtractWithCarryEngine *) base;
self->i += 1;
int *x = self->x, i = self->i, r = self->r;
if (UNLIKELY(i >= r))
{
int c = self->c, m = self->m, 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, int n)
{
for (int i = 0; i < n; i++)
SubtractWithCarryEngine_next(base);
}
void SubtractWithCarryEngine_delete(EngineBase *base)
{
SubtractWithCarryEngine *self = (SubtractWithCarryEngine *) base;
if (self != 0)
{
}
}
EngineBase *SubtractWithCarryEngine_new(int w, int s, int r)
{
assert(0 < w
&& w
< BITS
(int));
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(int, r);
self->c = 0;
self->i = r-1;
self->m = -1U >> (BITS(int) - w);
self->r = r;
self->s = s;
for (int i = 0; i < r; i++)
self
->x
[i
] = rand() & self
->m
; return (EngineBase *) self;
}
// Discard Block.
typedef struct {
EngineBase base, *owned;
int p, r, i;
} DiscardBlockEngine;
int 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, int n)
{
for (int i = 0; i < n; i++)
DiscardBlockEngine_next(base);
}
void DiscardBlockEngine_delete(EngineBase *base)
{
DiscardBlockEngine *self = (DiscardBlockEngine *) base;
if (self != 0)
{
EngineBase_delete(self->owned);
}
}
EngineBase *DiscardBlockEngine_new(EngineBase **unique, int p, int r)
{
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(void)
{
EngineBase *ranlux24_base = SubtractWithCarryEngine_new(24, 10, 24);
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)
{
EngineBase *ranlux24 = Ranlux24_new();
for (int i = 0; i < 24; i++)
{
int r = EngineBase_next(ranlux24);
}
int n = 10000000;
double t = clock_now();
EngineBase_discard(ranlux24, n);
printf("Elapsed: %.9fs\n", clock_now
() - t
);
EngineBase_delete(ranlux24);
return 0;
}
Ly8gUmFubHV4MjQgUFJORy4gKDEuMDMpCi8vIENoYW90aWMgd2l0aCBsb25nIHBlcmlvZCBidXQgc2xvdyBkdWUgdG8gZGlzY2FyZC4KCiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPGxpbWl0cy5oPgojaW5jbHVkZSA8YXNzZXJ0Lmg+CiNpbmNsdWRlIDxzdGRpby5oPgojaW5jbHVkZSA8dGltZS5oPgoKLy8gVXRpbGl0eS4KCiNkZWZpbmUgQklUUyh4KSBcCiAgICAoc2l6ZW9mKHgpICogQ0hBUl9CSVQpCgojZGVmaW5lIFVOTElLRUxZKGNvbmQpIFwKICAgIF9fYnVpbHRpbl9leHBlY3QoISEoY29uZCksIDApCgojZGVmaW5lIE5FV19UKFQsIG4pIFwKICAgICgoVCAqKSBhbGxvY2F0ZShuICogc2l6ZW9mKFQpKSkKCnZvaWQgKmFsbG9jYXRlKHNpemVfdCBuKQp7CiAgICB2b2lkICpyID0gbWFsbG9jKG4pOwogICAgYXNzZXJ0KHIgIT0gMCB8fCBuID09IDApOwogICAgcmV0dXJuIHI7Cn0KCi8vIEVuZ2luZSBCYXNlLgoKdHlwZWRlZiBzdHJ1Y3QgRW5naW5lQmFzZSB7CiAgICBpbnQgICgqbmV4dCkgICAgKHN0cnVjdCBFbmdpbmVCYXNlKik7CiAgICB2b2lkICgqZGlzY2FyZCkgKHN0cnVjdCBFbmdpbmVCYXNlKiwgaW50KTsKICAgIHZvaWQgKCpkZWxldGUpICAoc3RydWN0IEVuZ2luZUJhc2UqKTsKfSBFbmdpbmVCYXNlOwoKaW50IEVuZ2luZUJhc2VfbmV4dChFbmdpbmVCYXNlICpzZWxmKQp7CiAgICByZXR1cm4gc2VsZi0+bmV4dChzZWxmKTsKfQoKdm9pZCBFbmdpbmVCYXNlX2Rpc2NhcmQoRW5naW5lQmFzZSAqc2VsZiwgaW50IG4pCnsKICAgIHJldHVybiBzZWxmLT5kaXNjYXJkKHNlbGYsIG4pOwp9Cgp2b2lkIEVuZ2luZUJhc2VfZGVsZXRlKEVuZ2luZUJhc2UgKnNlbGYpCnsKICAgIHJldHVybiBzZWxmLT5kZWxldGUoc2VsZik7Cn0KCi8vIFN1YnRyYWN0IFdpdGggQ2FycnkuCgp0eXBlZGVmIHN0cnVjdCB7CiAgICBFbmdpbmVCYXNlIGJhc2U7CiAgICBpbnQgKngsIGMsIGksIG0sIHIsIHM7Cn0gU3VidHJhY3RXaXRoQ2FycnlFbmdpbmU7CgpzdGF0aWMgaW5saW5lIGludCBuZXh0KGludCAqeCwgaW50IG1hc2ssIGludCBpX3MsIGludCBpX3IsIGludCAqY2FycnkpCnsKICAgIGludCB5ID0geFtpX3NdIC0geFtpX3JdIC0gKmNhcnJ5OwogICAgKmNhcnJ5ID0gLSh5ID4+IChCSVRTKGludCkgLSAxKSk7CiAgICByZXR1cm4geSAmIG1hc2s7Cn0KCmludCBTdWJ0cmFjdFdpdGhDYXJyeUVuZ2luZV9uZXh0KEVuZ2luZUJhc2UgKmJhc2UpCnsKICAgIFN1YnRyYWN0V2l0aENhcnJ5RW5naW5lICpzZWxmID0gKFN1YnRyYWN0V2l0aENhcnJ5RW5naW5lICopIGJhc2U7CgogICAgc2VsZi0+aSArPSAxOwogICAgaW50ICp4ID0gc2VsZi0+eCwgaSA9IHNlbGYtPmksIHIgPSBzZWxmLT5yOwoKICAgIGlmIChVTkxJS0VMWShpID49IHIpKQogICAgewogICAgICAgIGludCBjID0gc2VsZi0+YywgbSA9IHNlbGYtPm0sIHMgPSBzZWxmLT5zLCB0ID0gciAtIHM7CgogICAgICAgIGZvciAoaSA9IDA7IGkgPCBzOyBpKyspCiAgICAgICAgICAgIHhbaV0gPSBuZXh0KHgsIG0sIGkgKyB0LCBpLCAmYyk7CiAgICAgICAgZm9yIChpID0gczsgaSA8IHI7IGkrKykKICAgICAgICAgICAgeFtpXSA9IG5leHQoeCwgbSwgaSAtIHMsIGksICZjKTsKICAgICAgICBzZWxmLT5jID0gYzsKICAgICAgICBzZWxmLT5pID0gaSA9IDA7CiAgICB9CiAgICByZXR1cm4geFtpXTsKfQoKdm9pZCBTdWJ0cmFjdFdpdGhDYXJyeUVuZ2luZV9kaXNjYXJkKEVuZ2luZUJhc2UgKmJhc2UsIGludCBuKQp7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IG47IGkrKykKICAgICAgICBTdWJ0cmFjdFdpdGhDYXJyeUVuZ2luZV9uZXh0KGJhc2UpOwp9Cgp2b2lkIFN1YnRyYWN0V2l0aENhcnJ5RW5naW5lX2RlbGV0ZShFbmdpbmVCYXNlICpiYXNlKQp7CiAgICBTdWJ0cmFjdFdpdGhDYXJyeUVuZ2luZSAqc2VsZiA9IChTdWJ0cmFjdFdpdGhDYXJyeUVuZ2luZSAqKSBiYXNlOwoKICAgIGlmIChzZWxmICE9IDApCiAgICB7CiAgICAgICAgZnJlZShzZWxmLT54KTsKICAgICAgICBmcmVlKHNlbGYpOwogICAgfQp9CgpFbmdpbmVCYXNlICpTdWJ0cmFjdFdpdGhDYXJyeUVuZ2luZV9uZXcoaW50IHcsIGludCBzLCBpbnQgcikKewogICAgYXNzZXJ0KDAgPCB3ICYmIHcgPCBCSVRTKGludCkpOwogICAgYXNzZXJ0KDAgPCBzICYmIHMgPCByKTsKCiAgICBTdWJ0cmFjdFdpdGhDYXJyeUVuZ2luZSAqc2VsZiA9IE5FV19UKFN1YnRyYWN0V2l0aENhcnJ5RW5naW5lLCAxKTsKCiAgICBzZWxmLT5iYXNlLm5leHQgPSBTdWJ0cmFjdFdpdGhDYXJyeUVuZ2luZV9uZXh0OwogICAgc2VsZi0+YmFzZS5kaXNjYXJkID0gU3VidHJhY3RXaXRoQ2FycnlFbmdpbmVfZGlzY2FyZDsKICAgIHNlbGYtPmJhc2UuZGVsZXRlID0gU3VidHJhY3RXaXRoQ2FycnlFbmdpbmVfZGVsZXRlOwogICAgc2VsZi0+eCA9IE5FV19UKGludCwgcik7CiAgICBzZWxmLT5jID0gMDsKICAgIHNlbGYtPmkgPSByLTE7CiAgICBzZWxmLT5tID0gLTFVID4+IChCSVRTKGludCkgLSB3KTsKICAgIHNlbGYtPnIgPSByOwogICAgc2VsZi0+cyA9IHM7CgogICAgZm9yIChpbnQgaSA9IDA7IGkgPCByOyBpKyspCiAgICAgICAgc2VsZi0+eFtpXSA9IHJhbmQoKSAmIHNlbGYtPm07CiAgICByZXR1cm4gKEVuZ2luZUJhc2UgKikgc2VsZjsKfQoKLy8gRGlzY2FyZCBCbG9jay4KCnR5cGVkZWYgc3RydWN0IHsKICAgIEVuZ2luZUJhc2UgYmFzZSwgKm93bmVkOwogICAgaW50IHAsIHIsIGk7Cn0gRGlzY2FyZEJsb2NrRW5naW5lOwoKaW50IERpc2NhcmRCbG9ja0VuZ2luZV9uZXh0KEVuZ2luZUJhc2UgKmJhc2UpCnsKICAgIERpc2NhcmRCbG9ja0VuZ2luZSAqc2VsZiA9IChEaXNjYXJkQmxvY2tFbmdpbmUgKikgYmFzZTsKCiAgICBpZiAoc2VsZi0+aSA9PSAwKQogICAgewogICAgICAgIEVuZ2luZUJhc2VfZGlzY2FyZChzZWxmLT5vd25lZCwgc2VsZi0+cCAtIHNlbGYtPnIpOwogICAgICAgIHNlbGYtPmkgPSBzZWxmLT5yOwogICAgfQogICAgc2VsZi0+aSAtPSAxOwogICAgcmV0dXJuIEVuZ2luZUJhc2VfbmV4dChzZWxmLT5vd25lZCk7Cn0KCnZvaWQgRGlzY2FyZEJsb2NrRW5naW5lX2Rpc2NhcmQoRW5naW5lQmFzZSAqYmFzZSwgaW50IG4pCnsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKQogICAgICAgIERpc2NhcmRCbG9ja0VuZ2luZV9uZXh0KGJhc2UpOwp9Cgp2b2lkIERpc2NhcmRCbG9ja0VuZ2luZV9kZWxldGUoRW5naW5lQmFzZSAqYmFzZSkKewogICAgRGlzY2FyZEJsb2NrRW5naW5lICpzZWxmID0gKERpc2NhcmRCbG9ja0VuZ2luZSAqKSBiYXNlOwoKICAgIGlmIChzZWxmICE9IDApCiAgICB7CiAgICAgICAgRW5naW5lQmFzZV9kZWxldGUoc2VsZi0+b3duZWQpOwogICAgICAgIGZyZWUoc2VsZik7CiAgICB9Cn0KCkVuZ2luZUJhc2UgKkRpc2NhcmRCbG9ja0VuZ2luZV9uZXcoRW5naW5lQmFzZSAqKnVuaXF1ZSwgaW50IHAsIGludCByKQp7CiAgICBhc3NlcnQoMCA8IHIgJiYgciA8PSBwKTsKCiAgICBEaXNjYXJkQmxvY2tFbmdpbmUgKnNlbGYgPSBORVdfVChEaXNjYXJkQmxvY2tFbmdpbmUsIDEpOwoKICAgIHNlbGYtPmJhc2UubmV4dCA9IERpc2NhcmRCbG9ja0VuZ2luZV9uZXh0OwogICAgc2VsZi0+YmFzZS5kaXNjYXJkID0gRGlzY2FyZEJsb2NrRW5naW5lX2Rpc2NhcmQ7CiAgICBzZWxmLT5iYXNlLmRlbGV0ZSA9IERpc2NhcmRCbG9ja0VuZ2luZV9kZWxldGU7CiAgICBzZWxmLT5vd25lZCA9ICp1bmlxdWU7ICp1bmlxdWUgPSAwOyAvLyBUcmFuc2ZlciBvd25lcnNoaXAuCiAgICBzZWxmLT5wID0gcDsKICAgIHNlbGYtPnIgPSByOwogICAgc2VsZi0+aSA9IHI7CiAgICByZXR1cm4gKEVuZ2luZUJhc2UgKikgc2VsZjsKfQoKLy8gUmFubHV4MjQuCgpFbmdpbmVCYXNlICpSYW5sdXgyNF9uZXcodm9pZCkKewogICAgRW5naW5lQmFzZSAqcmFubHV4MjRfYmFzZSA9IFN1YnRyYWN0V2l0aENhcnJ5RW5naW5lX25ldygyNCwgMTAsIDI0KTsKICAgIHJldHVybiBEaXNjYXJkQmxvY2tFbmdpbmVfbmV3KCZyYW5sdXgyNF9iYXNlLCAyMjMsIDIzKTsKfQoKLy8gTWFpbi4KCmRvdWJsZSBjbG9ja19ub3codm9pZCkKewogICAgc3RydWN0IHRpbWVzcGVjIG5vdzsKICAgIGNsb2NrX2dldHRpbWUoQ0xPQ0tfUFJPQ0VTU19DUFVUSU1FX0lELCAmbm93KTsKICAgIHJldHVybiBub3cudHZfc2VjICsgbm93LnR2X25zZWMgLyAxLjBFKzA5Owp9CgppbnQgbWFpbih2b2lkKQp7CiAgICBzcmFuZCh0aW1lKDApKTsKICAgIEVuZ2luZUJhc2UgKnJhbmx1eDI0ID0gUmFubHV4MjRfbmV3KCk7CgogICAgZm9yIChpbnQgaSA9IDA7IGkgPCAyNDsgaSsrKQogICAgewogICAgICAgIGludCByID0gRW5naW5lQmFzZV9uZXh0KHJhbmx1eDI0KTsKICAgICAgICBwcmludGYoIiVkXG4iLCByKTsKICAgIH0KCiAgICBpbnQgbiA9IDEwMDAwMDAwOwogICAgZG91YmxlIHQgPSBjbG9ja19ub3coKTsKICAgIEVuZ2luZUJhc2VfZGlzY2FyZChyYW5sdXgyNCwgbik7CiAgICBwcmludGYoIkVsYXBzZWQ6ICUuOWZzXG4iLCBjbG9ja19ub3coKSAtIHQpOwoKICAgIEVuZ2luZUJhc2VfZGVsZXRlKHJhbmx1eDI0KTsKICAgIHJldHVybiAwOwp9