Add new test: maj2random

maj2random is a simplified floating-point hash function derived from SHA-2,
retaining its high-quality entropy compression function modified to permute
entropy from a vec2 (designed for UV coordinates) returning float values
between 0.0 and 1.0.

Source: https://github.com/michaeljclark/maj2random

Author: jserv

README.md

@@ -149,6 +149,7 @@ In `rv32emu` repository, there are some prebuilt ELF files for testing purpose.
 * `dhrystone.elf` : See [rv8-bench](https://github.com/michaeljclark/rv8-bench)
 * `doom.elf` : See [sysprog21/doom_riscv](https://github.com/sysprog21/doom_riscv) [RV32M]
 * `ieee754.elf` : See [tests/ieee754.c](tests/ieee754.c) [RV32F]
+* `maj2random.elf` : See [tests/maj2random.c](tests/maj2random.c) [RV32F]
 * `mandelbrot.elf` : See [tests/mandelbrot.c](tests/mandelbrot.c)
 * `nqueens.elf` : See [tests/nqueens.c](tests/nqueens.c)
 * `pi.elf` : See [tests/pi.c](tests/pi.c) [RV32M]

build/maj2random.elf

@@ -0,0 +1,255 @@
+/*
+ * maj2random
+ *
+ * maj2random is a simplified floating point hash function derived from SHA-2,
+ * retaining its high quality entropy compression function modified to permute
+ * entropy from a vec2 (designed for UV coordinates) returning float values
+ * between 0.0 and 1.0. since maj2_random is a hash function it will return
+ * coherent noise. vector argument can be truncated prior to increase grain.
+ *
+ * maj2random is named after the maj mixing function called within the SHA-2
+ * round function and the 2 is becuase maj2random uses 2 rounds instead of 64.
+ * This seems to be sufficient to create visually diffuse noise with only 3%
+ * of the overhead of the full SHA-256 hashing function. maj2random also pays
+ * homage to arc4random due to the similarity of its name form.
+ *
+ * Source: https://github.com/michaeljclark/maj2random
+ */
+
+#include <math.h>
+typedef unsigned uint;
+typedef struct {
+    float a[2];
+} vec2;
+typedef struct {
+    uint a[2];
+} uvec2;
+
+/* first 8 rounds of the SHA-256 k constant */
+static uint sha256_k[8] = {
+    0x428a2f98u, 0x71374491u, 0xb5c0fbcfu, 0xe9b5dba5u,
+    0x3956c25bu, 0x59f111f1u, 0x923f82a4u, 0xab1c5ed5u,
+};
+
+static uint ror(uint x, int d)
+{
+    return (x >> d) | (x << (32 - d));
+}
+static uint sigma0(uint h1)
+{
+    return ror(h1, 2) ^ ror(h1, 13) ^ ror(h1, 22);
+}
+static uint sigma1(uint h4)
+{
+    return ror(h4, 6) ^ ror(h4, 11) ^ ror(h4, 25);
+}
+static uint ch(uint x, uint y, uint z)
+{
+    return z ^ (x & (y ^ z));
+}
+static uint maj(uint x, uint y, uint z)
+{
+    return (x & y) ^ ((x ^ y) & z);
+}
+static uint gamma0(uint a)
+{
+    return ror(a, 7) ^ ror(a, 18) ^ (a >> 3);
+}
+static uint gamma1(uint b)
+{
+    return ror(b, 17) ^ ror(b, 19) ^ (b >> 10);
+}
+
+static vec2 unorm(uvec2 n)
+{
+    vec2 r = {(float) (n.a[0] & ((1u << 23) - 1u)) / (float) ((1u << 23) - 1u),
+              (float) (n.a[1] & ((1u << 23) - 1u)) / (float) ((1u << 23) - 1u)};
+    return r;
+}
+static uvec2 sign(vec2 v)
+{
+    uvec2 r = {v.a[0] < 0.0f, v.a[1] < 0.0f};
+    return r;
+}
+
+static uvec2 maj_extract(vec2 uv)
+{
+    /*
+     * extract 48-bits of entropy from mantissas to create truncated
+     * two word initialization vector 'W' composed using the 48-bits
+     * of 'uv' entropy rotated and copied to keep the field equalized.
+     * the exponent is ignored because the inputs are expected to be
+     * normalized 'uv' values such as texture coordinates. it would be
+     * beneficial to include the exponent entropy but we can't depend
+     * on frexp or ilogb and log2 would be inaccurate.
+     */
+    uvec2 s = sign(uv);
+    uint x = (uint) (fabsf(uv.a[0]) * (float) (1u << 23)) | (s.a[0] << 23);
+    uint y = (uint) (fabsf(uv.a[1]) * (float) (1u << 23)) | (s.a[1] << 23);
+
+    uvec2 r = {(x) | (y << 24), (y >> 8) | (x << 16)};
+    return r;
+}
+
+static vec2 maj_random(vec2 uv, uint NROUNDS)
+{
+    uint H[8] = {0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u};
+    uint W[2];
+
+    uvec2 st = maj_extract(uv);
+
+    W[0] = st.a[0];
+    W[1] = st.a[1];
+
+    for (uint i = 0; i < NROUNDS; i++) {
+        W[i & 1] = gamma1(W[(i - 2) & 1]) + W[(i - 7) & 1] +
+                   gamma0(W[(i - 15) & 1]) + W[(i - 16) & 1];
+    }
+
+    /* we use N rounds instead of 64 and alternate 2 words of iv in W */
+    for (uint i = 0; i < NROUNDS; i++) {
+        uint T0 =
+            W[i & 1] + H[7] + sigma1(H[4]) + ch(H[4], H[5], H[6]) + sha256_k[i];
+        uint T1 = maj(H[0], H[1], H[2]) + sigma0(H[0]);
+        H[7] = H[6];
+        H[6] = H[5];
+        H[5] = H[4];
+        H[4] = H[3] + T0;
+        H[3] = H[2];
+        H[2] = H[1];
+        H[1] = H[0];
+        H[0] = T0 + T1;
+    }
+    uvec2 u = {H[0] ^ H[1] ^ H[2] ^ H[3], H[4] ^ H[5] ^ H[6] ^ H[7]};
+    return unorm(u);
+}
+
+#include <inttypes.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+typedef unsigned long long ullong;
+
+void sncatprintf(char *buf, size_t buflen, const char *fmt, int64_t n)
+{
+    size_t len = strlen(buf);
+    snprintf(buf + len, buflen - len, fmt, n);
+}
+
+char *format_comma(int64_t n)
+{
+    static char buf[65];
+
+    buf[0] = 0;
+    int n2 = 0;
+    int scale = 1;
+    if (n < 0) {
+        buf[0] = '-';
+        buf[1] = 0;
+        n = -n;
+    }
+    while (n >= 1000) {
+        n2 = n2 + scale * (n % 1000);
+        n /= 1000;
+        scale *= 1000;
+    }
+    sncatprintf(buf, sizeof(buf), "%d", n);
+    while (scale != 1) {
+        scale /= 1000;
+        n = n2 / scale;
+        n2 = n2 % scale;
+        sncatprintf(buf, sizeof(buf), ",%03d", n);
+    }
+    return buf;
+}
+
+char *format_rate(double rate)
+{
+    static char buf[65];
+    char unit = ' ';
+    double divisor = 1;
+    if (rate > 1e9) {
+        divisor = 1e9;
+        unit = 'G';
+    } else if (rate > 1e6) {
+        divisor = 1e6;
+        unit = 'M';
+    } else if (rate > 1e3) {
+        divisor = 1e3;
+        unit = 'K';
+    }
+    snprintf(buf, sizeof(buf), "%.1f %c", rate / divisor, unit);
+    return buf;
+}
+
+vec2 f1(ullong i, ullong j)
+{
+    return (vec2){0.0f, (float) i / (float) j};
+}
+vec2 f2(ullong i, ullong j)
+{
+    return (vec2){(float) i / (float) j, 0.0f};
+}
+vec2 f3(ullong i, ullong j)
+{
+    return (vec2){(float) i / (float) j, (float) i / (float) j};
+}
+
+void test_maj(const char *name,
+              int NROUNDS,
+              ullong count,
+              ullong range,
+              vec2 (*f)(ullong, ullong))
+{
+    vec2 sum = {0.f, 0.f};
+    vec2 var = {0.f, 0.f};
+
+    fflush(stdout);
+    for (ullong i = 0; i < count; i++) {
+        vec2 p = f(i, range);
+        vec2 q = maj_random(p, NROUNDS);
+        sum.a[0] += q.a[0];
+        sum.a[1] += q.a[1];
+        var.a[0] += (q.a[0] - .5f) * (q.a[0] - .5f);
+        var.a[1] += (q.a[1] - .5f) * (q.a[1] - .5f);
+    }
+
+    printf("%-32s%12s%12.5f%12.5f%12.5f%12.5f%12.5f%12.5f\n", name,
+           format_comma(count), sum.a[0] / count, sum.a[1] / count,
+           var.a[0] / count, var.a[1] / count, sqrt(var.a[0] / count),
+           sqrt(var.a[1] / count));
+}
+
+void test_header(const char *name)
+{
+    printf("%-32s%12s%12s%12s%12s%12s%12s%12s\n", name, "count", "mean(x)",
+           "mean(y)", "variance(x)", "variance(y)", "std-dev(x)", "std-dev(y)");
+}
+
+void run_all_tests(int NROUNDS, ullong i, ullong j)
+{
+    char name[32];
+    snprintf(name, sizeof(name), "maj_random (NROUNDS=%d)", NROUNDS);
+    printf("\n");
+    test_header(name);
+    test_maj("(             0, (0 - 1K)/8K )", NROUNDS, i, j, f1);
+    test_maj("(             0, (0 - 8K)/8K )", NROUNDS, j, j, f1);
+    test_maj("( (0 - 1K)/8K ),           0 )", NROUNDS, i, j, f2);
+    test_maj("( (0 - 8K)/8K ),           0 )", NROUNDS, j, j, f2);
+    test_maj("( (0 - 1K)/8K ), (0 - 1K)/8K )", NROUNDS, i, j, f3);
+    test_maj("( (0 - 8K)/8K ), (0 - 8K)/8K )", NROUNDS, j, j, f3);
+    printf("\n");
+}
+
+int main(int argc, char **argv)
+{
+    ullong i = 1000;
+    ullong j = 8000;
+    run_all_tests(2, i, j);
+    run_all_tests(4, i, j);
+    run_all_tests(6, i, j);
+    run_all_tests(8, i, j);
+    return 0;
+}