Benchmarking random float functions. Extended.

random_float_bench_ext.c

/*
    James William Fletcher (github.com/mrbid)
        October 2022

    compile: gcc random_float_bench_2022.c -Ofast -lm -o main

    Benchmarking random float functions, an update to the original in August 2021:
    https://gist.github.com/mrbid/51ed2963c88981452a5f87a3b072f8fb
    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 interval 6000000
// #define AVGITER 30000000

#define interval 1000000
#define AVGITER 3000000

// secure random
float rand_float1()
{
    static const float RECIP_FLOAT_UINT64_MAX = 1.f/(float)UINT64_MAX;
    int f = open("/dev/urandom", O_RDONLY | O_CLOEXEC);
    uint64_t s = 0;
    read(f, &s, sizeof(uint64_t));
    close(f);
    return ((float)s) * RECIP_FLOAT_UINT64_MAX;
}

// classic random
float rand_float2()
{
    static const float rmax = 1.f/(float)RAND_MAX;
    return ((float)rand()) * rmax;
}

// https://stackoverflow.com/questions/686353/random-float-number-generation
float rand_float3() // normal range 0-1
{
    uint_least32_t r = (rand() & 0xffff) + ((rand() & 0x00ff) << 16);
    return (float)r * 5.9604645E-8f;
}
float rand_float4() // normal range 0-1
{
    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()
{
    static __int64_t q = 74235;
    __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 q * 1.084202172e-19F;
}

// https://www.musicdsp.org/en/latest/Other/273-fast-float-random-numbers.html
// [email protected]
float rand_float6()
{
    static int srandfq = 74235;
    srandfq *= 16807;
    return (float)(srandfq) * 4.6566129e-010f;
}

// https://iquilezles.org/articles/sfrand/
// Inigo Quilez ([email protected])
float rand_float7()
{
    static int seed = 74235;
    float res;
    seed *= 16807;
    *((unsigned int *) &res) = ( ((unsigned int)seed)>>9 ) | 0x40000000;
    return( res-3.0f );
}

///

uint64_t microtime()
{
    struct timeval tv;
    struct timezone tz;
    memset(&tz, 0, sizeof(struct timezone));
    gettimeofday(&tv, &tz);
    return 1000000 * tv.tv_sec + tv.tv_usec;
}

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);
}

int main()
{
    secure_random_srand();

    printf("\n");
    setlocale(LC_NUMERIC, "");
    float ret = 0;
    unsigned long e = 0, ui = 0;
    uint64_t st = 0, et = 0, avg = 0;

    ////////

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

    e = 0;
    st = microtime();
    while(microtime() - st <= interval)
    {
        ret += rand_float1(-20, 20);
        e++;
    }
    ui = interval / 1000000;
    printf("Executions in %'lu seconds: %'lu\n", ui, e);
    printf("Executions per millisecond: %'lu\n", e/(1000*ui));
    printf("Executions per microsecond: %'lu\n", e/(1000000*ui));
    printf("~%'.8f executions every nanosecond\n\n", (float)e/(1000000000*ui));

    ////////

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

    e = 0;
    st = microtime();
    while(microtime() - st <= interval)
    {
        ret += rand_float2(-20, 20);
        e++;
    }
    ui = interval / 1000000;
    printf("Executions in %'lu seconds: %'lu\n", ui, e);
    printf("Executions per millisecond: %'lu\n", e/(1000*ui));
    printf("Executions per microsecond: %'lu\n", e/(1000000*ui));
    printf("~%'.8f executions every nanosecond\n\n", (float)e/(1000000000*ui));

    ////////

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

    e = 0;
    st = microtime();
    while(microtime() - st <= interval)
    {
        ret += rand_float3();
        e++;
    }
    ui = interval / 1000000;
    printf("Executions in %'lu seconds: %'lu\n", ui, e);
    printf("Executions per millisecond: %'lu\n", e/(1000*ui));
    printf("Executions per microsecond: %'lu\n", e/(1000000*ui));
    printf("~%'.8f executions every nanosecond\n\n", (float)e/(1000000000*ui));

    ////////

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

    e = 0;
    st = microtime();
    while(microtime() - st <= interval)
    {
        ret += rand_float4();
        e++;
    }
    ui = interval / 1000000;
    printf("Executions in %'lu seconds: %'lu\n", ui, e);
    printf("Executions per millisecond: %'lu\n", e/(1000*ui));
    printf("Executions per microsecond: %'lu\n", e/(1000000*ui));
    printf("~%'.8f executions every nanosecond\n\n", (float)e/(1000000000*ui));

    ////////

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

    e = 0;
    st = microtime();
    while(microtime() - st <= interval)
    {
        ret += rand_float5();
        e++;
    }
    ui = interval / 1000000;
    printf("Executions in %'lu seconds: %'lu\n", ui, e);
    printf("Executions per millisecond: %'lu\n", e/(1000*ui));
    printf("Executions per microsecond: %'lu\n", e/(1000000*ui));
    printf("~%'.8f executions every nanosecond\n\n", (float)e/(1000000000*ui));

    ////////

    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_float6();
        avg += __rdtsc()-st;
    }
    printf("rand_float6() AVG Cycles: %'lu\n", avg / AVGITER);
    printf("SAMPLE: %.3f\n", rand_float6());

    e = 0;
    st = microtime();
    while(microtime() - st <= interval)
    {
        ret += rand_float6();
        e++;
    }
    ui = interval / 1000000;
    printf("Executions in %'lu seconds: %'lu\n", ui, e);
    printf("Executions per millisecond: %'lu\n", e/(1000*ui));
    printf("Executions per microsecond: %'lu\n", e/(1000000*ui));
    printf("~%'.8f executions every nanosecond\n\n", (float)e/(1000000000*ui));

    ////////

    avg = 0;
    for(int i = 0; i < AVGITER; i++)
    {
        st = __rdtsc();
        ret += rand_float7();
        avg += __rdtsc()-st;
    }
    printf("rand_float7() AVG Cycles: %'lu\n", avg / AVGITER);
    printf("SAMPLE: %.3f\n", rand_float7());

    e = 0;
    st = microtime();
    while(microtime() - st <= interval)
    {
        ret += rand_float7();
        e++;
    }
    ui = interval / 1000000;
    printf("Executions in %'lu seconds: %'lu\n", ui, e);
    printf("Executions per millisecond: %'lu\n", e/(1000*ui));
    printf("Executions per microsecond: %'lu\n", e/(1000000*ui));
    printf("~%'.8f executions every nanosecond\n\n", (float)e/(1000000000*ui));

    ////////

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