//
// Cforall Version 1.0.0 Copyright (C) 2021 University of Waterloo
//
// PRNG.c -- high-perforamnce pseudo-random numbers
//
// The contents of this file are covered under the licence agreement in the
// file "LICENCE" distributed with Cforall.
//
// Author           : Peter A. Buhr
// Created On       : Wed Dec 29 09:38:12 2021
// Last Modified By : Peter A. Buhr
// Last Modified On : Tue Dec  5 08:14:57 2023
// Update Count     : 428
// 

#include <fstream.hfa>									// sout
#include <stdlib.hfa>									// PRNG
#include <clock.hfa>
#include <limits.hfa>									// MAX
#include <math.hfa>										// sqrt
#include <malloc.h>										// malloc_stats
#include <locale.h>										// setlocale
#include <thread.hfa>
#include <mutex_stmt.hfa>

#define xstr(s) str(s)
#define str(s) #s

//#define TIME

#ifdef TIME												// use -O2 -nodebug
#define STARTTIME start = timeHiRes()
#define ENDTIME( extra ) sout | wd(0,1, (timeHiRes() - start)`ms / 1000.) | extra "seconds"
enum { BUCKETS = 100_000, TRIALS = 100_000_000 };
#else
#define STARTTIME
#define ENDTIME( extra )
enum { BUCKETS = 100_000, TRIALS = 10_000_000 };
#endif // TIME

static void avgstd( size_t trials, size_t buckets[] ) {
	size_t min = MAX, max = 0;
	double sum = 0.0, diff;
	for ( i; BUCKETS ) {
		if ( buckets[i] < min ) min = buckets[i];
		if ( buckets[i] > max ) max = buckets[i];
		sum += buckets[i];
	} // for

	double avg = sum / BUCKETS;							// average
	sum = 0.0;
	for ( i; BUCKETS ) {								// sum squared differences from average
		diff = buckets[i] - avg;
		sum += diff * diff;
	} // for
	double std = sqrt( sum / BUCKETS );
	mutex( sout ) sout | "trials"  | trials | "buckets" | BUCKETS
		| "min" | min | "max" | max
		| "avg" | wd(0,1, avg) | "std" | wd(0,1, std) | "rstd" | wd(0,1, (avg == 0 ? 0.0 : std / avg * 100)) | "%";
} // avgstd


size_t seed = 1009;

thread T1 {};
void main( T1 & ) {
	size_t * buckets = calloc( BUCKETS );				// too big for task stack
	for ( TRIALS / 50 ) {
		buckets[rand() % BUCKETS] += 1;					// concurrent
	} // for
	avgstd( TRIALS / 50, buckets );
	free( buckets );
} // main

thread T2 {};
void main( T2 & ) {
	PRNG prng;
	if ( seed != 0 ) set_seed( prng, seed );
	size_t * buckets = calloc( BUCKETS );				// too big for task stack
	for ( TRIALS ) {
		buckets[prng( prng ) % BUCKETS] += 1;			// concurrent
	} // for
	avgstd( TRIALS, buckets );
	free( buckets );
} // main

thread T3 {};
void main( T3 & th ) {
	size_t * buckets = calloc( BUCKETS );				// too big for task stack
	for ( TRIALS / 5 ) {
		buckets[prng() % BUCKETS] += 1;					// concurrent
	} // for
	avgstd( TRIALS / 5, buckets );
	free( buckets );
} // main

thread T4 {};
void main( T4 & th ) {
	size_t * buckets = calloc( BUCKETS );				// too big for task stack
	for ( TRIALS ) {
		buckets[prng( th ) % BUCKETS] += 1;				// concurrent
	} // for
	avgstd( TRIALS, buckets );
	free( buckets );
} // main

// Compiler bug requires hiding declaration of th from the bucket access, otherwise the compiler thinks th is aliased
// and continually reloads it from memory, which doubles the cost.
static void dummy( thread$ & th ) __attribute__(( noinline ));
static void dummy( thread$ & th ) {
	size_t * buckets = (size_t *)calloc( BUCKETS, sizeof(size_t) ); // too big for task stack
	for ( size_t i = 0; i < TRIALS; i += 1 ) {
		buckets[prng( th ) % BUCKETS] += 1;				// sequential
	} // for
	avgstd( TRIALS, buckets );
	free( buckets );
} // dummy


