@@ -82,7 +82,6 @@ uint32_t GapsRng::uniform32(uint32_t a, uint32_t b)

     {

         return a;

     }

-

     uint32_t range = b + 1 - a;

     uint32_t x = uniform32();

     uint32_t iPart = std::numeric_limits<uint32_t>::max() / range;

@@ -90,7 +89,10 @@ uint32_t GapsRng::uniform32(uint32_t a, uint32_t b)

     {

         x = uniform32();

     }

-    return x / iPart + a;

+    uint32_t unif = x / iPart + a;

+    GAPS_ASSERT(unif >= a);

+    GAPS_ASSERT(unif <= b);

+    return unif;

 }

 uint64_t GapsRng::uniform64()

@@ -107,7 +109,6 @@ uint64_t GapsRng::uniform64(uint64_t a, uint64_t b)

     {

         return a;

     }

-

     uint64_t range = b + 1 - a;

     uint64_t x = uniform64();

     uint64_t iPart = std::numeric_limits<uint64_t>::max() / range;

@@ -115,7 +116,10 @@ uint64_t GapsRng::uniform64(uint64_t a, uint64_t b)

     {

         x = uniform64();

     }

-    return x / iPart + a;

+    uint64_t unif = x / iPart + a;

+    GAPS_ASSERT(unif >= a);

+    GAPS_ASSERT(unif <= b);

+    return unif;

 }

 int GapsRng::poisson(double lambda)

@@ -144,7 +148,6 @@ int GapsRng::poissonLarge(double lambda)

     double beta = gaps::pi_double / sqrt(3.0 * lambda);

     double alpha = beta * lambda;

     double k = std::log(c) - lambda - std::log(beta);

