#include <random>
#include <cmath>

const double pi = 3.14159265358979323846;

void rejection_sampling(int count)
{
	std::mt19937 rnd;
	std::uniform_real_distribution<double> dist(-1.0, 1.0);

	for(int i=0; i<count; ++i)
	{
		bool rejected = false;
		do
		{
			double x = dist(rnd);
			double y = dist(rnd);
			double z = dist(rnd);
			double r2 = x*x + y*y + z*z;
			if(r2 > 1.0 || r2 == 0)
				rejected = true;
			else
			{
				rejected = false;
				double r = sqrt(r2);
				printf("%f %f %f\n", x/r, y/r, z/r);
			}
		}
		while(rejected);
	}
}

void normal_deviate(int count)
{
	std::mt19937 rnd;
	std::normal_distribution<double> dist(0.0, 1.0);

	for(int i=0; i<count; ++i)
	{
		bool bad_luck = false;
		do
		{
			double x = dist(rnd);
			double y = dist(rnd);
			double z = dist(rnd);
			double r2 = x*x + y*y + z*z;
			if(r2 == 0)
				bad_luck = true;
			else
			{
				bad_luck = false;
				double r = sqrt(r2);
				printf("%f %f %f\n", x/r, y/r, z/r);
			}
		}
		while(bad_luck);
	}
}

void trig_method(int count)
{
	std::mt19937 rnd;
	std::uniform_real_distribution<double> dist(0.0, 1.0);

	for(int i=0; i<count; ++i)
	{
		double z = 2.0*dist(rnd) - 1.0;
		double t = 2.0*pi*dist(rnd);
		double r = sqrt(1.0-z*z);
		printf("%f %f %f\n", r*cos(t), r*sin(t), z);
	}
}

void coordinate_approach(int count)
{
	std::mt19937 rnd;
	std::uniform_real_distribution<double> dist(-1.0, 1.0);

	for(int i=0; i<count; ++i)
	{
		bool rejected = false;
		do
		{
			double u = dist(rnd);
			double v = dist(rnd);
			double s = u*u + v*v;
			if(s > 1.0)
				rejected = true;
			else
			{
				rejected = false;
				double a = 2.0*sqrt(1.0-s);
				printf("%f %f %f\n", a*u, a*v, 2.0*s-1.0);
			}
		}
		while(rejected);
	}
}

int main(int argc, char* argv[])
{
	printf("[rejection sampling]\n\n");
	rejection_sampling(500);
	
	printf("\n[normal deviate]\n\n");
	normal_deviate(500);
	
	printf("\n[trig method]\n\n");
	trig_method(500);
	
	printf("\n[coordinate approach]\n\n");
	coordinate_approach(500);

	return 0;
}