Fast Real-Time Safe Pseudorandom Number Generator + Benchmark

prng_test.c

#include <stdint.h>
#include <stdio.h>
#include <time.h>
#include <math.h>

// xorshift code taken from wikipedia
struct xorshift32_state {
    uint32_t a;
};

/* The state must be initialized to non-zero */
uint32_t xorshift32(struct xorshift32_state *state)
{
    /* Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" */
    uint32_t x = state->a;
    x ^= x << 13;
    x ^= x >> 17;
    x ^= x << 5;
    return state->a = x;
}


float prng(struct xorshift32_state *state) {
    uint32_t x = xorshift32(state);
    float out = (float) x / 0xFFFFFFFF;
    return out;
}

float central_limit_normal(int iters, struct xorshift32_state *state) {
    float sum = 0;
    for (int i = 0; i < iters; i++) {
        sum += prng(state);
    }
    return sum / iters;
}

struct pair {
    float x;
    float y;
};

struct pair box_muller_normal(struct xorshift32_state *state) {
    float u1 = prng(state);
    float u2 = prng(state);
    float r = sqrt(-2 * log(u1));
    float theta = 2 * M_PI * u2;
    struct pair p;
    p.x = r * cos(theta);
    p.y = r * sin(theta);
    return p;
}

// from juce::FastMathApproximations
static float fast_cos (float x) {
    float x2 = x * x;
    float numerator = -(-39251520 + x2 * (18471600 + x2 * (-1075032 + 14615 * x2)));
    float denominator = 39251520 + x2 * (1154160 + x2 * (16632 + x2 * 127));
    return numerator / denominator;
}

struct pair approx_box_muller_normal(struct xorshift32_state *state) {
    float u1 = prng(state);
    float u2 = prng(state);
    float r = sqrt(-2 * log(u1));
    float theta = 2 * M_PI * u2;
    struct pair p;
    p.x = r * fast_cos(theta);
    p.y = r * fast_cos(theta + M_PI_2);
    return p;
}

int main() {
    struct xorshift32_state state;
    state.a = 12345678u;
    int iters = 10000000;
    float out;

    printf("\n\n running tests with %i iterations \n\n", iters);

    float startTime = (float)clock()/CLOCKS_PER_SEC;
    for (int i = 0; i < iters; ++i) {
        out = central_limit_normal(32, &state);
    }
    float endTime = (float)clock()/CLOCKS_PER_SEC;
    float timeElapsed = endTime - startTime;
    printf("Central Limit (30 iters) elapsed: %0.3f seconds\n", timeElapsed);

    struct pair out2;

    startTime = (float)clock()/CLOCKS_PER_SEC;
    for (int i = 0; i < iters / 2.0; ++i) {
        out2 = box_muller_normal(&state);
    }
    endTime = (float)clock()/CLOCKS_PER_SEC;
    timeElapsed = endTime - startTime;
    printf("Box Muller elapsed: %0.3f seconds\n", timeElapsed);

    startTime = (float)clock()/CLOCKS_PER_SEC;
    for (int i = 0; i < iters / 2.0; ++i) {
        out2 = approx_box_muller_normal(&state);
    }
    endTime = (float)clock()/CLOCKS_PER_SEC;
    timeElapsed = endTime - startTime;
    printf("Fast Approx Box Muller Time elapsed: %0.3f seconds\n", timeElapsed);


    struct xorshift32_state state2;
    state2.a = state.a;
    float sum = 0;
    // get standard deviation of box muller vs. approx
    for (int i = 0; i < iters; ++i) {
        float diff = box_muller_normal(&state).x - approx_box_muller_normal(&state2).x;
        sum += diff * diff;
    }
    float std_dev = sqrt(sum / iters);
    printf("Standard Deviation of approx: %f\n", std_dev);

    printf("\n\n\n");
}