From f0e482ad78910e93be215154d1cf9aae796a59ac Mon Sep 17 00:00:00 2001
From: Ichthyostega <prg@ichthyostega.de>
Date: Mon, 31 Dec 2018 09:05:45 +0100
Subject: [PATCH] Yoshimi: investigation of the random number algorithm from
 the C standard lib

(ab)using the Lumiera tree here for research work on behalf of the Yoshimi project
For context, we stumbled over sonic changes due to using different random number algorighims,
in spite of all those algorithms producing mathematically sane numbers
---
 research/try.cpp | 159 ++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 138 insertions(+), 21 deletions(-)

diff --git a/research/try.cpp b/research/try.cpp
index 8f8108be1..a52a2ea15 100644
--- a/research/try.cpp
+++ b/research/try.cpp
@@ -40,13 +40,13 @@
 // 04/18 - investigate construction of static template members
 // 08/18 - Segfault when compiling some regular expressions for EventLog search
 // 10/18 - investigate insidious reinterpret cast
+// 12/18 - investigate the trinomial random number algorithm from the C standard lib
 
 
 /** @file try.cpp
- * Document an insidious wild cast, caused by the syntax `Type(arg)`.
- * I was under the wrong assumption this would be handled equivalent to a constructor invocation.
- * Seemingly it is rather handled as a C-style cast, i.e. equivalent to `(Type)arg`.
- * @see [Question on Stackoverflow](https://stackoverflow.com/q/52782967/444796)
+ * Investigate the trinomial random number algorithm from the C standard library (actually GLibc 2.28).
+ * Actually this is work for the yoshimi project; we try there to build an in-tree version of the PRNG,
+ * in order to reduce dependencies to external libraries, which might change the sound of existing synth patches.
  */
 
 typedef unsigned int uint;
@@ -56,9 +56,10 @@ typedef unsigned int uint;
 #include "lib/util.hpp"
 
 #include <string>
+#include <boost/lexical_cast.hpp>
 
 using std::string;
-using util::isSameObject;
+using boost::lexical_cast;
 
 
 
@@ -67,35 +68,151 @@ using util::isSameObject;
 #define SHOW_EXPR(_XX_) \
     cout << "Probe " << STRINGIFY(_XX_) << " ? = " << _XX_ <<endl;
 
-
+namespace {
   
+class StdlibPRNG
+{
+        char random_state[256];
+        struct random_data random_buf;
 
-class Wau
-  {
-    int i = -1;
-  };
+    public:
+        StdlibPRNG()
+        {
+            memset(&random_state, 0, sizeof(random_state));
+        }
 
-class Miau
-  {
-  public:
-    uint u = 1;
-  };
+        bool init(uint32_t seed)
+        {
+            memset(random_state, 0, sizeof(random_state));
+            memset(&random_buf, 0, sizeof(random_buf));
+            return 0 == initstate_r(seed, random_state, sizeof(random_state), &random_buf);
+        }
+
+        uint32_t prngval()
+        {
+            int32_t random_result;
+            random_r(&random_buf, &random_result);
+            // can not fail, since &random_buf can not be NULL
+            // random_result holds number 0...INT_MAX
+            return random_result;
+        }
+
+        float numRandom()
+        {
+            return prngval() / float(INT32_MAX);
+        }
+
+        // random number in the range 0...INT_MAX
+        uint32_t randomINT()
+        {
+            return prngval();
+        }
+};
+
+class TrinomialPRNG
+{
+        int32_t state[63];
+        int32_t *fptr;      /* Front pointer.  */
+        int32_t *rptr;      /* Rear pointer.  */
+
+    public:
+        TrinomialPRNG() : fptr(NULL), rptr(NULL) { }
+
+        bool init(uint32_t seed)
+        {
+            int kc = 63; /* random generation uses this trinomial: x**63 + x + 1.  */
+
+            /* We must make sure the seed is not 0.  Take arbitrarily 1 in this case.  */
+            if (seed == 0)
+              seed = 1;
+            state[0] = seed;
+
+            int32_t *dst = state;
+            int32_t word = seed;  // must be signed, see below
+            for (int i = 1; i < kc; ++i)
+            {
+                /* This does:
+                   state[i] = (16807 * state[i - 1]) % 2147483647;
+                   but avoids overflowing 31 bits. */
+                // Ichthyo 12/2018 : the above comment is only true for seed <= INT_MAX
+                //                   for INT_MAX < seed <= UINT_MAX the calculation diverges from correct
+                //                   modulus result, however, its values show a similar distribution pattern.
+                //                   Moreover the original code used long int for 'hi' and 'lo'.
+                //                   It behaves identical when using uint32_t, but not with int32_t
+                uint32_t hi = word / 127773;
+                uint32_t lo = word % 127773;
+                word = 16807 * lo - 2836 * hi;
+                if (word < 0)
+                    word += 2147483647;
+                *++dst = word;
+            }
+
+            fptr = &state[1];
+            rptr = &state[0];
+            kc *= 10;
+            while (--kc >= 0)
+                prngval();
+            return true;
+        }
 
 
+        uint32_t prngval()
+        {
+            uint32_t val = *fptr += uint32_t(*rptr);
+            uint32_t result = val >> 1;  // Chucking least random bit.
+                                         // Rationale: it has a less-then optimal repetition cycle.
+            int32_t *end = &state[62];
+            ++fptr;
+            if (fptr >= end)
+              {
+                fptr = state;
+                ++rptr;
+              }
+            else
+              {
+                ++rptr;
+                if (rptr >= end)
+                  rptr = state;
+              }
+            // random_result holds number 0...INT_MAX
+            return result;
+        }
 
 
+        float numRandom()
+        {
+            return prngval() / float(INT32_MAX);
+        }
 
+        // random number in the range 0...INT_MAX
+        uint32_t randomINT()
+        {
+            return prngval();
+        }
+};
+
+}
 
 int
 main (int, char**)
   {
-    Wau wau;
-    using ID = Miau &;
-    ID wuff = ID(wau);
+    StdlibPRNG oldGen;
+    TrinomialPRNG newGen;
+    
+    for (uint64_t seed=0; seed <= UINT32_MAX; ++seed)
+      {
+        oldGen.init(seed);
+        newGen.init(seed);
+        
+        for (uint i=0; i < 5*48000; ++i)
+          {
+            uint32_t oval = oldGen.prngval();
+            uint32_t nval = newGen.prngval();
+            if (oval != nval)
+              cout << "seed="<<seed << " i="<<i<< " \t "<<oval<< " != "<<nval<<endl;
+          }
+      }
     
-    cout << "Miau=" << wuff.u
-         << " ref to same object: " << isSameObject (wau, wuff)
-         << endl;
     
     cout <<  "\n.gulp.\n";