Random circle point benchmark

circle.c

// Random point in a circle benchmark
// Ref: https://youtu.be/4y_nmpv-9lI
// This is free and unencumbered software released into the public domain.
#define _POSIX_C_SOURCE 200112L
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <time.h>

static uint64_t
xoshiro256ss(uint64_t s[4])
{
    uint64_t x = s[1] * 5;
    uint64_t r = (x<<7 | x>>57) * 9;
    uint64_t t = s[1] << 17;
    s[2] ^= s[0]; s[3] ^= s[1];
    s[1] ^= s[2]; s[0] ^= s[3];
    s[2] ^= t;
    s[3] = s[3]<<45 | s[3]>>19;
    return r;
}

static double
uniform(uint64_t s[4])
{
    return xoshiro256ss(s) / 1.8446744073709552e19;
}

static void
rejection(uint64_t s[4], double c[2])
{
    do {
        c[0] = uniform(s)*2 - 1;
        c[1] = uniform(s)*2 - 1;
    } while (c[0]*c[0] + c[1]*c[1] > 1);
}

static void
sqrt_dist(uint64_t s[4], double c[2])
{
    double t = 2.0 * 3.141592653589793 * uniform(s);
    double r = sqrt(uniform(s));
    c[0] = r * cos(t);
    c[1] = r * sin(t);
}

static void
sum_dist(uint64_t s[4], double c[2])
{
    double t = 2.0 * 3.141592653589793 * uniform(s);
    double r = uniform(s) + uniform(s);
    if (r >= 1.0) {
        r = 2.0 - r;
    }
    c[0] = r * cos(t);
    c[1] = r * sin(t);
}

static void
max_dist(uint64_t s[4], double c[2])
{
    double t = 2.0 * 3.141592653589793 * uniform(s);
    double r = uniform(s);
    double x = uniform(s);
    if (x > r) {
        r = x;
    }
    c[0] = r * cos(t);
    c[1] = r * sin(t);
}

static double
now(void)
{
    struct timespec tp;
    clock_gettime(CLOCK_MONOTONIC, &tp);
    return tp.tv_sec + tp.tv_nsec/1e9;
}

int
main(void)
{
    volatile double sink;
    double start, c[2], sum[2];
    uint64_t s[] = {
        0x3243f6a8885a308d, 0x313198a2e0370734,
        0x4a4093822299f31d, 0x0082efa98ec4e6c8,
    };

    #define N (1L << 25)
    #define XSTR(s) STR(s)
    #define STR(s) #s
    #define BENCH(f) do { \
        start = now(); \
        sum[0] = sum[1] = 0.0; \
        for (long i = 0; i < N; i++) { \
            f(s, c); \
            sum[0] += c[0]; \
            sum[1] += c[1]; \
        } \
        sink = sum[0] + sum[1]; \
        printf("%-12s%.3f M-ops/s\n", STR(f), N / (now() - start) / 1e6); \
    } while (0)

    BENCH(rejection);
    BENCH(sqrt_dist);
    BENCH(sum_dist);
    BENCH(max_dist);
    (void)sink;
}