/*

This is the C program for full testing MT19937, or   "MERSENNE TWISTER"
algorithm for pseudo random number generation, with initialization improved,
by  MAKOTO MATSUMOTO  and  TAKUJI NISHIMURA, of  2002/1/26.

This test program was made and tested by Pablo Mariano Ronchi (2005-Sep-12)



COMMENTS ABOUT THIS C PROGRAM FOR TESTING THE  "Mersenne Twister"  ALGORITHM:

- All the statements made by the authors of the original algorithm and implementation,
  including but not limiting to those regarding the license, the usage without
  any warranties, and the conditions for distribution, apply to this Visual Basic
  program for testing the algorithm.


- USAGE:

  In your Visual Basic application (MS Access, MS Excel, VB compiler, etc):

  1) Compile and run this program. 
     SEE THE COMMENT THAT YOU WILL FIND BELOW IN THIS SOURCE FILE, JUST ABOVE THE LINE:

	 #include "mt19937ar.h"

  2) If you have run the Visual Basic version of this test program, you can compare
     the output files manually following this procedure:
     a) Open both text files.
     b) Maximize their windows (from keyboard: <left Alt> + <Space bar>, and then <X>).
     c) Quickly alternate between both files typing:
        <left Alt> + <Tab>, <left Alt> + <Tab>, <left Alt> + <Tab>, ...
        Any difference between the files (time values, for example) will be
        self evident.

  Note: If you cannot find the output of this C program, it is probably in the
        same folder where is the executable file (.exe).


  HOW TO ADAPT THE TEST TO YOUR PARTICULAR NEEDS:

  Within MtArFullTest() you will find sets of constants for each part of the test.
  Just change their values. Of course, change also the same values in the C code to
  run the test in exactly the same way.
  For instance, for the first part of the test, you will find:

     Const kRun1stPart As Boolean = True
     Const kMaxNr1 As Long = 1000000        '1000000
     Const kPrintEachNth1 As Long = 10000   '10000

  In this case, the 1st part of the test will be run, a million numbers will be generated,
  but just a hundred (1000000/10000) of them will be printed into the result test file.

  Note:
     The actual generation and printing is performed by MtArTest(), which is called as
     many times as needed from within MtArFullTest()



  Pablo Mariano Ronchi
  Buenos Aires, Argentina



End of comments for the C program for testing mt19937ar
*/




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

/*
To include the original C code of the Mersenne Twister, you have 2 possibilities. 
Either: 

1)	To include "mt19937ar.h"
	This is the default option, as you can see below. As of this writing, this file 
	and the corresponding mt19937ar.c can be downloaded from the following link:

	http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.sep.tgz

	The file is very small (16Kb) and it usually downloads very quickly even thru
	a dial-up connection.

Or:
2)	To comment out the following "mt19937ar.h", and then copy the code for the 
	Mersenne Twister algorithm in this file. Then, DELETE the main() function
	provided with the original code, because it will be replaced by the one in
	this source file.
*/
#include "mt19937ar.h"





/* 
ADDITIONAL CODE, NOT PRESENT IN THE ORIGINAL C CODE OF THE MERSENNE TWISTER ALGORITHM, 
BUT ADDED TO ITS VISUAL BASIC VERSION.
*/


/* 
The value of the following constants MUST MATCH the equivalent constants in
the module "mt19937ar.c":
*/
#define k2_31			(2147483648)		/* 2^31   ==  2147483648 == 80000000		*/
#define k2_31Neg		(0 - k2_31)			/* -2^31  == -2147483648 == 80000000		*/
#define k2_31b			(k2_31 - 1)			/* 2^31-1 ==  2147483647 == 7FFFFFFF		*/
#define k2_32			(4294967296)		/* 2^32-1 ==  4294967295 == FFFFFFFF == -1	*/
#define k2_32b			(4294967295)		/* 2^32-1 ==  4294967295 == FFFFFFFF == -1	*/
#define kDefaultSeed	(5489)


/*
The value of the following constants MUST MATCH the equivalent constants in
the Visual Basic version:
*/
#define kMT_Gap			(0.0000000000005)	/* 5.0E-13			*/
#define kMT_Gap2		(2.0*kMT_Gap)		/* 1.0E-12			*/
#define kMT_GapInterval	(1.0-kMT_Gap2)		/* 0.9999999999990	*/
#define kMT_2b			(kMT_GapInterval/k2_32b)
#define kMT_2c			(kMT_2b)
#define kMT_3b			(kMT_2b)
#define kMT_4b			(2.0/k2_32b)
#define kMT_5b			((2.0*kMT_GapInterval)/k2_32b)		/* 1.999999999998/k2_32b */




