/*

 
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);
}

//"fill" a quarter unit circle
REAL pimc(int N) {
    REAL pi,x,y;
    int i,count;
    unsigned int seed=23571;
    
    count=0;
    for(i=0;i<N;i++) {
        x=midsquare(seed); //getRN();
        y=midsquare(seed); //getRN();
        if(x*x+y*y<=1.0) count++;
    }
    pi=(1.0*count)/(1.0*N);
    return 4.0*pi;
}



//uses lim_{n->inf} (1-1/n)^n=1/e
//+ hit n targets randomly n times and count the empty ones
REAL emc(int N) {
    unsigned char *hits;
    int i,j,count;
    
    hits=new unsigned char[N];
    for(i=0;i<N;i++) hits[i]=0;
    for(i=0;i<N;i++) {
        j=(int) (N*getRN());if(j>=N) j=N-1;
        hits[j]=1;
    }
    count=0;
    for(i=0;i<N;i++) if(hits[i]==0) count++;
    delete[] hits;
    
    return (1.0*N)/(1.0*count);
}

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

int main(int argc,char *argv[])
{
    int N;
    REAL pi,euler;
    
    
    //welcome info
    
    printf("Monte Carlo pi & e calculation using ");
#ifdef SINGLE
    printf("single");
#else
    printf("double");
#endif
    printf(" precision arithmetics.\n");
    
    //default parameters
    
    N=100;

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


    //excute Monte Carlo calculations
    
    pi=pimc(N);
    euler=emc(N);
    
    printf("result for N=%d: pi=%.16le, e=%.16le\n",N,pi,euler);
    
    
    return 0;
    }

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