/*

 
Use "REAL" as floating point type instead of double or float
 
Compile with optional flag:
    -DSINGLE to use float instead of double
*/

#include <stdio.h>
#include <stdlib.h> /* for rand() */
#include <unistd.h> /* for getpid() */
#include <time.h> /* for time() */
#include <math.h>
#include <assert.h>

// construct REAL "type," depending on desired precision

#ifdef SINGLE
#define REAL float
#else
#define REAL double
#endif

// ******************************************************
// your code here

REAL midsquare(unsigned int &seed) {
    unsigned int sq;
 
    sq=seed*seed;
    seed=(sq>>8)&0xFFFF;
    
    return 1.0*seed/(1.0*0xFFFF);
}


REAL getRN() {
    return drand48();
    //return (1.0*rand())/(1.0*RAND_MAX);
}

// ******************************************************

int main(int argc,char *argv[])
{
    int N,NC,NB,i,j;
    REAL x,y,sum,sqsum,mean,variance,sdev;
    int *hist;
    REAL dx;
    unsigned int seed;
    
    //welcome info
    
    printf("# Random number demo using ");
#ifdef SINGLE
    printf("single");
#else
    printf("double");
#endif
    printf(" precision arithmetics.\n");
    
    //default parameters
    
    N=100;  //number of generated numbers
    NB=100; //bin number
    NC=2;   //order of CLT test

    // check if arguments are present and read them
    
    if(argc > 1 ) N = atoi(argv[1]);
    if(argc > 2 ) NB = atoi(argv[2]);
    if(argc > 3 ) NC = atoi(argv[3]);
    


    //excute  calculations
    
    printf("# parameters: N=%d, NB=%d, NC=%d\n",N,NB,NC);
    
    //scatter plot

    seed=34583;
    
    for(i=0;i<N;i++) {
        x=midsquare(seed);
        y=midsquare(seed);
        printf("%le %le\n",x,y);
    }

    
    //mean and stdev
/*
    sum=0.0;
    sqsum=0.0;
    for(i=0;i<N;i++) {
        x=getRN();
        sum+=x;
        sqsum+=x*x;
        mean=sum/(1.0*(i+1));
        variance=sqsum/(1.0*(i+1))-mean*mean;
        sdev=sqrt(variance);
        printf("%d %le %le %le %le\n",i+1,x,mean,variance,sdev);
    }
 */
    
    //statistics
/*
    hist=new int[NB];
    for(i=0;i<NB;i++) hist[i]=0;
    dx=1.0/(1.0*NB);
    
    for(i=0;i<N;i++) {
        x=getRN();
        j=(int) (x/dx);if(j>=NB) j=NB-1;
        hist[j]++;
    }
    for(i=0;i<NB;i++) {
        x=(i+0.5)*dx;
        y=1.0*hist[i]/(1.0*N);
        printf("%le %le\n",x,y);
    }
    delete[] hist;
*/
    
    //CLT
/*
    hist=new int[NB];
    for(i=0;i<NB;i++) hist[i]=0;
    dx=1.0/(1.0*NB);
    
    for(i=0;i<N;i++) {
        x=0.0;
        for(j=0;j<NC;j++) x+=getRN();
        x=x/(1.0*NC);
        j=(int) (x/dx);if(j>=NB) j=NB-1;
        hist[j]++;
    }
    for(i=0;i<NB;i++) {
        x=(i+0.5)*dx;
        y=1.0*hist[i]/(1.0*N);
        printf("%le %le\n",x,y);
    }
    delete[] hist;
 */
 
 
 
    return 0;
    }

// ******************************************************