Skip to content

Commit

Permalink
Merge branch 'master' into cornu/remove_regexp
Browse files Browse the repository at this point in the history
  • Loading branch information
@alkino
alkino authored Jan 26, 2024
2 parents f98ebed + f64b609 commit 3eef6d3
Show file tree
Hide file tree
Showing 3 changed files with 92 additions and 82 deletions.
8 changes: 6 additions & 2 deletions src/ivoc/ivocrand.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -337,8 +337,12 @@ static double r_nrnran123(void* r) {
id2 = (uint32_t) (chkarg(2, 0., dmaxuint));
if (ifarg(3))
id3 = (uint32_t) (chkarg(3, 0., dmaxuint));
NrnRandom123* r123 = new NrnRandom123(id1, id2, id3);
x->rand->generator(r123);
try {
NrnRandom123* r123 = new NrnRandom123(id1, id2, id3);
x->rand->generator(r123);
} catch (const std::bad_alloc& e) {
hoc_execerror("Bad allocation for 'NrnRandom123'", e.what());
}
delete x->gen;
x->gen = x->rand->generator();
x->type_ = 4;
Expand Down
117 changes: 57 additions & 60 deletions src/oc/nrnran123.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,97 +1,97 @@
#include <../../nrnconf.h>
#include <inttypes.h>
#include <nrnran123.h>
#include <hocdec.h>
#include <stdlib.h>
#include <math.h>
#include <cstdint>
#include "nrnran123.h"
#include <cstdlib>
#include <cmath>
#include <Random123/philox.h>

static const double SHIFT32 = 1.0 / 4294967297.0; /* 1/(2^32 + 1) */
using RNG = r123::Philox4x32;

static philox4x32_key_t k = {{0}};
static RNG::key_type k = {{0}};

struct nrnran123_State {
philox4x32_ctr_t c;
philox4x32_ctr_t r;
RNG::ctr_type c;
RNG::ctr_type r;
char which_;
};

void nrnran123_set_globalindex(uint32_t gix) {
k.v[0] = gix;
void nrnran123_set_globalindex(std::uint32_t gix) {
k[0] = gix;
}

/* if one sets the global, one should reset all the stream sequences. */
uint32_t nrnran123_get_globalindex() {
return k.v[0];
std::uint32_t nrnran123_get_globalindex() {
return k[0];
}

nrnran123_State* nrnran123_newstream(uint32_t id1, uint32_t id2) {
nrnran123_State* nrnran123_newstream(std::uint32_t id1, std::uint32_t id2) {
return nrnran123_newstream3(id1, id2, 0);
}
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_State* nrnran123_newstream3(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3) {
auto* s = new nrnran123_State;
s->c[1] = id3;
s->c[2] = id1;
s->c[3] = id2;
nrnran123_setseq(s, 0, 0);
return s;
}

void nrnran123_deletestream(nrnran123_State* s) {
free(s);
delete s;
}

void nrnran123_getseq(nrnran123_State* s, uint32_t* seq, char* which) {
*seq = s->c.v[0];
void nrnran123_getseq(nrnran123_State* s, std::uint32_t* seq, char* which) {
*seq = s->c[0];
*which = s->which_;
}

void nrnran123_setseq(nrnran123_State* s, uint32_t seq, char which) {
void nrnran123_setseq(nrnran123_State* s, std::uint32_t seq, char which) {
if (which > 3 || which < 0) {
s->which_ = 0;
} else {
s->which_ = which;
}
s->c.v[0] = seq;
s->c[0] = seq;
s->r = philox4x32(s->c, k);
}

void nrnran123_getids(nrnran123_State* s, uint32_t* id1, uint32_t* id2) {
*id1 = s->c.v[2];
*id2 = s->c.v[3];
void nrnran123_getids(nrnran123_State* s, std::uint32_t* id1, std::uint32_t* id2) {
*id1 = s->c[2];
*id2 = s->c[3];
}

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];
void nrnran123_getids3(nrnran123_State* s,
std::uint32_t* id1,
std::uint32_t* id2,
std::uint32_t* id3) {
*id3 = s->c[1];
*id1 = s->c[2];
*id2 = s->c[3];
}

uint32_t nrnran123_ipick(nrnran123_State* s) {
uint32_t rval;
std::uint32_t nrnran123_ipick(nrnran123_State* s) {
char which = s->which_;
assert(which < 4);
rval = s->r.v[which++];
std::uint32_t rval = s->r[which++];
if (which > 3) {
which = 0;
s->c.v[0]++;
s->c.incr();
s->r = philox4x32(s->c, k);
}
s->which_ = which;
return rval;
}

double nrnran123_dblpick(nrnran123_State* s) {
return nrnran123_uint2dbl(nrnran123_ipick(s));
static const double SHIFT32 = 1.0 / 4294967297.0; /* 1/(2^32 + 1) */
auto u = nrnran123_ipick(s);
return ((double) u + 1.0) * SHIFT32;
}

double nrnran123_negexp(nrnran123_State* s) {
/* min 2.3283064e-10 to max 22.18071 */
return -log(nrnran123_dblpick(s));
return -std::log(nrnran123_dblpick(s));
}

/* At cost of a cached value we could compute two at a time. */
/* 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 */
double nrnran123_normal(nrnran123_State* s) {
double w, x, y;
Expand All @@ -104,31 +104,28 @@ double nrnran123_normal(nrnran123_State* s) {
w = (u1 * u1) + (u2 * u2);
} while (w > 1);

y = sqrt((-2. * log(w)) / w);
y = std::sqrt((-2. * log(w)) / w);
x = u1 * y;
return x;
}

nrnran123_array4x32 nrnran123_iran(uint32_t seq, uint32_t id1, uint32_t id2) {
nrnran123_array4x32 nrnran123_iran(std::uint32_t seq, std::uint32_t id1, std::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 nrnran123_iran3(std::uint32_t seq,
std::uint32_t id1,
std::uint32_t id2,
std::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];
RNG::ctr_type c;
c[0] = seq;
c[1] = id3;
c[2] = id1;
c[3] = id2;
RNG::ctr_type r = philox4x32(c, k);
a.v[0] = r[0];
a.v[1] = r[1];
a.v[2] = r[2];
a.v[3] = r[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;
}
49 changes: 29 additions & 20 deletions src/oc/nrnran123.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,40 +23,49 @@ of the full distribution available from
http://www.deshawresearch.com/resources_random123.html
*/

#include <inttypes.h>
#include <cstdint>


typedef struct nrnran123_State nrnran123_State;
struct nrnran123_State;

typedef struct nrnran123_array4x32 {
uint32_t v[4];
} nrnran123_array4x32;
struct nrnran123_array4x32 {
std::uint32_t v[4];
};

/* global index. eg. run number */
/* all generator instances share this global index */
extern void nrnran123_set_globalindex(uint32_t gix);
extern uint32_t nrnran123_get_globalindex();
extern void nrnran123_set_globalindex(std::uint32_t gix);
extern std::uint32_t nrnran123_get_globalindex();

/* minimal data stream */
extern nrnran123_State* nrnran123_newstream(uint32_t id1, uint32_t id2);
extern nrnran123_State* nrnran123_newstream3(uint32_t id1, uint32_t id2, uint32_t id3);
extern nrnran123_State* nrnran123_newstream(std::uint32_t id1, std::uint32_t id2);
extern nrnran123_State* nrnran123_newstream3(std::uint32_t id1,
std::uint32_t id2,
std::uint32_t id3);
extern void nrnran123_deletestream(nrnran123_State*);
extern void nrnran123_getseq(nrnran123_State*, uint32_t* seq, char* which);
extern void nrnran123_setseq(nrnran123_State*, uint32_t seq, char which);
extern void nrnran123_getids(nrnran123_State*, uint32_t* id1, uint32_t* id2);
extern void nrnran123_getids3(nrnran123_State*, uint32_t* id1, uint32_t* id2, uint32_t* id3);
extern void nrnran123_getseq(nrnran123_State*, std::uint32_t* seq, char* which);
extern void nrnran123_setseq(nrnran123_State*, std::uint32_t seq, char which);
extern void nrnran123_getids(nrnran123_State*, std::uint32_t* id1, std::uint32_t* id2);
extern void nrnran123_getids3(nrnran123_State*,
std::uint32_t* id1,
std::uint32_t* id2,
std::uint32_t* id3);

extern double nrnran123_negexp(nrnran123_State*); /* mean 1.0 */
extern uint32_t nrnran123_ipick(nrnran123_State*); /* uniform 0 to 2^32-1 */
extern double nrnran123_dblpick(nrnran123_State*); /* uniform open interval (0,1)*/
/* nrnran123_dblpick minimum value is 2.3283064e-10 and max value is 1-min */
// Get a random uint32_t in [0, 2^32-1]
extern std::uint32_t nrnran123_ipick(nrnran123_State*);
// Get a random double on [0, 1]
// nrnran123_dblpick minimum value is 2.3283064e-10 and max value is 1-min
extern double nrnran123_dblpick(nrnran123_State*);

/* nrnran123_negexp min value is 2.3283064e-10, max is 22.18071 */
extern double nrnran123_negexp(nrnran123_State*); /* mean 1.0 */
extern double nrnran123_normal(nrnran123_State*); /* mean 0.0, std 1.0 */

/* more fundamental (stateless) (though the global index is still used) */
extern nrnran123_array4x32 nrnran123_iran(uint32_t seq, uint32_t id1, uint32_t id2);
extern nrnran123_array4x32 nrnran123_iran3(uint32_t seq, uint32_t id1, uint32_t id2, uint32_t id3);
extern double nrnran123_uint2dbl(uint32_t);
extern nrnran123_array4x32 nrnran123_iran(std::uint32_t seq, std::uint32_t id1, std::uint32_t id2);
extern nrnran123_array4x32 nrnran123_iran3(std::uint32_t seq,
std::uint32_t id1,
std::uint32_t id2,
std::uint32_t id3);

#endif

0 comments on commit 3eef6d3

Please sign in to comment.