-

     while(true)

     {

         double u = uniformd();

@@ -154,7 +157,6 @@ int GapsRng::poissonLarge(double lambda)

         {

             continue;

         }

-

         double v = uniformd();

         double y = alpha - beta * x;

         double w = 1.0 + std::exp(y);

@@ -1,11 +1,12 @@

 #include "Random.h"

 #include "Math.h"

+#include "../utils/Archive.h"

 #include "../utils/GapsAssert.h"

 #include <algorithm>

 #include <stdint.h>

-#include <limits>

 #include <cmath>

+#include <limits>

 const float maxU32AsFloat = static_cast<float>(std::numeric_limits<uint32_t>::max());

 const double maxU32AsDouble = static_cast<double>(std::numeric_limits<uint32_t>::max());

@@ -180,7 +181,6 @@ OptionalFloat GapsRng::truncNormal(float a, float b, float mean, float sd)

     {

         GAPS_ASSERT(pLower > 0.f);

         GAPS_ASSERT(pUpper < 1.f);

-

         float z = mRandState->q_norm_fast(uniform(pLower, pUpper), mean, sd); 

         z = gaps::max(a, gaps::min(z, b));

         return z;

@@ -192,7 +192,8 @@ OptionalFloat GapsRng::truncNormal(float a, float b, float mean, float sd)

 float GapsRng::truncGammaUpper(float b, float scale)

 {

     float upper = 1.f - std::exp(-b / scale) * (1.f + b / scale);

-    unsigned ndx = static_cast<unsigned>(uniform(0.f, upper * 5000.f));

+    const unsigned ndx = static_cast<unsigned>(uniform(0.f, upper * 5000.f));

+    GAPS_ASSERT(ndx < Q_GAMMA_LOOKUP_TABLE_SIZE);

     return mRandState->mQgammaLookupTable[ndx] * scale;

 }

@@ -219,9 +220,7 @@ Xoroshiro128plus::Xoroshiro128plus(uint64_t seed)

 {

     mState[0] = seed|1;

     mState[1] = seed|1;

-

-    // warmup

-    for (unsigned i = 0; i < 5000; ++i)

+    for (unsigned i = 0; i < 5000; ++i) // warmup

     {

         next();

     }

@@ -231,7 +230,6 @@ uint64_t Xoroshiro128plus::next()

 {

     mPreviousState[0] = mState[0];

     mPreviousState[1] = mState[1];

-

     const uint64_t s0 = mState[0];

     uint64_t s1 = mState[1];

     uint64_t result = s0 + s1;

@@ -269,29 +267,29 @@ GapsRandomState::GapsRandomState(unsigned seed) : mSeeder(seed)

 void GapsRandomState::initLookupTables()

 {

     // erf

-    for (unsigned i = 0; i < 3001; ++i)

+    for (unsigned i = 0; i < ERF_LOOKUP_TABLE_SIZE; ++i)

     {

         float x = static_cast<float>(i) / 1000.f;

         mErfLookupTable[i] = 2.f * gaps::p_norm(x * gaps::sqrt2, 0.f, 1.f) - 1.f;

     }

-    GAPS_ASSERT(mErfLookupTable[3000] < 1.f);

+    GAPS_ASSERT(mErfLookupTable[ERF_LOOKUP_TABLE_SIZE - 1] < 1.f);

     // erfinv

-    for (unsigned i = 0; i < 5000; ++i)

+    for (unsigned i = 0; i < ERF_INV_LOOKUP_TABLE_SIZE - 1; ++i)

     {

-        float x = static_cast<float>(i) / 5000.f;

+        float x = static_cast<float>(i) / static_cast<float>(ERF_INV_LOOKUP_TABLE_SIZE - 1);

         mErfinvLookupTable[i] = gaps::q_norm((1.f + x) / 2.f, 0.f, 1.f) / gaps::sqrt2;

     }

-    mErfinvLookupTable[5000] = gaps::q_norm(1.9998f / 2.f, 0.f, 1.f) / gaps::sqrt2;

+    mErfinvLookupTable[ERF_INV_LOOKUP_TABLE_SIZE - 1] = gaps::q_norm(1.9998f / 2.f, 0.f, 1.f) / gaps::sqrt2;

     // qgamma

     mQgammaLookupTable[0] = 0.f;

-    for (unsigned i = 1; i < 5000; ++i)

+    for (unsigned i = 1; i < Q_GAMMA_LOOKUP_TABLE_SIZE - 1; ++i)

     {

-        float x = static_cast<float>(i) / 5000.f;

+        float x = static_cast<float>(i) / static_cast<float>(Q_GAMMA_LOOKUP_TABLE_SIZE - 1);

         mQgammaLookupTable[i] = gaps::q_gamma(x, 2.f, 1.f);

     }

-    mQgammaLookupTable[5000] = gaps::q_gamma(0.9998f, 2.f, 1.f);

+    mQgammaLookupTable[Q_GAMMA_LOOKUP_TABLE_SIZE - 1] = gaps::q_gamma(0.9998f, 2.f, 1.f);

 }

 uint64_t GapsRandomState::nextSeed()

@@ -311,12 +309,16 @@ float GapsRandomState::p_norm_fast(float p, float mean, float sd) const

     if (term < 0.f)

     {

         term = gaps::max(term, -3.f);

-        erf = -mErfLookupTable[static_cast<unsigned>(-term * 1000.f)];

+        const unsigned ndx = static_cast<unsigned>(-term * 1000.f);

+        GAPS_ASSERT(ndx < ERF_LOOKUP_TABLE_SIZE);

+        erf = -mErfLookupTable[ndx];

     }

     else

     {

         term = gaps::min(term, 3.f);

-        erf = mErfLookupTable[static_cast<unsigned>(term * 1000.f)];

+        const unsigned ndx = static_cast<unsigned>(term * 1000.f);

+        GAPS_ASSERT(ndx < ERF_LOOKUP_TABLE_SIZE);

+        erf = mErfLookupTable[ndx];

     }

     return 0.5f * (1.f + erf);

 }

@@ -327,11 +329,15 @@ float GapsRandomState::q_norm_fast(float q, float mean, float sd) const

     float erfinv = 0.f;

     if (term < 0.f)

     {

-        erfinv = -mErfinvLookupTable[static_cast<unsigned>(-term * 5000.f)];

+        const unsigned ndx = static_cast<unsigned>(-term * (ERF_INV_LOOKUP_TABLE_SIZE - 1));

+        GAPS_ASSERT(ndx < ERF_INV_LOOKUP_TABLE_SIZE);

+        erfinv = -mErfinvLookupTable[ndx];

     }

     else

     {

-        erfinv = mErfinvLookupTable[static_cast<unsigned>(term * 5000.f)];

+        const unsigned ndx = static_cast<unsigned>(term * (ERF_INV_LOOKUP_TABLE_SIZE - 1));

+        GAPS_ASSERT(ndx < ERF_INV_LOOKUP_TABLE_SIZE);

+        erfinv = mErfinvLookupTable[ndx];

     }

     return mean + sd * gaps::sqrt2 * erfinv;

 }

@@ -4,6 +4,8 @@

 #include <algorithm>

 #include <stdint.h>

+#include <limits>

+#include <cmath>

 const float maxU32AsFloat = static_cast<float>(std::numeric_limits<uint32_t>::max());

 const double maxU32AsDouble = static_cast<double>(std::numeric_limits<uint32_t>::max());

old mode 100644

new mode 100755

@@ -42,7 +42,6 @@ uint32_t GapsRng::next()

 void GapsRng::advance()

 {

-    mPreviousState = mState;

     mState = mState * 6364136223846793005ull + (54u|1);

 }

@@ -28,8 +28,8 @@ bool OptionalFloat::hasValue() const

 GapsRng::GapsRng(GapsRandomState *randState)

 :

-mState(randState->nextSeed()),

-mRandState(randState)

+mRandState(randState),

+mState(randState->nextSeed())

 {

     advance();

 }

@@ -1,61 +1,13 @@

-// [[Rcpp::depends(BH)]]

-

 #include "Random.h"

-

 #include "Math.h"

 #include "../utils/GapsAssert.h"

-#include <boost/math/distributions/exponential.hpp>

-#include <boost/math/distributions/gamma.hpp>

-#include <boost/math/distributions/normal.hpp>

-

 #include <algorithm>

 #include <stdint.h>

-#define Q_GAMMA_THRESHOLD 0.000001f

-#define Q_GAMMA_MIN_VALUE 0.f

-

 const float maxU32AsFloat = static_cast<float>(std::numeric_limits<uint32_t>::max());

 const double maxU32AsDouble = static_cast<double>(std::numeric_limits<uint32_t>::max());

-static Xoroshiro128plus seeder;

-

-////////////////////////////// Lookup Tables ///////////////////////////////////

-

-static bool tables_are_initialized = false;

-static float erf_lookup_table[3001];

-static float erfinv_lookup_table[5001];

-static float qgamma_lookup_table[5001];

-

-static void initLookupTables()

-{

-    // erf

-    for (unsigned i = 0; i < 3001; ++i)

-    {

-        float x = static_cast<float>(i) / 1000.f;

-        erf_lookup_table[i] = 2.f * gaps::p_norm(x * gaps::sqrt2, 0.f, 1.f) - 1.f;

-    }

-

-    // erfinv

-    for (unsigned i = 0; i < 5000; ++i)

-    {

-        float x = static_cast<float>(i) / 5000.f;

-        erfinv_lookup_table[i] = gaps::q_norm((1.f + x) / 2.f, 0.f, 1.f) / gaps::sqrt2;

-    }

-    erfinv_lookup_table[5000] = gaps::q_norm(1.9998f / 2.f, 0.f, 1.f) / gaps::sqrt2;

-

-    // qgamma

-    qgamma_lookup_table[0] = 0.f;

-    for (unsigned i = 1; i < 5000; ++i)

-    {

-        float x = static_cast<float>(i) / 5000.f;

-        qgamma_lookup_table[i] = gaps::q_gamma(x, 2.f, 1.f);

-    }

-    qgamma_lookup_table[5000] = gaps::q_gamma(0.9998f, 2.f, 1.f);

-

-    GAPS_ASSERT(erf_lookup_table[3000] < 1.f);

-}

-

 /////////////////////////////// OptionalFloat //////////////////////////////////

 OptionalFloat::OptionalFloat() : mValue(0.f), mHasValue(false) {}

@@ -72,92 +24,14 @@ bool OptionalFloat::hasValue() const

     return mHasValue;

 }

-///////////////////////////// Xoroshiro128plus /////////////////////////////////

-

-static uint64_t rotl(const uint64_t x, int k)

-{

-    return (x << k) | (x >> (64 - k));

-}

-

-void Xoroshiro128plus::seed(uint64_t seed)

-{

-    mState[0] = seed|1;

-    mState[1] = seed|1;

-    warmup();

-}

-

-uint64_t Xoroshiro128plus::next()

-{

-    mPreviousState[0] = mState[0];

-    mPreviousState[1] = mState[1];

-

-    const uint64_t s0 = mState[0];

-    uint64_t s1 = mState[1];

-    uint64_t result = s0 + s1;

-    s1 ^= s0;

-    mState[0] = rotl(s0, 24) ^ s1 ^ (s1 << 16); // a, b

-    mState[1] = rotl(s1, 37); // c

-    return result;

-}

-

-void Xoroshiro128plus::rollBackOnce()

-{

-    mState[0] = mPreviousState[0];

-    mState[1] = mPreviousState[1];

-}

-

-void Xoroshiro128plus::warmup()

-{

-    for (unsigned i = 0; i < 5000; ++i)

-    {

-        next();

-    }

-}

-

-Archive& operator<<(Archive &ar, Xoroshiro128plus &gen)

-{

-    ar << gen.mState[0] << gen.mState[1];

-    return ar;

-}

-

-Archive& operator>>(Archive &ar, Xoroshiro128plus &gen)

-{

-    ar >> gen.mState[0] >> gen.mState[1];

-    return ar;

-}

-

-////////////////////////////////// GapsRng /////////////////////////////////////

-

-void GapsRng::setSeed(uint32_t sd)

-{

-    GAPS_ASSERT(!tables_are_initialized);

-    initLookupTables();

-    tables_are_initialized = true;

-    seeder.seed(sd);

-}

-

-Archive& GapsRng::save(Archive &ar)

-{

-    ar << seeder;

-    return ar;

-}

+//////////////////////////////// GapsRng ///////////////////////////////////////

-Archive& GapsRng::load(Archive &ar)

+GapsRng::GapsRng(GapsRandomState *randState)

+:

+mState(randState->nextSeed()),

+mRandState(randState)

 {

-    initLookupTables();

-    ar >> seeder;

-    return ar;

-}

-

-void GapsRng::rollBackOnce()

-{

-    seeder.rollBackOnce();

-}

-

-GapsRng::GapsRng()

-    : mState(seeder.next())

-{

-    next();

+    advance();

 }

 uint32_t GapsRng::next()

@@ -299,14 +173,14 @@ float GapsRng::exponential(float lambda)

 // fails if too far in tail

 OptionalFloat GapsRng::truncNormal(float a, float b, float mean, float sd)

 {

-    float pLower = gaps::p_norm_fast(a, mean, sd);

-    float pUpper = gaps::p_norm_fast(b, mean, sd);

+    float pLower = mRandState->p_norm_fast(a, mean, sd);

+    float pUpper = mRandState->p_norm_fast(b, mean, sd);

     if (!(pLower > 0.95f || pUpper < 0.05f)) // too far in tail

     {

         GAPS_ASSERT(pLower > 0.f);

         GAPS_ASSERT(pUpper < 1.f);

-        float z = gaps::q_norm_fast(uniform(pLower, pUpper), mean, sd); 

+        float z = mRandState->q_norm_fast(uniform(pLower, pUpper), mean, sd); 

         z = gaps::max(a, gaps::min(z, b));

         return z;

     }

@@ -318,10 +192,10 @@ float GapsRng::truncGammaUpper(float b, float scale)

 {

     float upper = 1.f - std::exp(-b / scale) * (1.f + b / scale);

     unsigned ndx = static_cast<unsigned>(uniform(0.f, upper * 5000.f));

-    return qgamma_lookup_table[ndx] * scale;

+    return mRandState->mQgammaLookupTable[ndx] * scale;

 }

-Archive& operator<<(Archive &ar, GapsRng &gen)

+Archive& operator<<(Archive &ar, const GapsRng &gen)

 {

     ar << gen.mState;

     return ar;

@@ -333,81 +207,142 @@ Archive& operator>>(Archive &ar, GapsRng &gen)

     return ar;

 }

-////////////////////////// Distribution Calculations ///////////////////////////

+///////////////////////////// Xoroshiro128plus /////////////////////////////////

-float gaps::d_gamma(float d, float shape, float scale)

+static uint64_t rotl(const uint64_t x, int k)

 {

-    boost::math::gamma_distribution<> gam(shape, scale);

-    return pdf(gam, d);

+    return (x << k) | (x >> (64 - k));

 }

-float gaps::p_gamma(float p, float shape, float scale)

+Xoroshiro128plus::Xoroshiro128plus(uint64_t seed)

 {

-    boost::math::gamma_distribution<> gam(shape, scale);

-    return cdf(gam, p);

-}

+    mState[0] = seed|1;

+    mState[1] = seed|1;

-float gaps::q_gamma(float q, float shape, float scale)

-{

-    if (q < Q_GAMMA_THRESHOLD)

+    // warmup

+    for (unsigned i = 0; i < 5000; ++i)

     {

-        return Q_GAMMA_MIN_VALUE;

+        next();

     }

-    boost::math::gamma_distribution<> gam(shape, scale);

-    return quantile(gam, q);

 }

-float gaps::d_norm(float d, float mean, float sd)

+uint64_t Xoroshiro128plus::next()

 {

-    boost::math::normal_distribution<> norm(mean, sd);

-    return pdf(norm, d);

+    mPreviousState[0] = mState[0];

+    mPreviousState[1] = mState[1];

+

+    const uint64_t s0 = mState[0];

+    uint64_t s1 = mState[1];

+    uint64_t result = s0 + s1;

+    s1 ^= s0;

+    mState[0] = rotl(s0, 24) ^ s1 ^ (s1 << 16); // a, b

+    mState[1] = rotl(s1, 37); // c

+    return result;

 }

-float gaps::q_norm(float q, float mean, float sd)

+void Xoroshiro128plus::rollBackOnce()

 {

-    boost::math::normal_distribution<> norm(mean, sd);

-    return quantile(norm, q);

+    mState[0] = mPreviousState[0];

+    mState[1] = mPreviousState[1];

 }

-float gaps::p_norm(float p, float mean, float sd)

+Archive& operator<<(Archive &ar, const Xoroshiro128plus &gen)

 {

-    boost::math::normal_distribution<> norm(mean, sd);

-    return cdf(norm, p);

+    ar << gen.mState[0] << gen.mState[1];

+    return ar;

 }

-float gaps::p_norm_fast(float p, float mean, float sd)

+Archive& operator>>(Archive &ar, Xoroshiro128plus &gen)

+{

+    ar >> gen.mState[0] >> gen.mState[1];

+    return ar;

+}

+

+///////////////////////////// GapsRandomState //////////////////////////////////

+

+GapsRandomState::GapsRandomState(unsigned seed) : mSeeder(seed)

+{

+    initLookupTables();

+}

+

+void GapsRandomState::initLookupTables()

+{

+    // erf

+    for (unsigned i = 0; i < 3001; ++i)

+    {

+        float x = static_cast<float>(i) / 1000.f;

+        mErfLookupTable[i] = 2.f * gaps::p_norm(x * gaps::sqrt2, 0.f, 1.f) - 1.f;

+    }

+    GAPS_ASSERT(mErfLookupTable[3000] < 1.f);

+

+    // erfinv

+    for (unsigned i = 0; i < 5000; ++i)

+    {

+        float x = static_cast<float>(i) / 5000.f;

+        mErfinvLookupTable[i] = gaps::q_norm((1.f + x) / 2.f, 0.f, 1.f) / gaps::sqrt2;

+    }

+    mErfinvLookupTable[5000] = gaps::q_norm(1.9998f / 2.f, 0.f, 1.f) / gaps::sqrt2;

+

+    // qgamma

+    mQgammaLookupTable[0] = 0.f;

+    for (unsigned i = 1; i < 5000; ++i)

+    {

+        float x = static_cast<float>(i) / 5000.f;

+        mQgammaLookupTable[i] = gaps::q_gamma(x, 2.f, 1.f);

+    }

+    mQgammaLookupTable[5000] = gaps::q_gamma(0.9998f, 2.f, 1.f);

+}

+

+uint64_t GapsRandomState::nextSeed()

+{

+    return mSeeder.next();

+}

+

+void GapsRandomState::rollBackOnce()

+{

+    mSeeder.rollBackOnce();

+}

+

+float GapsRandomState::p_norm_fast(float p, float mean, float sd) const

 {

     float term = (p - mean) / (sd * gaps::sqrt2);

     float erf = 0.f;

     if (term < 0.f)

     {

         term = gaps::max(term, -3.f);

-        erf = -erf_lookup_table[static_cast<unsigned>(-term * 1000.f)];

+        erf = -mErfLookupTable[static_cast<unsigned>(-term * 1000.f)];

     }

     else

     {

         term = gaps::min(term, 3.f);

-        erf = erf_lookup_table[static_cast<unsigned>(term * 1000.f)];

+        erf = mErfLookupTable[static_cast<unsigned>(term * 1000.f)];

     }

     return 0.5f * (1.f + erf);

 }

-float gaps::q_norm_fast(float q, float mean, float sd)

+float GapsRandomState::q_norm_fast(float q, float mean, float sd) const

 {

     float term = 2.f * q - 1.f;

     float erfinv = 0.f;

     if (term < 0.f)

     {

-        erfinv = -erfinv_lookup_table[static_cast<unsigned>(-term * 5000.f)];

+        erfinv = -mErfinvLookupTable[static_cast<unsigned>(-term * 5000.f)];

     }

     else

     {

-        erfinv = erfinv_lookup_table[static_cast<unsigned>(term * 5000.f)];

+        erfinv = mErfinvLookupTable[static_cast<unsigned>(term * 5000.f)];

     }

     return mean + sd * gaps::sqrt2 * erfinv;

 }

-double gaps::lgamma(double x)

+Archive& operator<<(Archive &ar, const GapsRandomState &s)

 {

-    return boost::math::lgamma(x);

+    ar << s.mSeeder;

+    return ar;

+}

+

+Archive& operator>>(Archive &ar, GapsRandomState &s)

+{

+    ar >> s.mSeeder;

+    return ar;

 }

\ No newline at end of file

@@ -22,6 +22,7 @@ static Xoroshiro128plus seeder;

 ////////////////////////////// Lookup Tables ///////////////////////////////////

+static bool tables_are_initialized = false;

 static float erf_lookup_table[3001];

 static float erfinv_lookup_table[5001];

 static float qgamma_lookup_table[5001];

@@ -129,7 +130,9 @@ Archive& operator>>(Archive &ar, Xoroshiro128plus &gen)

 void GapsRng::setSeed(uint32_t sd)

 {

+    GAPS_ASSERT(!tables_are_initialized);

     initLookupTables();

+    tables_are_initialized = true;

     seeder.seed(sd);

 }

@@ -2,7 +2,7 @@

 #include "Random.h"

-#include "Algorithms.h"

+#include "Math.h"

 #include "../utils/GapsAssert.h"

 #include <boost/math/distributions/exponential.hpp>

@@ -20,6 +20,41 @@ const double maxU32AsDouble = static_cast<double>(std::numeric_limits<uint32_t>:

 static Xoroshiro128plus seeder;

+////////////////////////////// Lookup Tables ///////////////////////////////////

+

+static float erf_lookup_table[3001];

+static float erfinv_lookup_table[5001];

+static float qgamma_lookup_table[5001];

+

+static void initLookupTables()

+{

+    // erf

+    for (unsigned i = 0; i < 3001; ++i)

+    {

+        float x = static_cast<float>(i) / 1000.f;

+        erf_lookup_table[i] = 2.f * gaps::p_norm(x * gaps::sqrt2, 0.f, 1.f) - 1.f;

+    }

+

+    // erfinv

+    for (unsigned i = 0; i < 5000; ++i)

+    {

+        float x = static_cast<float>(i) / 5000.f;

+        erfinv_lookup_table[i] = gaps::q_norm((1.f + x) / 2.f, 0.f, 1.f) / gaps::sqrt2;

+    }

+    erfinv_lookup_table[5000] = gaps::q_norm(1.9998f / 2.f, 0.f, 1.f) / gaps::sqrt2;

+

+    // qgamma

+    qgamma_lookup_table[0] = 0.f;

+    for (unsigned i = 1; i < 5000; ++i)

+    {

+        float x = static_cast<float>(i) / 5000.f;

+        qgamma_lookup_table[i] = gaps::q_gamma(x, 2.f, 1.f);

+    }

+    qgamma_lookup_table[5000] = gaps::q_gamma(0.9998f, 2.f, 1.f);

+

+    GAPS_ASSERT(erf_lookup_table[3000] < 1.f);

+}

+

 /////////////////////////////// OptionalFloat //////////////////////////////////

 OptionalFloat::OptionalFloat() : mValue(0.f), mHasValue(false) {}

@@ -94,6 +129,7 @@ Archive& operator>>(Archive &ar, Xoroshiro128plus &gen)

 void GapsRng::setSeed(uint32_t sd)

 {

+    initLookupTables();

     seeder.seed(sd);

 }

@@ -105,6 +141,7 @@ Archive& GapsRng::save(Archive &ar)

 Archive& GapsRng::load(Archive &ar)

 {

+    initLookupTables();

     ar >> seeder;

     return ar;

 }

@@ -204,14 +241,7 @@ uint64_t GapsRng::uniform64(uint64_t a, uint64_t b)

 int GapsRng::poisson(double lambda)

 {

-    if (lambda <= 5.0)

-    {

-        return poissonSmall(lambda);

-    }

-    else

-    {

-        return poissonLarge(lambda);

-    }

+    return lambda <= 5.0 ? poissonSmall(lambda) : poissonLarge(lambda);

 }

 // lambda <= 5

@@ -263,33 +293,29 @@ float GapsRng::exponential(float lambda)

     return -1.f * std::log(uniform()) / lambda;

 }

+// fails if too far in tail

 OptionalFloat GapsRng::truncNormal(float a, float b, float mean, float sd)

 {

-    float pLower = gaps::p_norm(a, mean, sd);

-    float pUpper = gaps::p_norm(b, mean, sd);

+    float pLower = gaps::p_norm_fast(a, mean, sd);

+    float pUpper = gaps::p_norm_fast(b, mean, sd);

     if (!(pLower > 0.95f || pUpper < 0.05f)) // too far in tail

     {

-        float u = uniform(pLower, pUpper);

-        while (u == 0.f || u == 1.f)

-        {

-            u = uniform(pLower, pUpper);

-        }

-        float m = gaps::q_norm(u, mean, sd); 

-        m = gaps::max(a, gaps::min(m, b));

-        return m;

+        GAPS_ASSERT(pLower > 0.f);

+        GAPS_ASSERT(pUpper < 1.f);

+

+        float z = gaps::q_norm_fast(uniform(pLower, pUpper), mean, sd); 

+        z = gaps::max(a, gaps::min(z, b));

+        return z;

     }

     return OptionalFloat();

 }

-float GapsRng::truncGammaUpper(float b, float shape, float scale)

+// shape is hardcoded to 2 since it never changes

+float GapsRng::truncGammaUpper(float b, float scale)

 {

-    float upper = gaps::p_gamma(b, shape, scale);

-    float u = uniform(0.f, upper);

-    while (u == 0.f || u == 1.f)

-    {

-        u = uniform(0.f, upper);

-    }

-    return gaps::q_gamma(u, shape, scale);

+    float upper = 1.f - std::exp(-b / scale) * (1.f + b / scale);

+    unsigned ndx = static_cast<unsigned>(uniform(0.f, upper * 5000.f));

+    return qgamma_lookup_table[ndx] * scale;

 }

 Archive& operator<<(Archive &ar, GapsRng &gen)

@@ -346,13 +372,39 @@ float gaps::p_norm(float p, float mean, float sd)

     return cdf(norm, p);

 }

-double gaps::lgamma(double x)

+float gaps::p_norm_fast(float p, float mean, float sd)

 {

-    return boost::math::lgamma(x);

+    float term = (p - mean) / (sd * gaps::sqrt2);

+    float erf = 0.f;

+    if (term < 0.f)

+    {

+        term = gaps::max(term, -3.f);

+        erf = -erf_lookup_table[static_cast<unsigned>(-term * 1000.f)];

+    }

+    else

+    {

+        term = gaps::min(term, 3.f);

+        erf = erf_lookup_table[static_cast<unsigned>(term * 1000.f)];

+    }

+    return 0.5f * (1.f + erf);

 }

-float gaps::d_norm_fast(float d, float mean, float sd)

+float gaps::q_norm_fast(float q, float mean, float sd)

 {

-    return std::exp((d - mean) * (d - mean) / (-2.f * sd * sd))

-        / std::sqrt(2.f * gaps::pi * sd * sd);

+    float term = 2.f * q - 1.f;

+    float erfinv = 0.f;

+    if (term < 0.f)

+    {

+        erfinv = -erfinv_lookup_table[static_cast<unsigned>(-term * 5000.f)];

+    }

+    else

+    {

+        erfinv = erfinv_lookup_table[static_cast<unsigned>(term * 5000.f)];

+    }

+    return mean + sd * gaps::sqrt2 * erfinv;

+}

+

+double gaps::lgamma(double x)

+{

+    return boost::math::lgamma(x);

 }

\ No newline at end of file

@@ -263,14 +263,22 @@ float GapsRng::exponential(float lambda)

     return -1.f * std::log(uniform()) / lambda;

 }

