isaac64.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#if HAVE_SYS_TYPES_H
#include <sys/types.h>
#endif
 
/*
------------------------------------------------------------------------------
isaac64.cpp: A Fast cryptographic random number generator
           for 32-bit and 64-bit machines.
By Bob Jenkins, 1996.  Public Domain.
 
Modified for modularity by Tom Bartol and Rex Kerr
------------------------------------------------------------------------------
*/
#include "isaac64.h"
 
 
#define ind(mm, x) (*(ub8*) ((ub1*) (mm) + ((x) & ((RANDSIZ - 1) << 3))))
 
#define rngstep(mix, a, b, mm, m, m2, r, x)      \
    {                                            \
        x = *m;                                  \
        a = (mix) + *(m2++);                     \
        *(m++) = y = ind(mm, x) + a + b;         \
        *(r++) = b = ind(mm, y >> RANDSIZL) + x; \
    }
 
#define mix(a, b, c, d, e, f, g, h) \
    {                               \
        a -= e;                     \
        f ^= h >> 9;                \
        h += a;                     \
        b -= f;                     \
        g ^= a << 9;                \
        a += b;                     \
        c -= g;                     \
        h ^= b >> 23;               \
        b += c;                     \
        d -= h;                     \
        a ^= c << 15;               \
        c += d;                     \
        e -= a;                     \
        b ^= d >> 14;               \
        d += e;                     \
        f -= b;                     \
        c ^= e << 20;               \
        e += f;                     \
        g -= c;                     \
        d ^= f >> 17;               \
        f += g;                     \
        h -= d;                     \
        e ^= g << 14;               \
        g += h;                     \
    }
 
 
void isaac64_generate(struct isaac64_state* rng) {
    ub8 a, b, x, y, *m, *m2, *r, *mend;
 
    m = rng->mm;
    r = rng->randrsl;
    a = rng->aa;
    b = rng->bb + (++rng->cc);
    for (m = rng->mm, mend = m2 = m + (RANDSIZ / 2); m < mend;) {
        rngstep(~(a ^ (a << 21)), a, b, rng->mm, m, m2, r, x);
        rngstep(a ^ (a >> 5), a, b, rng->mm, m, m2, r, x);
        rngstep(a ^ (a << 12), a, b, rng->mm, m, m2, r, x);
        rngstep(a ^ (a >> 33), a, b, rng->mm, m, m2, r, x);
    }
    for (m2 = rng->mm; m2 < mend;) {
        rngstep(~(a ^ (a << 21)), a, b, rng->mm, m, m2, r, x);
        rngstep(a ^ (a >> 5), a, b, rng->mm, m, m2, r, x);
        rngstep(a ^ (a << 12), a, b, rng->mm, m, m2, r, x);
        rngstep(a ^ (a >> 33), a, b, rng->mm, m, m2, r, x);
    }
    rng->bb = b;
    rng->aa = a;
}
 
 
void isaac64_init(struct isaac64_state* rng, ub4 seed) {
    ub8 *r, *m;
    ub8 a, b, c, d, e, f, g, h;
    ub4 i;
 
    rng->aa = (ub8) 0;
    rng->bb = (ub8) 0;
    rng->cc = (ub8) 0;
 
    a = b = c = d = e = f = g = h = 0x9e3779b97f4a7c13LL; /* the golden ratio */
 
    r = rng->randrsl;
    m = rng->mm;
 
    for (i = 0; i < RANDSIZ; ++i)
        r[i] = (ub8) 0;
 
    r[0] = seed;
 
    for (i = 0; i < 4; ++i) /* scramble it */
    {
        mix(a, b, c, d, e, f, g, h);
    }
 
    for (i = 0; i < RANDSIZ; i += 8) /* fill in m[] with messy stuff */
    {
        /* use all the information in the seed */
        a += r[i];
        b += r[i + 1];
        c += r[i + 2];
        d += r[i + 3];
        e += r[i + 4];
        f += r[i + 5];
        g += r[i + 6];
        h += r[i + 7];
        mix(a, b, c, d, e, f, g, h);
        m[i] = a;
        m[i + 1] = b;
        m[i + 2] = c;
        m[i + 3] = d;
        m[i + 4] = e;
        m[i + 5] = f;
        m[i + 6] = g;
        m[i + 7] = h;
    }
 
    /* do a second pass to make all of the seed affect all of m[] */
    for (i = 0; i < RANDSIZ; i += 8) {
        a += m[i];
        b += m[i + 1];
        c += m[i + 2];
        d += m[i + 3];
        e += m[i + 4];
        f += m[i + 5];
        g += m[i + 6];
        h += m[i + 7];
        mix(a, b, c, d, e, f, g, h);
        m[i] = a;
        m[i + 1] = b;
        m[i + 2] = c;
        m[i + 3] = d;
        m[i + 4] = e;
        m[i + 5] = f;
        m[i + 6] = g;
        m[i + 7] = h;
    }
 
    isaac64_generate(rng);  /* fill in the first set of results */
    rng->randcnt = RANDMAX; /* prepare to use the first set of results */
}