#include <float.h>

#ifdef __cplusplus
extern "C" {
#endif
    void ranlux_(float* f,const int* ngen);
    void rluxgo_(const int* lux,const int* seed,const int* k1,const int* k2);
#ifdef __cplusplus
}
#endif


/*          Joel Heinrich  March 18, 2005

C and C++ interface to CERNLIB's RANLUX (Fortran) uniform random
number generator. Must be linked with CERNLIB's libmathlib.a and the
fortran library; e.g. /usr/lib/libg2c.a on Linux.

*/



/* rlx returns a uniform random deviate */

#define NV 32

double rlx() {
  static const double scl = 1.0 + 0.375*FLT_EPSILON;
  static float v[NV];
  static const int nv = NV;
  static const float* const end = &v[NV];
  static const float* p = &v[NV];
  if(p==end) {
    p=v;
    ranlux_(v,&nv);
  }
  return *p++ * scl ;
  /* Multiplication by scl guarantees that the result is not
     any simple fraction like 0.5, while still maintaining
     0 < rlx() < 1.  This reduces the chance of problems like
     division by zero occurring in the users code; it has nothing
     to do with increasing randomness.  The value used for scl is
     specifically tailored for ranlux.f and may fail with
     other generators.       */
}


#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <assert.h>

static int pid(); /* defined below */

/*          Joel Heinrich  June 13, 2005

   rlxrandomize initializes RANLUX with seed taken from, in this
    order of precedence:

   1. Environmental variable RLXSEED (assumed to be decimal)
   2. System device /dev/urandom + proccess id
   3. System time (number of seconds since Jan 1 1970)

   (lux=3 is ranlux's default luxury level)

   /dev/urandom has a bug in some SMP kernel versions: different processes
   that read /dev/urandom simultaneously on a hyperthreaded CPU
   will read the same bytes.  Adding the process id is a workaround.

*/

void rlxrandomize() {
  int seed=0;
  const char* const env = getenv("RLXSEED");

  if(env)
    seed = abs((int)strtol(env,NULL,10));

  if(seed==0) {
    FILE* const rdev = fopen("/dev/urandom","r");
    if(rdev) {
      const size_t nr =  fread(&seed,sizeof(seed),1,rdev);
      seed = (nr) ? abs(seed+pid()) : 0;
      fclose(rdev);
    }
  }
    
  if(seed==0) {
    const time_t systime = time(NULL), bad = (time_t) -1 ;
    assert(systime != bad);
    seed = abs((int)systime);
  }
  
  assert(seed != 0);

  {
    float v[24];
    const int lux=3, zero=0, nv = sizeof(v)/sizeof(v[0]);
    int i;
    rluxgo_(&lux,&seed,&zero,&zero);

    /* exercise RANLUX and perform simple sanity check */
    memset(v,0,sizeof(v));
    ranlux_(v,&nv);
    for(i=0;i<nv;++i)
      assert(v[i] > 0.0 && v[i] < 1.0); 
  }

  return;
}

/* Unix/Linux getpid returns id number of current process.
   Isolate call in separate function because getpid
   is not part of standard C library.                         */

#include<unistd.h>
static int pid() {
  return getpid();
}