// ------------------------------------------------------------------------------
// np_rnd.c 
// Functions to generate random numbers with uniform/gaussian probability distributions
// 
// (c) 2024 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy
// ------------------------------------------------------------------------------
//
//    This source is free software: you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation, either version 3 of the License, or
//    (at your option) any later version.
//    This file is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with this source distribution.  
//    If not, see <http://www.gnu.org/licenses/>.

#include "np_rnd.h"

#if _WIN32 // note the underscore: without it, it's not msdn official!
	// Windows (x64 and x86)
	#include <windows.h>   // required only for GetTickCount(...)
	#define K_RAND_MAX UINT_MAX
#elif _SVID_SOURCE || _XOPEN_SOURCE || __unix__ || (defined (__APPLE__) && defined(__MACH__)) /* POSIX or Unix or Apple */
	#include <stdlib.h>
	#define rand_s(x) (*x)=(unsigned int)lrand48()	// returns unsigned integers in the range 0..0x7FFFFFFF
	#define K_RAND_MAX 0x7FFFFFFF					// that's the max number
																				// generated by lrand48
#else
  #error "No good quality PRNG found"
#endif


void np_normrnd_real(float *dst, int nitems, float mean, float stdev)
{
	// for odd or even nitems
	
	unsigned int r;
	float phi=0, u=0;
	int set = 0;
	float scalephi = (M_2PI / (1.0f + K_RAND_MAX));
	float scaleu = (1.0f / (1.0f + K_RAND_MAX));


	while (nitems--) 
		if (set==1) {
			// generate a second normally distributed number
			*dst++ = sinf(phi)*u*stdev+mean;
			set = 0;
			}
		else {
			// generate a uniform distributed phase phi in the range [0,2*PI)
			rand_s((unsigned int*)&r); phi = scalephi * r;

			// generate a uniform distributed random number u in the range (0,1]
			rand_s((unsigned int*)&r); u = scaleu * (1.0f + r);

			// generate normally distributed number
			u = (float)sqrt(-2.0f * log(u));
			*dst++ = cosf(phi) * u * stdev + mean;

			set=1;
			}
}

void np_normrnd_cpx(float* dst, int nitems, float mean, float stdev)
{
	// as qpc_normrnd_real, but generates always nitems complex numbers (2*nitems reals)

	unsigned int r;
	float phi = 0, u = 0;
	// JHT	int set = 0;
	float scalephi = (M_2PI / (1.0f + K_RAND_MAX));
	float scaleu   = (1.0f / (1.0f + K_RAND_MAX));

	while (nitems--) {
		// generate a uniform distributed phase phi in the range [0,2*PI)
		rand_s((unsigned int*)&r); phi = scalephi * r;

		// generate a uniform distributed random number u in the range (0,1]
		rand_s((unsigned int*)&r); u = scaleu * (1.0f + r);
			
		// generate normally distributed real and imaginary parts
		u = (float)sqrt(-2.0f * log(u));
		*dst++ = cosf(phi) * u * stdev + mean;
		*dst++ = sinf(phi) * u * stdev + mean;
		}
}

void np_unidrnd(unsigned int* dst, int nitems, unsigned int nsetsize)
{
	// generate uniform unsigned int random numbers in the range [0..nsetsize)

	unsigned int r;
	float u;
	float scaleu = (1.0f / (1.0f + K_RAND_MAX));

	while (nitems--) {
		rand_s((unsigned int*)&r); u = scaleu * nsetsize * r;
		*dst++ = (unsigned int)floorf(u);
	}

}
void np_unidrnd_uc(unsigned char* dst, int nitems, unsigned char nsetsize)
{
	// generate uniform unsigned char random numbers in the range [0..nsetsize)

	unsigned int r;
	float u;
	float scaleu = (1.0f / (1.0f + K_RAND_MAX));

	while (nitems--) {
		rand_s((unsigned int*)&r); u = scaleu * nsetsize * r;
		*dst++ = (unsigned char)floorf(u);
	}

}

void np_unifrnd(float* dst, int nitems, float fmax)
{
	// generate uniform float random numbers in the range [0..fmax)

	unsigned int r;
	float u;
	float scaleu = (1.0f / (1.0f + K_RAND_MAX));

	while (nitems--) {
		rand_s((unsigned int*)&r); u = scaleu * fmax * r;
		*dst++ = u;
	}

}