Move random number related code to src/gnu (#2689)

* This patch moves most of the code that is not HOC to the directory gnu.
   * `Isaac`, `MCellRan4`, `Random123`, `Rand`
* In the future this directory will be cleaned to remove gnu stuff and
   renamed something like 'src/random'.
* This way nrn will have an external random library, and a clean hoc_interface.

cmake/CompilerHelper.cmake

@@ -51,7 +51,7 @@ if(CMAKE_C_COMPILER_ID MATCHES "PGI" OR CMAKE_C_COMPILER_ID MATCHES "NVHPC")
     # Random123 does not play nicely with NVHPC 21.11+'s detection of ABM features, see:
     # https://github.com/BlueBrain/CoreNeuron/issues/724 and
     # https://github.com/DEShawResearch/random123/issues/6.
-    list(APPEND NRN_COMPILE_DEFS R123_USE_INTRIN_H=0)
+    list(APPEND NRN_R123_COMPILE_DEFS R123_USE_INTRIN_H=0)
   endif()
 else()
   set(NRN_HAVE_NVHPC_COMPILER OFF)

cmake/NeuronFileLists.cmake

@@ -9,6 +9,9 @@ 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/nrnisaac.h
+    gnu/nrnran123.h
     nrniv/backtrace_utils.h
     nrniv/bbsavestate.h
     nrnmpi/nrnmpidec.h
@@ -37,15 +40,12 @@ set(HEADER_FILES_TO_INSTALL
     oc/hocgetsym.h
     oc/hoclist.h
     oc/hocparse.h
-    oc/mcran4.h
     oc/mech_api.h
     oc/memory.hpp
     oc/nrnapi.h
     oc/nrnassrt.h
-    oc/nrnisaac.h
     oc/nrnmpi.h
     oc/nrnrandom.h
-    oc/nrnran123.h
     oc/oc_ansi.h
     oc/ocfunc.h
     oc/ocmisc.h
@@ -77,16 +77,10 @@ set(HEADER_FILES_TO_INSTALL
 # =============================================================================
 set(NRN_HEADERS_INCLUDE_LIST)
 
-# =============================================================================
-# Lists of random number related files
-# =============================================================================
-set(RAN_FILE_LIST isaac64.cpp mcran4.cpp nrnisaac.cpp nrnran123.cpp)
-
 # =============================================================================
 # Files in oc directory
 # =============================================================================
 set(OC_FILE_LIST
-    ${RAN_FILE_LIST}
     audit.cpp
     axis.cpp
     code.cpp
@@ -102,6 +96,7 @@ set(OC_FILE_LIST
     hoc_oop.cpp
     list.cpp
     math.cpp
+    oc_mcran4.cpp
     memory.cpp
     mswinprt.cpp
     nonlin.cpp

src/gnu/CMakeLists.txt

@@ -6,14 +6,23 @@ add_library(
   Erlang.cpp
   Geom.cpp
   HypGeom.cpp
+  isaac64.cpp
+  Isaac64RNG.cpp
   LogNorm.cpp
+  MCellRan4RNG.cpp
+  mcran4.cpp
   MLCG.cpp
   NegExp.cpp
   Normal.cpp
+  nrnisaac.cpp
+  nrnran123.cpp
   Poisson.cpp
+  Rand.cpp
   Random.cpp
   RndInt.cpp
   RNG.cpp
   Uniform.cpp
   Weibull.cpp)
 set_property(TARGET nrngnu PROPERTY POSITION_INDEPENDENT_CODE ON)
+target_include_directories(nrngnu SYSTEM PRIVATE "${PROJECT_SOURCE_DIR}/external/Random123/include")
+target_compile_definitions(nrngnu PRIVATE ${NRN_R123_COMPILE_DEFS})

src/gnu/Isaac64RNG.cpp

@@ -0,0 +1,20 @@
+#include "Isaac64RNG.hpp"
+
+uint32_t Isaac64::cnt_ = 0;
+
+Isaac64::Isaac64(std::uint32_t seed) {
+    if (cnt_ == 0) {
+        cnt_ = 0xffffffff;
+    }
+    --cnt_;
+    seed_ = seed;
+    if (seed_ == 0) {
+        seed_ = cnt_;
+    }
+    rng_ = nrnisaac_new();
+    reset();
+}
+
+Isaac64::~Isaac64() {
+    nrnisaac_delete(rng_);
+}

src/gnu/Isaac64RNG.hpp

@@ -0,0 +1,33 @@
+#pragma once
+
+#include <cstdint>
+
+#include "RNG.h"
+#include "nrnisaac.h"
+
+class Isaac64: public RNG {
+  public:
+    Isaac64(std::uint32_t seed = 0);
+    ~Isaac64();
+    std::uint32_t asLong() {
+        return nrnisaac_uint32_pick(rng_);
+    }
+    void reset() {
+        nrnisaac_init(rng_, seed_);
+    }
+    double asDouble() {
+        return nrnisaac_dbl_pick(rng_);
+    }
+    std::uint32_t seed() {
+        return seed_;
+    }
+    void seed(std::uint32_t s) {
+        seed_ = s;
+        reset();
+    }
+
+  private:
+    std::uint32_t seed_;
+    void* rng_;
+    static std::uint32_t cnt_;
+};

src/gnu/MCellRan4RNG.cpp

@@ -0,0 +1,15 @@
+#include "MCellRan4RNG.hpp"
+
+MCellRan4::MCellRan4(std::uint32_t ihigh, std::uint32_t ilow) {
+    ++cnt_;
+    ilow_ = ilow;
+    ihigh_ = ihigh;
+    if (ihigh_ == 0) {
+        ihigh_ = cnt_;
+        ihigh_ = (std::uint32_t) asLong();
+    }
+    orig_ = ihigh_;
+}
+MCellRan4::~MCellRan4() {}
+
+std::uint32_t MCellRan4::cnt_ = 0;

src/gnu/MCellRan4RNG.hpp

@@ -0,0 +1,34 @@
+#pragma once
+
+#include <cstdint>
+
+#include "RNG.h"
+#include "mcran4.h"
+
+// The decision that has to be made is whether each generator instance
+// should have its own seed or only one seed for all. We choose separate
+// seed for each but if arg not present or 0 then seed chosen by system.
+
+// the addition of ilow > 0 means that value is used for the lowindex
+// instead of the mcell_ran4_init global 32 bit lowindex.
+
+class MCellRan4: public RNG {
+  public:
+    MCellRan4(std::uint32_t ihigh = 0, std::uint32_t ilow = 0);
+    virtual ~MCellRan4();
+    virtual std::uint32_t asLong() {
+        return (std::uint32_t) (ilow_ == 0 ? mcell_iran4(&ihigh_) : nrnRan4int(&ihigh_, ilow_));
+    }
+    virtual void reset() {
+        ihigh_ = orig_;
+    }
+    virtual double asDouble() {
+        return (ilow_ == 0 ? mcell_ran4a(&ihigh_) : nrnRan4dbl(&ihigh_, ilow_));
+    }
+    std::uint32_t ihigh_;
+    std::uint32_t orig_;
+    std::uint32_t ilow_;
+
+  private:
+    static std::uint32_t cnt_;
+};

src/gnu/NrnRandom123RNG.hpp

@@ -0,0 +1,28 @@
+#pragma once
+
+#include <cstdint>
+
+#include "nrnran123.h"
+#include "RNG.h"
+
+class NrnRandom123: public RNG {
+  public:
+    NrnRandom123(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3 = 0);
+    ~NrnRandom123();
+    std::uint32_t asLong() {
+        return nrnran123_ipick(s_);
+    }
+    double asDouble() {
+        return nrnran123_dblpick(s_);
+    }
+    void reset() {
+        nrnran123_setseq(s_, 0, 0);
+    }
+    nrnran123_State* s_;
+};
+NrnRandom123::NrnRandom123(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3) {
+    s_ = nrnran123_newstream3(id1, id2, id3);
+}
+NrnRandom123::~NrnRandom123() {
+    nrnran123_deletestream(s_);
+}

src/gnu/Rand.cpp

@@ -0,0 +1,19 @@
+#include "Rand.hpp"
+
+#include "ACG.h"
+#include "Normal.h"
+
+Rand::Rand(unsigned long seed, int size, Object* obj) {
+    // printf("Rand\n");
+    gen = new ACG(seed, size);
+    rand = new Normal(0., 1., gen);
+    type_ = 0;
+    obj_ = obj;
+}
+
+Rand::~Rand() {
+    // printf("~Rand\n");
+    delete gen;
+    delete rand;
+}
+

src/ivoc/random1.h → src/gnu/Rand.hpp

@@ -1,20 +1,24 @@
-#ifndef random1_h
-#define random1_h
+#pragma once
 
 #include "RNG.h"
 #include "Random.h"
 
 struct Object;
 
+/* type_:
+ * 0: ACG
+ * 1: MLCG
+ * 2: MCellRan4
+ * 3: Isaac64
+ * 4: Random123
+ */
 class Rand {
   public:
-    Rand(unsigned long seed = 0, int size = 55, Object* obj = NULL);
+    Rand(unsigned long seed = 0, int size = 55, Object* obj = nullptr);
     ~Rand();
     RNG* gen;
     Random* rand;
     int type_;  // can do special things with some kinds of RNG
     // double* looks like random variable that gets changed on every fadvance
     Object* obj_;
 };
-
-#endif

src/oc/isaac64.cpp → src/gnu/isaac64.cpp

@@ -1,8 +1,3 @@
-#include <../../nrnconf.h>
-#if HAVE_SYS_TYPES_H
-#include <sys/types.h>
-#endif
-
 /*
 ------------------------------------------------------------------------------
 isaac64.cpp: A Fast cryptographic random number generator

src/oc/isaac64.h → src/gnu/isaac64.h

@@ -1,3 +1,5 @@
+#pragma once
+
 /*
 ------------------------------------------------------------------------------
 isaac64.h: Definitions for a fast cryptographic random number generator
@@ -12,10 +14,7 @@ Jenkins, R.J. (1996) ISAAC, in Fast Software Encryption, vol. 1039,
 
 ------------------------------------------------------------------------------
 */
-#ifndef ISAAC64_H
-#define ISAAC64_H
 
-#include <nrnconf.h>
 #include <stdint.h>
 
 #define RANDSIZL (4) /* I recommend 8 for crypto, 4 for simulations */
@@ -87,5 +86,3 @@ Macros to get individual random numbers
                          rng->randcnt = RANDMAX - 2,                                          \
                          DBL64 * (*((ub8*) (((ub4*) (rng->randrsl)) + rng->randcnt)))))
 
-
-#endif /* ISAAC64_H */

src/oc/mcran4.cpp → src/gnu/mcran4.cpp

@@ -39,17 +39,19 @@ contained the header:
     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */
 
-#include <../../nrnconf.h>
 #include <stdio.h>
 #include <stddef.h>
 #include <math.h>
 #include <float.h>
 #include <stdlib.h>
-#include <mcran4.h>
-#include "hocdec.h"
+#include "mcran4.h"
 
 static uint32_t lowindex = 0;
 
+double mcell_lowindex() {
+    return static_cast<double>(lowindex);
+}
+
 void mcell_ran4_init(uint32_t low) {
     lowindex = low;
 }
@@ -70,40 +72,6 @@ uint32_t mcell_iran4(uint32_t* high) {
     return nrnRan4int(high, lowindex);
 }
 
-/* Hoc interface */
-extern double chkarg();
-extern int use_mcell_ran4_;
-
-
-void hoc_mcran4() {
-    uint32_t idx;
-    double* xidx;
-    double x;
-    xidx = hoc_pgetarg(1);
-    idx = (uint32_t) (*xidx);
-    x = mcell_ran4a(&idx);
-    *xidx = idx;
-    hoc_ret();
-    hoc_pushx(x);
-}
-void hoc_mcran4init() {
-    double prev = (double) lowindex;
-    if (ifarg(1)) {
-        uint32_t idx = (uint32_t) chkarg(1, 0., 4294967295.);
-        mcell_ran4_init(idx);
-    }
-    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);
-}
-
 uint32_t nrnRan4int(uint32_t* idx1, uint32_t idx2) {
     uint32_t u, v, w, m, n;
     /* 64-bit hash */

src/oc/mcran4.h → src/gnu/mcran4.h

@@ -1,6 +1,7 @@
 #pragma once
 #include <cstdint>
 
+double mcell_lowindex();
 void mcell_ran4_init(uint32_t);
 double mcell_ran4(uint32_t* idx1, double* x, unsigned int n, double range);
 double mcell_ran4a(uint32_t* idx1);

src/gnu/nrnisaac.cpp

@@ -0,0 +1,25 @@
+#include <cstdint>
+#include "nrnisaac.h"
+#include "isaac64.h"
+
+using RNG = struct isaac64_state;
+
+void* nrnisaac_new() {
+    return new RNG;
+}
+
+void nrnisaac_delete(void* v) {
+    delete static_cast<RNG*>(v);
+}
+
+void nrnisaac_init(void* v, unsigned long int seed) {
+    isaac64_init(static_cast<RNG*>(v), seed);
+}
+
+double nrnisaac_dbl_pick(void* v) {
+    return isaac64_dbl32(static_cast<RNG*>(v));
+}
+
+std::uint32_t nrnisaac_uint32_pick(void* v) {
+    return isaac64_uint32(static_cast<RNG*>(v));
+}