circle_sampling.c

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

/*
 * Test different methods of uniformly sampling a point on the unit circle.
 * Inspired by https://www.youtube.com/watch?v=4y_nmpv-9lI
 * Sample outputs:
 *
 * $ gcc circle_sampling.c -o circle_sampling -lm
 * $ ./circle_sampling
 * warmup (reflection): 85.358 ns/sample
 * rejection:           30.625 ns/sample
 * transformation:      67.108 ns/sample
 * reflection:          82.702 ns/sample
 * maximum:             78.505 ns/sample
 *
 * $ gcc circle_sampling.c -o circle_sampling -lm -O3
 * $ ./circle_sampling
 * warmup (reflection): 60.714 ns/sample
 * rejection:           25.570 ns/sample
 * transformation:      44.802 ns/sample
 * reflection:          57.651 ns/sample
 * maximum:             56.620 ns/sample
 *
 */

void rejection_sample(float* target)
{
	while (1) {
		float x = (((float) rand()) / RAND_MAX - .5)*2;
		float y = (((float) rand()) / RAND_MAX - .5)*2;
		if (x*x + y*y <= 1) {
			*target = x;
			*(target+1) = y;
			return;
		}
	}
}

void transformation_sample(float* target)
{
	float phi = (((float) rand()) / RAND_MAX)*2*M_PI;
	float r = sqrt(((float) rand()) / RAND_MAX);
	float x = r * cos(phi);
	float y = r * sin(phi);
	*target = x;
	*(target+1) = y;
	return;
}

void reflection_sample(float* target)
{
	float phi = (((float) rand()) / RAND_MAX)*2*M_PI;
	float r1 = (((float) rand()) / RAND_MAX);
	float r2 = (((float) rand()) / RAND_MAX);

	float r = r1 + r2;
	if (r > 1) r = 2 - r;

	float x = r * cos(phi);
	float y = r * sin(phi);
	*target = x;
	*(target+1) = y;
	return;
}

void max_sample(float* target)
{
	float phi = (((float) rand()) / RAND_MAX)*2*M_PI;
	float r1 = (((float) rand()) / RAND_MAX);
	float r2 = (((float) rand()) / RAND_MAX);

	float r = r1 > r2 ? r1 : r2;

	float x = r * cos(phi);
	float y = r * sin(phi);
	*target = x;
	*(target+1) = y;
	return;
}

int main(int argc, char* argv)
{
	srand(31415926);
	float* ptr = (float*) malloc(2 * sizeof(float));

	int REPS = 10000000;
	clock_t t;

	t = clock();
	for (int i=0; i<REPS; i++) {
		reflection_sample(ptr);
	}
	t = clock() - t;
	double time_taken = (1e9/REPS)*((double)t)/CLOCKS_PER_SEC;
	printf("warmup (reflection): %.3f ns/sample\n", time_taken);

	t = clock();
	for (int i=0; i<REPS; i++) {
		rejection_sample(ptr);
	}
	t = clock() - t;
	time_taken = (1e9/REPS)*((double)t)/CLOCKS_PER_SEC;
	printf("rejection:           %.3f ns/sample\n", time_taken);

	t = clock();
	for (int i=0; i<REPS; i++) {
		transformation_sample(ptr);
	}
	t = clock() - t;
	time_taken = (1e9/REPS)*((double)t)/CLOCKS_PER_SEC;
	printf("transformation:      %.3f ns/sample\n", time_taken);

	t = clock();
	for (int i=0; i<REPS; i++) {
		reflection_sample(ptr);
	}
	t = clock() - t;
	time_taken = (1e9/REPS)*((double)t)/CLOCKS_PER_SEC;
	printf("reflection:          %.3f ns/sample\n", time_taken);

	t = clock();
	for (int i=0; i<REPS; i++) {
		max_sample(ptr);
	}
	t = clock() - t;
	time_taken = (1e9/REPS)*((double)t)/CLOCKS_PER_SEC;
	printf("maximum:             %.3f ns/sample\n", time_taken);

	return 0;
}