/*==========================================================================*/
/*                                                                          */
/*       Test of the random number routines to generate uniformly           */
/*       distributed random numbers                                         */
/*                                                                          */
/*           Three types of random generators can be called                 */
/*                -  modulo generator                                       */
/*                -  lagged Fibonacci                                       */
/*                -  shift register generator                               */
/*                                                                          */
/*       PROGRAM NAME:      random                                          */
/*                                                                          */
/*       CALLING SEQUENCE:                                                  */
/*                                                                          */
/*         Main program without command line arguments.                     */
/*                                                                          */
/*       RETURN:                                                            */
/*                                                                          */
/*         Exit code    -1  Unable to allocate space in main memory for     */
/*                          desired amount of random numbers                */
/*                                                                          */
/*   Version 1.00                                                           */
/*                                                                          */
/*   (c) Satoru Ozawa and Dieter W. Heermann                                */
/*                                                                          */
/*==========================================================================*/

/*--------------------------------------------------------------------------*/
/*                                                                          */
/*  Required files for this program:                                        */
/*                                                                          */
/*        mod.c           Basic random number generator for uniform pseudo- */
/*                        random numbers in the interval (0,1).             */
/*        fibgen.c        Fibonacci random number generator to generate     */
/*                        uniformly distributed random numbers in the       */
/*                        interval [0,1]                                    */
/*        r250.c          Basic random number generator for uniform pseudo- */
/*                        random numbers in the interval (0,1) using the    */
/*                        shift-register method                             */
/*                                                                          */
/*--------------------------------------------------------------------------*/

/* ==== Include files for this program ==== */

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

# include "sgl.h"               /* include for the graphics library         */
# include "mod.c"               /* modulo random number generator           */
# include "r250.c"              /* shift-register random number generator   */
# include "fibgen.c"            /* Fibonacci random number generator        */

/*--------------------------------------------------------------------------*/
/*                         Main program                                     */
/*--------------------------------------------------------------------------*/
#ifdef  __STDC__
void    main(int argc, char *argv[])
#else
void    main(argc,argv)
        int     argc;
        char    *argv[];