signed long genrand_int32SignedLong() 
{
/* generates a random number on [0,0xffffffff]-interval
BUT RETURNS IT AS A (signed long) in the range [-2^31, 2^31-1]
*/
double tmp=genrand_int32();

if (tmp>k2_31b)
	return (signed long)(tmp-k2_32);
else
    return (signed long)tmp;

}    /* genrand_int32SignedLong */



double genrand_real2b()
{
/*
Returns results in the range [0,1) == [0, 1-kMT_Gap2]
Its lowest value is : 0.0
Its highest value is: 0.9999999999990
*/
return (genrand_int32()*kMT_2b);
}		/* genrand_real2b */



double genrand_real2c()
{
/*
Returns results in the range (0,1] == [0+kMT_Gap2, 1.0]
Its lowest value is : 0.0000000000010  (1E-12)
Its highest value is: 1.0
*/
return (kMT_Gap2+(genrand_int32()*kMT_2c));	/* ==kMT_Gap2+genrand_real2b() */
}		/* genrand_real2c */


double genrand_real3b()
{
/*
Returns results in the range (0,1) == [0+kMT_Gap, 1-kMT_Gap]
Its lowest value is : 0.0000000000005  (5E-13)
Its highest value is: 0.9999999999995
*/
return (kMT_Gap+(genrand_int32()*kMT_3b));
}		/* genrand_real3b */


double genrand_real4b()
{
/*
Returns results in the range [-1,1] == [-1.0, 1.0]
Its lowest value is : -1.0
Its highest value is: 1.0
*/
return ((genrand_int32()*kMT_4b)-1.0);
}		/* genrand_real4b */


double genrand_real5b()
{
/*
Returns results in the range (-1,1) == [-kMT_GapInterval, kMT_GapInterval]
Its lowest value is : -0.9999999999990
Its highest value is: 0.9999999999990
*/
return (kMT_Gap2+((genrand_int32()*kMT_5b)-1.0));
}		/* genrand_real5b */



/* END OF ADDITIONAL CODE */





/* Definitions for this tester: */

#define kgenrand_int32				(1)
#define kgenrand_int31				(2)
#define kgenrand_real1				(3)
#define kgenrand_real2				(4)
#define kgenrand_real3				(5)
#define kgenrand_res53				(6)
#define kgenrand_int32SignedLong    (7)
#define kgenrand_real2b				(8)
#define kgenrand_real2c				(9)
#define kgenrand_real3b				(10)
#define kgenrand_real4b				(11)
#define kgenrand_real5b				(12)


#define MAXSTR	(256)
#define TRUE	(-1)
#define FALSE	(0)

#define TIMER() ((long)time(NULL))


double TestN, TestP;		/* for testing */
FILE *ff;







