Benchmarking random float functions with box muller transformation.

rand_float_normal_bench.c

/*
    James William Fletcher (github.com/mrbid)
        August 2021

    Benchmarking random float functions using
    the box muller transformation.
    
    https://james-william-fletcher.medium.com/benchmarking-random-float-functions-in-c-13b9febb3d5b
*/
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <string.h>
#include <locale.h>
#include <time.h>
#include <sys/time.h>
#include <unistd.h>
#include <fcntl.h>
#include <x86intrin.h>

#pragma GCC diagnostic ignored "-Wunused-result"

#define TEST_SQRTPS 1
#define AVGITER 3000000
#define rand_min -20
#define rand_max 20

//////// RANDOM FLOATS

// classic random
float rand_float1()
{
    static const float rmax = (float)RAND_MAX;
    return (((float)rand())+1e-7f) / rmax;
}

float rand_float2()
{
    static const float FLOAT_UINT64_MAX = (float)UINT64_MAX;
    int f = open("/dev/urandom", O_RDONLY | O_CLOEXEC);
    uint64_t s = 0;
    ssize_t result = read(f, &s, sizeof(uint64_t));
    close(f);
    return (((float)s)+1e-7f) / FLOAT_UINT64_MAX;
}

// https://stackoverflow.com/questions/686353/random-float-number-generation
float rand_float3()
{
    uint_least32_t r = (rand() & 0xffff) + ((rand() & 0x00ff) << 16);
    return (float)r * 5.9604645E-8f;
}
float rand_float4()
{
    uint32_t pattern = 0x3f800000;
    uint32_t random23 = 0x7fffff & (rand() << 8 ^ rand());
    pattern |= random23;
    char buffer[sizeof(float)];
    memcpy(buffer, &pattern, sizeof(float));
    float f;
    memcpy(&f, buffer, sizeof(float));
    return f - 1.0f;
}

// adapted from ogre3d asm_math.h
// https://www.flipcode.com/archives/07-15-2002.shtml
// https://www.cs.cmu.edu/afs/andrew/scs/cs/oldfiles/15-494-sp09/dst/A/sw/ogre-1.6.4/OgreMain/include/asm_math.h
// https://gist.github.com/mrbid/9a050ee747a9188bc0aa849385bef865#file-rand_float_normal_bench-c-L63
float rand_float5(const __int64_t seed)
{
    static const float FLOAT_MAX = 9223372036854775807.0f;

    static __int64_t q = 8008135;
    if(seed != 0)
        q = seed;

    __m64 mm0 = _mm_cvtsi64_m64(q);
    __m64 mm1 = _m_pshufw(mm0, 0x1E);
    mm0 = _mm_add_pi32(mm0, mm1);
    q = _m_to_int64(mm0);

    _m_empty();
    return ( fabsf(q+1e-7f) / FLOAT_MAX );
}

//////// HELPERS
static inline float scale(const float normal, const float min, const float max)
{
    return ( normal * (max-min) ) + min;
}

static inline float rsqrtss(float f)
{
    return _mm_cvtss_f32(_mm_rsqrt_ss(_mm_set_ss(f)));
}

static inline float sqrtps(float f)
{
    return _mm_cvtss_f32(_mm_sqrt_ps(_mm_set_ss(f)));
}

void secure_random_srand()
{
    unsigned int s = 0;
    int f = open("/dev/urandom", O_RDONLY | O_CLOEXEC);
    read(f, &s, sizeof(unsigned int));
    close(f);
    srand(s);
}

uint64_t secure_random_seed()
{
    uint64_t s = 0;
    int f = open("/dev/urandom", O_RDONLY | O_CLOEXEC);
    read(f, &s, sizeof(uint64_t));
    close(f);
    return s;
}

//////// BOX MULLERS

float rand_normal_float1()
{
    float u = rand_float1() * 2.f - 1.f;
    float v = rand_float1() * 2.f - 1.f;
    float r = u * u + v * v;
    while(r == 0.f || r > 1.f)
    {
        u = rand_float1() * 2.f - 1.f;
        v = rand_float1() * 2.f - 1.f;
        r = u * u + v * v;
    }
    return u * sqrtf(-2.f * logf(r) / r);
}