#endif
 {
    /*---------------------------------------------------------- */
    /*                     Declarations                          */
    /*---------------------------------------------------------- */

    /* ====  First the areas which we need  ==== */

    float     *x;               /* will take the produced random numbers */

    /* ==== dcl's in general ==== */

    int    i, sweep, max_sweeps;
    int    L;
    int    return_code;
    int    modul,multi,inc;
    int    seed;
    int    LagP, LagQ;
    int    mf[251];
    int    InputGenerator = 0;
    int    GeneratorType;

    float dx;
    float s1,s2,s3,s4;
    float us1,us2,us3,us4;
    float d1,d2,d3,d4;
    float umean,uvariance,ustandard_dev,uskew,uexcess;

    int FiboGenerator();
    int ModGenerator();

    /*---------------------------------------------------------- */
    /*                End of declares                            */
    /*---------------------------------------------------------- */

    /*---------------------------------------------------------- */
    /*      Start with setting up procedures:                    */
    /*      Read in the data for the test run                    */
    /*---------------------------------------------------------- */

    /* ====  Print the header for this test  ==== */

    printf("=======================================================\n");
    printf("======>  Test of the random number generator   <=======\n");
    printf("=======================================================\n");
    printf("\n\n");
    printf("(c) Satoru Ozawa and Dieter W. Heermann \n");
    printf("Version 1.00 (1994) \n");
    printf("\n\n");

    printf("The following random number generators can be tested:\n");
    while ( !InputGenerator ) {
        printf("Please input  1 for the modulo generator\n");
        printf("              2 for the Fibonacci generator\n");
        printf("              3 for the R250 generator\n");
        printf("GeneratorType :");
        scanf("%d",&GeneratorType);
        if ( (GeneratorType > 0) && (GeneratorType < 4)) InputGenerator = 1;
    }

    /* ==== Now get the parameters needed for the type of generator ==== */
    switch ( GeneratorType ) {
        case 1 :   /* Modulo Generator */
                   printf("Modulus (M)    :");
                   scanf("%d",&modul);
                   printf("Multiplier (a) :");
                   scanf("%d",&multi);
                   printf("Increment (b)  :");
                   scanf("%d",&inc);
                   break;
        case 2 :   /* Fibonacci Generator */
                   printf("Lag p          :");
                   scanf("%d",&LagP);
                   printf("Lag q          :");
                   scanf("%d",&LagQ);
                   break;
        default:
                break;
    }

    printf("Seed                     :");
    scanf("%d",&seed);
    printf("Sample size for this run :");
    scanf("%d",&max_sweeps);

    printf("\n\n");

    /* ==== Grep the storage for the random numbers  ==== */

    x = (float *) calloc(max_sweeps, sizeof(float));
    if ( x == NULL ) {
        printf("Out of memory: exit code -1\n");
        printf("Resize your request\n");
        exit(-1);
    }

    /*---------------------------------------------------------- */
    /*                      T e s t                              */
    /*---------------------------------------------------------- */

    printf("sampling uniform distribution ... ");
    switch ( GeneratorType ) {
        case 1 : /* Modulo Generator */
                 return_code = ModGenerator(modul,multi,inc,seed,max_sweeps,x);
                 break;
        case 2 : /* Fibonacci Generator */
                 return_code = FiboGenerator(LagP,LagQ,seed,max_sweeps,x);
        case 3 : /* R250 */
                 return_code = init_r250(seed,mf);
                 return_code = r250( max_sweeps,x,mf);
        default:
                break;
    }
    printf("done %d\n", return_code);

    printf("computing the moment of the distribution ... ");
    us1 = 0;
    us2 = 0;
    us3 = 0;
    us4 = 0;

    printf("loop 1");
    for(sweep=0; sweep < max_sweeps; sweep++) {
        us1  += x[sweep];
    }
    umean = us1 / max_sweeps;
    printf(" done\n");

    printf("loop 2 ...");

    for(sweep=0; sweep < max_sweeps; sweep++) {
        d1  = x[sweep] - umean;
        d2  = d1 * d1;
        d3  = d2 * d1;
        d4  = d2 * d2;

        us2 += d2;
        us3 += d3;
        us4 += d4;
    }

    printf("done\n");
    us2 /= max_sweeps;
    us3 /= max_sweeps;
    us4 /= max_sweeps;

    uvariance     = us2;
    ustandard_dev = sqrt(us2);
    uskew         = us3 / sqrt ( us2 * us2 * us2 );
    uexcess       = us4 / ( us2 * us2) - 3.0;
    printf("done\n");

    printf("visual analysis of the distribution\n");
    printf("type in any character to continue after you see ->finish-<\n\n");

    L = (int) sqrt( (float) max_sweeps / 2.0);
    dx = 1.0 / (float) L;
    L *= 10;
    if ( L > 700 ) {
         L = 700;
         dx = 0.05;
    }
    sgl_winsize(L,L);
    sgl_init();
    for(i=0;i<max_sweeps-1;i=i+2) {
       sgl_rectarea(x[i],x[i+1],dx,dx);
    }
    sgl_display();

    printf("->finish-<");
    i = getchar();
    i = getchar();

    printf("\n");
    printf("Results for Sampling a Uniform Distribution:\n");
    printf("-----------------------------------------------------------\n");
    printf("mean value           : %f \n", umean);
    printf("standard deviation   : %f \n", ustandard_dev);
    printf("variance             : %f \n", uvariance);
    printf("sample skew          : %f \n", uskew);
    printf("sample excess        : %f \n", uexcess);
    printf("\n");
    printf("Second moment        : %f \n", us2);
    printf("Third moment         : %f \n", us3);
    printf("Fourth moment        : %f \n", us4);
    printf("\n");

    /*---------------------------------------------------------- */
    /*               F I N I S H E D                             */
    /*---------------------------------------------------------- */

    /* ==== Clean up === */

    printf("finished\n");

    free(x);

    exit(0);    /* Normal return */
}