#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include <time.h>

#define unlikely(x)  __builtin_expect (!!(x), 0)

uint32_t state[2] = { 0, 0 };

static const double dummy = (double)(1 << 20);

inline double my_rand(){
 uint32_t temp;
 uint64_t temp64;
 double *temp64_d = (double *) &temp64;
 uint64_t two_20 = * ( uint64_t * ) &dummy;

 if ( unlikely ( state[0] == 0 ) ) {
  state[0] = random();
  state[1] = random();
 }

  // Mix the bits.  Never replaces state[i] with 0 if it is nonzero.
  state[0] = 18273 * (state[0] & 0xFFFF) + (state[0] >> 16);
  state[1] = 36969 * (state[1] & 0xFFFF) + (state[1] >> 16);

  temp = (state[0] << 14) + (state[1] & 0x3FFFF);

  temp64 = temp;
  temp64 ^= two_20;
  return *temp64_d - dummy;
}

int trial() {
    int count = 0;
    double sum = 0;
    while (sum <= 1.0) {
        sum += rand() / (double) RAND_MAX;
        count++;
    }
    return count;
}

int trial2() {
    int count = 0;
    double sum = 0;
    while (sum <= 1.0) {
        sum += my_rand();
        count++;
    }
    return count;
}

double monteCarlo(int trial_type, int n) {
    int i, total = 0;
    if ( trial_type == 1 ) { // Use rand()
      printf("Using original rand() implementation.\n");
      for (i = 0; i < n; i++) {
          total += trial();
      }
    } else { // Use my_rand()
      printf("Using V8-style rand function.\n");
      for (i = 0; i < n; i++) {
          total += trial2();
      }
    }
    return total / (double) n;
}

int main(int argc, char *argv[]) {
    int i = 0;
    srandom(time(NULL));
    if ( argv[1][0] == '1' )
      printf("%f\n", monteCarlo(1,100000000));
    else
      printf("%f\n", monteCarlo(2,100000000));
    return 0;
}