void MtArTest(unsigned long kMaxNr, unsigned long kPrintEachNth,
                     int fcod, int withtime)
{
unsigned long ii, cp;
int intfn;
double tmp;
char s[MAXSTR], fn[MAXSTR];
double sec1, sec2;

/* (before calling the TIMER() following code generates and prints the 1st number: */
TestN+=1;
TestP+=1;
cp   = 0;
sec1 = TIMER();

/* print the first number: */
switch (fcod) {
    case kgenrand_int32:
        tmp = genrand_int32();
        strcpy(fn,"genrand_int32()");
		break;
    case kgenrand_int31:
        tmp = genrand_int31();
        strcpy(fn,"genrand_int31()");
		break;
    case kgenrand_real1:
        tmp = genrand_real1();
        strcpy(fn,"genrand_real1()");
		break;
    case kgenrand_real2:
        tmp = genrand_real2();
        strcpy(fn,"genrand_real2()");
		break;
    case kgenrand_real3:
        tmp = genrand_real3();
        strcpy(fn,"genrand_real3()");
		break;
    case kgenrand_res53:
        tmp = genrand_res53();
        strcpy(fn,"genrand_res53()");
		break;

	/* Not defined in the original C code: */

    case kgenrand_int32SignedLong:
        tmp = genrand_int32SignedLong();
        strcpy(fn,"genrand_int32SignedLong()");
		break;
    case kgenrand_real2b:
        tmp = genrand_real2b();
        strcpy(fn,"genrand_real2b()");
		break;
    case kgenrand_real2c:
        tmp = genrand_real2c();
        strcpy(fn,"genrand_real2c()");
		break;
    case kgenrand_real3b:
        tmp = genrand_real3b();
        strcpy(fn,"genrand_real3b()");
		break;
    case kgenrand_real4b:
        tmp = genrand_real4b();
        strcpy(fn,"genrand_real4b()");
		break;
    case kgenrand_real5b:
        tmp = genrand_real5b();
        strcpy(fn,"genrand_real5b()");
		break;

	default: 
		tmp = 0.0; 
        strcpy(fn,"NO FUNCTION");
	}

fprintf(ff,"%d of %d outputs of %s:\n",((kMaxNr/kPrintEachNth)+1),kMaxNr,fn);
fprintf(ff,"\n");

intfn=(fcod==2 || fcod==1 || fcod==7);

if (intfn)
    sprintf(s,"%11.0f",tmp);
else
    sprintf(s,"%10.8f",tmp);

fprintf(ff,"%s%s <-- this is the first number; then each ",s,((intfn)?" ":""));
if (kPrintEachNth>1) sprintf(s,"%dth number",kPrintEachNth); else strcpy(s,"one");
fprintf(ff,"%s is printed:\n",s);


for (ii=2; ii<=kMaxNr; ii++) {

    switch (fcod) {
        case kgenrand_int32: tmp = genrand_int32(); break;
        case kgenrand_int31: tmp = genrand_int31(); break;
        case kgenrand_real1: tmp = genrand_real1(); break;
        case kgenrand_real2: tmp = genrand_real2(); break;
        case kgenrand_real3: tmp = genrand_real3(); break;
        case kgenrand_res53: tmp = genrand_res53(); break;

		/* Not defined in the original C code: */

        case kgenrand_int32SignedLong: tmp = genrand_int32SignedLong(); break;
        case kgenrand_real2b: tmp = genrand_real2b(); break;
        case kgenrand_real2c: tmp = genrand_real2c(); break;
        case kgenrand_real3b: tmp = genrand_real3b(); break;
        case kgenrand_real4b: tmp = genrand_real4b(); break;
        case kgenrand_real5b: tmp = genrand_real5b(); break;

		default: 
			 tmp = 0.0; 
	}

    if (ii % kPrintEachNth==0) {
        cp++;
        
		if (intfn)
			sprintf(s,"%11.0f",tmp);
		else
			sprintf(s,"%10.8f",tmp);
        
        fprintf(ff,s);
		TestP+=1;
        
        if (cp==5) {
            cp=0; 
            fprintf(ff," \n");
		}
        else
            fprintf(ff," ");
    
	}

}

sec2 = TIMER();
TestN+=(kMaxNr-1);


if (withtime) {
    fprintf(ff,"\n");
    fprintf(ff,"Elapsed time in seconds for generating %d numbers and printing %d of them: %.2f\n",
			   kMaxNr,((kMaxNr/kPrintEachNth)+1),(double)(sec2-sec1));
}

fprintf(ff,"\n");
fprintf(ff,"\n");
}		/* MtArTest */






