Imported Upstream version 3.0+20111124upstream/3.0+20111124upstream

5 files changed, 1038 insertions, 0 deletions

diff --git a/src/common/ChiSquare.c b/src/common/ChiSquare.c
new file mode 100644
index 0000000..baee413
--- /dev/null
+++ b/src/common/ChiSquare.c
@@ -0,0 +1,67 @@
+/*
+	ChiSquare.c
+		the chi-square cdf
+		after W.J. Kennedy and J.E. Gentle,
+		Statistical computing, p. 116
+		last modified 9 Feb 05 th
+*/
+
+#ifdef HAVE_ERF
+#define Erf erf
+#else
+#include "Erf.c"
+#endif
+
+static inline real Normal(creal x)
+{
+  return .5*Erf(x/1.414213562373095048801689) + .5;
+}
+
+/*********************************************************************/
+
+static real ChiSquare(creal x, cint df)
+{
+  real y;
+
+  if( df <= 0 ) return -999;
+
+  if( x <= 0 ) return 0;
+  if( x > 1000*df ) return 1;
+
+  if( df > 1000 ) {
+    if( x < 2 ) return 0;
+    y = 2./(9*df);
+    y = (pow(x/df, 1/3.) - (1 - y))/sqrt(y);
+    if( y > 5 ) return 1;
+    if( y < -18.8055 ) return 0;
+    return Normal(y);
+  }
+
+  y = .5*x;
+
+  if( df & 1 ) {
+    creal sqrty = sqrt(y);
+    real h = Erf(sqrty);
+    count i;
+
+    if( df == 1 ) return h;
+
+    y = sqrty*exp(-y)/.8862269254527579825931;
+    for( i = 3; i < df; i += 2 ) {
+      h -= y;
+      y *= x/i;
+    }
+    y = h - y;
+  }
+  else {
+    real term = exp(-y), sum = term;
+    count i;
+
+    for( i = 1; i < df/2; ++i )
+      sum += term *= y/i;
+    y = 1 - sum;
+  }
+
+  return Max(0., y);
+}
+
diff --git a/src/common/DoSample.c b/src/common/DoSample.c
new file mode 100644
index 0000000..a637674
--- /dev/null
+++ b/src/common/DoSample.c
@@ -0,0 +1,352 @@
+/*
+	DoSample.c
+		the actual sampling routine, serial and parallel,
+		for the C versions of the Cuba routines
+		by Thomas Hahn
+		last modified 24 Nov 11 th
+*/
+
+#define MINSLICE 10
+#define MINCORES 1
+//#define MINCORES 2
+
+#if defined(VEGAS) || defined(SUAVE)
+#define VEG_ONLY(...) __VA_ARGS__
+#else
+#define VEG_ONLY(...)
+#endif
+
+#ifdef DIVONNE
+#define DIV_ONLY(...) __VA_ARGS__
+#define LDX(ldx) ldx
+#else
+#define DIV_ONLY(...)
+#define LDX(ldx) t->ndim
+#endif
+
+typedef struct {
+  real *f;
+  number n;
+  VEG_ONLY(count iter;)
+  DIV_ONLY(number neval_opt, neval_cut;
+           count ldx, phase, iregion;)
+#define NREGIONS ldx
+#define NEVAL n
+#define RETVAL phase
+} Slice;
+
+/*********************************************************************/
+
+#ifndef MSG_WAITALL
+/* Windows */
+#define MSG_WAITALL 0
+#endif
+
+static inline int readsock(int fd, void *data, size_t n)
+{
+  ssize_t got;
+  size_t remain = n;
+  do got = recv(fd, data, remain, MSG_WAITALL);
+  while( got > 0 && (data += got, remain -= got) > 0 );
+  return got;
+}
+
+static inline int writesock(int fd, const void *data, size_t n)
+{
+  ssize_t got;
+  size_t remain = n;
+  do got = send(fd, data, remain, MSG_WAITALL);
+  while( got > 0 && (data += got, remain -= got) > 0 );
+  return got;
+}
+
+/*********************************************************************/
+
+static inline bool SampleSerial(cThis *t, number n, creal *x, real *f
+  VEG_ONLY(, creal *w, ccount iter)
+  DIV_ONLY(, ccount ldx))
+{
+  while( n-- ) {
+    if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata
+          VEG_ONLY(, w++, &iter)
+          DIV_ONLY(, &t->phase)) == ABORT ) return true;
+    x += LDX(ldx);
+    f += t->ncomp;
+  }
+  return false;
+}
+
+/*********************************************************************/
+
+static void DoSample(This *t, number n, creal *x, real *f
+  VEG_ONLY(, creal *w, ccount iter)
+  DIV_ONLY(, ccount ldx))
+{
+  char s[128];
+  Slice slice;
+  int ncores;
+
+  t->neval += n;
+
+  ncores = IMin(t->ncores, n/MINSLICE);
+
+  if( ncores < MINCORES ) {
+    if( VERBOSE > 2 ) {
+      sprintf(s, "sampling " NUMBER " points serially", n);
+      Print(s);
+    }
+
+    if( SampleSerial(t, n, x, f
+          VEG_ONLY(, w, iter)
+          DIV_ONLY(, ldx)) ) longjmp(t->abort, -99);
+  }
+  else {
+    int core, abort;
+
+    slice.n = (n + ncores - 1)/ncores;
+
+    if( VERBOSE > 2 ) {
+      sprintf(s, "sampling " NUMBER " points each on %d cores",
+        slice.n, ncores);
+      Print(s);
+    }
+
+    slice.f = f;
+    VEG_ONLY(slice.iter = iter;)
+    DIV_ONLY(slice.ldx = ldx;)
+    DIV_ONLY(slice.phase = t->phase;)
+
+    for( core = 0; core < ncores; ++core ) {
+      cint fd = t->child[core];
+      writesock(fd, &slice, sizeof slice);
+      VEG_ONLY(writesock(fd, w, slice.n*sizeof *w);)
+      writesock(fd, x, slice.n*LDX(ldx)*sizeof *x);
+
+      VEG_ONLY(w += n;)
+      x += slice.n*LDX(ldx);
+      slice.f += slice.n*t->ncomp;
+      n -= slice.n;
+      slice.n = IMin(slice.n, n);
+    }
+
+    abort = 0;
+    for( core = ncores; --core >= 0; ) {
+      cint fd = t->child[core];
+      readsock(fd, &slice, sizeof slice);
+      if( slice.n == 0 ) abort = 1;
+      else readsock(fd, slice.f, slice.n*t->ncomp*sizeof *f);
+    }
+    if( abort ) longjmp(t->abort, -99);
+  }
+}
+
+/*********************************************************************/
+
+#ifdef DIVONNE
+static inline int ReadyCore(cThis *t)
+{
+  int core;
+  fd_set ready;
+
+  memcpy(&ready, &t->children, sizeof ready);
+  select(t->nchildren, &ready, NULL, NULL, NULL);
+
+  for( core = 0; core < t->ncores; ++core )
+    if( FD_ISSET(t->child[core], &ready) ) break;
+
+  return core;
+}
+
+/*********************************************************************/
+
+static int ExploreParent(This *t, cint iregion)
+{
+  TYPEDEFREGION;
+  Region *region;
+  Slice slice;
+  int ireg = iregion, core = t->running;
+
+  if( t->ncores < MINCORES ) return Explore(t, iregion);
+
+  if( t->running >= ((iregion < 0) ? 1 : t->ncores) ) {
+    Totals totals[t->ncomp];
+    count comp, succ;
+    cint fd = t->child[core = ReadyCore(t)];
+
+    --t->running;
+    readsock(fd, &slice, sizeof slice);
+//DEBSLICE("parent read", fd, slice);
+    ireg = slice.iregion;
+    region = RegionPtr(ireg);
+    succ = ireg + region->next;
+    readsock(fd, region, sizeof(Region));
+    if( --slice.NREGIONS > 0 ) {
+      region->next = t->nregions - ireg;
+      EnlargeRegions(t, slice.NREGIONS);
+      readsock(fd, RegionPtr(t->nregions), slice.NREGIONS*sizeof(Region));
+      t->nregions += slice.NREGIONS;
+      RegionPtr(t->nregions-1)->next = succ - t->nregions + 1;
+    }
+
+    readsock(fd, totals, sizeof totals);
+    for( comp = 0; comp < t->ncomp; ++comp )
+      t->totals[comp].secondspread =
+        Max(t->totals[comp].secondspread, totals[comp].secondspread);
+
+    t->neval += slice.NEVAL;
+    t->neval_opt += slice.neval_opt;
+    t->neval_cut += slice.neval_cut;
+
+    if( slice.RETVAL == -1 ) return -1;
+  }
+
+  if( iregion >= 0 ) {
+    region = RegionPtr(iregion);
+    cint fd = t->child[core];
+    slice.n = 0;
+    slice.phase = t->phase;
+    slice.iregion = iregion;
+//DEBSLICE("  parent write", fd, slice);
+    writesock(fd, &slice, sizeof slice);
+    writesock(fd, &t->samples[region->isamples], sizeof(Samples));
+    writesock(fd, region, sizeof *region);
+    writesock(fd, t->totals, sizeof *t->totals);
+    region->depth = 0;
+    ++t->running;
+  }
+
+  return ireg;
+}
+#endif
+
+/*********************************************************************/
+
+static inline void DoChild(This *t, cint fd)
+{
+  Slice slice;
+
+#ifdef DIVONNE
+  TYPEDEFREGION;
+  Totals totals[t->ncomp];
+
+  t->totals = totals;
+  t->ncores = 0;	/* no recursive forks */
+  AllocRegions(t);
+  SamplesIni(&t->samples[0]);
+  t->samples[0].n = 0;
+  SamplesIni(&t->samples[1]);
+  t->samples[1].n = 0;
+  SamplesIni(&t->samples[2]);
+  t->samples[2].n = 0;
+#endif
+
+  while( readsock(fd, &slice, sizeof slice) ) {
+    number n = slice.n;
+    DIV_ONLY(t->phase = slice.phase;)
+//DEBSLICE("  child read", fd, slice);
+    if( n > 0 ) {
+      VEG_ONLY(real w[n];)
+      real x[n*LDX(slice.ldx)];
+      real f[n*t->ncomp];
+
+      VEG_ONLY(readsock(fd, w, sizeof w);)
+      readsock(fd, x, sizeof x);
+
+      if( SampleSerial(t, n, x, f
+            VEG_ONLY(, w, slice.iter)
+            DIV_ONLY(, slice.ldx)) ) slice.n = 0;
+      writesock(fd, &slice, sizeof slice);
+      if( slice.n ) writesock(fd, f, sizeof f);
+    }
+#ifdef DIVONNE
+    else {
+      Samples *samples, psamples;
+
+      readsock(fd, &psamples, sizeof psamples);
+      readsock(fd, RegionPtr(0), sizeof(Region));
+      readsock(fd, totals, sizeof totals);
+      t->nregions = 1;
+      t->neval = t->neval_opt = t->neval_cut = 0;
+
+      samples = &t->samples[RegionPtr(0)->isamples];
+      if( psamples.n != samples->n ) {
+        SamplesFree(samples);
+        *samples = psamples;
+        SamplesAlloc(t, samples);
+      }
+
+      slice.RETVAL = Explore(t, 0);
+      slice.NREGIONS = t->nregions;
+      slice.NEVAL = t->neval;
+      slice.neval_opt = t->neval_opt;
+      slice.neval_cut = t->neval_cut;
+//DEBSLICE("child write", fd, slice);
+      writesock(fd, &slice, sizeof slice);
+      writesock(fd, RegionPtr(0), t->nregions*sizeof(Region));
+      writesock(fd, totals, sizeof totals);
+    }
+#endif
+  }
+
+  exit(0);
+}
+
+/*********************************************************************/
+
+static inline void ForkCores(This *t)
+{
+  int core;
+  cchar *env = getenv("CUBACORES");
+
+  t->ncores = env ? atoi(env) : sysconf(_SC_NPROCESSORS_ONLN);
+#ifdef HAVE_GETLOADAVG
+  if( env == NULL || t->ncores < 0 ) {
+    double load = 0;
+    getloadavg(&load, 1);
+    t->ncores = abs(t->ncores) - floor(load);
+  }
+#endif
+
+#ifdef DIVONNE
+  t->nchildren = t->running = 0;
+#endif
+
+  if( t->ncores < MINCORES ) return;
+  if( VERBOSE ) printf("using %d cores\n", t->ncores);
+  fflush(stdout);
+
+  Alloc(t->child, t->ncores);
+  for( core = 0; core < t->ncores; ++core ) {
+    int fd[2];
+    pid_t pid;
+    assert(
+      socketpair(AF_LOCAL, SOCK_STREAM, 0, fd) != -1 &&
+      (pid = fork()) != -1 );
+    if( pid == 0 ) {
+      close(fd[0]);
+      DoChild(t, fd[1]);
+    }
+    close(fd[1]);
+    t->child[core] = fd[0];
+#ifdef DIVONNE
+    FD_SET(fd[0], &t->children);
+    t->nchildren = IMax(t->nchildren, fd[0] + 1);
+#endif
+  }
+}
+
+/*********************************************************************/
+
+static inline void WaitCores(cThis *t)
+{
+  if( t->ncores >= MINCORES ) {
+    int core;
+    pid_t pid;
+    for( core = 0; core < t->ncores; ++core )
+      close(t->child[core]);
+    free(t->child);
+    for( core = 0; core < t->ncores; ++core )
+      wait(&pid);
+  }
+}
+
diff --git a/src/common/Erf.c b/src/common/Erf.c
new file mode 100644
index 0000000..635b6e5
--- /dev/null
+++ b/src/common/Erf.c
@@ -0,0 +1,51 @@
+/*
+	Erf.c
+		Gaussian error function
+		= 2/Sqrt[Pi] Integrate[Exp[-t^2], {t, 0, x}]
+		Code from Takuya Ooura's gamerf2a.f
+		http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html
+		last modified 8 Feb 05 th
+*/
+
+
+static real Erfc(creal x)
+{
+  static creal c[] = {
+    2.96316885199227378e-01, 6.12158644495538758e-02,
+    1.81581125134637070e-01, 5.50942780056002085e-01,
+    6.81866451424939493e-02, 1.53039662058770397e+00,
+    1.56907543161966709e-02, 2.99957952311300634e+00,
+    2.21290116681517573e-03, 4.95867777128246701e+00,
+    1.91395813098742864e-04, 7.41471251099335407e+00,
+    9.71013284010551623e-06, 1.04765104356545238e+01,
+    1.66642447174307753e-07, 1.48455557345597957e+01,
+    6.10399733098688199e+00, 1.26974899965115684e+01 };
+  real y = x*x;
+  y = exp(-y)*x*(
+    c[0]/(y + c[1]) + c[2]/(y + c[3]) +
+    c[4]/(y + c[5]) + c[6]/(y + c[7]) +
+    c[8]/(y + c[9]) + c[10]/(y + c[11]) +
+    c[12]/(y + c[13]) + c[14]/(y + c[15]) );
+  if( x < c[16] ) y += 2/(exp(c[17]*x) + 1);
+  return y;
+}
+
+
+static real Erf(creal x)
+{
+  static creal c[] = {
+    1.12837916709551257e+00,
+   -3.76126389031833602e-01,
+    1.12837916706621301e-01,
+   -2.68661698447642378e-02,
+    5.22387877685618101e-03,
+   -8.49202435186918470e-04 };
+  real y = fabs(x);
+  if( y > .125 ) {
+    y = 1 - Erfc(y);
+    return (x > 0) ? y : -y;
+  }
+  y *= y;
+  return x*(c[0] + y*(c[1] + y*(c[2] +
+    y*(c[3] + y*(c[4] + y*c[5])))));
+}
diff --git a/src/common/Random.c b/src/common/Random.c
new file mode 100644
index 0000000..66b180e
--- /dev/null
+++ b/src/common/Random.c
@@ -0,0 +1,345 @@
+/*
+	Random.c
+		quasi- and pseudo-random-number generation
+		last modified 8 Jun 10 th
+*/
+
+
+/*
+	PART 1: Sobol quasi-random-number generator
+	adapted from ACM TOMS algorithm 659
+*/
+
+static void SobolGet(This *t, real *x)
+{
+  number seq = t->rng.sobol.seq++;
+  count zerobit = 0, dim;
+
+  while( seq & 1 ) {
+    ++zerobit;
+    seq >>= 1;
+  }
+
+  for( dim = 0; dim < t->ndim; ++dim ) {
+    t->rng.sobol.prev[dim] ^= t->rng.sobol.v[dim][zerobit];
+    x[dim] = t->rng.sobol.prev[dim]*t->rng.sobol.norm;
+  }
+}
+
+
+static void SobolSkip(This *t, number n)
+{
+  while( n-- ) {
+    number seq = t->rng.sobol.seq++;
+    count zerobit = 0, dim;
+
+    while( seq & 1 ) {
+      ++zerobit;
+      seq >>= 1;
+    }
+
+    for( dim = 0; dim < t->ndim; ++dim )
+      t->rng.sobol.prev[dim] ^= t->rng.sobol.v[dim][zerobit];
+  }
+}
+
+
+static inline void SobolIni(This *t)
+{
+  static number ini[9*40] = {
+      3,   1,   0,   0,   0,   0,   0,   0,   0,
+      7,   1,   1,   0,   0,   0,   0,   0,   0,
+     11,   1,   3,   7,   0,   0,   0,   0,   0,
+     13,   1,   1,   5,   0,   0,   0,   0,   0,
+     19,   1,   3,   1,   1,   0,   0,   0,   0,
+     25,   1,   1,   3,   7,   0,   0,   0,   0,
+     37,   1,   3,   3,   9,   9,   0,   0,   0,
+     59,   1,   3,   7,  13,   3,   0,   0,   0,
+     47,   1,   1,   5,  11,  27,   0,   0,   0,
+     61,   1,   3,   5,   1,  15,   0,   0,   0,
+     55,   1,   1,   7,   3,  29,   0,   0,   0,
+     41,   1,   3,   7,   7,  21,   0,   0,   0,
+     67,   1,   1,   1,   9,  23,  37,   0,   0,
+     97,   1,   3,   3,   5,  19,  33,   0,   0,
+     91,   1,   1,   3,  13,  11,   7,   0,   0,
+    109,   1,   1,   7,  13,  25,   5,   0,   0,
+    103,   1,   3,   5,  11,   7,  11,   0,   0,
+    115,   1,   1,   1,   3,  13,  39,   0,   0,
+    131,   1,   3,   1,  15,  17,  63,  13,   0,
+    193,   1,   1,   5,   5,   1,  27,  33,   0,
+    137,   1,   3,   3,   3,  25,  17, 115,   0,
+    145,   1,   1,   3,  15,  29,  15,  41,   0,
+    143,   1,   3,   1,   7,   3,  23,  79,   0,
+    241,   1,   3,   7,   9,  31,  29,  17,   0,
+    157,   1,   1,   5,  13,  11,   3,  29,   0,
+    185,   1,   3,   1,   9,   5,  21, 119,   0,
+    167,   1,   1,   3,   1,  23,  13,  75,   0,
+    229,   1,   3,   3,  11,  27,  31,  73,   0,
+    171,   1,   1,   7,   7,  19,  25, 105,   0,
+    213,   1,   3,   5,   5,  21,   9,   7,   0,
+    191,   1,   1,   1,  15,   5,  49,  59,   0,
+    253,   1,   1,   1,   1,   1,  33,  65,   0,
+    203,   1,   3,   5,  15,  17,  19,  21,   0,
+    211,   1,   1,   7,  11,  13,  29,   3,   0,
+    239,   1,   3,   7,   5,   7,  11, 113,   0,
+    247,   1,   1,   5,   3,  15,  19,  61,   0,
+    285,   1,   3,   1,   1,   9,  27,  89,   7,
+    369,   1,   1,   3,   7,  31,  15,  45,  23,
+    299,   1,   3,   3,   9,   9,  25, 107,  39 };
+
+  count dim, bit, nbits;
+  number max, *pini = ini;
+  cnumber nmax = 2*t->maxeval;
+
+  for( nbits = 0, max = 1; max <= nmax; max <<= 1 ) ++nbits;
+  t->rng.sobol.norm = 1./max;
+
+  for( bit = 0; bit < nbits; ++bit )
+    t->rng.sobol.v[0][bit] = (max >>= 1);
+
+  for( dim = 1; dim < t->ndim; ++dim ) {
+    number *pv = t->rng.sobol.v[dim], *pvv = pv;
+    number powers = *pini++, j;
+    int inibits = -1, bit;
+    for( j = powers; j; j >>= 1 ) ++inibits;
+
+    memcpy(pv, pini, inibits*sizeof(*pini));
+    pini += 8;
+
+    for( bit = inibits; bit < nbits; ++bit ) {
+      number newv = *pvv, j = powers;
+      int b;
+      for( b = 0; b < inibits; ++b ) {
+        if( j & 1 ) newv ^= pvv[b] << (inibits - b);
+        j >>= 1;
+      }
+      pvv[inibits] = newv;
+      ++pvv;
+    }
+
+    for( bit = 0; bit < nbits - 1; ++bit )
+      pv[bit] <<= nbits - bit - 1;
+  }
+
+  t->rng.sobol.seq = 0;
+  VecClear(t->rng.sobol.prev);
+
+  t->rng.getrandom = SobolGet;
+  t->rng.skiprandom = SobolSkip;
+}
+
+
+/*
+	PART 2: Mersenne Twister pseudo-random-number generator
+	adapted from T. Nishimura's and M. Matsumoto's C code at
+	http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+*/
+
+/* 32 or 53 random bits */
+#define RANDOM_BITS 32
+
+
+static inline state_t Twist(state_t a, state_t b)
+{
+  state_t mixbits = (a & 0x80000000) | (b & 0x7fffffff);
+  state_t matrixA = (-(b & 1)) & 0x9908b0df;
+  return (mixbits >> 1) ^ matrixA;
+}
+
+
+static inline void MersenneReload(state_t *state)
+{
+  state_t *s = state;
+  int j;
+
+  for( j = MERSENNE_N - MERSENNE_M + 1; --j; ++s )
+    *s = s[MERSENNE_M] ^ Twist(s[0], s[1]);
+  for( j = MERSENNE_M; --j; ++s )
+    *s = s[MERSENNE_M - MERSENNE_N] ^ Twist(s[0], s[1]);
+  *s = s[MERSENNE_M - MERSENNE_N] ^ Twist(s[0], state[0]);
+}
+
+
+static inline state_t MersenneInt(state_t s)
+{
+  s ^= s >> 11;
+  s ^= (s << 7) & 0x9d2c5680;
+  s ^= (s << 15) & 0xefc60000;
+  return s ^ (s >> 18);
+}
+
+
+static void MersenneGet(This *t, real *x)
+{
+  count next = t->rng.mersenne.next, dim;
+
+  for( dim = 0; dim < t->ndim; ++dim ) {
+#if RANDOM_BITS == 53
+    state_t a, b;
+#endif
+
+    if( next >= MERSENNE_N ) {
+      MersenneReload(t->rng.mersenne.state);
+      next = 0;
+    }
+
+#if RANDOM_BITS == 53
+    a = MersenneInt(t->rng.mersenne.state[next++]) >> 5;
+    b = MersenneInt(t->rng.mersenne.state[next++]) >> 6;
+    x[dim] = (67108864.*a + b)/9007199254740992.;
+#else
+    x[dim] = MersenneInt(t->rng.mersenne.state[next++])/4294967296.;
+#endif
+  }
+
+  t->rng.mersenne.next = next;
+}
+
+
+static void MersenneSkip(This *t, number n)
+{
+#if RANDOM_BITS == 53
+  n = 2*n*t->ndim + t->rng.mersenne.next;
+#else
+  n = n*t->ndim + t->rng.mersenne.next;
+#endif
+  t->rng.mersenne.next = n % MERSENNE_N;
+  n /= MERSENNE_N;
+  while( n-- ) MersenneReload(t->rng.mersenne.state);
+}
+
+
+static inline void MersenneIni(This *t)
+{
+  state_t seed = t->seed;
+  state_t *next = t->rng.mersenne.state;
+  count j;
+
+  for( j = 1; j <= MERSENNE_N; ++j ) {
+    *next++ = seed;
+    seed = 0x6c078965*(seed ^ (seed >> 30)) + j;
+    /* see Knuth TAOCP Vol 2, 3rd Ed, p. 106 for multiplier */
+  }
+
+  MersenneReload(t->rng.mersenne.state);
+  t->rng.mersenne.next = 0;
+
+  t->rng.getrandom = MersenneGet;
+  t->rng.skiprandom = MersenneSkip;
+}
+
+
+/*
+	PART 3: Ranlux subtract-and-borrow random-number generator 
+	proposed by Marsaglia and Zaman, implemented by F. James with 
+	the name RCARRY in 1991, and later improved by Martin Luescher 
+	in 1993 to produce "Luxury Pseudorandom Numbers".
+	Adapted from the CERNlib Fortran 77 code by F. James, 1993.
+
+	The available luxury levels are:
+
+	level 0  (p = 24): equivalent to the original RCARRY of Marsaglia
+	         and Zaman, very long period, but fails many tests.
+	level 1  (p = 48): considerable improvement in quality over level 0,
+	         now passes the gap test, but still fails spectral test.
+	level 2  (p = 97): passes all known tests, but theoretically still
+	         defective.
+	level 3  (p = 223): DEFAULT VALUE.  Any theoretically possible
+	         correlations have very small chance of being observed.
+	level 4  (p = 389): highest possible luxury, all 24 bits chaotic.
+*/
+
+
+static inline int RanluxInt(This *t, count n)
+{
+  int s = 0;
+
+  while( n-- ) {
+    s = t->rng.ranlux.state[t->rng.ranlux.j24] -
+        t->rng.ranlux.state[t->rng.ranlux.i24] + t->rng.ranlux.carry;
+    s += (t->rng.ranlux.carry = NegQ(s)) & (1 << 24);
+    t->rng.ranlux.state[t->rng.ranlux.i24] = s;
+    --t->rng.ranlux.i24;
+    t->rng.ranlux.i24 += NegQ(t->rng.ranlux.i24) & 24;
+    --t->rng.ranlux.j24;
+    t->rng.ranlux.j24 += NegQ(t->rng.ranlux.j24) & 24;
+  }
+
+  return s;
+}
+
+
+static void RanluxGet(This *t, real *x)
+{
+/* The Generator proper: "Subtract-with-borrow",
+   as proposed by Marsaglia and Zaman, FSU, March 1989 */
+
+  count dim;
+
+  for( dim = 0; dim < t->ndim; ++dim ) {
+    cint nskip = (--t->rng.ranlux.n24 >= 0) ? 0 :
+      (t->rng.ranlux.n24 = 24, t->rng.ranlux.nskip);
+    cint s = RanluxInt(t, 1 + nskip);
+    x[dim] = s*0x1p-24;
+/* small numbers (with less than 12 significant bits) are "padded" */
+    if( s < (1 << 12) )
+      x[dim] += t->rng.ranlux.state[t->rng.ranlux.j24]*0x1p-48;
+  }
+}
+
+
+static void RanluxSkip(This *t, cnumber n)
+{
+  RanluxInt(t, n + t->rng.ranlux.nskip*(n/24));
+  t->rng.ranlux.n24 = 24 - n % 24;
+}
+
+
+static inline void RanluxIni(This *t)
+{
+  cint skip[] = {24, 48, 97, 223, 389,
+    223, 223, 223, 223, 223, 223, 223, 223, 223, 223,
+    223, 223, 223, 223, 223, 223, 223, 223, 223, 223};
+  state_t seed = t->seed;
+  state_t level = RNG;
+  count i;
+
+  if( level < sizeof skip ) level = skip[level];
+  t->rng.ranlux.nskip = level - 24;
+
+  t->rng.ranlux.i24 = 23;
+  t->rng.ranlux.j24 = 9;
+  t->rng.ranlux.n24 = 24;
+
+  for( i = 0; i < 24; ++i ) {
+    cint k = seed/53668;
+    seed = 40014*(seed - k*53668) - k*12211;
+    seed += NegQ(seed) & 2147483563;
+    t->rng.ranlux.state[i] = seed & ((1 << 24) - 1);
+  }
+
+  t->rng.ranlux.carry = ~TrueQ(t->rng.ranlux.state[23]) & (1 << 24);
+
+  t->rng.getrandom = RanluxGet;
+  t->rng.skiprandom = RanluxSkip;
+}
+
+
+/*
+	PART 4: User routines:
+
+	- IniRandom sets up the random-number generator to produce a
+	  sequence of at least n ndim-dimensional random vectors.
+
+	- GetRandom retrieves one random vector.
+
+	- SkipRandom skips over n random vectors.
+*/
+
+static inline void IniRandom(This *t)
+{
+  if( t->seed == 0 ) SobolIni(t);
+  else if( RNG == 0 ) MersenneIni(t);
+  else RanluxIni(t);
+}
+
diff --git a/src/common/stddecl.h b/src/common/stddecl.h
new file mode 100644
index 0000000..6bf3c27
--- /dev/null
+++ b/src/common/stddecl.h
@@ -0,0 +1,223 @@
+/*
+	stddecl.h
+		Type declarations common to all Cuba routines
+		last modified 28 Sep 11 th
+*/
+
+
+#ifndef _stddecl_h_
+#define _stddecl_h_
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <float.h>
+#include <limits.h>
+#include <unistd.h>
+#include <assert.h>
+#include <fcntl.h>
+#include <setjmp.h>
+#include <sys/wait.h>
+#include <sys/stat.h>
+#include <sys/types.h>
+#include <sys/socket.h>
+
+
+#ifndef NDIM
+#define NDIM t->ndim
+#endif
+#ifndef NCOMP
+#define NCOMP t->ncomp
+#endif
+
+
+#define VERBOSE (t->flags & 3)
+#define LAST (t->flags & 4)
+#define SHARPEDGES (t->flags & 8)
+#define REGIONS (t->flags & 128)
+#define RNG (t->flags >> 8)
+
+#define INFTY DBL_MAX
+
+#define NOTZERO 0x1p-104
+
+#define ABORT -999
+
+#define Elements(x) (sizeof(x)/sizeof(*x))
+
+#define Copy(d, s, n) memcpy(d, s, (n)*sizeof(*(d)))
+
+#define VecCopy(d, s) Copy(d, s, t->ndim)
+
+#define ResCopy(d, s) Copy(d, s, t->ncomp)
+
+#define Clear(d, n) memset(d, 0, (n)*sizeof(*(d)))
+
+#define VecClear(d) Clear(d, t->ndim)
+
+#define ResClear(d) Clear(d, t->ncomp)
+
+#define Zap(d) memset(d, 0, sizeof(d))
+
+#define MaxErr(avg) Max(t->epsrel*fabs(avg), t->epsabs)
+
+#ifdef __cplusplus
+#define mallocset(p, n) (*(void **)&p = malloc(n))
+#define reallocset(p, n) (*(void **)&p = realloc(p, n))
+#else
+#define mallocset(p, n) (p = malloc(n))
+#define reallocset(p, n) (p = realloc(p, n))
+#endif
+
+#define ChkAlloc(r) if( r == NULL ) { \
+  fprintf(stderr, "Out of memory in " __FILE__ " line %d.\n", __LINE__); \
+  exit(1); \
+}
+
+#define Alloc(p, n) MemAlloc(p, (n)*sizeof(*p))
+#define MemAlloc(p, n) ChkAlloc(mallocset(p, n))
+#define ReAlloc(p, n) ChkAlloc(reallocset(p, n))
+
+
+#ifdef __cplusplus
+#define Extern extern "C"
+#else
+#define Extern extern
+typedef enum { false, true } bool;
+#endif
+
+typedef const char cchar;
+
+typedef const bool cbool;
+
+typedef const int cint;
+
+typedef const long clong;
+
+#define COUNT "%d"
+typedef /*unsigned*/ int count;
+typedef const count ccount;
+
+#ifdef LONGLONGINT
+#define PREFIX(s) ll##s
+#define NUMBER "%lld"
+#define NUMBER7 "%7lld"
+typedef long long int number;
+#else
+#define PREFIX(s) s
+#define NUMBER "%d"
+#define NUMBER7 "%7d"
+typedef int number;
+#endif
+typedef const number cnumber;
+
+#define REAL "%g"
+#define REALF "%f"
+typedef /*long*/ double real;
+	/* Switching to long double is not as trivial as it
+	   might seem here.  sqrt, erf, exp, pow need to be
+	   replaced by their long double versions (sqrtl, ...),
+	   printf formats need to be updated similarly, and
+	   ferrying long doubles to Mathematica is of course
+	   quite another matter, too. */
+
+typedef const real creal;
+
+
+struct _this;
+
+typedef unsigned int state_t;
+
+#define SOBOL_MINDIM 1
+#define SOBOL_MAXDIM 40
+
+/* length of state vector */
+#define MERSENNE_N 624
+
+/* period parameter */
+#define MERSENNE_M 397
+
+typedef struct {
+  void (*getrandom)(struct _this *t, real *x);
+  void (*skiprandom)(struct _this *t, cnumber n);
+  union {
+    struct {
+      real norm;
+      number v[SOBOL_MAXDIM][30], prev[SOBOL_MAXDIM];
+      number seq;
+    } sobol;
+    struct {
+      state_t state[MERSENNE_N];
+      count next;
+    } mersenne;
+    struct {
+      count n24, i24, j24, nskip;
+      int carry, state[24];
+    } ranlux;
+  };
+} RNGState;
+
+
+#if NOUNDERSCORE
+#define SUFFIX(s) s
+#else
+#define SUFFIX(s) s##_
+#endif
+
+#define EXPORT(s) EXPORT_(PREFIX(s))
+#define EXPORT_(s) SUFFIX(s)
+
+
+static inline real Sq(creal x)
+{
+  return x*x;
+}
+
+static inline real Min(creal a, creal b)
+{
+  return (a < b) ? a : b;
+}
+
+static inline real Max(creal a, creal b)
+{
+  return (a > b) ? a : b;
+}
+
+static inline real Weight(creal sum, creal sqsum, cnumber n)
+{
+  creal w = sqrt(sqsum*n);
+  return (n - 1)/Max((w + sum)*(w - sum), NOTZERO);
+}
+
+
+/* (a < 0) ? -1 : 0 */
+#define NegQ(a) ((a) >> (sizeof(a)*8 - 1))
+
+/* (a < 0) ? -1 : 1 */
+#define Sign(a) (1 + 2*NegQ(a))
+
+/* (a < 0) ? 0 : a */
+#define IDim(a) ((a) & NegQ(-(a)))
+
+/* (a < b) ? a : b */
+#define IMin(a, b) ((a) - IDim((a) - (b)))
+
+/* (a > b) ? a : b */
+#define IMax(a, b) ((b) + IDim((a) - (b)))
+
+/* (a == 0) ? 0 : -1 */
+#define TrueQ(a) NegQ((a) | (-a))
+
+/* a + (a == 0) */
+#define Min1(a) ((a) + 1 + TrueQ(a))
+
+/* abs(a) + (a == 0) */
+#define Abs1(a) (((a) ^ NegQ(a)) - NegQ((a) - 1))
+
+#endif
+