-float GapsRng::inverseNormSample(float a, float b, float mean, float sd)

+OptionalFloat GapsRng::truncNormal(float a, float b, float mean, float sd)

 {

-    float u = uniform(a, b);

-    while (u == 0.f || u == 1.f)

+    float pLower = gaps::p_norm(a, mean, sd);

+    float pUpper = gaps::p_norm(b, mean, sd);

+    if (!(pLower > 0.95f || pUpper < 0.05f)) // too far in tail

     {

-        u = uniform(a, b);

+        float u = uniform(pLower, pUpper);

+        while (u == 0.f || u == 1.f)

+        {

+            u = uniform(pLower, pUpper);

+        }

+        float m = gaps::q_norm(u, mean, sd); 

+        m = gaps::max(a, gaps::min(m, b));

+        return m;

     }

-    return gaps::q_norm(u, mean, sd);

+    return OptionalFloat();

 }

 float GapsRng::truncGammaUpper(float b, float shape, float scale)

@@ -22,9 +22,9 @@ static Xoroshiro128plus seeder;

 /////////////////////////////// OptionalFloat //////////////////////////////////

-OptionalFloat::OptionalFloat() : mHasValue(false), mValue(0.f) {}

+OptionalFloat::OptionalFloat() : mValue(0.f), mHasValue(false) {}