int main() {
	// setlocale( LC_NUMERIC, getenv( "LANG" ) );		// causes leaked storage message

	// only works on the current pthread thread
	// locale_t loc = newlocale( LC_NUMERIC_MASK, getenv( "LANG" ), (locale_t)0p );
	// if ( loc == (locale_t)0p ) abort( "newlocale" );
	// uselocale( loc );

	enum { TASKS = 4 };
	Time start;

#ifdef TIME												// too slow for test and generates non-repeatable results
#if 1
	sout | "glib rand" | nl | nl;

	size_t rseed;
	if ( seed != 0 ) rseed = seed;
	else rseed = rdtscl();
	srand( rseed );

	sout | sepOff;
	sout | nl | wd(26, "rand()" ) | wd(12, "rand(5)") | wd(12, "rand(0,5)" );
	for ( 20 ) {
		sout | wd(26, rand()) | nonl;
		sout | wd(12, rand() % 5) | nonl;
		sout | wd(12, rand() % (5 - 0 + 1) + 0);
	} // for
	sout | sepOn;
	sout | "seed" | rseed;

	sout | nl | "Sequential";
	STARTTIME;
	{
		size_t * buckets = calloc( BUCKETS );			// too big for task stack
		for ( i; TRIALS / 5 ) {
			buckets[rand() % BUCKETS] += 1;				// sequential
		} // for
		avgstd( TRIALS / 5, buckets );
		free( buckets );
	}
	ENDTIME( " x 5 " );

	sout | nl | "Concurrent";
	STARTTIME;
	{
		processor p[TASKS - 1];							// already 1 processor
		{
			T1 t[TASKS];
		} // wait for threads to complete
	}
	ENDTIME( " x 50 " );
#endif // 0
#endif // TIME

	sout | nl | "CFA " xstr(PRNG_NAME);

#if 1
	PRNG prng;

	if ( seed != 0 ) set_seed( prng, seed );

	sout | sepOff;
	sout | nl | wd(26, "PRNG()" ) | wd(12, "PRNG(5)") | wd(12, "PRNG(0,5)" );
	for ( 20 ) {
		sout | wd(26, prng( prng )) | nonl;				// cascading => side-effect functions called in arbitary order
		sout | wd(12, prng( prng, 5 )) | nonl;
		sout | wd(12, prng( prng, 0, 5 ));
	} // for
	sout | sepOn;
	sout | "seed" | get_seed( prng );

	sout | nl | "Sequential";
	STARTTIME;
	{
		size_t * buckets = calloc( BUCKETS );			// too big for task stack
		for ( TRIALS ) {
			buckets[prng( prng ) % BUCKETS] += 1;		// sequential
		} // for
		avgstd( TRIALS, buckets );
		free( buckets );
	}
	ENDTIME();

	sout | nl | "Concurrent";
	STARTTIME;
	{
		processor p[TASKS - 1];							// already 1 processor
		{
			T2 t[TASKS];
		} // wait for threads to complete
	}
	ENDTIME();
#endif // 0
#if 1
	if ( seed != 0 ) set_seed( seed );

	sout | sepOff;
	sout | nl | wd(26, "prng()" ) | wd(12, "prng(5)") | wd(12, "prng(0,5)" );
	for ( 20 ) {
		sout | wd(26, prng()) | nonl;					// cascading => side-effect functions called in arbitary order
		sout | wd(12, prng( 5 )) | nonl;
		sout | wd(12, prng( 0, 5 ));
	} // for
	sout | sepOn;
	sout | "seed" | get_seed( prng );

	sout | nl | "Sequential";
	STARTTIME;
	{
		size_t * buckets = calloc( BUCKETS );			// too big for task stack
		for ( TRIALS / 5 ) {
			buckets[prng() % BUCKETS] += 1;
		} // for
		avgstd( TRIALS / 5, buckets );
		free( buckets );
	}
	ENDTIME( " x 5 " );

	sout | nl | "Concurrent";
	STARTTIME;
	{
		processor p[TASKS - 1];							// already 1 processor
		{
			T3 t[TASKS];
		} // wait for threads to complete
	}
	ENDTIME( " x 5 " );
#endif // 0
#if 1
	if ( seed != 0 ) set_seed( seed );
	thread$ & th = *active_thread();

	sout | sepOff;
	sout | nl | wd(26, "prng(t)" ) | wd(12, "prng(t,5)") | wd(12, "prng(t,0,5)" );
	for ( 20 ) {
		sout | wd(26, prng( th )) | nonl;				// cascading => side-effect functions called in arbitary order
		sout | wd(12, prng( th, 5 )) | nonl;
		sout | wd(12, prng( th, 0, 5 ));
	} // for
	sout | sepOn;
	sout | "seed" | get_seed( prng );

	sout | nl | "Sequential";
	STARTTIME;
	{
		dummy( th );
	}
	ENDTIME();

	sout | nl | "Concurrent";
	STARTTIME;
	{
		processor p[TASKS - 1];							// already 1 processor
		{
			T4 t[TASKS];
		} // wait for threads to complete
	}
	ENDTIME();
#endif // 0
//	malloc_stats();
	// freelocale( loc );
} // main


// Local Variables: //
// compile-command: "cfa -DTIME -O2 -nodebug PRNG.cfa" //
// End: //