// --------------------------------------------------------------
// --------------------------------------------------------------

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

#include "noise.h"

// --------------------------------------------------------------
// internal funcs
// --------------------------------------------------------------
#define PM_IA 16807
#define PM_IM 2147483647
#define PM_AM (1.0/PM_IM)
#define PM_IQ 127773
#define PM_IR 2836
#define PM_MASK 123459876

//	'Minimal' random number generator of Park and Miller. Returns a uniform random deviate
//	between 0.0 and 1.0. Set or resetidumto any integer value (except the unlikely value MASK)
//	to initialize the sequence; idummust not be altered between calls for successive deviates in
//	a sequence.
double ran0(long *idum)
{
	long k;
	double ans;
	
	*idum ^= PM_MASK;				// XORing with MASKallows use of zero and other
	k=(*idum)/PM_IQ;				// simple bit patterns for idum.
	*idum=PM_IA*(*idum-k*PM_IQ)-PM_IR *k;		// Compute idum=(IA*idum) % IM without over-
							// flows by Schrage’s method.
	if (*idum < 0) 
		*idum += PM_IM;
		
	ans=PM_AM*(*idum);				// Convert idumto a floating result.
	*idum ^= PM_MASK;				// Unmask before return.
	
	return ans;
}

// --------------------------------------------------------------
// Exported functions
// --------------------------------------------------------------
void Noise_Init(noise_gen_t *pObj, long seed)
{
	pObj->state = seed;
}

void Noise_Free(noise_gen_t *pObj)
{
}

double Noise_Uniform(noise_gen_t *pObj, double mu)
{
	return ran0(&pObj->state) + mu -0.5;

}

double Noise_Gaussian(noise_gen_t *pObj, double mu, double sigma)
{
	double U1, U2, V1, V2, S, Y;

	do
	{
		U1 = Noise_Uniform(pObj, 0.5);                // U1=[0,1]
		U2 = Noise_Uniform(pObj, 0.5);                // U2=[0,1]
		V1 = 2 * U1 - 1;          // V1=[-1,1]
		V2 = 2 * U2 - 1;          // V2=[-1,1]
		S = V1 * V1 + V2 * V2;

	} while (S >= 1);

	// X = sqrt(-2 * log(S) / S) * V1;
	Y = sqrt(-2 * log(S) / S) * V2;
	Y = mu + sigma*Y;

	return Y;
}