float rand_normal_float2()
{
    float u = rand_float2() * 2.f - 1.f;
    float v = rand_float2() * 2.f - 1.f;
    float r = u * u + v * v;
    while(r == 0.f || r > 1.f)
    {
        u = rand_float2() * 2.f - 1.f;
        v = rand_float2() * 2.f - 1.f;
        r = u * u + v * v;
    }
    return u * sqrtf(-2.f * logf(r) / r);
}

float rand_normal_float3()
{
    float u = rand_float3() * 2.f - 1.f;
    float v = rand_float3() * 2.f - 1.f;
    float r = u * u + v * v;
    while(r == 0.f || r > 1.f)
    {
        u = rand_float3() * 2.f - 1.f;
        v = rand_float3() * 2.f - 1.f;
        r = u * u + v * v;
    }
    return u * sqrtf(-2.f * logf(r) / r);
}

float rand_normal_float4()
{
    float u = rand_float4() * 2.f - 1.f;
    float v = rand_float4() * 2.f - 1.f;
    float r = u * u + v * v;
    while(r == 0.f || r > 1.f)
    {
        u = rand_float4() * 2.f - 1.f;
        v = rand_float4() * 2.f - 1.f;
        r = u * u + v * v;
    }
    return u * sqrtf(-2.f * logf(r) / r);
}

float rand_normal_float5()
{
    float u = rand_float5(0) * 2.f - 1.f;
    float v = rand_float5(0) * 2.f - 1.f;
    float r = u * u + v * v;
    while(r == 0.f || r > 1.f)
    {
        u = rand_float5(0) * 2.f - 1.f;
        v = rand_float5(0) * 2.f - 1.f;
        r = u * u + v * v;
    }
    return u * sqrtf(-2.f * logf(r) / r);
}

//////// BOX MULLERS -- SQRTPS

float rand_normal_float1_sqrtps()
{
    float u = rand_float1() * 2.f - 1.f;
    float v = rand_float1() * 2.f - 1.f;
    float r = u * u + v * v;
    while(r == 0.f || r > 1.f)
    {
        u = rand_float1() * 2.f - 1.f;
        v = rand_float1() * 2.f - 1.f;
        r = u * u + v * v;
    }
    return u * sqrtps(-2.f * logf(r) / r);
}

float rand_normal_float2_sqrtps()
{
    float u = rand_float2() * 2.f - 1.f;
    float v = rand_float2() * 2.f - 1.f;
    float r = u * u + v * v;
    while(r == 0.f || r > 1.f)
    {
        u = rand_float2() * 2.f - 1.f;
        v = rand_float2() * 2.f - 1.f;
        r = u * u + v * v;
    }
    return u * sqrtps(-2.f * logf(r) / r);
}

float rand_normal_float3_sqrtps()
{
    float u = rand_float3() * 2.f - 1.f;
    float v = rand_float3() * 2.f - 1.f;
    float r = u * u + v * v;
    while(r == 0.f || r > 1.f)
    {
        u = rand_float3() * 2.f - 1.f;
        v = rand_float3() * 2.f - 1.f;
        r = u * u + v * v;
    }
    return u * sqrtps(-2.f * logf(r) / r);
}

float rand_normal_float4_sqrtps()
{
    float u = rand_float4() * 2.f - 1.f;
    float v = rand_float4() * 2.f - 1.f;
    float r = u * u + v * v;
    while(r == 0.f || r > 1.f)
    {
        u = rand_float4() * 2.f - 1.f;
        v = rand_float4() * 2.f - 1.f;
        r = u * u + v * v;
    }
    return u * sqrtps(-2.f * logf(r) / r);
}

float rand_normal_float5_sqrtps()
{
    float u = rand_float5(0) * 2.f - 1.f;
    float v = rand_float5(0) * 2.f - 1.f;
    float r = u * u + v * v;
    while(r == 0.f || r > 1.f)
    {
        u = rand_float5(0) * 2.f - 1.f;
        v = rand_float5(0) * 2.f - 1.f;
        r = u * u + v * v;
    }
    return u * sqrtps(-2.f * logf(r) / r);
}

