nrnran123.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#include <inttypes.h>
#include <nrnran123.h>
#include <hocdec.h>
#include <stdlib.h>
#include <math.h>
#include <Random123/philox.h>
 
static const double SHIFT32 = 1.0 / 4294967297.0; /* 1/(2^32 + 1) */
 
static philox4x32_key_t k = {{0}};
 
struct nrnran123_State {
    philox4x32_ctr_t c;
    philox4x32_ctr_t r;
    char which_;
};
 
void nrnran123_set_globalindex(uint32_t gix) {
    k.v[0] = gix;
}
 
/* if one sets the global, one should reset all the stream sequences. */
uint32_t nrnran123_get_globalindex() {
    return k.v[0];
}
 
extern "C" nrnran123_State* nrnran123_newstream(uint32_t id1, uint32_t id2) {
    return nrnran123_newstream3(id1, id2, 0);
}
extern "C" nrnran123_State* nrnran123_newstream3(uint32_t id1, uint32_t id2, uint32_t id3) {
    nrnran123_State* s;
    s = (nrnran123_State*) ecalloc(sizeof(nrnran123_State), 1);
    s->c.v[1] = id3;
    s->c.v[2] = id1;
    s->c.v[3] = id2;
    nrnran123_setseq(s, 0, 0);
    return s;
}
 
extern "C" void nrnran123_deletestream(nrnran123_State* s) {
    free(s);
}
 
extern "C" void nrnran123_getseq(nrnran123_State* s, uint32_t* seq, char* which) {
    *seq = s->c.v[0];
    *which = s->which_;
}
 
extern "C" void nrnran123_setseq(nrnran123_State* s, uint32_t seq, char which) {
    if (which > 3 || which < 0) {
        s->which_ = 0;
    } else {
        s->which_ = which;
    }
    s->c.v[0] = seq;
    s->r = philox4x32(s->c, k);
}
 
extern "C" void nrnran123_getids(nrnran123_State* s, uint32_t* id1, uint32_t* id2) {
    *id1 = s->c.v[2];
    *id2 = s->c.v[3];
}
 
extern "C" void nrnran123_getids3(nrnran123_State* s, uint32_t* id1, uint32_t* id2, uint32_t* id3) {
    *id3 = s->c.v[1];
    *id1 = s->c.v[2];
    *id2 = s->c.v[3];
}
 
extern "C" uint32_t nrnran123_ipick(nrnran123_State* s) {
    uint32_t rval;
    char which = s->which_;
    assert(which < 4);
    rval = s->r.v[which++];
    if (which > 3) {
        which = 0;
        s->c.v[0]++;
        s->r = philox4x32(s->c, k);
    }
    s->which_ = which;
    return rval;
}
 
extern "C" double nrnran123_dblpick(nrnran123_State* s) {
    return nrnran123_uint2dbl(nrnran123_ipick(s));
}
 
extern "C" double nrnran123_negexp(nrnran123_State* s) {
    /* min 2.3283064e-10 to max 22.18071 */
    return -log(nrnran123_dblpick(s));
}
 
/* At cost of a cached  value we could compute two at a time. */
/* But that would make it difficult to transfer to coreneuron for t > 0 */
extern "C" double nrnran123_normal(nrnran123_State* s) {
    double w, x, y;
    double u1, u2;
    do {
        u1 = nrnran123_dblpick(s);
        u2 = nrnran123_dblpick(s);
        u1 = 2. * u1 - 1.;
        u2 = 2. * u2 - 1.;
        w = (u1 * u1) + (u2 * u2);
    } while (w > 1);
 
    y = sqrt((-2. * log(w)) / w);
    x = u1 * y;
    return x;
}
 
nrnran123_array4x32 nrnran123_iran(uint32_t seq, uint32_t id1, uint32_t id2) {
    return nrnran123_iran3(seq, id1, id2, 0);
}
nrnran123_array4x32 nrnran123_iran3(uint32_t seq, uint32_t id1, uint32_t id2, uint32_t id3) {
    nrnran123_array4x32 a;
    philox4x32_ctr_t c;
    c.v[0] = seq;
    c.v[1] = id3;
    c.v[2] = id1;
    c.v[3] = id2;
    philox4x32_ctr_t r = philox4x32(c, k);
    a.v[0] = r.v[0];
    a.v[1] = r.v[1];
    a.v[2] = r.v[2];
    a.v[3] = r.v[3];
    return a;
}
 
double nrnran123_uint2dbl(uint32_t u) {
    /* 0 to 2^32-1 transforms to double value in open (0,1) interval */
    /* min 2.3283064e-10 to max (1 - 2.3283064e-10) */
    return ((double) u + 1.0) * SHIFT32;
}