/////////////////////////////////////////////////////////////////////////////////
//                                                                             //
// Interface to different RANDOM NUMBER GENERATORS                             //
// Make generators interchangeable                                             //                 
//                                                                             //
// Burkhard Militzer                                    Urbana 1998            //
//                                                                             //
/////////////////////////////////////////////////////////////////////////////////

#ifndef _RANDOM_
#define _RANDOM_

#include "Standard.h"

#include <stdlib.h>
#include <unistd.h> 
#include <sys/times.h>
#include <sys/time.h>

inline void InitRandom(const bool randomSeedFlag, const int stream=0) {
  if (stream>0) error("No implemented for DRAND48");

  // added to run properly on G5 using gcc 3.3
  unsigned short seed16v[3]={0,0,0};

  if (randomSeedFlag) {
    struct timeval val;
    struct timezone zone;
    gettimeofday(&val,&zone);
    long t = val.tv_sec*1000000 + val.tv_usec; // 'long' should have 8 bytes, 'short' should have 2 bytes
    seed16v[0] = (unsigned short)((t >> 32) & 0xFFFF) ;
    seed16v[1] = (unsigned short)((t >> 16) & 0xFFFF) ;
    seed16v[2] = (unsigned short)( t        & 0xFFFF) ;
    /*
    unsigned short* tt = (unsigned short*) &t;
    seed16v[0] = tt[0];
    seed16v[1] = tt[1];
    seed16v[2] = tt[2];
    */
    //    srand (time (NULL));
  }

  //  unsigned short *j;
  //  j = seed48(seed16v);
  seed48(seed16v);
}

inline double Random() { 
  return drand48();
}

////////////////////////////////////////////////////////////////////////////////////////

inline double RandomRange(const double d) { // return -d ... +d
  return (2.0*Random()-1.0)*d;
}
inline double RandomRange(const double a, const double b) { // return a ... b
  return (b-a)*Random()+a;
}
inline double RandomRangeLogScale(const double a, const double b) {
  return a*exp(log(b/a)*Random());
}

inline double BoxMuller() {
  double x1 = Random();
  double x2 = Random();
  double r = sqrt(-2.0*log(x1)) * cos(2.0*pi*x2);
  return r;
}

inline void BoxMuller(double & r1, double & r2) {
  double x1 = Random();
  double x2 = Random();
  double f  = sqrt(-2.0*log(x1));
  double a  = 2.0*pi*x2;
  r1 = f * cos(a);
  r2 = f * sin(a);
}

#endif // _RANDOM_