/* -*- Mode: C; indent-tabs-mode: t; c-basic-offset: 4; tab-width: 4 -*-  */
/*
 * main.c
 * Copyright (C) 2018 TikalonLLC <gualtieri@ieee.org>
 * 
 * KISS 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.
 * 
 * KISS 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 program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
This program is based on the code contained in the article, "64-bit KISS RNGs,"
by George Marsaglia (Feb 28, 2009).  As explained in the article, "This 64-bit 
KISS RNG has three components, each nearly good enough to serve alone. The 
components are:
Multiply-With-Carry (MWC), period (2^121+2^63-1)
Xorshift (XSH), period 2^64-1
Congruential (CNG), period 2^64

see https://www.thecodingforums.com/threads/64-bit-kiss-rngs.673657/
*/

#include <stdio.h>
#include <string.h>
#include <math.h>

#define MWC (t=(x<<58)+c, c=(x>>6), x+=t, c+=(x<t), x)
#define XSH ( y^=(y<<13), y^=(y>>17), y^=(y<<43) )
#define CNG ( z=6906969069LL*z+1234567 )
#define KISS (MWC+XSH+CNG)
#define limit 100000000

int i,j;
char str[256]="";
long int bins[16];
double average,stdev;
static unsigned long long
x=1234567890987654321ULL,c=123456123456123456ULL,
y=362436362436362436ULL,z=1066149217761810ULL,t;

int main(void)
{
//initialize the digit counting bins
for(j=0;j<16;j++)
{
	bins[j]=0;
}

for(i=0;i<limit;i++)
{
//t, and its string representation, str, are the random numbers
//you could write them to a file, but this is not recommended
//when the limit is 100 million!
t=KISS;
sprintf(str, "%llX", t); //Coverts hex number into the equivalent string

for(j=0;j<strlen(str);j++)
{
	if(str[j]<58)
	{
		bins[str[j]-48]++;
	}
	else
	{
		bins[str[j]-55]++;
	}

}

//remember that some of these hex numbers are smaller than 16 digits
//need to adjust for some leading zeros

bins[0] = bins[0]+(16-strlen(str));

//show that the program is executing by printing an asterisk every million
//calculations

if((i%1000000) == 0)
{
printf("*");
//wont show unless we flush the buffer while running in a loop
fflush(stdout);
}

}

printf("\n");
//print results and get average
//the average should be the limit of the loop variable,
//but it's good to check for errors
average = 0;
for(j=0;j<16;j++)
{
	printf("[%X] = %ld\n",j,bins[j]);
	average = average+bins[j];
}

average = average/16;

//get standard deviation of results
stdev = 0;
for(j=0;j<16;j++)
{
	stdev= stdev + ((bins[j]-average)*(bins[j]-average));
}
stdev =stdev/15;
stdev = sqrt(stdev);

printf("\naverage = %f\tstdev = %f\tratio = %f\n",average,stdev,stdev/average);

printf("\nDone.\n");

return 0;
}