Added xorshift generator and longest runs test

rng/rng/.vscode/c_cpp_properties.json

@@ -0,0 +1,16 @@
+{
+    "configurations": [
+        {
+            "name": "Mac",
+            "includePath": [
+                "${workspaceFolder}/**"
+            ],
+            "defines": [],
+            "compilerPath": "/usr/bin/clang",
+            "cStandard": "c17",
+            "cppStandard": "c++20",
+            "intelliSenseMode": "macos-clang-arm64"
+        }
+    ],
+    "version": 4
+}

rng/rng/.vscode/settings.json

@@ -1,63 +0,0 @@
-{
-    "files.associations": {
-        "*.ts": "typescriptreact",
-        "array": "cpp",
-        "iostream": "cpp",
-        "__bit_reference": "cpp",
-        "__hash_table": "cpp",
-        "__locale": "cpp",
-        "__node_handle": "cpp",
-        "__split_buffer": "cpp",
-        "__threading_support": "cpp",
-        "__verbose_abort": "cpp",
-        "bitset": "cpp",
-        "cctype": "cpp",
-        "clocale": "cpp",
-        "cmath": "cpp",
-        "complex": "cpp",
-        "cstdarg": "cpp",
-        "cstddef": "cpp",
-        "cstdint": "cpp",
-        "cstdio": "cpp",
-        "cstdlib": "cpp",
-        "cstring": "cpp",
-        "ctime": "cpp",
-        "cwchar": "cpp",
-        "cwctype": "cpp",
-        "deque": "cpp",
-        "execution": "cpp",
-        "memory": "cpp",
-        "fstream": "cpp",
-        "initializer_list": "cpp",
-        "iomanip": "cpp",
-        "ios": "cpp",
-        "iosfwd": "cpp",
-        "istream": "cpp",
-        "limits": "cpp",
-        "locale": "cpp",
-        "mutex": "cpp",
-        "new": "cpp",
-        "optional": "cpp",
-        "ostream": "cpp",
-        "print": "cpp",
-        "queue": "cpp",
-        "ratio": "cpp",
-        "sstream": "cpp",
-        "stack": "cpp",
-        "stdexcept": "cpp",
-        "streambuf": "cpp",
-        "string": "cpp",
-        "string_view": "cpp",
-        "tuple": "cpp",
-        "typeinfo": "cpp",
-        "unordered_map": "cpp",
-        "variant": "cpp",
-        "vector": "cpp",
-        "algorithm": "cpp",
-        "__tree": "cpp",
-        "csignal": "cpp",
-        "map": "cpp",
-        "set": "cpp",
-        "unordered_set": "cpp"
-    }
-}

rng/rng/generators/LCG/lehmer.cpp

@@ -6,38 +6,26 @@
 //
 
 #include "../../rng.h"
+#include "../generator.h"
 
 // parameters recommended by Nakazawa & Nakazawa
 // https://en.wikipedia.org/wiki/Lehmer_random_number_generator
 
-uint64_t lehmer_minstd(uint64_t prev){
+namespace splat {
-    uint64_t a = 7759097958782935LL;
+    class lehmer_generator : public PRNG {
-    uint64_t m = 18055400005099021LL;
+        public:
-    
+            lehmer_generator(uint32_t genSeed) : PRNG(genSeed) {
-    uint64_t result = (a * prev) % m;
+                seed = genSeed;
-    
+            }
-    return result;
+            uint32_t generate() override {
-}
+                seed = (a * seed) % m;
-
+                return seed >> 32;
-class lehmer_generator {
+            }
-    uint64_t seed;
+            std::string getName() override {
-    public: 
+                return "lehmer";
-        uint32_t generate(){
+            }
-            seed = lehmer_minstd(seed);
+        private:
-            return seed;
+            uint64_t a = 7759097958782935LL;
-        }
+            uint64_t m = 18055400005099021LL;
-        lehmer_generator(uint64_t genSeed){
+    };
-            seed = genSeed;
+}
-        }
-};
-
-// std::vector<std::bitset<32>> lehmer_generate(int seed, int amount){
-//     std::vector<std::bitset<32>> bits(amount);
-    
-//     for(int i=0; i<amount; i++){
-//         seed = lehmer_minstd(seed);
-//         bits[i]=seed;
-//     }
-    
-//     return bits;
-// }

