readme.md

Sampling Uniformly in a Sphere

• Method 1: rejection sampling: sample uniformly in a cube, then reject points outside the unit sphere.

• Method 2: polar sampling: use polar coordinates to directly generate points uniformly in the sphere

There are two implementations here. One in JavaScript, one in C. Which method is faster is different depending on the language!

Results (JavaScript)

Using Node v14.16.1.

rejection x 14,358,214 ops/sec ±0.73% (91 runs sampled)
polar x 16,253,126 ops/sec ±0.63% (90 runs sampled)

Converted to timing: rejection is 69.6 ns/sample, polar is 61.5 ns/sample.

Rejection sampling takes 14% more time.

Results (C)

Using GCC 10.2.1, with -O3 -ffast-math.

rejection: 21.1 ns/sample, ±4.3 ns/sample
polar: 32.1 ns/sample, ±0.7 ns/sample

Rejection sampling takes 34% less time.

Analysis

The JavaScript version does everything at double precision, which may explain the differences.

samplec.c

#define _DEFAULT_SOURCE 1

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

struct random_state {
    uint64_t state;
    uint64_t inc;
};

static float random_next(struct random_state *st) {
    uint64_t s = st->state;
    st->state = s * 6364136223846793005ull + st->inc;
    uint32_t x = ((s >> 18) ^ s) >> 27;
    uint32_t r = s >> 59;
    uint32_t v = (x >> r) | (x << ((-r) & 31));
    return (float)v * (1.0f / 4294967296.0f);
}

// Generate a point in the unit sphere, uniformly.
void rejection_sample(struct random_state *st, int n, float *arr) {
    for (int i = 0; i < n;) {
        float x = random_next(st) * 2.0f - 1.0f;
        float y = random_next(st) * 2.0f - 1.0f;
        float z = random_next(st) * 2.0f - 1.0f;
        if (x * x + y * y + z * z <= 1) {
            arr[3 * i] = x;
            arr[3 * i + 1] = y;
            arr[3 * i + 2] = z;
            i++;
        }
    }
}

// Generate a point in the unit sphere, uniformly.
void polar_sample(struct random_state *st, int n, float *arr) {
    for (int i = 0; i < n; i++) {
        float z = random_next(st) * 2.0f - 1.0f;
        float angle = (8.0f * atanf(1.0f)) * random_next(st);
        float z_scale = cbrtf(random_next(st));
        float xy_scale = z_scale * sqrtf(1.0f - z * z);
        arr[3 * i] = cosf(angle) * xy_scale;
        arr[3 * i + 1] = cosf(angle) * xy_scale;
        arr[3 * i + 2] = z * z_scale;
    }
}

typedef void (*sample_func)(struct random_state *st, int n, float *arr);

void benchmark(const char *name, int runs, int iter, int n, sample_func func) {
    float *arr = malloc(n * 3 * sizeof(float));
    if (arr == NULL) {
        abort();
    }
    double *times = malloc(runs * sizeof(double));
    if (times == NULL) {
        abort();
    }
    struct random_state st = {.state = 1, .inc = 1};
    func(&st, n, arr);
    for (int i = 0; i < runs; i++) {
        struct timespec t0, t1;
        clock_gettime(CLOCK_MONOTONIC, &t0);
        for (int j = 0; j < iter; j++) {
            func(&st, n, arr);
        }
        clock_gettime(CLOCK_MONOTONIC, &t1);
        times[i] = (double)(t1.tv_sec - t0.tv_sec) +
                   1.0e-9 * (double)(t1.tv_nsec - t0.tv_nsec);
    }
    double sum = 0.0;
    for (int i = 0; i < runs; i++) {
        sum += times[i];
    }
    double mean = sum / (double)runs;
    sum = 0.0;
    for (int i = 0; i < runs; i++) {
        double delta = times[i] - mean;
        sum += delta * delta;
    }
    double var = sum / (double)runs;
    printf("%s: %.1f ns/sample, ±%.1f ns/sample\n", name,
           mean * 1.0e9 / (double)(iter * n),
           sqrt(var) * 1.0e9 / (double)(iter * n));
    free(arr);
    free(times);
}

int main(int argc, char**argv) {
    (void)argc;
    (void)argv;
    benchmark("rejection", 1000, 100, 1000, rejection_sample);
    benchmark("polar", 1000, 100, 1000, polar_sample);
    return 0;
}

samplejs.js

const benchmark = require('benchmark');
const suite = new benchmark.Suite();

// Generate a point in the unit sphere, uniformly.
function rejectionSample() {
  while (true) {
    let x = Math.random() * 2 - 1;
    let y = Math.random() * 2 - 1;
    let z = Math.random() * 2 - 1;
    if (x * x + y * y + z * z <= 1) {
      return { x, y, z };
    }
  }
}

// Generate a point in the unit sphere, uniformly.
function polarSample() {
  let z = Math.random() * 2 - 1;
  let angle = 2 * Math.PI * Math.random();
  let zScale = Math.cbrt(Math.random());
  let xyScale = zScale * Math.sqrt(1 - z * z);
  return {
    x: Math.cos(angle) * xyScale,
    y: Math.sin(angle) * xyScale,
    z: z * zScale,
  };
}

suite.add('rejection', rejectionSample);
suite.add('polar', polarSample);
suite.on('cycle', event => console.log(String(event.target)));
suite.run();