Remove usage of MCellRan4

cmake/NeuronFileLists.cmake

@@ -9,7 +9,6 @@ set(STRUCTURED_HEADER_FILES_TO_INSTALL
     neuron/container/generic_data_handle.hpp neuron/container/non_owning_soa_identifier.hpp
     neuron/model_data_fwd.hpp)
 set(HEADER_FILES_TO_INSTALL
-    gnu/mcran4.h
     gnu/nrnran123.h
     nrniv/backtrace_utils.h
     nrniv/bbsavestate.h

src/gnu/CMakeLists.txt

@@ -6,7 +6,6 @@ add_library(
   Geom.cpp
   HypGeom.cpp
   LogNorm.cpp
-  mcran4.cpp
   NegExp.cpp
   Normal.cpp
   NrnRandom123RNG.cpp

src/nrniv/kschan.cpp

@@ -279,9 +279,12 @@ static double ks_vres(void* v) {
 
 static double ks_rseed(void* v) {
     if (ifarg(1)) {
-        KSSingle::idum_ = (unsigned int) chkarg(1, 0, 1e9);
+        nrnran123_setseq(KSSingle::rand_state,  static_cast<std::uint32_t>(chkarg(1, 0, 1e9)), 0);
     }
-    return (double) KSSingle::idum_;
+    std::uint32_t seq;
+    char which;
+    nrnran123_getseq(KSSingle::rand_state, &seq, &which);
+    return (double) seq;
 }
 
 static double ks_usetable(void* v) {
@@ -804,7 +807,7 @@ void KSChan_reg() {
     ksgate_sym = hoc_lookup("KSGate");
     kstrans_sym = hoc_lookup("KSTrans");
     KSSingle::vres_ = 0.1;
-    KSSingle::idum_ = 0;
+    KSSingle::rand_state = nrnran123_newstream(0, 0);
 }
 
 // param is gmax, g, i --- if change then change numbers below

src/nrniv/kssingle.cpp

@@ -130,7 +130,7 @@
 extern NetCvode* net_cvode_instance;
 
 double KSSingle::vres_;
-unsigned int KSSingle::idum_;
+nrnran123_State* KSSingle::rand_state;
 
 KSSingle::KSSingle(KSChan* c) {
     // implemenation assumes one ks gate complex with power 1

src/nrniv/kssingle.h

@@ -5,8 +5,7 @@
 #include <mymath.h>
 #include <kschan.h>
 #include <netcon.h>
-
-#include <mcran4.h>
+#include <nrnran123.h>
 
 class KSSingleTrans;
 class KSSingleState;
@@ -71,10 +70,10 @@ class KSSingle {
         return MyMath::eq(x, y, vres_);
     }
     double exprand() {
-        return -log(mcell_ran4a(&idum_));
+        return nrnran123_negexp(rand_state);
     }
     double unifrand(double range) {
-        return mcell_ran4a(&idum_) * range;
+        return nrnran123_dblpick(rand_state) * range;
     }
     int rvalrand(int);
 
@@ -84,7 +83,7 @@ class KSSingle {
     double* rval_;
     bool uses_ligands_;
     static double vres_;
-    static unsigned int idum_;
+    static nrnran123_State* rand_state;
 
     static unsigned long singleevent_deliver_;
     static unsigned long singleevent_move_;

src/oc/hoc_init.cpp

@@ -255,7 +255,6 @@ void hoc_init(void) /* install constants and built-ins table */
         }
     }
 
-    set_use_mcran4(false);
     nrn_xopen_broadcast_ = 255;
     extern void hoc_init_space(void);
     hoc_init_space();

src/oc/oc_mcran4.cpp

@@ -1,41 +1,33 @@
 #include "hocdec.h"
-#include "mcran4.h"
+#include "nrnran123.h"
 
-int use_mcell_ran4_;
-
-void set_use_mcran4(bool value) {
-    use_mcell_ran4_ = value ? 1 : 0;
-}
-
-bool use_mcran4() {
-    return use_mcell_ran4_ != 0;
-}
+static const nrnran123_State* state = nrnran123_newstream(0, 0);
 
 void hoc_mcran4() {
     uint32_t idx;
     double* xidx;
     double x;
     xidx = hoc_pgetarg(1);
     idx = (uint32_t) (*xidx);
-    x = mcell_ran4a(&idx);
-    *xidx = idx;
+    nrnran123_setseq(state, idx, 0);
+    x = nrnran123_dblpick(state);
+    *xidx = idx + 1;
     hoc_ret();
     hoc_pushx(x);
 }
 void hoc_mcran4init() {
-    double prev = mcell_lowindex();
+    std::uint32_t seq;
+    char which;
+    nrnran123_getseq(state, &seq, &which);
+    double prev = static_cast<double>(seq);
     if (ifarg(1)) {
         uint32_t idx = (uint32_t) chkarg(1, 0., 4294967295.);
-        mcell_ran4_init(idx);
+        nrnran123_setseq(state, idx, 0);
     }
     hoc_ret();
     hoc_pushx(prev);
 }
 void hoc_usemcran4() {
-    double prev = (double) use_mcell_ran4_;
-    if (ifarg(1)) {
-        use_mcell_ran4_ = (int) chkarg(1, 0., 1.);
-    }
     hoc_ret();
-    hoc_pushx(prev);
+    hoc_pushx(0);
 }

src/oc/oc_mcran4.hpp

@@ -1,5 +1,3 @@
-void set_use_mcran4(bool value);
-bool use_mcran4();
 void hoc_mcran4();
 void hoc_mcran4init();
 void hoc_usemcran4();

src/oc/scoprand.cpp

@@ -23,8 +23,7 @@ static char RCSid[] = "random.cpp,v 1.4 1999/01/04 12:46:49 hines Exp";
 #endif
 
 #include <math.h>
-#include "oc_mcran4.hpp"
-#include "mcran4.h"
+#include "nrnran123.h"
 #include "scoplib.h"
 static uint32_t value = 1;
 
@@ -59,17 +58,8 @@ static uint32_t value = 1;
  *--------------------------------------------------------------------------- */
 
 double scop_random(void) {
-    if (use_mcran4()) {
-        /*perhaps 4 times slower but much higher quality*/
-        return mcell_ran4a(&value);
-    } else {
-        uint32_t a = 2147437301, c = 453816981,
-                 /* m = 2^32 - 1, the largest long int value that can be represented */
-            /*m = 0xFFFFFFFF;*/ /* limited to 32 bit integers*/
-            m = ~0;
-        value = a * value + c;
-        return (fabs((double) value / (double) m));
-    }
+    auto *rng = nrnran123_newstream(value, 0);
+    return nrnran123_dblpick(rng);
 }
 
 /*-----------------------------------------------------------------------------