diff --git a/src/ivoc/ivocrand.cpp b/src/ivoc/ivocrand.cpp index d3501913db..fa45403dcc 100644 --- a/src/ivoc/ivocrand.cpp +++ b/src/ivoc/ivocrand.cpp @@ -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; diff --git a/src/oc/nrnran123.cpp b/src/oc/nrnran123.cpp index d9db49351c..eeb2768115 100644 --- a/src/oc/nrnran123.cpp +++ b/src/oc/nrnran123.cpp @@ -1,81 +1,79 @@ -#include <../../nrnconf.h> -#include -#include -#include -#include -#include +#include +#include "nrnran123.h" +#include +#include #include -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; @@ -83,15 +81,17 @@ uint32_t nrnran123_ipick(nrnran123_State* s) { } 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; @@ -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; -} diff --git a/src/oc/nrnran123.h b/src/oc/nrnran123.h index 0a0acb6a71..a4e4844bf6 100644 --- a/src/oc/nrnran123.h +++ b/src/oc/nrnran123.h @@ -23,40 +23,49 @@ of the full distribution available from http://www.deshawresearch.com/resources_random123.html */ -#include +#include -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