From 378b5f1c3e80db0240a9692fb554d926d85f6530 Mon Sep 17 00:00:00 2001
From: Justin Ethier <justin.ethier@gmail.com>
Date: Mon, 20 Jun 2016 22:50:35 -0400
Subject: [PATCH] Experimenting with MRG32k3a

---
 include/cyclone/runtime.h |  1 +
 runtime.c                 | 55 +++++++++++++++++++++++++++++++++++++++
 srfi/27.sld               | 15 +++++++++--
 3 files changed, 69 insertions(+), 2 deletions(-)

diff --git a/include/cyclone/runtime.h b/include/cyclone/runtime.h
index 44d852bc..ce55e66e 100644
--- a/include/cyclone/runtime.h
+++ b/include/cyclone/runtime.h
@@ -441,4 +441,5 @@ void Cyc_rt_raise(void *data, object err);
 void Cyc_rt_raise2(void *data, const char *msg, object err);
 void Cyc_rt_raise_msg(void *data, const char *err);
 
+double MRG32k3a (double seed);
 #endif                          /* CYCLONE_RUNTIME_H */
diff --git a/runtime.c b/runtime.c
index 37663352..cd5a771e 100644
--- a/runtime.c
+++ b/runtime.c
@@ -4231,3 +4231,58 @@ object copy2heap(void *data, object obj)
   return gc_alloc(Cyc_heap, gc_allocated_bytes(obj, NULL, NULL), obj, data,
                   &on_stack);
 }
+
+/* RNG section */
+#define norm 2.328306549295728e-10
+#define m1   4294967087.0
+#define m2   4294944443.0
+#define a12     1403580.0
+#define a13n     810728.0
+#define a21      527612.0
+#define a23n    1370589.0
+
+/***
+The seeds for s10, s11, s12 must be integers in [0, m1 - 1] and not all 0. 
+The seeds for s20, s21, s22 must be integers in [0, m2 - 1] and not all 0. 
+***/
+
+//#define SEED 12345
+
+// JAE TODO: OK not to have these static?
+//static double s10 = SEED, s11 = SEED, s12 = SEED,
+//              s20 = SEED, s21 = SEED, s22 = SEED;
+
+
+double MRG32k3a (double seed)
+{
+   double s10 = seed, s11 = seed, s12 = seed,
+          s20 = seed, s21 = seed, s22 = seed;
+   long k;
+   double p1, p2;
+   /* Component 1 */
+   p1 = a12 * s11 - a13n * s10;
+   k = p1 / m1;
+   p1 -= k * m1;
+   if (p1 < 0.0)
+      p1 += m1;
+   s10 = s11;
+   s11 = s12;
+   s12 = p1;
+
+   /* Component 2 */
+   p2 = a21 * s22 - a23n * s20;
+   k = p2 / m2;
+   p2 -= k * m2;
+   if (p2 < 0.0)
+      p2 += m2;
+   s20 = s21;
+   s21 = s22;
+   s22 = p2;
+
+   /* Combination */
+   if (p1 <= p2)
+      return ((p1 - p2 + m1) * norm);
+   else
+      return ((p1 - p2) * norm);
+}
+/* END RNG */
diff --git a/srfi/27.sld b/srfi/27.sld
index 1c3a76d4..e439eda9 100644
--- a/srfi/27.sld
+++ b/srfi/27.sld
@@ -11,17 +11,18 @@
           (scheme case-lambda)
           (scheme time))
   (export random-integer random-real default-random-source
+          next-mrg32k3a ;; TODO: only here for testing
           make-random-source random-source?
           random-source-state-ref random-source-state-set!
           random-source-randomize! random-source-pseudo-randomize!
           random-source-make-integers random-source-make-reals)
   (begin
     ;; Numbers taken from bsd random
-    (define mult 1103515245)
+    ;(define mult 1103515245)
     (define incr 12345)
     (define m 536870912)
     ;; Cutting off seems like a good idea
-    (define cutoff 100)
+    ;(define cutoff 100)
 
     (define-c next-lcg
       "(void *data, int argc, closure _, object k, object seed)"
@@ -34,6 +35,16 @@
        result = obj_int2obj(next_seed);
        return_closcall1(data, k, result);")
 
+    ;; Testing this out
+    ;; TODO: handle ints, too. of course that also adds overhead...
+    (define-c next-mrg32k3a
+      "(void *data, int argc, closure _, object k, object seed)"
+      "double dval = MRG32k3a( double_value(seed) );
+       {
+        make_double(result, dval);
+        return_closcall1(data, k, &result);
+       }")
+
     (define-record-type <random-source>
       (raw-random-source n)
       random-souce?