rng/rng/generators/generator.h

@@ -0,0 +1,18 @@
+#include "../rng.h"
+
+#pragma once
+
+namespace splat {
+  class PRNG {
+    public:
+      virtual uint32_t generate() {return -1;};
+      virtual ~PRNG() = default;
+      virtual std::string getName() {return "N/A";};
+      PRNG(uint32_t genseed){
+        seed = genseed;
+      }
+    protected:
+      uint32_t seed;
+  };
+}
+

rng/rng/generators/mt19937.cpp

@@ -1,4 +1,5 @@
 #include "../rng.h"
+#include "./generator.h"
 
 // more specifically this will be mt19937 - 32 bit
 
@@ -63,33 +64,33 @@ std::array<uint32_t,624> mt19937_temper(std::array<uint32_t,624> state){
     return tempered;
 }
 
-class mt19937_generator {
+namespace splat {
-    public: 
+    class mt19937_generator : public PRNG {
-        uint32_t seed;
+        public:
-        uint32_t generate() {
+            mt19937_generator(uint32_t genSeed) : PRNG(genSeed) {
-            uint32_t generated = random_values[position];
+                state = mt19937_init(seed);
-            position++;
-            if(position>=624){
                 nextblock();
             }
-            return generated;
+            uint32_t generate() override {
-        }
+                uint32_t generated = random_values[position];
-        double generate_01(){
+                position++;
-            return ((double)generate())/(4294967295);
+                if(position>=624){
-        }
+                    nextblock();
-        mt19937_generator(uint32_t genSeed){
+                }
-            seed = genSeed;
+                return generated;
-            state = mt19937_init(seed);
+            }
-            nextblock();
+            std::string getName() override {
-        }
+                return "mt19937-32";
-    private:
+            }
-        std::array<uint32_t,624> state;
+        private:
-        std::array<uint32_t,624> random_values;
+            std::array<uint32_t,624> state;
-        int position;
+            std::array<uint32_t,624> random_values;
-        // goes to next block of 624 values
+            int position;
-        void nextblock() {
+            // goes to next block of 624 values
-            state = mt19937_twist(state);
+            void nextblock() {
-            random_values = mt19937_temper(state);
+                state = mt19937_twist(state);
-            position = 0;
+                random_values = mt19937_temper(state);
-        }
+                position = 0;
-};
+            }
+    };
+}

rng/rng/generators/xorshift.cpp

@@ -0,0 +1,25 @@
+#include "../rng.h"
+#include "./generator.h"
+
+// parameters recommended by Nakazawa & Nakazawa
+// https://en.wikipedia.org/wiki/Lehmer_random_number_generator
+
+namespace splat {
+    class xorshift32 : public PRNG {
+        public:
+            xorshift32(uint32_t genSeed) : PRNG(genSeed) {
+                state = genSeed;
+            }
+            uint32_t generate() override {
+                state ^= state << 13;
+                state ^= state >> 17;
+                state ^= state << 5;
+                return state;
+            }
+            std::string getName() override {
+                return "xorshift32";
+            }
+        private:
+            uint32_t state;
+    };
+}

rng/rng/main.cpp

@@ -8,45 +8,78 @@
 #include "rng.h"
 #include "generators/LCG/lehmer.cpp"
 #include "generators/mt19937.cpp"
+#include "generators/xorshift.cpp"
 #include "randomness_tests/frequency_monobit.cpp"
 #include "randomness_tests/frequency_block.cpp"
 #include "randomness_tests/runs.cpp"
+#include "randomness_tests/runs_ones.cpp"
+
+
+
+#include "generators/generator.h"
+#include "randomness_tests/rngtest.h"
+
+namespace splat {
+  template <typename T>
+  void addGen(std::vector<std::unique_ptr<PRNG>> &generators, uint32_t &seed){
+    generators.push_back(std::make_unique<T>(seed));
+  }
+
+  std::vector<std::unique_ptr<PRNG>> getAllGenerators(uint32_t seed){
+    std::vector<std::unique_ptr<PRNG>> generators;
+
+    // add generators here
+    addGen<mt19937_generator>(generators, seed);
+    addGen<lehmer_generator>(generators, seed);
+    addGen<xorshift32>(generators,seed);
+    
+    return generators;
+  }
+
+  template <typename T>
+  void addTest(std::vector<std::unique_ptr<RNGTEST>> &tests, std::vector<std::bitset<32>> &data){
+    tests.push_back(std::make_unique<T>(data));
+  }
+
+  std::vector<std::unique_ptr<RNGTEST>> getAllTests(std::vector<std::bitset<32>> &data){
+    std::vector<std::unique_ptr<RNGTEST>> tests;
+
+    // add generators here
+    addTest<frequency_monobit_test>(tests, data);
+    addTest<frequency_block_test>(tests, data);
+    addTest<runs_test>(tests, data);
+    addTest<runs_ones_test>(tests, data);
+    
+    return tests;
+  }
+
+}
 
 int main(int argc, const char * argv[]) {
 
-    int seed = 8898;
+    int seed = 235211213;
     int blocksgenerated = 1000000;
 
-    mt19937_generator mt2 = mt19937_generator(seed);
+    std::vector<std::unique_ptr<splat::PRNG>> generators = splat::getAllGenerators(1238124);
-    lehmer_generator lg = lehmer_generator(seed);
+  
+    for(const std::unique_ptr<splat::PRNG> &generator : generators){
+        std::cout << generator->getName() << ": \n";
+        std::vector<std::bitset<32>> bits;
+        for(int i=0; i<blocksgenerated; i++){
+            bits.push_back(generator->generate());
+        }
 
-    std::vector<std::bitset<32>> mtbits;
+        std::vector<std::unique_ptr<splat::RNGTEST>> tests = splat::getAllTests(bits);
-    std::vector<std::bitset<32>> lgbits;
 
-    for(int i=0; i<blocksgenerated; i++){
+        for(const std::unique_ptr<splat::RNGTEST> &test : tests){
-        mtbits.push_back(mt2.generate());
+          std::string namestr = test->getName();
-        lgbits.push_back(lg.generate());
+          while(namestr.length() < 30){
+            namestr.append(" ");
+          }
+          std::cout << " - " << namestr << " = " << (test->passed()?"PASS":"FAIL") << " ("<< std::fixed << std::setprecision(8) << test->value() <<")\n";
+        }
     }
 
-    // for(auto bs : lgbits){
-    //     for(int i=0; i<32; i++){
-    //         std::cout << bs[i];
-    //     }
-    // }
-
-    int freqblocksize = (blocksgenerated * 32) / 120;
-
-    std::cout << "Linear Congruential Generator:\n";
-    std::cout << " - FreqMonobit       = "<< std::fixed << std::setprecision(10) << test_frequency_monobit(lgbits) << "\n";
-    std::cout << " - FreqBlock (n/120) = "<< std::fixed << std::setprecision(10) << test_frequency_block(lgbits, freqblocksize) << "\n";
-    std::cout << " - Runs              = "<< std::fixed << std::setprecision(10) << test_runs(lgbits) << "\n";
-    
-    std::cout << "Mersenne Twister mt19937-32:\n";
-    std::cout << " - FreqMonobit       = "<< std::fixed << std::setprecision(10) << test_frequency_monobit(mtbits) << "\n";
-    std::cout << " - FreqBlock (n/120) = "<< std::fixed << std::setprecision(10) << test_frequency_block(mtbits, freqblocksize) << "\n";
-    std::cout << " - Runs              = "<< std::fixed << std::setprecision(10) << test_runs(mtbits) << "\n";
-
-
 
     return 0;
-}
+};

rng/rng/plot_data.py

@@ -1,20 +0,0 @@
-import matplotlib.pyplot as plt
-
-x_coords = []
-y_coords = []
-
-with open("data.txt", "r") as file:
-    for line in file:
-        x_coords.append(float(line.split(",")[0]))
-        y_coords.append(float(line.split(",")[1]))
-
-plt.figure(figsize=(10, 8))
-plt.scatter(x_coords, y_coords, alpha=0.6, s=1)
-plt.grid(True, alpha=0.3)
-
-plt.axis('equal')
-plt.xlim(-0.05, 1.05)
-plt.ylim(-0.05, 1.05)
-
-plt.tight_layout()
-plt.show()

rng/rng/randomness_tests/frequency_block.cpp

@@ -8,44 +8,50 @@
 #include "../rng.h"
 #include "../math/incomplete_gamma.cpp"
 
-/**
+#include "./rngtest.h"
- * Returns the p value from the block frequency test
- */
-double test_frequency_block(std::vector<std::bitset<32>> data, u_long block_size){
-    
-    // std::cout << "BLOCK SIZE: "<< block_size<<"\n";
 
-    long long bitcount = data.size() * 32;
+namespace splat {
-    long long chunks = bitcount / block_size;
+    class frequency_block_test : public RNGTEST {
+        public: 
+            frequency_block_test(std::vector<std::bitset<32>> &testData) : RNGTEST(testData) {
+                block_size = data.size() * 32 / 150;
+                testPValue = runTest(data);
+                testPassed = testPValue > 0.01;
+            }
+            std::string getName() override {
+                return std::format("Frequency Block [{}]", block_size);
+            }
+            double runTest(std::vector<std::bitset<32>> &data) override {
 
-    double x2stat = 0;
+                // std::cout << "DEBUG: "<<data.size() << "\n";
-    
-    for(int chunkIndex = 0; chunkIndex < chunks; chunkIndex++){
-        // std::cout << "NEWCHUNK: ";
-        double onecount = 0;
-        for(int i=0; i<block_size; i++){
-            int dataIndex = (((chunkIndex * block_size)+i) / 32);
-            int bitIndex = (((chunkIndex * block_size)+i) % 32);
-            onecount += data[dataIndex][bitIndex];
-            // std::cout << data[dataIndex][bitIndex];
-        }
-        double oneproportion = onecount / (double)block_size;
 
-        x2stat += 4*block_size*std::pow((oneproportion - 0.5),2);
+                long long bitcount = data.size() * 32;
+                long long chunks = bitcount / block_size;
 
-        // std::cout << "\n";
+                double x2stat = 0;
-    }
+                
+                for(int chunkIndex = 0; chunkIndex < chunks; chunkIndex++){
+                    // std::cout << "NEWCHUNK: ";
+                    double onecount = 0;
+                    for(int i=0; i<block_size; i++){
+                        int dataIndex = (((chunkIndex * block_size)+i) / 32);
+                        int bitIndex = (((chunkIndex * block_size)+i) % 32);
+                        onecount += data[dataIndex][bitIndex];
+                        // std::cout << data[dataIndex][bitIndex];
+                    }
+                    double oneproportion = onecount / (double)block_size;
 
-    // std::cout << "[debug] calling igamc with "<<(chunks/2) <<','<<(x2stat/2) << "\n";
+                    x2stat += 4*block_size*std::pow((oneproportion - 0.5),2);
-    double p = igam(chunks/2, x2stat/2);
-    return p;
-}
 
-/**
+                    // std::cout << "\n";
- * It is recommended that each sequence to be tested consist of a minimum of 100 bits (i.e., n ≥ 100). Note
+                }
-that n ≥ MN. The block size M should be selected such that M ≥ 20, M > .01n and N < 100.
+
- */
+                // std::cout << "[debug] calling igamc with "<<(chunks/2) <<','<<(x2stat/2) << "\n";
-bool pass_frequency_monobit(std::vector<std::bitset<32>> data, int chunkSize){
+                double p = igam(chunks/2, x2stat/2);
-    double pv = test_frequency_block(data, chunkSize);
+                return p;
-    return pv >= 0.01;
+            }
+        private:
+            long long block_size;
+        
+    };
 }

rng/rng/randomness_tests/frequency_monobit.cpp

@@ -6,30 +6,28 @@
 //
 
 #include "../rng.h"
+#include "./rngtest.h"
 
-
+namespace splat {
-double test_frequency_monobit(std::vector<std::bitset<32>> data){
+    class frequency_monobit_test : public RNGTEST {
-    
+        public: 
-    long long s = 0;
+        frequency_monobit_test(std::vector<std::bitset<32>> &testData) : RNGTEST(testData) {
-    
+            testPValue = runTest(data);
-    for(auto bitset : data){
+            testPassed = testPValue > 0.01;
-        for(int i=0; i<32; i++){
-            s += (2*bitset[i])-1;
         }
-    }
+        std::string getName() override {
-    
+            return "Frequency Monobit";
-    // std::cout << "  [sval: "<<s<<"]  ";
+        }
-
+        double runTest(std::vector<std::bitset<32>> &data) override {
-    double sobs = ((double)std::abs(s))/std::sqrt(data.size()*32);
+            long long s = 0;
-
+            for(auto bitset : data){
-
+                for(int i=0; i<32; i++){
-
+                    s += (2*bitset[i])-1;
-    double pvalue = std::erfc(sobs/std::sqrt(2));
+                }
-
+            }
-    return pvalue;
+            double sobs = ((double)std::abs(s))/std::sqrt(data.size()*32);
-}
+            double pvalue = std::erfc(sobs/std::sqrt(2));
-
+            return pvalue;
-bool pass_frequency_monobit(std::vector<std::bitset<32>> data){
+        }
-    double pv = test_frequency_monobit(data);
+    };
-    return pv >= 0.01;
 }

rng/rng/randomness_tests/rngtest.h

@@ -0,0 +1,25 @@
+#include "../rng.h"
+
+#pragma once
+
+namespace splat {
+  class RNGTEST {
+    public:
+      bool passed(){
+        return testPassed;
+      }
+      double value(){
+        return testPValue;
+      }
+      virtual ~RNGTEST() = default;
+      virtual std::string getName() {return "N/A";};
+      RNGTEST(std::vector<std::bitset<32>> &testData) : data(testData){
+      };
+    protected:
+      std::vector<std::bitset<32>> &data;
+      bool testPassed;
+      double testPValue;
+      virtual double runTest(std::vector<std::bitset<32>> &data) {return -1.0;};
+  };
+}
+

rng/rng/randomness_tests/runs.cpp

@@ -1,26 +1,38 @@
 #include "../rng.h"
+#include "./rngtest.h"
 
-double test_runs(std::vector<std::bitset<32>> data){
+namespace splat {
-
+    class runs_test : public RNGTEST {
-    long long totalSize = data.size() * 32;
+        public: 
-
+        runs_test(std::vector<std::bitset<32>> &testData) : RNGTEST(testData) {
-    double vnobs = 1;
+            testPValue = runTest(data);
-
+            testPassed = testPValue > 0.01;
-    double onesproportion = 0;
-
-    for(int i=0; i<totalSize-1; i++){
-
-        if(data[i/32][i%32]!=data[(i+1)/32][(i+1)%32]){
-            vnobs++;
         }
-        onesproportion+=data[i/32][i%32];
+        std::string getName() override {
-    }
+            return "Runs";
+        }
+        double runTest(std::vector<std::bitset<32>> &data) override {
+            long long totalSize = data.size() * 32;
 
-    onesproportion/=totalSize;
+            double vnobs = 1;
 
-    return std::erfc(
+            double onesproportion = 0;
-        std::abs(vnobs-(2*totalSize*onesproportion*(1-onesproportion)))
+
-        /
+            for(int i=0; i<totalSize-1; i++){
-        (2*std::sqrt(2*totalSize)*onesproportion*(1-onesproportion))
+
-    );
+                if(data[i/32][i%32]!=data[(i+1)/32][(i+1)%32]){
+                    vnobs++;
+                }
+                onesproportion+=data[i/32][i%32];
+            }
+
+            onesproportion/=totalSize;
+
+            return std::erfc(
+                std::abs(vnobs-(2*totalSize*onesproportion*(1-onesproportion)))
+                /
+                (2*std::sqrt(2*totalSize)*onesproportion*(1-onesproportion))
+            );
+        }
+    };
 }

rng/rng/randomness_tests/runs_ones.cpp

@@ -0,0 +1,78 @@
+#include "../rng.h"
+#include "./rngtest.h"
+#include "../math/incomplete_gamma.cpp"
+#define vecdebug(list) for(auto vecdebugitem : list){ std::cout << vecdebugitem <<","; }std::cout << "\n";
+
+
+
+namespace splat {
+    class runs_ones_test : public RNGTEST {
+        public: 
+            runs_ones_test(std::vector<std::bitset<32>> &testData) : RNGTEST(testData) {
+                testPValue = runTest(data);
+                testPassed = testPValue > 0.01;
+            }
+            std::string getName() override {
+                return "Runs of 1s";
+            }
+            double runTest(std::vector<std::bitset<32>> &data) override {
+                long long totalSize = data.size() * 32;
+
+                if(totalSize < 750000){
+                    std::cout << "Error, runs of ones test requires more than 750K bits of input\n";
+                    return 0;
+                }
+                // we do it with M=10^4 here
+
+                std::vector<double> counts(7,0);
+
+                double chunks = totalSize / 10000;
+
+                for(int chunkIndex=0; chunkIndex<chunks; chunkIndex++){
+
+                    int longestRun = 0;
+
+                    int ones = 0;
+
+                    for(int i=0; i<10000; i++){
+                        long long index = chunkIndex*10000 + i;
+                        int dataIndex = (index / 32);
+                        int bitIndex = (index % 32);
+
+                        if(data[dataIndex][bitIndex]){
+                            ones++;
+                        }else{
+                            if(ones>longestRun) longestRun = ones;
+                            ones = 0;
+                        }
+                    }
+
+                    if(longestRun <=10){
+                        counts[0]++;
+                    }else if(longestRun>=16){
+                        counts[6]++;
+                    }else{
+                        counts[longestRun-10]++;
+                    }
+                }
+                
+                double x2 = 0;
+
+                std::vector<double> probabilities = {
+                    0.0882,
+                    0.2092,
+                    0.2483,
+                    0.1933,
+                    0.1208,
+                    0.0675,
+                    0.0727
+                };
+
+                for(int i=0; i<7; i++){
+                    x2 += std::pow((counts[i] - (chunks * probabilities[i])),2.0) /(chunks * probabilities[i]);
+
+                }
+                return igamc(3, x2/2.0);
+            }
+    };
+}

rng/rng/readme.md

@@ -0,0 +1,2 @@
+🫟 Splatter
+c++ rng library with generators and tests

rng/rng/rng.h

@@ -14,3 +14,4 @@
 #include <vector>
 #include <cmath>
 #include <iomanip>
+#include <format>

rng/rng/run.sh

@@ -1,3 +1,2 @@
-g++ -std=c++11 main.cpp -o main.o
+g++ -std=c++20 main.cpp -o main.o
-./main.o
+./main.o
-rm main.o