-OptionalFloat::OptionalFloat(float f) : mHasValue(true), mValue(f) {}

+OptionalFloat::OptionalFloat(float f) : mValue(f), mHasValue(true) {}

 float OptionalFloat::value()

 {

@@ -116,7 +116,9 @@ void GapsRng::rollBackOnce()

 GapsRng::GapsRng()

     : mState(seeder.next())

-{}

+{

+    next();

+}

 uint32_t GapsRng::next()

 {

@@ -52,6 +52,9 @@ void Xoroshiro128plus::seed(uint64_t seed)

 uint64_t Xoroshiro128plus::next()

 {

+    mPreviousState[0] = mState[0];

+    mPreviousState[1] = mState[1];

+

     const uint64_t s0 = mState[0];

     uint64_t s1 = mState[1];

     uint64_t result = s0 + s1;

@@ -61,6 +64,12 @@ uint64_t Xoroshiro128plus::next()

     return result;

 }

+void Xoroshiro128plus::rollBackOnce()

+{

+    mState[0] = mPreviousState[0];

+    mState[1] = mPreviousState[1];

+}

+

 void Xoroshiro128plus::warmup()

 {

     for (unsigned i = 0; i < 5000; ++i)

@@ -100,15 +109,15 @@ Archive& GapsRng::load(Archive &ar)

     return ar;

 }

-GapsRng::GapsRng()

-    : mState(0)

-{}

-

-void GapsRng::init()

+void GapsRng::rollBackOnce()

 {

-    mState = seeder.next();

+    seeder.rollBackOnce();

 }