void MtArFullTest()
{
/*
This test consists of many parts, each part generating kMaxNr numbers but printing
into the output file just each kPrintEachNth of them. The output file could be compared
with a similar VBA test also provided.
*/

/* The following constants are for testing: */
#define kRun1stPart  (TRUE)
#define kMaxNr1 (1000000)         /* 1000000 */
#define kPrintEachNth1 (10000)    /* 10000   */

#define kRun2ndPart  (TRUE)
/* #define kMaxNr2 (100000000)       /* 100000000==1e8 */
/* #define kMaxNr2 (10000000)        /* 10000000==1e7 */
#define kMaxNr2 (1000000)         /* 1000000==1e6 */
#define kPrintEachNth2 (0)        /* not used  */

#define kRun3rdPart  (TRUE)
#define kMaxNr3 (1000)            /* 1000 */
#define kPrintEachNth3 (100)      /* 100  */

#define kRun4thPart  (TRUE)
#define kMaxNr4 (1000)            /* 1000 */
#define kPrintEachNth4 (100)      /* 100  */





unsigned long init[4] = {0x123, 0x234, 0x345, 0x456}, length=4;
/* the following array differs from the previous one in just the last bit: */
unsigned long init2[4]= {0x123, 0x234, 0x345, 0x457};
unsigned long init3[4]= {0x0, 0x0, 0x0, 0x0};
unsigned long init4[4]= {0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff};

unsigned long ii, seed, aa, bb;
char s[MAXSTR];
long sec1, sec2, secA, secB;


ff=fopen("MT19937arCFullTest.txt","w");
TestN=0;
TestP=0;

strcpy(s,"==========================================================\n");
fprintf(ff,s);
fprintf(ff,"Full Testing of the Mersenne Twister functions (MT19937ar)\n");
fprintf(ff,"C version\n");
fprintf(ff,s);
fprintf(ff,"\n");
fprintf(ff,"\n");

secA=TIMER(); 


/*
First part of test: generation and printing:
===================
*/
if (kRun1stPart) {
    strcpy(s,"-------------------------------------------\n");
    fprintf(ff,s);
    fprintf(ff,"First part of test: generation and printing\n");
    fprintf(ff,s);
    fprintf(ff,"\n");
    
    init_by_array(init,length);
    sprintf(s,"SEED:  initialization by array of %d elements "
			  "(same as in the original C code)\n",length);
    fprintf(ff,s);
	MtArTest(kMaxNr1, kPrintEachNth1, kgenrand_int32, TRUE);
    
    fprintf(ff,"\n");
}


/*
Second part of test: time needed just for generation
====================
*/
if (kRun2ndPart) {
    strcpy(s,"----------------------------------------------------\n");
    fprintf(ff,s);
    fprintf(ff,"Second part of test: time needed just for generation\n");
    fprintf(ff,s);
    fprintf(ff,"\n");
    
    init_by_array(init,length);
    
	aa=genrand_int32();
    sec1 = TIMER();
    for (ii=0; ii<kMaxNr2; ii++) bb=genrand_int32();
    sec2 = TIMER();
	TestN+=kMaxNr2+1;
    
    fprintf(ff,"Elapsed time in seconds for generating %d numbers and printing 2 of them: %.2f\n",
			kMaxNr2+1,(double)(sec2-sec1));
	fprintf(ff,"Control numbers:   A=%u   B=%u   C=%u\n", aa,bb,(kMaxNr2+1));
    fprintf(ff,"     A should be equal to 1067595299\n");
	if (kMaxNr2==100000000)
		fprintf(ff,"     B should be equal to 2891345831 if C is 100000001");
	else if (kMaxNr2==10000000)
		fprintf(ff,"     B should be equal to 2982969678 if C is 10000001");
	else {
		fprintf(ff,"     B should be equal to 3661023188 if C is 1000001");
		if (kMaxNr2!=1000000) fprintf(ff,", or usually different otherwise");
	}
	fprintf(ff,"\n");
    fprintf(ff,"     (A is the first and B is the last of C pseudorandom numbers generated)\n");
	TestP+=2;	/* for A and B */
    
	fprintf(ff,"\n");
    fprintf(ff,"\n");
    fprintf(ff,"\n");
}


/*
Third part of test: Generate numbers with different seeds
===================
*/

if (kRun3rdPart) {
    strcpy(s,"--------------------------------------------------------------\n");
    fprintf(ff,s);
    fprintf(ff,"Third part of test: Generation of numbers with different seeds\n");
    fprintf(ff,s);
    fprintf(ff,"\n");
    
    /* do not call init_genrand() */
    
    seed=kDefaultSeed; init_genrand(seed); fprintf(ff,"SEED:  default = %d\n",seed);
	MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    seed=(unsigned long)(k2_31Neg); init_genrand(seed); fprintf(ff,"SEED:  -2^31 = %d\n",seed);
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    seed=(unsigned long)(k2_31Neg+1); init_genrand(seed); fprintf(ff,"SEED:  -2^31+1 = %d\n",seed);
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    seed=(unsigned long)(-2); init_genrand(seed); fprintf(ff,"SEED:  %d\n",seed);
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    seed=(unsigned long)(-1); init_genrand(seed); fprintf(ff,"SEED:  %d\n",seed);
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    seed=0; init_genrand(seed); fprintf(ff,"SEED:  %d\n",seed);
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    seed=1; init_genrand(seed); fprintf(ff,"SEED:  %d\n",seed);
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    seed=2; init_genrand(seed); fprintf(ff,"SEED:  %d\n",seed);
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    seed=k2_31b - 1; init_genrand(seed); fprintf(ff,"SEED:  2^31-2 = %d\n",seed);
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    seed = k2_31b; init_genrand(seed); fprintf(ff,"SEED:  2^31-1 = %d\n",seed);
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);

    init_by_array(init,length);
    fprintf(ff,"SEED:  initialization by array of %d elements, equal to { 0x123, 0x234, 0x345, 0x456 } .\n",length);
    fprintf(ff,"\n");
    fprintf(ff,"       IT IS THE SAME TEST AS THE ONE IN THE ORIGINAL C CODE, \n");
    fprintf(ff,"       but prints only ");
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    init_by_array(init2,length);
    fprintf(ff,"SEED:  initialization by array of %d elements, equal to { 0x123, 0x234, 0x345, 0x457 } .\n",length);
    fprintf(ff,"\n");
    fprintf(ff,"       The initialization array differs from the one in "
			   "the original C code, \n");
    fprintf(ff,"       in JUST THE LAST BIT. Prints only ");
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);

    init_by_array(init3,length);
    fprintf(ff,"SEED:  initialization by array of %d elements,\n",length);
    fprintf(ff,"       each one equal to 0 .\n");
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    init_by_array(init4,length);
    fprintf(ff,"SEED:  initialization by array of %d elements,\n",length);
    fprintf(ff,"       each one equal to 0xFFFFFFFF == 2^32-1 == 4294967295 .\n");
    MtArTest(kMaxNr3, kPrintEachNth3, kgenrand_int32, FALSE);
    
    fprintf(ff,"\n");
}