int main()
{
    secure_random_srand();
    rand_float5(secure_random_seed());

    printf("\n");
    setlocale(LC_NUMERIC, "");
    float ret = 0.f;
    uint64_t st = 0, et = 0, avg = 0;

    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_normal_float1();
        avg += __rdtsc()-st;
    }
    printf("rand_normal_float1() Cycles: %'lu\n", avg / AVGITER);
    printf("0-1:  %.20f\n", rand_float1());
    printf("norm: %.20f\n\n", rand_normal_float1());

#if TEST_SQRTPS == 1
    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_normal_float1_sqrtps();
        avg += __rdtsc()-st;
    }
    printf("rand_normal_float1_sqrtps() Cycles: %'lu\n", avg / AVGITER);
    printf("0-1:  %.20f\n", rand_float1());
    printf("norm: %.20f\n\n", rand_normal_float1_sqrtps());

    printf("\n--------\n\n");
#endif
    ////////

    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_normal_float2();
        avg += __rdtsc()-st;
    }
    printf("rand_normal_float2() Cycles: %'lu\n", avg / AVGITER);
    printf("0-1:  %.20f\n", rand_float2());
    printf("norm: %.20f\n\n", rand_normal_float2());

#if TEST_SQRTPS == 1
    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_normal_float2_sqrtps();
        avg += __rdtsc()-st;
    }
    printf("rand_normal_float2_sqrtps() Cycles: %'lu\n", avg / AVGITER);
    printf("0-1:  %.20f\n", rand_float2());
    printf("norm: %.20f\n\n", rand_normal_float2_sqrtps());

    printf("\n--------\n\n");
#endif
    ////////

    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_normal_float3();
        avg += __rdtsc()-st;
    }
    printf("rand_normal_float3() Cycles: %'lu\n", avg / AVGITER);
    printf("0-1:  %.20f\n", rand_float3());
    printf("norm: %.20f\n\n", rand_normal_float3());

#if TEST_SQRTPS == 1
    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_normal_float3_sqrtps();
        avg += __rdtsc()-st;
    }
    printf("rand_normal_float3_sqrtps() Cycles: %'lu\n", avg / AVGITER);
    printf("0-1:  %.20f\n", rand_float3());
    printf("norm: %.20f\n\n", rand_normal_float3_sqrtps());

    printf("\n--------\n\n");
#endif
    ////////

    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_normal_float4();
        avg += __rdtsc()-st;
    }
    printf("rand_normal_float4() Cycles: %'lu\n", avg / AVGITER);
    printf("0-1:  %.20f\n", rand_float4());
    printf("norm: %.20f\n\n", rand_normal_float4());

#if TEST_SQRTPS == 1
    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_normal_float4_sqrtps();
        avg += __rdtsc()-st;
    }
    printf("rand_normal_float4_sqrtps() Cycles: %'lu\n", avg / AVGITER);
    printf("0-1:  %.20f\n", rand_float4());
    printf("norm: %.20f\n\n", rand_normal_float4_sqrtps());

    printf("\n--------\n\n");
#endif
    ////////

    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_normal_float5();
        avg += __rdtsc()-st;
    }
    printf("rand_normal_float5() Cycles: %'lu\n", avg / AVGITER);
    printf("0-1:  %.20f\n", rand_float5(0));
    printf("norm: %.20f\n\n", rand_normal_float5());

#if TEST_SQRTPS == 1
    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_normal_float5_sqrtps();
        avg += __rdtsc()-st;
    }
    printf("rand_normal_float5_sqrtps() Cycles: %'lu\n", avg / AVGITER);
    printf("0-1:  %.20f\n", rand_float5(0));
    printf("norm: %.20f\n\n", rand_normal_float5_sqrtps());
#endif

    ////////

    // done
    printf("%c\n", (char)ret);  // forces the compiler to not disregard the functions we are testing
    return 0;
}