diff --git a/apps/devApps/random explorer/addons.make b/apps/devApps/random explorer/addons.make
new file mode 100644
index 00000000000..f47512dc0e8
--- /dev/null
+++ b/apps/devApps/random explorer/addons.make	
@@ -0,0 +1 @@
+ofxGui
diff --git a/apps/devApps/random explorer/src/Dist.hpp b/apps/devApps/random explorer/src/Dist.hpp
new file mode 100644
index 00000000000..86d1b15468d
--- /dev/null
+++ b/apps/devApps/random explorer/src/Dist.hpp	
@@ -0,0 +1,202 @@
+#ifndef DIST_HPP
+#define DIST_HPP
+
+#include "ofMain.h"
+
+
+struct Dist {
+    
+    inline static string base_url_{ "https://en.cppreference.com/w/cpp/numeric/random/" };
+    string url_;
+	ofParameterGroup parameters_;
+    ofColor color_ {128,128,128};
+    std::string info_;
+    glm::vec2 range_ { };
+    bool discrete_ { false };
+    size_t underflow_;
+    size_t overflow_;
+    std::size_t max_ { 0 };
+    float cost_ { 0.0f };
+    
+    virtual auto gen() -> void = 0;
+    virtual auto clear() -> void = 0;
+    virtual auto compile() -> void = 0;
+    virtual auto draw(float x, float y, float w, float h) -> void = 0;
+    virtual ~Dist() = default;
+
+    Dist() {};
+};
+
+template <typename T>
+struct ConcreteDist: public Dist {
+    
+    std::vector<T> data_;
+    std::function<T()> gen_;
+    std::vector<T> bins_;
+    int autorot_ { 0 };
+	
+	ofParameter<void> url_button_;
+
+    ConcreteDist(std::string label,
+                 std::string info,
+                 std::string url,
+                 std::vector<ofAbstractParameter *> params,
+                 std::function<T()> gen,
+                 std::size_t num_bins = 101,
+                 glm::vec2 range = {0, 100},
+                 bool discrete = false )
+    : gen_(gen)
+    {
+        url_ = base_url_+url;
+        info_ = info;
+        range_ = range;
+        discrete_ = discrete;
+        bins_.resize(num_bins);
+        parameters_.setName(label);
+        for (auto & p: params) parameters_.add(*p);
+		url_button_.set("click for reference");
+		url_button_.addListener(this, &ConcreteDist::open_url);
+		parameters_.add(url_button_);
+		
+		//panel_.getControl("seed")->setTextColor(ofColor::green);
+
+
+    }
+    
+	auto open_url() {
+		ofLaunchBrowser(url_);
+	}
+	
+    auto gen() -> void override {
+        data_.push_back(gen_());
+    }
+    
+    auto clear() -> void override {
+        overflow_ = 0;
+        underflow_ = 0;
+        data_.clear();
+        std::fill(bins_.begin(), bins_.end(), T{0});
+    }
+    
+    auto compile() -> void override {
+        // histograms non-vecs only
+        if constexpr (std::is_arithmetic_v<T>) {
+            float divisor = (range_.y-range_.x)/float(bins_.size()-1);
+            for (auto & v: data_) {
+                v /= divisor;
+                if (v < range_.x) {
+                    underflow_++;
+                } else if (v >= bins_.size()) {
+                    overflow_++;
+                } else {
+                    bins_.at(v) = bins_.at(v)+1;
+                }
+            }
+
+			max_ = 0.0;
+			for (size_t i=0;i<bins_.size(); i++) {
+				if (bins_.at(i) > max_) max_ = bins_.at(i);
+			}
+		} else {
+			// no histograms for vecN
+		}
+    }
+    
+    auto draw(float x, float y, float w, float h) -> void override {
+        ofPushStyle();
+        ofPushMatrix();
+        {
+            ofTranslate(x,y);
+            ofSetColor(color_);
+            ofDrawRectangle(0,0,w,h);
+
+            ofSetColor(ofColor::darkRed);
+            if (underflow_) ofDrawBitmapString("undershoot: " + ofToString(underflow_), w+5, 58);
+            if (overflow_) ofDrawBitmapString("overshoot: " + ofToString(overflow_), w+5, 74);
+
+            ofSetColor(192,192,192,255);
+            ofDrawBitmapString(info_, w+5, 35);
+            ofSetColor(255,255,255,255);
+            ofDrawBitmapStringHighlight(parameters_.getName() + " " + ofToString(cost_*1000, 2, 5)+"ms", w+5, 12);
+    
+            if constexpr (std::is_arithmetic_v<T>) {
+
+                auto p = 0.0f;
+                double incr = w/bins_.size();
+                auto fact = h/max_;
+                if (discrete_) {
+                    
+                    // line bars for discrete
+                    ofTranslate(incr/2, 0);
+                    for (auto y: bins_) {
+                        ofDrawLine(0, h, 0, h-float(y)*fact);
+                        if (y==0) {
+                            ofNoFill();
+                            ofDrawCircle(0, h-float(y)*fact, 2.5);
+                            ofFill();
+                        } else {
+                            ofDrawCircle(0, h-float(y)*fact, 3);
+                        }
+                        ofTranslate(int(incr), 0);
+                    }
+                } else {
+
+                    // integral for reals
+                    ofPolyline line;
+                    line.addVertex(0,h-bins_[0]*fact);
+                    for (auto y: bins_) line.lineTo(p+=incr, h-float(y)*fact);
+                    line.draw();
+                }
+                
+            } else if constexpr (std::is_same_v<T, glm::vec2>) {
+                
+                ofSetColor(255,255,255,96);
+                for (const auto & d: data_) ofDrawCircle(d, .5);
+                
+            } else if constexpr (std::is_same_v<T, glm::vec3>) {
+                
+                ofSetColor(255,255,255,32);
+                of3dPrimitive prim;
+                prim.getMesh().getVertices() = data_;
+                prim.getMesh().setMode(OF_PRIMITIVE_POINTS);
+                prim.rotateDeg(70,{ 0.2, 0.3, 0.5 }); // just some perspective
+                
+                ofPushMatrix();
+                {
+                    ofTranslate(w * 0.2, h * 0.2);
+                    prim.drawWireframe();
+                    prim.drawAxes(w * 0.5);
+                }
+                ofPopMatrix();
+			} else {
+				ofDrawBitmapString("unsupported visualisation", 10,10);
+			}
+        }
+        ofPopMatrix();
+        ofPopStyle();
+    }
+};
+
+struct DistGroup {
+    std::vector<std::shared_ptr<Dist>> dists_;
+    ofxPanel panel_;
+    DistGroup(std::vector<std::shared_ptr<Dist>> dists): dists_(dists) { }
+    
+    auto draw(std::string label, int square, int gap) {
+        panel_.draw();
+        ofPushMatrix();
+        {
+            ofTranslate(panel_.getPosition());
+            ofDrawBitmapString(label,0,-10);
+            ofTranslate(panel_.getWidth()+20,0);
+            for (const auto & dist: dists_) {
+                dist->draw(0,0,square,square);
+                ofTranslate(0,square + gap);
+            }
+        }
+        ofPopMatrix();
+    }
+};
+
+
+#endif /* DIST_HPP */
diff --git a/apps/devApps/random explorer/src/OKColor.h b/apps/devApps/random explorer/src/OKColor.h
new file mode 100644
index 00000000000..93345bc0f7c
--- /dev/null
+++ b/apps/devApps/random explorer/src/OKColor.h	
@@ -0,0 +1,690 @@
+#ifndef OK_COLOR_H
+#define OK_COLOR_H
+
+// Copyright(c) 2021 Björn Ottosson
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy of
+// this softwareand associated documentation files(the "Software"), to deal in
+// the Software without restriction, including without limitation the rights to
+// use, copy, modify, merge, publish, distribute, sublicense, and /or sell copies
+// of the Software, and to permit persons to whom the Software is furnished to do
+// so, subject to the following conditions :
+// The above copyright noticeand this permission notice shall be included in all
+// copies or substantial portions of the Software.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+#include <cmath>
+#include <cfloat>
+
+class OKColor
+{
+
+struct Lab { float L; float a; float b; };
+struct RGB { float r; float g; float b; };
+struct HSV { float h; float s; float v; };
+struct HSL { float h; float s; float l; };
+struct LC { float L; float C; };
+
+// Alternative representation of (L_cusp, C_cusp)
+// Encoded so S = C_cusp/L_cusp and T = C_cusp/(1-L_cusp)
+// The maximum value for C in the triangle is then found as fmin(S*L, T*(1-L)), for a given L
+struct ST { float S; float T; };
+
+static constexpr float pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062f;
+
+static float clamp(float x, float min, float max)
+{
+    if (x < min)
+        return min;
+    if (x > max)
+        return max;
+
+    return x;
+}
+
+static float sgn(float x)
+{
+    return (float)(0.f < x) - (float)(x < 0.f);
+}
+
+static float srgb_transfer_function(float a)
+{
+    return .0031308f >= a ? 12.92f * a : 1.055f * powf(a, .4166666666666667f) - .055f;
+}
+
+static float srgb_transfer_function_inv(float a)
+{
+    return .04045f < a ? powf((a + .055f) / 1.055f, 2.4f) : a / 12.92f;
+}
+
+static Lab linear_srgb_to_oklab(RGB c)
+{
+    float l = 0.4122214708f * c.r + 0.5363325363f * c.g + 0.0514459929f * c.b;
+    float m = 0.2119034982f * c.r + 0.6806995451f * c.g + 0.1073969566f * c.b;
+    float s = 0.0883024619f * c.r + 0.2817188376f * c.g + 0.6299787005f * c.b;
+
+    float l_ = cbrtf(l);
+    float m_ = cbrtf(m);
+    float s_ = cbrtf(s);
+
+    return {
+        0.2104542553f * l_ + 0.7936177850f * m_ - 0.0040720468f * s_,
+        1.9779984951f * l_ - 2.4285922050f * m_ + 0.4505937099f * s_,
+        0.0259040371f * l_ + 0.7827717662f * m_ - 0.8086757660f * s_,
+    };
+}
+
+static RGB oklab_to_linear_srgb(Lab c)
+{
+    float l_ = c.L + 0.3963377774f * c.a + 0.2158037573f * c.b;
+    float m_ = c.L - 0.1055613458f * c.a - 0.0638541728f * c.b;
+    float s_ = c.L - 0.0894841775f * c.a - 1.2914855480f * c.b;
+
+    float l = l_ * l_ * l_;
+    float m = m_ * m_ * m_;
+    float s = s_ * s_ * s_;
+
+    return {
+        +4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s,
+        -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s,
+        -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s,
+    };
+}
+
+// Finds the maximum saturation possible for a given hue that fits in sRGB
+// Saturation here is defined as S = C/L
+// a and b must be normalized so a^2 + b^2 == 1
+static float compute_max_saturation(float a, float b)
+{
+    // Max saturation will be when one of r, g or b goes below zero.
+
+    // Select different coefficients depending on which component goes below zero first
+    float k0, k1, k2, k3, k4, wl, wm, ws;
+
+    if (-1.88170328f * a - 0.80936493f * b > 1)
+    {
+        // Red component
+        k0 = +1.19086277f; k1 = +1.76576728f; k2 = +0.59662641f; k3 = +0.75515197f; k4 = +0.56771245f;
+        wl = +4.0767416621f; wm = -3.3077115913f; ws = +0.2309699292f;
+    }
+    else if (1.81444104f * a - 1.19445276f * b > 1)
+    {
+        // Green component
+        k0 = +0.73956515f; k1 = -0.45954404f; k2 = +0.08285427f; k3 = +0.12541070f; k4 = +0.14503204f;
+        wl = -1.2684380046f; wm = +2.6097574011f; ws = -0.3413193965f;
+    }
+    else
+    {
+        // Blue component
+        k0 = +1.35733652f; k1 = -0.00915799f; k2 = -1.15130210f; k3 = -0.50559606f; k4 = +0.00692167f;
+        wl = -0.0041960863f; wm = -0.7034186147f; ws = +1.7076147010f;
+    }
+
+    // Approximate max saturation using a polynomial:
+    float S = k0 + k1 * a + k2 * b + k3 * a * a + k4 * a * b;
+
+    // Do one step Halley's method to get closer
+    // this gives an error less than 10e6, except for some blue hues where the dS/dh is close to infinite
+    // this should be sufficient for most applications, otherwise do two/three steps
+
+    float k_l = +0.3963377774f * a + 0.2158037573f * b;
+    float k_m = -0.1055613458f * a - 0.0638541728f * b;
+    float k_s = -0.0894841775f * a - 1.2914855480f * b;
+
+    {
+        float l_ = 1.f + S * k_l;
+        float m_ = 1.f + S * k_m;
+        float s_ = 1.f + S * k_s;
+
+        float l = l_ * l_ * l_;
+        float m = m_ * m_ * m_;
+        float s = s_ * s_ * s_;
+
+        float l_dS = 3.f * k_l * l_ * l_;
+        float m_dS = 3.f * k_m * m_ * m_;
+        float s_dS = 3.f * k_s * s_ * s_;
+
+        float l_dS2 = 6.f * k_l * k_l * l_;
+        float m_dS2 = 6.f * k_m * k_m * m_;
+        float s_dS2 = 6.f * k_s * k_s * s_;
+
+        float f = wl * l + wm * m + ws * s;
+        float f1 = wl * l_dS + wm * m_dS + ws * s_dS;
+        float f2 = wl * l_dS2 + wm * m_dS2 + ws * s_dS2;
+
+        S = S - f * f1 / (f1 * f1 - 0.5f * f * f2);
+    }
+
+    return S;
+}
+
+// finds L_cusp and C_cusp for a given hue
+// a and b must be normalized so a^2 + b^2 == 1
+static LC find_cusp(float a, float b)
+{
+    // First, find the maximum saturation (saturation S = C/L)
+    float S_cusp = compute_max_saturation(a, b);
+
+    // Convert to linear sRGB to find the first point where at least one of r,g or b >= 1:
+    RGB rgb_at_max = oklab_to_linear_srgb({ 1, S_cusp * a, S_cusp * b });
+    float L_cusp = cbrtf(1.f / fmax(fmax(rgb_at_max.r, rgb_at_max.g), rgb_at_max.b));
+    float C_cusp = L_cusp * S_cusp;
+
+    return { L_cusp , C_cusp };
+}
+
+// Finds intersection of the line defined by
+// L = L0 * (1 - t) + t * L1;
+// C = t * C1;
+// a and b must be normalized so a^2 + b^2 == 1
+static float find_gamut_intersection(float a, float b, float L1, float C1, float L0, LC cusp)
+{
+    // Find the intersection for upper and lower half seprately
+    float t;
+    if (((L1 - L0) * cusp.C - (cusp.L - L0) * C1) <= 0.f)
+    {
+        // Lower half
+
+        t = cusp.C * L0 / (C1 * cusp.L + cusp.C * (L0 - L1));
+    }
+    else
+    {
+        // Upper half
+
+        // First intersect with triangle
+        t = cusp.C * (L0 - 1.f) / (C1 * (cusp.L - 1.f) + cusp.C * (L0 - L1));
+
+        // Then one step Halley's method
+        {
+            float dL = L1 - L0;
+            float dC = C1;
+
+            float k_l = +0.3963377774f * a + 0.2158037573f * b;
+            float k_m = -0.1055613458f * a - 0.0638541728f * b;
+            float k_s = -0.0894841775f * a - 1.2914855480f * b;
+
+            float l_dt = dL + dC * k_l;
+            float m_dt = dL + dC * k_m;
+            float s_dt = dL + dC * k_s;
+
+
+            // If higher accuracy is required, 2 or 3 iterations of the following block can be used:
+            {
+                float L = L0 * (1.f - t) + t * L1;
+                float C = t * C1;
+
+                float l_ = L + C * k_l;
+                float m_ = L + C * k_m;
+                float s_ = L + C * k_s;
+
+                float l = l_ * l_ * l_;
+                float m = m_ * m_ * m_;
+                float s = s_ * s_ * s_;
+
+                float ldt = 3 * l_dt * l_ * l_;
+                float mdt = 3 * m_dt * m_ * m_;
+                float sdt = 3 * s_dt * s_ * s_;
+
+                float ldt2 = 6 * l_dt * l_dt * l_;
+                float mdt2 = 6 * m_dt * m_dt * m_;
+                float sdt2 = 6 * s_dt * s_dt * s_;
+
+                float r = 4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s - 1;
+                float r1 = 4.0767416621f * ldt - 3.3077115913f * mdt + 0.2309699292f * sdt;
+                float r2 = 4.0767416621f * ldt2 - 3.3077115913f * mdt2 + 0.2309699292f * sdt2;
+
+                float u_r = r1 / (r1 * r1 - 0.5f * r * r2);
+                float t_r = -r * u_r;
+
+                float g = -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s - 1;
+                float g1 = -1.2684380046f * ldt + 2.6097574011f * mdt - 0.3413193965f * sdt;
+                float g2 = -1.2684380046f * ldt2 + 2.6097574011f * mdt2 - 0.3413193965f * sdt2;
+
+                float u_g = g1 / (g1 * g1 - 0.5f * g * g2);
+                float t_g = -g * u_g;
+
+                float b = -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s - 1;
+                float b1 = -0.0041960863f * ldt - 0.7034186147f * mdt + 1.7076147010f * sdt;
+                float b2 = -0.0041960863f * ldt2 - 0.7034186147f * mdt2 + 1.7076147010f * sdt2;
+
+                float u_b = b1 / (b1 * b1 - 0.5f * b * b2);
+                float t_b = -b * u_b;
+
+                t_r = u_r >= 0.f ? t_r : FLT_MAX;
+                t_g = u_g >= 0.f ? t_g : FLT_MAX;
+                t_b = u_b >= 0.f ? t_b : FLT_MAX;
+
+                t += fmin(t_r, fmin(t_g, t_b));
+            }
+        }
+    }
+
+    return t;
+}
+
+static float find_gamut_intersection(float a, float b, float L1, float C1, float L0)
+{
+    // Find the cusp of the gamut triangle
+    LC cusp = find_cusp(a, b);
+
+    return find_gamut_intersection(a, b, L1, C1, L0, cusp);
+}
+
+static RGB gamut_clip_preserve_chroma(RGB rgb)
+{
+    if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
+        return rgb;
+
+    Lab lab = linear_srgb_to_oklab(rgb);
+
+    float L = lab.L;
+    float eps = 0.00001f;
+    float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
+    float a_ = lab.a / C;
+    float b_ = lab.b / C;
+
+    float L0 = clamp(L, 0, 1);
+
+    float t = find_gamut_intersection(a_, b_, L, C, L0);
+    float L_clipped = L0 * (1 - t) + t * L;
+    float C_clipped = t * C;
+
+    return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
+}
+
+static RGB gamut_clip_project_to_0_5(RGB rgb)
+{
+    if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
+        return rgb;
+
+    Lab lab = linear_srgb_to_oklab(rgb);
+
+    float L = lab.L;
+    float eps = 0.00001f;
+    float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
+    float a_ = lab.a / C;
+    float b_ = lab.b / C;
+
+    float L0 = 0.5;
+
+    float t = find_gamut_intersection(a_, b_, L, C, L0);
+    float L_clipped = L0 * (1 - t) + t * L;
+    float C_clipped = t * C;
+
+    return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
+}
+
+static RGB gamut_clip_project_to_L_cusp(RGB rgb)
+{
+    if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
+        return rgb;
+
+    Lab lab = linear_srgb_to_oklab(rgb);
+
+    float L = lab.L;
+    float eps = 0.00001f;
+    float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
+    float a_ = lab.a / C;
+    float b_ = lab.b / C;
+
+    // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once.
+    LC cusp = find_cusp(a_, b_);
+
+    float L0 = cusp.L;
+
+    float t = find_gamut_intersection(a_, b_, L, C, L0);
+
+    float L_clipped = L0 * (1 - t) + t * L;
+    float C_clipped = t * C;
+
+    return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
+}
+
+static RGB gamut_clip_adaptive_L0_0_5(RGB rgb, float alpha = 0.05f)
+{
+    if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
+        return rgb;
+
+    Lab lab = linear_srgb_to_oklab(rgb);
+
+    float L = lab.L;
+    float eps = 0.00001f;
+    float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
+    float a_ = lab.a / C;
+    float b_ = lab.b / C;
+
+    float Ld = L - 0.5f;
+    float e1 = 0.5f + fabs(Ld) + alpha * C;
+    float L0 = 0.5f * (1.f + sgn(Ld) * (e1 - sqrtf(e1 * e1 - 2.f * fabs(Ld))));
+
+    float t = find_gamut_intersection(a_, b_, L, C, L0);
+    float L_clipped = L0 * (1.f - t) + t * L;
+    float C_clipped = t * C;
+
+    return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
+}
+
+static RGB gamut_clip_adaptive_L0_L_cusp(RGB rgb, float alpha = 0.05f)
+{
+    if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
+        return rgb;
+
+    Lab lab = linear_srgb_to_oklab(rgb);
+
+    float L = lab.L;
+    float eps = 0.00001f;
+    float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
+    float a_ = lab.a / C;
+    float b_ = lab.b / C;
+
+    // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once.
+    LC cusp = find_cusp(a_, b_);
+
+    float Ld = L - cusp.L;
+    float k = 2.f * (Ld > 0 ? 1.f - cusp.L : cusp.L);
+
+    float e1 = 0.5f * k + fabs(Ld) + alpha * C / k;
+    float L0 = cusp.L + 0.5f * (sgn(Ld) * (e1 - sqrtf(e1 * e1 - 2.f * k * fabs(Ld))));
+
+    float t = find_gamut_intersection(a_, b_, L, C, L0);
+    float L_clipped = L0 * (1.f - t) + t * L;
+    float C_clipped = t * C;
+
+    return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
+}
+
+static float toe(float x)
+{
+    constexpr float k_1 = 0.206f;
+    constexpr float k_2 = 0.03f;
+    constexpr float k_3 = (1.f + k_1) / (1.f + k_2);
+    return 0.5f * (k_3 * x - k_1 + sqrtf((k_3 * x - k_1) * (k_3 * x - k_1) + 4 * k_2 * k_3 * x));
+}
+
+static float toe_inv(float x)
+{
+    constexpr float k_1 = 0.206f;
+    constexpr float k_2 = 0.03f;
+    constexpr float k_3 = (1.f + k_1) / (1.f + k_2);
+    return (x * x + k_1 * x) / (k_3 * (x + k_2));
+}
+
+static ST to_ST(LC cusp)
+{
+    float L = cusp.L;
+    float C = cusp.C;
+    return { C / L, C / (1 - L) };
+}
+
+// Returns a smooth approximation of the location of the cusp
+// This polynomial was created by an optimization process
+// It has been designed so that S_mid < S_max and T_mid < T_max
+static ST get_ST_mid(float a_, float b_)
+{
+    float S = 0.11516993f + 1.f / (
+        +7.44778970f + 4.15901240f * b_
+        + a_ * (-2.19557347f + 1.75198401f * b_
+            + a_ * (-2.13704948f - 10.02301043f * b_
+                + a_ * (-4.24894561f + 5.38770819f * b_ + 4.69891013f * a_
+                    )))
+        );
+
+    float T = 0.11239642f + 1.f / (
+        +1.61320320f - 0.68124379f * b_
+        + a_ * (+0.40370612f + 0.90148123f * b_
+            + a_ * (-0.27087943f + 0.61223990f * b_
+                + a_ * (+0.00299215f - 0.45399568f * b_ - 0.14661872f * a_
+                    )))
+        );
+
+    return { S, T };
+}
+
+struct Cs { float C_0; float C_mid; float C_max; };
+static Cs get_Cs(float L, float a_, float b_)
+{
+    LC cusp = find_cusp(a_, b_);
+
+    float C_max = find_gamut_intersection(a_, b_, L, 1, L, cusp);
+    ST ST_max = to_ST(cusp);
+
+    // Scale factor to compensate for the curved part of gamut shape:
+    float k = C_max / fmin((L * ST_max.S), (1 - L) * ST_max.T);
+
+    float C_mid;
+    {
+        ST ST_mid = get_ST_mid(a_, b_);
+
+        // Use a soft minimum function, instead of a sharp triangle shape to get a smooth value for chroma.
+        float C_a = L * ST_mid.S;
+        float C_b = (1.f - L) * ST_mid.T;
+        C_mid = 0.9f * k * sqrtf(sqrtf(1.f / (1.f / (C_a * C_a * C_a * C_a) + 1.f / (C_b * C_b * C_b * C_b))));
+    }
+
+    float C_0;
+    {
+        // for C_0, the shape is independent of hue, so ST are constant. Values picked to roughly be the average values of ST.
+        float C_a = L * 0.4f;
+        float C_b = (1.f - L) * 0.8f;
+
+        // Use a soft minimum function, instead of a sharp triangle shape to get a smooth value for chroma.
+        C_0 = sqrtf(1.f / (1.f / (C_a * C_a) + 1.f / (C_b * C_b)));
+    }
+
+    return { C_0, C_mid, C_max };
+}
+
+public:
+static RGB okhsl_to_srgb(HSL hsl)
+{
+    float h = hsl.h;
+    float s = hsl.s;
+    float l = hsl.l;
+
+    if (l == 1.0f)
+    {
+        return { 1.f, 1.f, 1.f };
+    }
+
+    else if (l == 0.f)
+    {
+        return { 0.f, 0.f, 0.f };
+    }
+
+    float a_ = cosf(2.f * pi * h);
+    float b_ = sinf(2.f * pi * h);
+    float L = toe_inv(l);
+
+    Cs cs = get_Cs(L, a_, b_);
+    float C_0 = cs.C_0;
+    float C_mid = cs.C_mid;
+    float C_max = cs.C_max;
+
+    float mid = 0.8f;
+    float mid_inv = 1.25f;
+
+    float C, t, k_0, k_1, k_2;
+
+    if (s < mid)
+    {
+        t = mid_inv * s;
+
+        k_1 = mid * C_0;
+        k_2 = (1.f - k_1 / C_mid);
+
+        C = t * k_1 / (1.f - k_2 * t);
+    }
+    else
+    {
+        t = (s - mid)/ (1 - mid);
+
+        k_0 = C_mid;
+        k_1 = (1.f - mid) * C_mid * C_mid * mid_inv * mid_inv / C_0;
+        k_2 = (1.f - (k_1) / (C_max - C_mid));
+
+        C = k_0 + t * k_1 / (1.f - k_2 * t);
+    }
+
+    RGB rgb = oklab_to_linear_srgb({ L, C * a_, C * b_ });
+    return {
+        srgb_transfer_function(rgb.r),
+        srgb_transfer_function(rgb.g),
+        srgb_transfer_function(rgb.b),
+    };
+}
+
+static HSL srgb_to_okhsl(RGB rgb)
+{
+    Lab lab = linear_srgb_to_oklab({
+        srgb_transfer_function_inv(rgb.r),
+        srgb_transfer_function_inv(rgb.g),
+        srgb_transfer_function_inv(rgb.b)
+        });
+
+    float C = sqrtf(lab.a * lab.a + lab.b * lab.b);
+    float a_ = lab.a / C;
+    float b_ = lab.b / C;
+
+    float L = lab.L;
+    float h = 0.5f + 0.5f * atan2f(-lab.b, -lab.a) / pi;
+
+    Cs cs = get_Cs(L, a_, b_);
+    float C_0 = cs.C_0;
+    float C_mid = cs.C_mid;
+    float C_max = cs.C_max;
+
+    // Inverse of the interpolation in okhsl_to_srgb:
+
+    float mid = 0.8f;
+    float mid_inv = 1.25f;
+
+    float s;
+    if (C < C_mid)
+    {
+        float k_1 = mid * C_0;
+        float k_2 = (1.f - k_1 / C_mid);
+
+        float t = C / (k_1 + k_2 * C);
+        s = t * mid;
+    }
+    else
+    {
+        float k_0 = C_mid;
+        float k_1 = (1.f - mid) * C_mid * C_mid * mid_inv * mid_inv / C_0;
+        float k_2 = (1.f - (k_1) / (C_max - C_mid));
+
+        float t = (C - k_0) / (k_1 + k_2 * (C - k_0));
+        s = mid + (1.f - mid) * t;
+    }
+
+    float l = toe(L);
+    return { h, s, l };
+}
+
+
+static RGB okhsv_to_srgb(HSV hsv)
+{
+    float h = hsv.h;
+    float s = hsv.s;
+    float v = hsv.v;
+
+    float a_ = cosf(2.f * pi * h);
+    float b_ = sinf(2.f * pi * h);
+
+    LC cusp = find_cusp(a_, b_);
+    ST ST_max = to_ST(cusp);
+    float S_max = ST_max.S;
+    float T_max = ST_max.T;
+    float S_0 = 0.5f;
+    float k = 1 - S_0 / S_max;
+
+    // first we compute L and V as if the gamut is a perfect triangle:
+
+    // L, C when v==1:
+    float L_v = 1     - s * S_0 / (S_0 + T_max - T_max * k * s);
+    float C_v = s * T_max * S_0 / (S_0 + T_max - T_max * k * s);
+
+    float L = v * L_v;
+    float C = v * C_v;
+
+    // then we compensate for both toe and the curved top part of the triangle:
+    float L_vt = toe_inv(L_v);
+    float C_vt = C_v * L_vt / L_v;
+
+    float L_new = toe_inv(L);
+    C = C * L_new / L;
+    L = L_new;
+
+    RGB rgb_scale = oklab_to_linear_srgb({ L_vt, a_ * C_vt, b_ * C_vt });
+    float scale_L = cbrtf(1.f / fmax(fmax(rgb_scale.r, rgb_scale.g), fmax(rgb_scale.b, 0.f)));
+
+    L = L * scale_L;
+    C = C * scale_L;
+
+    RGB rgb = oklab_to_linear_srgb({ L, C * a_, C * b_ });
+    return {
+        srgb_transfer_function(rgb.r),
+        srgb_transfer_function(rgb.g),
+        srgb_transfer_function(rgb.b),
+    };
+}
+
+static HSV srgb_to_okhsv(RGB rgb)
+{
+    Lab lab = linear_srgb_to_oklab({
+        srgb_transfer_function_inv(rgb.r),
+        srgb_transfer_function_inv(rgb.g),
+        srgb_transfer_function_inv(rgb.b)
+        });
+
+    float C = sqrtf(lab.a * lab.a + lab.b * lab.b);
+    float a_ = lab.a / C;
+    float b_ = lab.b / C;
+
+    float L = lab.L;
+    float h = 0.5f + 0.5f * atan2f(-lab.b, -lab.a) / pi;
+
+    LC cusp = find_cusp(a_, b_);
+    ST ST_max = to_ST(cusp);
+    float S_max = ST_max.S;
+    float T_max = ST_max.T;
+    float S_0 = 0.5f;
+    float k = 1 - S_0 / S_max;
+
+    // first we find L_v, C_v, L_vt and C_vt
+
+    float t = T_max / (C + L * T_max);
+    float L_v = t * L;
+    float C_v = t * C;
+
+    float L_vt = toe_inv(L_v);
+    float C_vt = C_v * L_vt / L_v;
+
+    // we can then use these to invert the step that compensates for the toe and the curved top part of the triangle:
+    RGB rgb_scale = oklab_to_linear_srgb({ L_vt, a_ * C_vt, b_ * C_vt });
+    float scale_L = cbrtf(1.f / fmax(fmax(rgb_scale.r, rgb_scale.g), fmax(rgb_scale.b, 0.f)));
+
+    L = L / scale_L;
+    C = C / scale_L;
+
+    C = C * toe(L) / L;
+    L = toe(L);
+
+    // we can now compute v and s:
+
+    float v = L / L_v;
+    float s = (S_0 + T_max) * C_v / ((T_max * S_0) + T_max * k * C_v);
+
+    return { h, s, v };
+}
+
+
+}; // namespace ok_color
+
+#endif
diff --git a/apps/devApps/random explorer/src/main.cpp b/apps/devApps/random explorer/src/main.cpp
new file mode 100644
index 00000000000..55bb12ac3fc
--- /dev/null
+++ b/apps/devApps/random explorer/src/main.cpp	
@@ -0,0 +1,16 @@
+#include "ofMain.h"
+#include "ofApp.h"
+
+int main( ){
+
+	//Use ofGLFWWindowSettings for more options like multi-monitor fullscreen
+	ofGLWindowSettings settings;
+	settings.setSize(1024, 768);
+	settings.windowMode = OF_WINDOW; //can also be OF_FULLSCREEN
+
+	auto window = ofCreateWindow(settings);
+
+	ofRunApp(window, make_shared<ofApp>());
+	ofRunMainLoop();
+
+}
diff --git a/apps/devApps/random explorer/src/ofApp.cpp b/apps/devApps/random explorer/src/ofApp.cpp
new file mode 100644
index 00000000000..37b81e37e2e
--- /dev/null
+++ b/apps/devApps/random explorer/src/ofApp.cpp	
@@ -0,0 +1,87 @@
+#include "ofApp.h"
+#include "ofRandomEngine.h"
+#include "ofRandomDistributions.h"
+
+void ofApp::setup(){
+    
+    ofSetWindowShape(1920,1000);
+    ofEnableBlendMode(OF_BLENDMODE_ALPHA);
+    ofSetWindowTitle("Random Explorer");
+    
+    panel_.setup("Global state");
+    panel_.add(size_);
+    panel_.add(seed_);
+	panel_.add(reinit_);
+    
+//	panel_.add(ok_color_);
+//	panel_.add(saturation_);
+//	panel_.add(value_);
+//	panel_.add(offset_);
+	
+	for (auto & [name, group]: dists_) {
+		group->panel_.setup();
+		//        group->panel_.disableHeader();
+		for (auto & dist: group->dists_)  {
+			group->panel_.add(dist->parameters_);
+			group->panel_.getGroup(dist->parameters_.getName()).getButton("click for reference").setTextColor(ofColor{60,60,120});
+		}
+	}
+	
+	colorize();
+    
+    auto y = panel_.getPosition().y + panel_.getHeight() + 120;
+
+    dists_["core"]->panel_.setPosition(10,y);
+    dists_["special"]->panel_.setPosition(10+col_w_,y);
+    dists_["of"]->panel_.setPosition(10+col_w_*2,y);
+    dists_["old"]->panel_.setPosition(10+col_w_*2,y+((square_+gap_)*(dists_["of"]->dists_.size()+1)));
+
+    seed_.addListener(this, &ofApp::seed);
+	reinit_.addListener(this, &ofApp::reinit);
+
+	ok_color_.addListener(this, &ofApp::colorize_bool);
+	saturation_.addListener(this, &ofApp::colorize_float);
+	value_.addListener(this, &ofApp::colorize_float);
+	offset_.addListener(this, &ofApp::colorize_float);
+	
+    reinit(); // not required from the random perspective, but keeps the demo together
+    
+}
+
+void ofApp::update(){
+    
+    // DISTRIBUTIONS
+    auto num_samples = pow(size_.get(),3);
+	size_string_ = "samples per distribution per frame: "s+ofToString(num_samples)+" total random calls per frame: "+ofToString(num_samples*20)+" @ "+ofToString(ofGetFrameRate())+"fps";
+    for (auto & [name, group]: dists_) {
+        for (const auto & d: group->dists_) {
+            d->clear();
+            float start = ofGetElapsedTimef();
+            for (std::size_t i=0; i < num_samples; i++) d->gen();
+            d->cost_ = ofGetElapsedTimef() - start;
+            d->compile();
+        }
+    }
+}
+
+void ofApp::draw(){
+    
+    auto x =  220;
+	ofDrawBitmapStringHighlight(size_string_, x, 40);
+    ofDrawBitmapStringHighlight(dna_string_, x, 80);
+	ofDrawBitmapStringHighlight(shuffle_string_, x, 100);
+	ofDrawBitmapStringHighlight("other random functions (like std::shuffle) can feed from the same thread-safe engine, ensuring coherence", x, 140);
+
+    if (of::random::Engine::instance()->is_deterministic()) {
+        ofDrawBitmapStringHighlight("engine is deterministic (seeded)", x, 20, ofColor::black, ofColor::green);
+    } else {
+        ofDrawBitmapStringHighlight("engine is non-deterministic", x, 20, ofColor::black, ofColor::white);
+    }
+
+    panel_.draw();
+    dists_["core"]->draw("C++ fundamental distributions", square_, gap_);
+    dists_["special"]->draw("more specialized distributions", square_, gap_);
+    dists_["of"]->draw("OF/art-centric wrappers/utils", square_, gap_);
+    dists_["old"]->draw("Previous implementation (reference)", square_, gap_);
+    
+}
diff --git a/apps/devApps/random explorer/src/ofApp.h b/apps/devApps/random explorer/src/ofApp.h
new file mode 100644
index 00000000000..2a953f67909
--- /dev/null
+++ b/apps/devApps/random explorer/src/ofApp.h	
@@ -0,0 +1,197 @@
+#pragma once
+
+#include "ofMain.h"
+#include "ofxGui.h"
+#include "Dist.hpp"
+#include "OKColor.h"
+
+class ofApp : public ofBaseApp {
+
+    ofxPanel panel_;
+    std::vector<std::shared_ptr<ofxPanel>> group_panels_;
+    ofParameter<size_t> size_{"size (cube root)", 25, 1, 50};
+    ofParameter<unsigned long> seed_{"seed", 0, 0, 1000};
+    ofParameter<void> reinit_{"re-init engine"};
+	ofParameter<bool> ok_color_{"ok_color"};
+	ofParameter<float> saturation_{"saturation", 0.9};
+	ofParameter<float> value_{"value", .25};
+	ofParameter<float> offset_{"offset", 0};
+
+    size_t col_w_ = 640;
+    size_t square_ = 110;
+    size_t gap_ = 20;
+
+    std::array<uint8_t, 8> dna_;
+    std::vector<uint8_t> shuffle_{1,2,3,4,5,6,7,8};
+    std::string dna_string_;
+	std::string shuffle_string_;
+	std::string size_string_;
+
+    ofParameter<float> 		rand_min_		{"min", 0, 0, 100};
+    ofParameter<float> 		rand_max_		{"max", 100, 0, 100};
+    ofParameter<float> 		uni_min_		{"min", 0, 0, 100};
+    ofParameter<float> 		uni_max_		{"max", 100, 0, 100};
+    ofParameter<int> 		uni_int_min_	{"min", 1, 0, 9};
+    ofParameter<int> 		uni_int_max_	{"max", 6, 1, 10};
+    ofParameter<float>      norm_mean_  	{"mean", 100, 0, 200};
+    ofParameter<float>      norm_dev_   	{"stddev", 20, 0, 50};
+    ofParameter<float>      log_mean_   	{"mean", 1.6, 0, 5};
+    ofParameter<float>      log_dev_    	{"stddev", 0.25, 0, 10};
+    ofParameter<float> 		bn_min_			{"min", 35, 0, 200};
+    ofParameter<float> 		bn_max_			{"max", 165, 0, 200};
+    ofParameter<float> 		bn_focus_		{"focus", 1.75, 0.01, 20};
+    ofParameter<float> 		gamma_alpha_	{"alpha", 3, 1, 20};
+    ofParameter<float> 		gamma_beta_		{"beta", 5, 0, 50};
+    ofParameter<float> 		poiss_mean_		{"mean", 4, 0, 20};
+    ofParameter<float> 		exp_lambda_		{"lambda", .1, 0.01, 0.5};
+    ofParameter<glm::vec2> 	vec_min_		{"min", {40,10}, {0,0}, {200,200}};
+    ofParameter<glm::vec2> 	vec_max_		{"max", {60,90}, {0,0}, {200,200}};
+    ofParameter<glm::vec2> 	vec_mean_		{"mean", {50,50}, {0,0}, {100,100}};
+    ofParameter<glm::vec2> 	vec_dev_		{"stddev", {50,15}, {0,0}, {50,50}};
+    ofParameter<float> 		vec_gamma_a_	{"alpha", 1.5, 1, 20};
+    ofParameter<float>      vec_gamma_b_	{"beta", 10, 0, 50};
+    ofParameter<float>      yes_        	{"yes", .33, 0, 1};
+    ofParameter<int>        bin_p_      	{"p", 5, 0, 40};
+    ofParameter<float>      bin_t_      	{"t", .5, 0, .99};
+    ofParameter<float>      chi_n_      	{"freedom", 4, 0, 10};
+
+    std::map<std::string, std::shared_ptr<DistGroup>> dists_ {
+        {
+            "old", std::make_shared<DistGroup>(std::vector<std::shared_ptr<Dist>>  {
+                
+                std::make_shared<ConcreteDist<float>>("rand()", "old-skool (thread risky)\nbased on srand()", "binomial_distribution",
+                                                      std::vector<ofAbstractParameter *>{&rand_min_, &rand_max_}, [&]
+                                                      { return ofRandom(rand_min_, rand_max_); } ) } )
+        }, {
+            "core", std::make_shared<DistGroup>(std::vector<std::shared_ptr<Dist>>  {
+                
+                std::make_shared<ConcreteDist<float>>("uniform <real>", "modern, thread safe\nright bound improbable", "uniform_real_distribution",
+                                                      std::vector<ofAbstractParameter *>{&uni_min_, &uni_max_}, [&]
+                                                      { return of::random::uniform<float>(uni_min_, uni_max_); }),
+                
+                std::make_shared<ConcreteDist<int>>("uniform <int>", "unambiguously\nincludes bounds", "uniform_int_distribution",
+                                                    std::vector<ofAbstractParameter *>{&uni_int_min_, &uni_int_max_}, [&]
+                                                    { return of::random::uniform<int>(uni_int_min_, uni_int_max_); }, 11, glm::vec2{0,  10}, true),
+                
+                std::make_shared<ConcreteDist<double>>("normal", "also aliased to random::gaussian\for stochastic familiarity", "normal_distribution",
+                                                      std::vector<ofAbstractParameter *>{&norm_mean_, &norm_dev_}, [&]
+                                                      { return of::random::normal<double>(norm_mean_, norm_dev_); }, 101, glm::vec2{0, 200} ),
+                
+                std::make_shared<ConcreteDist<float>>("exponential", "all in the title", "exponential_distribution",
+                                                      std::vector<ofAbstractParameter *>{&exp_lambda_}, [&]
+                                                      { return of::random::exponential<float>(exp_lambda_);}),
+                
+               std::make_shared<ConcreteDist<int>>("poisson", "for sparse, \npositive discrete\n(if mean > 12: ~= normal)", "poisson_distribution",
+                                                    std::vector<ofAbstractParameter *>{&poiss_mean_}, [&]
+                                                    { return of::random::poisson<int>(poiss_mean_); }, 21, glm::vec2{0,  20}, true),
+                
+                // note int is used here instead of bool because of compiler complaints about "bit values" whatever
+                std::make_shared<ConcreteDist<int>>("bernoulli", "binary prob\ne.g. toss of a coin\naliased to of::random::yes",  "bernoulli_distribution",
+                                                    std::vector<ofAbstractParameter *>{&yes_}, [&]
+                                                    { return of::random::bernoulli(yes_); }, 2, glm::vec2{0,1}, true) } )
+        }, {
+            "special", std::make_shared<DistGroup>(std::vector<std::shared_ptr<Dist>>  {
+                
+                std::make_shared<ConcreteDist<int>>("lognormal", "like normal but log", "lognormal_distribution",
+                                                    std::vector<ofAbstractParameter *>{&log_mean_, &log_dev_},[&]
+                                                    { return of::random::lognormal<int>(log_mean_, log_dev_); }, 21, glm::vec2{0, 20}, true),
+                
+                std::make_shared<ConcreteDist<float>>("gamma", "for special purposes", "gamma_distribution",
+                                                      std::vector<ofAbstractParameter *>{&gamma_alpha_, &gamma_beta_}, [&]
+                                                      { return of::random::gamma<int>(gamma_alpha_, gamma_beta_); } ),
+                
+                std::make_shared<ConcreteDist<float>>("chi-squared", "cool energy bulge", "chi_squared_distribution",
+                                                      std::vector<ofAbstractParameter *>{&chi_n_}, [&]
+                                                      { return of::random::chi_squared<float>(chi_n_);}, 101, glm::vec2{0,  20} ),
+                
+                std::make_shared<ConcreteDist<int>>("binomial", "~normal for ints\nsquashes nicely on zero", "binomial_distribution",
+                                                    std::vector<ofAbstractParameter *>{&bin_p_, &bin_t_}, [&]
+                                                    { return of::random::binomial(bin_p_,bin_t_); }, 21, glm::vec2{0,  20}, true),
+                
+                std::make_shared<ConcreteDist<int>>("geometric", "~expon for sparse ints\nwith no parameter", "geometric_distribution",
+                                                    std::vector<ofAbstractParameter *>{}, [&]
+                                                    { return of::random::geometric<int>(); }, 21, glm::vec2{0,  20}, true) } )
+        }, {
+            "of", std::make_shared<DistGroup>(std::vector<std::shared_ptr<Dist>>  {
+                
+                std::make_shared<ConcreteDist<float>>("bound normal", "practical min/max\nwith enforced limits\n(allows to \"focus\")\nintuitiver", "normal_distribution",
+                                                      std::vector<ofAbstractParameter *>{&bn_min_, &bn_max_, &bn_focus_}, [&]
+                                                      { return ofRandomBoundNormal<int>(bn_min_, bn_max_, bn_focus_); }, 101, glm::vec2{0, 200}),
+                
+                std::make_shared<ConcreteDist<glm::vec2>>("vec2 uniform/axis", "all generators can\nbe typed vec2, 3 or 4", "uniform_int_distribution",
+                                                          std::vector<ofAbstractParameter *>{&vec_min_, &vec_max_}, [&]
+                                                          { return of::random::uniform<glm::vec2>(vec_min_, vec_max_); }, 101, glm::vec2{0, 200}),
+                
+                std::make_shared<ConcreteDist<glm::vec2>>("vec2 normal/axis", "vecN version can be\nrefined to vecN axis", "normal_distribution",
+                                                          std::vector<ofAbstractParameter *>{&vec_mean_, &vec_dev_}, [&]
+                                                          { return of::random::gaussian<glm::vec2>(vec_mean_, vec_dev_); }, 101, glm::vec2{0, 200}),
+                
+                std::make_shared<ConcreteDist<glm::vec3>>("vec3 gamma", "", "gamma_distribution",
+                                                          std::vector<ofAbstractParameter *>{&vec_gamma_a_, &vec_gamma_b_}, [&]
+                                                          { return of::random::gamma<glm::vec3>(vec_gamma_a_, vec_gamma_b_); }, 101, glm::vec2{0, 200} ) } )
+        }
+    };
+    
+public:
+    auto rebuild_signatures() {
+        dna_string_= "1. first 8 numbers: ";
+        for (size_t i = 0; i < dna_.size(); i++) {
+            dna_[i] = of::random::uniform<uint8_t>(0,255);
+            dna_string_ += ofToUpper(ofToHex(dna_[i])) + " ";
+        }
+        
+        shuffle_string_= "2. shuffled vector: ";
+        shuffle_ = { 1,2,3,4,5,6,7,8 }; // initial order
+        of::random::shuffle(shuffle_);
+        for (size_t i = 0; i < shuffle_.size(); i++) {
+            shuffle_string_ += ofToUpper(ofToHex(shuffle_[i])) + " ";
+        }
+    }
+    
+    auto seed (unsigned long & seed) {
+        of::random::engine()->seed(seed);
+        rebuild_signatures();
+        panel_.getControl("seed")->setTextColor(ofColor::green);
+    }
+    
+    auto reinit() {
+        seed_ = 0;
+        of::random::Engine::destruct();
+        of::random::Engine::construct(); // from scratch
+        rebuild_signatures();
+        panel_.getControl("seed")->setTextColor(ofColor::gray);
+    }
+    
+	auto colorize() {
+		
+		size_t sum = 0;
+		for (auto& group : {dists_["core"], dists_["special"], dists_["of"]}) sum+=group->dists_.size();
+		auto chunk = 255.0/sum;
+		
+		size_t i = 0; // to spread the hue; order below matters
+		for (auto & group:  {dists_["core"], dists_["special"], dists_["of"], dists_["old"]} ) {
+			for (auto & dist: group->dists_) {
+				if (i>=sum) {
+					dist->color_ = {48,48,48,255}; // sad
+				} else {
+					if (ok_color_) {
+						auto color = OKColor::okhsv_to_srgb({float(fmod((((chunk * i++)/256.0f)+offset_),1)), saturation_.get(), value_.get()});
+						dist->color_=ofColor(color.r*255.0, color.g*255.0, color.b*255.0);
+					} else {
+						dist->color_.setHsb(int((chunk * i++) + (offset_*255))%255, saturation_*255.0, value_*255.0);
+					}
+				}
+				group->panel_.getGroup(dist->parameters_.getName()).setHeaderBackgroundColor(dist->color_);
+			}
+		}
+	}
+	
+	auto colorize_float(float & v) { colorize(); } // for non-bool ofParameters
+	auto colorize_bool(bool & v) { colorize(); } // for non-bool for ofParameters
+
+	
+    void setup() override;
+    void update() override;
+    void draw() override;
+    
+};
diff --git a/examples/math/randomExample/README.md b/examples/math/randomExample/README.md
new file mode 100644
index 00000000000..8fe39bc507d
--- /dev/null
+++ b/examples/math/randomExample/README.md
@@ -0,0 +1,27 @@
+#level1_randomExample
+--
+### Learning Objectives
+
+This openFrameworks Example is designed to demonstrate the OF random Engine and it's seeding. 
+
+During the ofBaseApp initializations (i.e. before the setup() call), OF instanciates a thread-safe random engine based on Mersenne-Twister, with a non-deterministic random seed sequence. This provides a centralized source of unpredictable randomness (which is generally what is expected by a random function).
+
+However if desired, you can supply a determinist seed. This makes it possible to generate deterministic works (or test things with repeatable random).
+
+This example demonstrates the Engine coupled with the ofShuffle function.
+
+### Expected Behavior
+
+When launching this app, you should see:
+
+* numbers 1-8, in an unpredictably shuffled order.
+
+### Instructions for use:
+
+* Drag the seed value to activate determinism and notice how the sequence is shuffled differently. For a given seed, this random order will be the same, on all C++17 platforms (until the Mersenne-Twister implementation or std::uniform_int_distributions are modified or deprecated).
+
+* Click reset to destroy and re-instanciate the Engine, with a fresh random seed sequence.
+
+### Other classes used in this file
+
+This Example uses ofxGui. 
diff --git a/examples/math/randomExample/addons.make b/examples/math/randomExample/addons.make
new file mode 100644
index 00000000000..f47512dc0e8
--- /dev/null
+++ b/examples/math/randomExample/addons.make
@@ -0,0 +1 @@
+ofxGui
diff --git a/examples/math/randomExample/bin/data/.gitkeep b/examples/math/randomExample/bin/data/.gitkeep
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/examples/math/randomExample/src/main.cpp b/examples/math/randomExample/src/main.cpp
new file mode 100644
index 00000000000..b14f12ffb8e
--- /dev/null
+++ b/examples/math/randomExample/src/main.cpp
@@ -0,0 +1,17 @@
+#include "ofMain.h"
+#include "ofApp.h"
+
+//========================================================================
+int main( ){
+
+	//Use ofGLFWWindowSettings for more options like multi-monitor fullscreen
+	ofGLWindowSettings settings;
+	settings.setSize(1024, 768);
+	settings.windowMode = OF_WINDOW; //can also be OF_FULLSCREEN
+
+	auto window = ofCreateWindow(settings);
+
+	ofRunApp(window, make_shared<ofApp>());
+	ofRunMainLoop();
+
+}
diff --git a/examples/math/randomExample/src/ofApp.cpp b/examples/math/randomExample/src/ofApp.cpp
new file mode 100644
index 00000000000..4d4373f8a0a
--- /dev/null
+++ b/examples/math/randomExample/src/ofApp.cpp
@@ -0,0 +1,25 @@
+#include "ofApp.h"
+
+void ofApp::setup(){
+        
+    panel_.setup("Engine Seeding");
+    panel_.add(seed_);
+    panel_.add(reconstruct_);
+    
+    seed_.addListener(this, &ofApp::seed);
+    reconstruct_.addListener(this, &ofApp::reconstruct);
+    
+    reconstruct();
+}
+
+void ofApp::draw(){
+    
+    panel_.draw();
+    
+    ofDrawBitmapStringHighlight(shuffled_string_, 230, 80);
+    if (of::random::engine()->is_deterministic()) {
+        ofDrawBitmapStringHighlight("engine is deterministic (seeded)", 230, 40, ofColor::black, ofColor::green);
+    } else {
+        ofDrawBitmapStringHighlight("engine is non-deterministic", 230, 40, ofColor::black, ofColor::white);
+    }
+}
diff --git a/examples/math/randomExample/src/ofApp.h b/examples/math/randomExample/src/ofApp.h
new file mode 100644
index 00000000000..6996ef47198
--- /dev/null
+++ b/examples/math/randomExample/src/ofApp.h
@@ -0,0 +1,55 @@
+#pragma once
+
+#include "ofMain.h"
+#include "ofxGui.h"
+
+class ofApp : public ofBaseApp{
+    
+public:
+    
+    ofxPanel panel_;
+    std::string shuffled_string_;
+    ofParameter<unsigned long> seed_{"seed", 0, 0, 1000};
+    ofParameter<void> reconstruct_{"re-construct engine"};
+
+    void setup() override;
+    void draw() override;
+    
+    auto perform_shuffle() {
+        std::vector<uint8_t> values = { 1, 2, 3, 4, 5, 6, 7, 8 }; // initial order
+        ofShuffle(values);
+        shuffled_string_ = "shuffled values:";
+        for (const auto v: values) {
+            shuffled_string_ += " " + ofToUpper(ofToHex(v));
+        }
+    }
+    
+    auto seed (unsigned long & seed) {
+        of::random::seed(seed);
+        perform_shuffle();
+        panel_.getControl("seed")->setTextColor(ofColor::green);
+    }
+    
+    void reconstruct() {
+        seed_ = 0;
+        of::random::Engine::destruct();
+        of::random::Engine::construct(); // from scratch
+        perform_shuffle();
+        panel_.getControl("seed")->setTextColor(ofColor::gray);
+    }
+    
+    //  MARK: - ofBaseApp interface boilerplate:
+    //  void update() override;
+    //	void keyPressed(int key) override;
+    //	void keyReleased(int key) override;
+    //	void mouseMoved(int x, int y ) override;
+    //	void mouseDragged(int x, int y, int button) override;
+    //	void mousePressed(int x, int y, int button) override;
+    //	void mouseReleased(int x, int y, int button) override;
+    //	void mouseEntered(int x, int y) override;
+    //	void mouseExited(int x, int y) override;
+    //	void windowResized(int w, int h) override;
+    //	void dragEvent(ofDragInfo dragInfo) override;
+    //	void gotMessage(ofMessage msg) override;
+    
+};
diff --git a/libs/openFrameworks/math/ofMath.cpp b/libs/openFrameworks/math/ofMath.cpp
index 2b5ef740381..768aaf4a50b 100644
--- a/libs/openFrameworks/math/ofMath.cpp
+++ b/libs/openFrameworks/math/ofMath.cpp
@@ -18,9 +18,6 @@ int ofNextPow2(int a){
 //--------------------------------------------------
 void ofSeedRandom() {
 
-	// good info here:
-	// http://stackoverflow.com/questions/322938/recommended-way-to-initialize-srand
-
 	#ifdef TARGET_WIN32
 		srand(GetTickCount());
 	#elif !defined(TARGET_EMSCRIPTEN)
@@ -38,8 +35,12 @@ void ofSeedRandom() {
 }
 
 //--------------------------------------------------
+void ofSetRandomSeed(unsigned long new_seed) {
+    of::random::seed(new_seed);
+}
+
 void ofSeedRandom(int val) {
-	srand((long) val);
+    ofSetRandomSeed(val);
 }
 
 //--------------------------------------------------
diff --git a/libs/openFrameworks/math/ofMath.h b/libs/openFrameworks/math/ofMath.h
index 5a8e9c19afa..7291df07da0 100644
--- a/libs/openFrameworks/math/ofMath.h
+++ b/libs/openFrameworks/math/ofMath.h
@@ -82,11 +82,19 @@ float ofRandomWidth();
 /// \returns a random number between 0 and ofGetHeight().
 float ofRandomHeight();
 
-/// \brief Seed the seeds the random number generator with a unique value.
+/// \brief Seed the random number generator.
+///
+/// This passes a seed value to the random engine;
+/// see of::random::Engine for details
+///
+/// \param val The value with which to seed the generator.
+void ofSetRandomSeed(unsigned long new_seed);
+
+/// \brief Seeds the random number generator with a unique value.
 ///
-/// This seeds the random number generator with an acceptably random value, 
+/// This seeds the old-school srand-based random number generator with an acceptably random value, 
 /// generated from clock time and the PID.
-void ofSeedRandom();
+[[deprecated("to get a new non-deterministic state, destroy and re-construct the of::random::Engine")]] void ofSeedRandom();
 
 /// \brief Seed the random number generator.
 ///
@@ -95,7 +103,7 @@ void ofSeedRandom();
 /// setup.  This can be useful for debugging and testing.
 ///
 /// \param val The value with which to seed the generator.
-void ofSeedRandom(int val);
+[[deprecated("use ofSetRandomSeed() or of::random::seed() instead")]] void ofSeedRandom(int val);
 
 /// \}
 
diff --git a/libs/openFrameworks/utils/ofRandomDistributions.h b/libs/openFrameworks/utils/ofRandomDistributions.h
new file mode 100644
index 00000000000..18d88115209
--- /dev/null
+++ b/libs/openFrameworks/utils/ofRandomDistributions.h
@@ -0,0 +1,667 @@
+#ifndef OF_RANDOM_DISTRIBUTIONS_H_
+#define OF_RANDOM_DISTRIBUTIONS_H_
+
+#include <random>
+#include <glm/glm.hpp>
+
+// https://gist.github.com/imneme/540829265469e673d045
+// https://github.com/effolkronium/random/tree/master
+// https://www.theanalysisfactor.com/differences-between-normal-and-poisson-distributions
+// https://blogs.sas.com/content/iml/2019/07/22/extreme-value-normal-data.html
+// https://rovdownloads.com/blog/quick-overview-of-probability-distributions/
+
+namespace of::random
+{
+
+// MARK: - UNIFORM
+
+// special case of distinguishing real and integer for uniform distribution
+
+template<typename T = float, typename ... Args>
+std::enable_if_t<std::is_floating_point_v<T>, T>
+uniform(Args&&... args) {
+    std::uniform_real_distribution<T> distr(std::forward<Args>(args)...);
+    return distr(of::random::gen());
+}
+
+template<typename T, typename ... Args>
+std::enable_if_t<std::is_integral_v<T>, T>
+uniform(Args&&... args) {
+    std::uniform_int_distribution<T> distr(std::forward<Args>(args)...);
+    return distr(of::random::gen());
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+uniform(double min, double max) {
+    std::uniform_real_distribution<float> dist(min,max);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+uniform(T min, T max) {
+    return {
+        std::uniform_real_distribution<float>{min.x, max.x}(of::random::gen()),
+        std::uniform_real_distribution<float>{min.y, max.y}(of::random::gen())
+    };
+}
+
+template <typename T, typename ... Args>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+uniform(Args&&... args) {
+    std::uniform_real_distribution<float> dist(std::forward<Args>(args)...);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+uniform(T min, T max) {
+    return {
+        std::uniform_real_distribution<float>{min.x, max.x}(of::random::gen()),
+        std::uniform_real_distribution<float>{min.y, max.y}(of::random::gen()),
+        std::uniform_real_distribution<float>{min.z, max.z}(of::random::gen())
+    };
+}
+
+// MARK: NORMAL (gaussian)
+
+// better to cast params as double (otherwise casual ints produce unexpected results)
+
+template<typename T = float>
+std::enable_if_t<std::is_arithmetic_v<T>, T>
+normal(T mean, T stddev) {
+    return std::normal_distribution<T>{mean, stddev}(of::random::gen());
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+normal(float mean, float stddev) {
+    std::normal_distribution<float> dist(mean, stddev);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+normal(T mean, T stddev) {
+    return {
+        std::normal_distribution<float>{mean.x, stddev.x}(of::random::gen()),
+        std::normal_distribution<float>{mean.y, stddev.y}(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+normal(float mean, float stddev) {
+    std::normal_distribution<float> dist(mean, stddev);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+normal(T mean, T stddev) {
+    return {
+        std::normal_distribution<float>{mean.x, stddev.x}(of::random::gen()),
+        std::normal_distribution<float>{mean.y, stddev.y}(of::random::gen()),
+        std::normal_distribution<float>{mean.z, stddev.z}(of::random::gen())
+    };
+}
+
+// MARK: LOGNORMAL
+
+template<typename T = float>
+std::enable_if_t<std::is_arithmetic_v<T>, T>
+lognormal(double mean, double stddev) {
+    return std::lognormal_distribution{mean, stddev}(of::random::gen());
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+lognormal(double mean, double stddev) {
+    std::lognormal_distribution dist(mean, stddev);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+lognormal(T mean, T stddev) {
+    return {
+        std::lognormal_distribution{mean.x, stddev.x}(of::random::gen()),
+        std::lognormal_distribution{mean.y, stddev.y}(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+lognormal(double mean, double stddev) {
+    std::lognormal_distribution dist(mean, stddev);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+lognormal(T mean, T stddev) {
+    return {
+        std::lognormal_distribution{mean.x, stddev.x}(of::random::gen()),
+        std::lognormal_distribution{mean.y, stddev.y}(of::random::gen()),
+        std::lognormal_distribution{mean.z, stddev.z}(of::random::gen())
+    };
+}
+
+// MARK: GAMMA
+
+template<typename T = float>
+std::enable_if_t<std::is_arithmetic_v<T>, T>
+gamma(double alpha, double beta) {
+    return std::gamma_distribution{alpha, beta}(of::random::gen());
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+gamma(double alpha, double beta) {
+    std::gamma_distribution dist(alpha, beta);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+gamma(T alpha, T beta) {
+    return {
+        std::gamma_distribution{alpha.x, beta.x}(of::random::gen()),
+        std::gamma_distribution{alpha.y, beta.y}(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+gamma(double alpha, double beta) {
+    std::gamma_distribution dist(alpha, beta);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+gamma(T alpha, T beta) {
+    return {
+        std::gamma_distribution{alpha.x, beta.x}(of::random::gen()),
+        std::gamma_distribution{alpha.y, beta.y}(of::random::gen()),
+        std::gamma_distribution{alpha.z, beta.z}(of::random::gen())
+    };
+}
+// MARK: POISSON
+
+template<typename T = float>
+std::enable_if_t<std::is_arithmetic_v<T>, T>
+poisson(double mean) {
+    return std::poisson_distribution<T>{mean}(of::random::gen());
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+poisson(double mean) {
+    std::poisson_distribution<T> dist(mean);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+poisson(T mean) {
+    return {
+        std::poisson_distribution<T>{mean.x}(of::random::gen()),
+        std::poisson_distribution<T>{mean.y}(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+poisson(double mean) {
+    std::poisson_distribution<T> dist(mean);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+poisson(T mean) {
+    return {
+        std::poisson_distribution<T>{mean.x}(of::random::gen()),
+        std::poisson_distribution<T>{mean.y}(of::random::gen()),
+        std::poisson_distribution<T>{mean.z}(of::random::gen())
+    };
+}
+
+// MARK: EXPONENTIAL
+
+template<typename T = float>
+std::enable_if_t<std::is_arithmetic_v<T>, T>
+exponential(T lambda) {
+    return std::exponential_distribution<T>{lambda}(of::random::gen());
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+exponential(double lambda) {
+    std::exponential_distribution<T> dist(lambda);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+exponential(T lambda) {
+    return {
+        std::exponential_distribution<T>{lambda.x}(of::random::gen()),
+        std::exponential_distribution<T>{lambda.y}(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+exponential(double lambda) {
+    std::exponential_distribution<T> dist(lambda);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+exponential(T lambda) {
+    return {
+        std::exponential_distribution<T>{lambda.x}(of::random::gen()),
+        std::exponential_distribution<T>{lambda.y}(of::random::gen()),
+        std::exponential_distribution<T>{lambda.z}(of::random::gen())
+    };
+}
+
+// MARK: CHI_SQUARED
+
+template<typename T = float>
+std::enable_if_t<std::is_arithmetic_v<T>, T>
+chi_squared(T n) {
+    return std::chi_squared_distribution<T>{n}(of::random::gen());
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+chi_squared(double n) {
+    std::chi_squared_distribution<T> dist(n);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+chi_squared(T n) {
+    return {
+        std::chi_squared_distribution<T>{n.x}(of::random::gen()),
+        std::chi_squared_distribution<T>{n.y}(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+chi_squared(double n) {
+    std::chi_squared_distribution<T> dist(n);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+chi_squared(T n) {
+    return {
+        std::chi_squared_distribution<T>{n.x}(of::random::gen()),
+        std::chi_squared_distribution<T>{n.y}(of::random::gen()),
+        std::chi_squared_distribution<T>{n.z}(of::random::gen())
+    };
+}
+
+// MARK: BINOMIAL
+
+template<typename T = int>
+std::enable_if_t<std::is_arithmetic_v<T>, T>
+binomial(int p, double t) {
+    if (t>=1) {
+        std::cout << "of::random::binomial(): t must be < 1.0\n";
+        return 0;
+    } else {
+        std::binomial_distribution<T> dist(p,t);
+        return dist(of::random::gen());
+    }
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+binomial(int p, double t) {
+    std::binomial_distribution<int> dist(p,t);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+binomial(T p,T t) {
+    return {
+        std::binomial_distribution<int>{p.x,t.x}(of::random::gen()),
+        std::binomial_distribution<int>{p.y,t.y}(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+binomial(int p, double t) {
+    std::binomial_distribution<int> dist(p,t);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+binomial(T p,T t) {
+    return {
+        std::binomial_distribution{p.x,t.x}(of::random::gen()),
+        std::binomial_distribution{p.y,t.y}(of::random::gen()),
+        std::binomial_distribution{p.z,t.z}(of::random::gen())
+    };
+}
+
+// MARK: GEOMETRIC
+
+template<typename T = int>
+std::enable_if_t<std::is_arithmetic_v<T>, T>
+geometric() {
+    return std::geometric_distribution{}(of::random::gen());
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+geometric() {
+    std::geometric_distribution<int> dist;
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+geometric() {
+    std::geometric_distribution<int> dist;
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+// MARK: BERNOUILLI
+
+template <typename T = bool>
+auto bernoulli(double p) {
+    return std::bernoulli_distribution(p)(of::random::gen());
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+bernoulli(double p) {
+    std::bernoulli_distribution dist(p);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+bernoulli(T p) {
+    return {
+        std::bernoulli_distribution{p.x}(of::random::gen()),
+        std::bernoulli_distribution{p.y}(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+bernoulli(double p) {
+    std::bernoulli_distribution dist(p);
+    return {
+        dist(of::random::gen()),
+        dist(of::random::gen()),
+        dist(of::random::gen())
+    };
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+bernoulli(T p) {
+    return {
+        std::bernoulli_distribution{p.x}(of::random::gen()),
+        std::bernoulli_distribution{p.y}(of::random::gen()),
+        std::bernoulli_distribution{p.z}(of::random::gen())
+    };
+}
+// MARK: - CONVENIENT FORWARDING ALIASES
+
+template <typename... Args>
+auto gaussian(Args&&... args) -> decltype(normal(std::forward<Args>(args)...)) {
+  return normal(std::forward<Args>(args)...);
+}
+
+template <typename... Args>
+auto yes(Args&&... args) -> decltype(bernoulli(std::forward<Args>(args)...)) {
+  return bernoulli(std::forward<Args>(args)...);
+}
+
+// whilst nerdy enough, the above do not work with {initializer list} for glm::vec parameters such as:
+// of::random::gaussian<glm::vec2>({0,10},{1,100});
+// so this additional explicit nerdiness is required to wrap them:
+
+template<class T>
+auto gaussian(T min,  T max) -> decltype(normal<T>(min, max)) {
+  return normal(min, max);
+}
+template<class T>
+auto yes(T p) -> decltype(bernoulli(p)) {
+  return bernoulli(p);
+}
+
+// MARK: - OTHER FUNCTIONS
+
+template<class T = float> // works for all non-refined
+T bound_normal(float min, float max, float focus = 4.0f) {
+    if (min >= max) {
+        std::cout << "ofRandomNormalLimits()" << "max must be > than min\n";
+        return {};
+    } else {
+        if (focus <= .0099) {
+            std::cout << "ofRandomNormalLimits()" << "focus must be at least .01\n";
+            return {};
+        } else {
+            T v;
+            do { v = of::random::normal<T>((max+min)/2.0f, (max-min)/(2*focus)); } while (v < min || v > max);
+            return v;
+        }
+    }
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+bound_normal(T min, T max, T focus = {4.0f, 4.0f}) {
+    if (!min.x < max.x || !min.y < max.y) {
+        std::cout << "ofRandomNormalLimits()" << "max must be > than min\n";
+        return {};
+    } else {
+        if (focus.x < 1 || focus.y < 1 ) {
+            std::cout << "ofRandomNormalLimits()" << "focus must be at least 1\n";
+            return {};
+        } else {
+            T v;
+            do { v.x = of::random::normal<float>((max.x+min.x)/2.0f, (max.x-min.x)/(2*focus.x)); } while (v.x < min || v.x > max);
+            do { v.y = of::random::normal<float>((max.y+min.y)/2.0f, (max.y-min.y)/(2*focus.y)); } while (v.y < min || v.y > max);
+            return v;
+        }
+    }
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+bound_normal(T min, T max, T focus = {4.0f, 4.0f, 4.0f}) {
+    if (!min.x < max.x || !min.y < max.y || !min.z < max.z) {
+        std::cout << "ofRandomNormalLimits()" << "max must be > than min\n";
+        return {};
+    } else {
+        if (focus.x < 1 || focus.y < 1 || focus.z < 1) {
+            std::cout << "ofRandomNormalLimits()" << "focus must be at least 1\n";
+            return {};
+        } else {
+            T v;
+            do { v.x = of::random::normal<float>((max.x+min.x)/2.0f, (max.x-min.x)/(2*focus.x)); } while (v.x < min || v.x > max);
+            do { v.y = of::random::normal<float>((max.y+min.y)/2.0f, (max.y-min.y)/(2*focus.y)); } while (v.y < min || v.y > max);
+            do { v.z = of::random::normal<float>((max.z+min.z)/2.0f, (max.z-min.z)/(2*focus.z)); } while (v.z < min || v.z > max);
+            return v;
+        }
+    }
+}
+
+} // end namespace of::random
+
+// MARK: - ANONYMOUS NAMESPACE FUNCTIONS
+
+namespace {
+
+// again, 2 templates per function, for the convenience of implicit brace init of glm::vec parameters
+
+template<class T = float, typename ... Args>
+T ofRandomUniform(Args&&... args) { return of::random::uniform<T>(std::forward<Args>(args)...); }
+
+template<class T = float>
+T ofRandomUniform(T min, T max) { return of::random::uniform<T>(min, max); }
+
+template<class T, typename ... Args>
+T ofRandomNormal(Args&&... args) { return of::random::normal<T>(std::forward<Args>(args)...); }
+
+template<class T>
+T ofRandomNormal(T min, T max) { return of::random::normal<T>(min, max); }
+
+template<class T, typename ... Args>
+T ofRandomGaussian(Args&&... args) { return of::random::gaussian<T>(std::forward<Args>(args)...); }
+
+template<class T>
+T ofRandomGaussian(T mean, T stddev) { return of::random::gaussian<T>(mean, stddev); }
+
+template<class T, typename ... Args>
+T ofRandomBernoulli(Args&&... args) { return of::random::bernoulli<T>(std::forward<Args>(args)...); }
+
+template<class T>
+T ofRandomBernoulli(T prob) { return of::random::bernoulli<T>(prob); }
+
+template<class T, typename ... Args>
+T ofRandomYes(Args&&... args) { return of::random::yes<T>(std::forward<Args>(args)...); }
+
+template<class T>
+T ofRandomYes(T prob) { return of::random::yes<T>(prob); }
+
+template<class T, typename ... Args>
+T ofRandomPoisson(Args&&... args) { return of::random::poisson<T>(std::forward<Args>(args)...); }
+
+template<class T>
+T ofRandomPoisson(T mean) { return of::random::poisson<T>(mean); }
+
+template<class T, typename ... Args>
+T ofRandomExponential(Args&&... args) { return of::random::exponential<T>(std::forward<Args>(args)...); }
+
+template<class T>
+T ofRandomExponential(T lambda) { return of::random::exponential<T>(lambda); }
+
+template<class T, typename ... Args>
+T ofRandomChiSquared(Args&&... args) { return of::random::chi_squared<T>(std::forward<Args>(args)...); }
+
+template<class T>
+T ofRandomChiSquared(T freedom) { return of::random::chi_squared<T>(freedom); }
+
+template<class T, typename ... Args>
+T ofRandomGamma(Args&&... args) { return of::random::gamma<T>(std::forward<Args>(args)...); }
+
+template<class T>
+T ofRandomGamma(T a, T b) { return of::random::gamma<T>(a, b); }
+
+template<class T, typename ... Args>
+T ofRandomLogNormal(Args&&... args) { return of::random::lognormal<T>(std::forward<Args>(args)...); }
+
+template<class T>
+T ofRandomLogNormal(T mean, T stddev) { return of::random::lognormal<T>(mean, stddev); }
+
+template<class T, typename ... Args>
+T ofRandomGeometric(Args&&... args) { return of::random::geometric<T>(std::forward<Args>(args)...); }
+
+template<class T>
+T ofRandomGeometric() { return of::random::geometric<T>(); }
+
+template<class T = float>
+T ofRandomBoundNormal(float min, float max, float focus = 4.0f) {
+    return of::random::bound_normal(min, max, focus);
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec2>, T>
+ofRandomBoundNormal(T min, T max, T focus = {4.0f, 4.0f}) {
+    return of::random::bound_normal(min, max, focus);
+}
+
+template <typename T>
+std::enable_if_t<std::is_same_v<T, glm::vec3>, T>
+ofRandomBoundNormal(T min, T max, T focus = {4.0f, 4.0f, 4.0f}) {
+    return of::random::bound_normal(min, max, focus);
+}
+
+} // end anonymous namespace
+#endif // OF_RANDOM_DISTRIBUTIONS_H_
diff --git a/libs/openFrameworks/utils/ofRandomEngine.h b/libs/openFrameworks/utils/ofRandomEngine.h
new file mode 100644
index 00000000000..ce4d5c9a9fa
--- /dev/null
+++ b/libs/openFrameworks/utils/ofRandomEngine.h
@@ -0,0 +1,74 @@
+#ifndef OF_RANDOM_HPP_
+#define OF_RANDOM_HPP_
+
+#include <random>
+#include <glm/glm.hpp>
+
+#include "ofSingleton.hpp"
+#include "ofMath.h"
+
+namespace of::random
+{
+
+// https://stackoverflow.com/questions/25360241/using-random-number-generator-multiple-instances-or-singleton-approach
+// https://simplecxx.github.io/2018/11/03/seed-mt19937.html
+
+/// \class of::random::Engine
+///
+/// An mt19937 instance wrapped in a singleton, with default non-deterministic seeding
+/// A balance is aimed between ease-of-use, and versatility, and minimize reliance on advanced concepts
+/// (e.g. the lib is templated, but default usage does not require template parameters)
+/// The goal is to have a centralized, thread-safe source of randomness that can be deterministic or not.
+///
+class Engine: public of::utils::Singleton<Engine> {
+    
+	std::random_device rd_{ };
+	std::seed_seq seq_{ rd_(), rd_(), rd_(), rd_() }; // 4 is considered fine for non-cryptographic needs
+	std::mt19937 gen_{ seq_ };
+	bool deterministic_{ false }; // by default the degine is non-deterministic (unpredictable)
+    
+public:
+    
+    Engine() {
+        ofSeedRandom(); // called to maintain "parallelism" until old-school srand() is phase-out of OF
+    }
+
+    /// return the generator for use in random distributions or functions
+    ///
+    /// \returns a reference to the mt19937 generator
+    auto & gen() { return gen_; }
+    
+    /// passes a value to seed the mt19937 generator
+	void seed(unsigned long new_seed) {
+		deterministic_ = true;
+        gen_.seed(new_seed);
+	}
+    
+    /// \returns true or fall depending if the engine is deterministic
+    auto is_deterministic() const { return deterministic_; }
+};
+
+/// \returns a reference to the engine singleton instance
+inline auto engine() {
+    return of::random::Engine::instance();
+}
+
+/// \returns a reference to the generator within the engine instance
+inline auto & gen() {
+    return of::random::Engine::instance()->gen();
+}
+
+/// Passes a value to seed the mt19937 generator within the ending instance
+inline void seed(unsigned long seed) {
+	 of::random::Engine::instance()->seed(seed);
+}
+
+/// Shuffles the order of the elements within the passed container, using the centralized random engine
+template<class T>
+void shuffle(T & values) {
+    std::shuffle(values.begin(), values.end(), of::random::gen());
+}
+
+} // end namespace of::random
+
+#endif // OF_RANDOM_HPP_
diff --git a/libs/openFrameworks/utils/ofSingleton.hpp b/libs/openFrameworks/utils/ofSingleton.hpp
new file mode 100644
index 00000000000..a8c8965c1a9
--- /dev/null
+++ b/libs/openFrameworks/utils/ofSingleton.hpp
@@ -0,0 +1,58 @@
+#ifndef OF_SINGLETON_HPP_
+#define OF_SINGLETON_HPP_
+
+// atomic C++17 DCLP CRTP singleton adapted by burton@artificiel.org from
+// https://github.com/jimmy-park/singleton/blob/main/include/singleton_dclp.hpp (1d26f91)
+
+#include <cassert>
+#include <mutex>
+#include <atomic>
+
+namespace of::utils
+{
+
+template <typename Derived>
+class Singleton {
+public:
+	template <typename... Args>
+	static void construct(Args&&... args) {
+		struct Dummy : public Derived {
+			using Derived::Derived;
+			void prohibit_construct_from_derived() const override { }
+		};
+
+		if (!instance_.load(std::memory_order_acquire)) {
+			if (std::lock_guard lock { mutex_ }; !instance_.load(std::memory_order_relaxed)) {
+				instance_.store(new Dummy { std::forward<Args>(args)... }, std::memory_order_release);
+			}
+		}
+	}
+
+	static Derived * instance() {
+		auto * the_instance = instance_.load(std::memory_order_acquire);
+		assert(the_instance);
+		return the_instance;
+	}
+
+	static void destruct() {
+		if (auto * the_instance = instance_.exchange(nullptr, std::memory_order_acq_rel)) {
+			delete the_instance;
+		}
+	}
+
+protected:
+	Singleton() = default;
+	Singleton(const Singleton&) = delete;
+	Singleton(Singleton&&) noexcept = delete;
+	Singleton& operator=(const Singleton&) = delete;
+	Singleton& operator=(Singleton&&) noexcept = delete;
+	virtual ~Singleton() = default;
+
+private:
+	virtual void prohibit_construct_from_derived() const = 0;
+	inline static std::atomic<Derived*> instance_ { nullptr };
+	inline static std::mutex mutex_;
+};
+
+} // end namespace of::utils
+#endif // OF_SINGLETON_HPP_
diff --git a/libs/openFrameworks/utils/ofUtils.cpp b/libs/openFrameworks/utils/ofUtils.cpp
index d9f1c91261b..12bdd34cb7c 100644
--- a/libs/openFrameworks/utils/ofUtils.cpp
+++ b/libs/openFrameworks/utils/ofUtils.cpp
@@ -61,7 +61,7 @@ namespace of{
 namespace priv{
 	void initutils(){
         ofResetElapsedTimeCounter();
-        ofSeedRandom();
+        of::random::Engine::construct();
     }
 
 	void endutils(){
diff --git a/libs/openFrameworks/utils/ofUtils.h b/libs/openFrameworks/utils/ofUtils.h
index efce1e5032c..e3c4169c052 100644
--- a/libs/openFrameworks/utils/ofUtils.h
+++ b/libs/openFrameworks/utils/ofUtils.h
@@ -15,6 +15,9 @@
 #include <type_traits>
 #include <random>
 
+#include "ofRandomEngine.h"
+#include "ofRandomDistributions.h"
+
 /// \section Elapsed Time
 /// \brief Reset the elapsed time counter.
 ///
@@ -219,15 +222,24 @@ int ofGetDay();
 /// \returns the current weekday [0-6].
 int ofGetWeekday();
 
+/// \section Containers
+/// \brief Randomly reorder the values in a container.
+/// \tparam T Any container that meets std::shuffle's requirements
+/// which are: ValueSwappable and LegacyRandomAccessIterator.
+
+template<typename ... Args>
+void ofShuffle(Args&&... args) {
+    of::random::shuffle(std::forward<Args>(args)...);
+}
+
 /// \section Vectors
 /// \brief Randomly reorder the values in a vector.
 /// \tparam T the type contained by the vector.
 /// \param values The vector of values to modify.
-/// \sa http://www.cplusplus.com/reference/algorithm/random_shuffle/
+
 template<class T>
 void ofRandomize(std::vector<T>& values) {
-    //switch from random_shuffle ( removed in some C++17 impl )
-    std::shuffle(values.begin(), values.end(), std::default_random_engine(0));
+    of::random::shuffle(values);
 }
 
 /// \brief Conditionally remove values from a vector.