/*
Fourth part of test: Generate real numbers
====================
*/

if (kRun4thPart) {
    strcpy(s,"-------------------------------------------------------\n");
    fprintf(ff,s);
    fprintf(ff,"Fourth part of test: Generation using all the functions\n");
    fprintf(ff,s);
    fprintf(ff,"\n");
    
    fprintf(ff,"\n");
    fprintf(ff,"** FUNCTIONS THAT RETURN AN INTEGER VALUE:\n");
    fprintf(ff,"\n");
    fprintf(ff,"\n");
    
    seed=kDefaultSeed; init_genrand(seed); 
	fprintf(ff,"SEED:  default = %d\n",seed);
    fprintf(ff,"\n");
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_int32, FALSE);
    
    seed=kDefaultSeed; init_genrand(seed); 
	fprintf(ff,"SEED:  default = %d (same seed as before, to compare the "
			   "different values returned)\n",seed);
    fprintf(ff,"\n");
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_int31, FALSE);
    
    
    fprintf(ff,"\n");
    fprintf(ff,"** FUNCTIONS THAT RETURN A REAL VALUE:\n");
    fprintf(ff,"\n");
    fprintf(ff,"\n");
    
    fprintf(ff,"SEED:  NO NEW SEEDS ARE USED in the remaining tests of this 4th section. \n");
    fprintf(ff,"       The last one (default=%d) remains valid and is not reset "
			   "between consecutive tests.\n",seed);
    fprintf(ff,"\n");
    fprintf(ff,"\n");
    
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_real1, FALSE);
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_real2, FALSE);
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_real3, FALSE);
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_res53, FALSE);
    
    fprintf(ff,"\n");
    fprintf(ff,"** FUNCTIONS THAT ARE PRESENT ONLY IN THE VISUAL BASIC VERSION \n");
    fprintf(ff,"   WITH NO COUNTERPART IN THE ORIGINAL C VERSION:\n");
    fprintf(ff,"\n");
    fprintf(ff,"\n");
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_int32SignedLong, FALSE);
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_real2b, FALSE);
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_real2c, FALSE);
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_real3b, FALSE);
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_real4b, FALSE);
    MtArTest(kMaxNr4, kPrintEachNth4, kgenrand_real5b, FALSE);

    fprintf(ff,"\n");
}

secB=TIMER();


strcpy(s,"==========================================================\n");
fprintf(ff,s);
fprintf(ff,"Full Testing of the Mersenne Twister functions (MT19937ar)\n");
fprintf(ff,"C version\n");
fprintf(ff,"\n");
fprintf(ff,"TOTAL AMOUNT of pseudorandom numbers generated\n");
fprintf(ff,"                 in all the previous tests: %.0f\n",TestN);
fprintf(ff,"                                   printed: %.0f\n",TestP);
fprintf(ff,"ELAPSED TIME in seconds for the full test : %.2f\n",(double)(secB-secA));
fprintf(ff,s);

fclose(ff);    /* close the output file */

}		/* MtArFullTest */




int main(void)
{
MtArFullTest();
return 0;
}		/* main */
