Support exponential draws as well (closes #1)

Author: eddelbuettel

ChangeLog

@@ -1,3 +1,16 @@
+2025-01-02  Dirk Eddelbuettel  <[email protected]>
+
+	* inst/include/Ziggurat.h (kiss): Use an inline function avoiding a
+	compiler nag relative which KISS macro triggers
+	(rexp): Added, using inline kiss()
+	(efix): Added as well supporting rexp()
+	(init): Also initialize tables for exponential RNG
+
+	* src/ziggurat.cpp (zrexp): Added, calling rexp()
+	* man/ziggurat.Rd: Add alias for zrexp()
+	* src/RcppExports.cpp: Regenerated
+	* R/RcppExports.R: Idem
+
 2024-12-30  Dirk Eddelbuettel  <[email protected]>
 
 	* R/standardTest.R (.safeMaxCores): New helper function to provide

R/RcppExports.R

@@ -41,6 +41,10 @@ zrnorm <- function(n) {
     .Call(`_RcppZiggurat_zrnorm`, n)
 }
 
+zrexp <- function(n) {
+    .Call(`_RcppZiggurat_zrexp`, n)
+}
+
 zrnormVec <- function(x) {
     .Call(`_RcppZiggurat_zrnormVec`, x)
 }

inst/include/Ziggurat.h

@@ -1,8 +1,6 @@
-// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
-//
 // Ziggurat.h:  Marsaglia + Tsang / Leong, Zhang et al Ziggurat generator
 //
-// Copyright (C) 2013  Dirk Eddelbuettel
+// Copyright (C) 2013 - 2025  Dirk Eddelbuettel
 //
 // This file is part of RcppZiggurat.
 //
@@ -66,7 +64,8 @@ namespace Ziggurat {
 
 #define UNI  (0.5 + (int32_t) KISS * 0.2328306e-09)
 #define IUNI KISS
-#define RNOR (hz = KISS, iz = hz & 127, ( abs ( hz ) < kn[iz] ) ? hz * wn[iz] : nfix())
+//#define RNOR (hz = KISS, iz = hz & 127, ( abs ( hz ) < kn[iz] ) ? hz * wn[iz] : nfix())
+//#define REXP (jz = KISS, iz = jz & 255, (       jz   < ke[iz] ) ? jz * we[iz] : efix())
 
     class Ziggurat : public Zigg {
     public:
@@ -88,8 +87,14 @@ namespace Ziggurat {
         uint32_t getSeed() { 
             return jsr; 
         }
+        inline uint32_t kiss() {
+            return ((MWC ^ CONG ) + SHR3);
+        }
         inline double norm() {
-            return RNOR;
+            return (hz = kiss(), iz = hz & 127, ( abs ( hz ) < kn[iz] ) ? hz * wn[iz] : nfix());
+        }
+        inline double rexp() {
+            return (jz = kiss(), iz = jz & 255, (       jz   < ke[iz] ) ? jz * we[iz] : efix());
         }
         std::vector<uint32_t> getPars() {
             //C++11: std::vector<uint32_t> pars{ jsr, z, w, jcong };
@@ -109,23 +114,25 @@ namespace Ziggurat {
 
     private:
         double fn[128];
+        double fe[256];
         int32_t hz;
         uint32_t iz;
         uint32_t jcong;
         uint32_t jsr;
         uint32_t jz;
         uint32_t kn[128];
+        uint32_t ke[256];
         uint32_t w;
         double wn[128];
+        double we[256];
         uint32_t z;
 
         void init() {               // called from ctor, could be in ctor
-            double dn = 3.442619855899;
+            double dn = 3.442619855899, de=7.697117470131487;
             int i;
-            const double m1 = 2147483648.0;
-            double q;
-            double tn = 3.442619855899;
-            const double vn = 9.91256303526217E-03;
+            const double m1 = 2147483648.0, m2 = 4294967296.0;
+            double tn = dn, te = de, q;
+            const double vn = 9.91256303526217E-03, ve=3.949659822581572e-3;
 
             //  Set up the tables for the normal random number generator.
             q = vn / exp (- 0.5 * dn * dn);
@@ -145,6 +152,25 @@ namespace Ziggurat {
                 fn[i] = (double) (exp( - 0.5 * dn * dn));
                 wn[i] = (double) (dn / m1);
             }
+
+            //  Set up the tables for the exponential random number generator.
+            q = ve/exp(-de);
+            ke[0] = (uint32_t) (de/q)*m2;
+            ke[1]=0;
+
+            we[0] = (double) q/m2;
+            we[255] = (double) de/m2;
+
+            fe[0] = 1.0;
+            fe[255] = (double) exp(-de);
+
+            for (i=254; i>=1; i--) {
+                de = -log(ve/de+exp(-de));
+                ke[i+1] = (de/te) * m2;
+                te = de;
+                fe[i] = (double) exp(-de);
+                we[i] = (double) de/m2;
+            }
             return;
         }
 
@@ -177,6 +203,26 @@ namespace Ziggurat {
                 }
             }
         }
+        inline double efix(void) {
+            double x;
+            for ( ; ; ) {
+                //  IZ = 0.
+                if (iz == 0) {
+                    return (7.69711 - log (UNI));
+                }
+                x = jz * we[iz];
+                if (fe[iz] + UNI * (fe[iz-1] - fe[iz]) < exp (- x)) {
+                    return x;
+                }
+                //  Initiate, try to exit the loop.
+                jz = SHR3;      // note we still use SHR3; using KISS leads to degradation
+                iz = (jz & 255);
+                if (jz < ke[iz]) {
+                    return ((double) (jz * we[iz]));
+                }
+            }
+        }
+
     };
 
 #undef znew

man/ziggurat.Rd

@@ -1,5 +1,6 @@
 \name{zrnorm}
 \alias{zrnorm}
+\alias{zrexp}
 \alias{zrnormLZLLV}
 \alias{zrnormMT}
 \alias{zrnormV1}
@@ -36,6 +37,7 @@
 }
 \usage{
   zrnorm(n)
+  zrexp(n)
   zrnormLZLLV(n)
   zrnormMT(n)
   zrnormV1(n)

src/RcppExports.cpp

@@ -6,6 +6,11 @@
 
 using namespace Rcpp;
 
+#ifdef RCPP_USE_GLOBAL_ROSTREAM
+Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
+Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
+#endif
+
 // zrnormMT
 Rcpp::NumericVector zrnormMT(int n);
 RcppExport SEXP _RcppZiggurat_zrnormMT(SEXP nSEXP) {
@@ -112,6 +117,17 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// zrexp
+Rcpp::NumericVector zrexp(int n);
+RcppExport SEXP _RcppZiggurat_zrexp(SEXP nSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< int >::type n(nSEXP);
+    rcpp_result_gen = Rcpp::wrap(zrexp(n));
+    return rcpp_result_gen;
+END_RCPP
+}
 // zrnormVec
 Rcpp::NumericVector zrnormVec(Rcpp::NumericVector x);
 RcppExport SEXP _RcppZiggurat_zrnormVec(SEXP xSEXP) {
@@ -315,6 +331,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_RcppZiggurat_zsetseedV1", (DL_FUNC) &_RcppZiggurat_zsetseedV1, 1},
     {"_RcppZiggurat_zgetseedV1", (DL_FUNC) &_RcppZiggurat_zgetseedV1, 0},
     {"_RcppZiggurat_zrnorm", (DL_FUNC) &_RcppZiggurat_zrnorm, 1},
+    {"_RcppZiggurat_zrexp", (DL_FUNC) &_RcppZiggurat_zrexp, 1},
     {"_RcppZiggurat_zrnormVec", (DL_FUNC) &_RcppZiggurat_zrnormVec, 1},
     {"_RcppZiggurat_zrnormStl", (DL_FUNC) &_RcppZiggurat_zrnormStl, 1},
     {"_RcppZiggurat_zsetseed", (DL_FUNC) &_RcppZiggurat_zsetseed, 1},

src/ziggurat.cpp

@@ -130,6 +130,15 @@ Rcpp::NumericVector zrnorm(int n) {
     return x;
 }
 
+// [[Rcpp::export]]
+Rcpp::NumericVector zrexp(int n) {
+    Rcpp::NumericVector x(n);
+    for (int i=0; i<n; i++) {
+        x[i] = zigg.rexp();
+    }
+    return x;
+}
+
 // [[Rcpp::export]]
 Rcpp::NumericVector zrnormVec(Rcpp::NumericVector x) {
     int n = x.size();