diff --git a/configure.ac b/configure.ac
index 39ca5ca6..8b3cb01c 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1,7 +1,7 @@
 AC_REVISION([m4_esyscmd_s([git describe --always])])
 
 AC_PREREQ(2.63)
-AC_INIT([musrfit],[1.2.0],[andreas.suter@psi.ch])
+AC_INIT([musrfit],[1.2.1],[andreas.suter@psi.ch])
 AC_CONFIG_AUX_DIR(admin)
 AC_CANONICAL_HOST
 #AC_MSG_RESULT([${host} ${host_cpu} ${host_vendor} ${host_os}])
@@ -36,7 +36,7 @@ dnl -----------------------------------------------
 #release versioning
 MUSR_MAJOR_VERSION=1
 MUSR_MINOR_VERSION=2
-MUSR_MICRO_VERSION=0
+MUSR_MICRO_VERSION=1
 
 #release versioning
 MUSR_ROOT_MAJOR_VERSION=1
@@ -69,7 +69,7 @@ PLUGIN_MINOR_VERSION=0
 PLUGIN_MICRO_VERSION=0
 
 #release versioning
-CUBA_MAJOR_VERSION=3
+CUBA_MAJOR_VERSION=4
 CUBA_MINOR_VERSION=2
 CUBA_MICRO_VERSION=0
 
diff --git a/src/external/BMWtools/BMWIntegrator.cpp b/src/external/BMWtools/BMWIntegrator.cpp
index 63ae7a3c..96ce4154 100644
--- a/src/external/BMWtools/BMWIntegrator.cpp
+++ b/src/external/BMWtools/BMWIntegrator.cpp
@@ -31,6 +31,7 @@
 #include "cuba.h"
 
 #define USERDATA NULL
+#define SPIN NULL
 #define SEED 0
 #define STATEFILE NULL
 
@@ -45,6 +46,7 @@ std::vector<double> TDWaveGapIntegralCuhre::fPar;
 double TDWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -57,9 +59,9 @@ double TDWaveGapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -96,6 +98,7 @@ std::vector<double> TCosSqDWaveGapIntegralCuhre::fPar;
 double TCosSqDWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-8);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -108,9 +111,9 @@ double TCosSqDWaveGapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -147,6 +150,7 @@ std::vector<double> TSinSqDWaveGapIntegralCuhre::fPar;
 double TSinSqDWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-8);
   const double EPSABS (1e-10);
   const unsigned int VERBOSE (0);
@@ -159,9 +163,9 @@ double TSinSqDWaveGapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -198,6 +202,7 @@ std::vector<double> TAnSWaveGapIntegralCuhre::fPar;
 double TAnSWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -210,9 +215,9 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -249,6 +254,7 @@ std::vector<double> TAnSWaveGapIntegralDivonne::fPar;
 double TAnSWaveGapIntegralDivonne::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -268,10 +274,10 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Divonne(fNDim, NCOMP, Integrand, USERDATA,
+  Divonne(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE, SEED, MINEVAL, MAXEVAL,
     KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION,
-    NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL, STATEFILE,
+    NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -308,6 +314,7 @@ std::vector<double> TAnSWaveGapIntegralSuave::fPar;
 double TAnSWaveGapIntegralSuave::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -316,14 +323,15 @@ double TAnSWaveGapIntegralSuave::IntegrateFunc()
   const unsigned int MAXEVAL (1000000);
 
   const unsigned int NNEW (1000);
+  const unsigned int NMIN (2);
   const double FLATNESS (25.);
 
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Suave(fNDim, NCOMP, Integrand, USERDATA,
+  Suave(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, SEED, MINEVAL, MAXEVAL,
-    NNEW, FLATNESS, STATEFILE,
+    NNEW, NMIN, FLATNESS, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -360,6 +368,7 @@ std::vector<double> TNonMonDWave1GapIntegralCuhre::fPar;
 double TNonMonDWave1GapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -372,9 +381,9 @@ double TNonMonDWave1GapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -411,6 +420,7 @@ std::vector<double> TNonMonDWave2GapIntegralCuhre::fPar;
 double TNonMonDWave2GapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -423,9 +433,9 @@ double TNonMonDWave2GapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
diff --git a/src/external/libCuba/src/Makefile.am b/src/external/libCuba/src/Makefile.am
index 4a7d8f43..9771bde9 100644
--- a/src/external/libCuba/src/Makefile.am
+++ b/src/external/libCuba/src/Makefile.am
@@ -11,7 +11,7 @@ lib_LTLIBRARIES = libcuba.la
 
 libcuba_la_SOURCES =
 
-libcuba_la_LIBADD = cuhre/libcuhre.la divonne/libdivonne.la suave/libsuave.la vegas/libvegas.la -lm
+libcuba_la_LIBADD = common/libcommon.la cuhre/libcuhre.la divonne/libdivonne.la suave/libsuave.la vegas/libvegas.la -lm
 libcuba_la_LDFLAGS = -version-info $(CUBA_LIBRARY_VERSION) -release $(CUBA_RELEASE) $(AM_LDFLAGS)
 
 pkgconfigdir = $(libdir)/pkgconfig
diff --git a/src/external/libCuba/src/common/CSample.c b/src/external/libCuba/src/common/CSample.c
index 5705be06..dfdbae1e 100644
--- a/src/external/libCuba/src/common/CSample.c
+++ b/src/external/libCuba/src/common/CSample.c
@@ -3,19 +3,27 @@
 		the serial sampling routine
 		for the C versions of the Cuba routines
 		by Thomas Hahn
-		last modified 19 Dec 11 th
+		last modified 9 Oct 14 th
 */
 
 
-static inline number SampleRaw(cThis *t, number n, creal *x, real *f
-  VES_ONLY(, creal *w, ccount iter))
+coreinit cubafun_;
+extern int cubaverb_;
+extern corespec cubaworkers_;
+
+
+static inline number SampleRaw(This *t, number n, creal *x, real *f,
+  cint core VES_ONLY(, creal *w, ccount iter))
 {
-  for( ; n; --n ) {
-    if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata
-          VES_ONLY(, w++, &iter)
+  number nvec;
+  for( nvec = t->nvec; n > 0; n -= nvec ) {
+    nvec = IMin(n, nvec);
+    if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata, &nvec, &core
+          VES_ONLY(, w, &iter)
           DIV_ONLY(, &t->phase)) == ABORT ) return -1;
-    x += t->ndim;
-    f += t->ncomp;
+    VES_ONLY(w += nvec;)
+    x += nvec*t->ndim;
+    f += nvec*t->ncomp;
   }
   return 0;
 }
@@ -25,8 +33,9 @@ static inline number SampleRaw(cThis *t, number n, creal *x, real *f
 static inline void DoSampleSerial(This *t, cnumber n, creal *x, real *f
   VES_ONLY(, creal *w, ccount iter))
 {
+  MasterInit();
   t->neval += n;
-  if( SampleRaw(t, n, x, f VES_ONLY(, w, iter)) ) 
+  if( SampleRaw(t, n, x, f, -1 VES_ONLY(, w, iter)) ) 
     longjmp(t->abort, -99);
 }
 
@@ -43,7 +52,13 @@ DIV_ONLY(static int Explore(This *t, cint iregion);)
 #define DoSample DoSampleSerial
 #define Explore ExploreSerial
 #define ForkCores(t)
-#define WaitCores(t)
+
+static inline void WaitCores(This *t, Spin **pspin)
+{
+  if( Invalid(pspin) ) MasterExit();
+}
+
+#define WaitCores(t, pspin)
 
 #endif
 
@@ -51,7 +66,7 @@ DIV_ONLY(static int Explore(This *t, cint iregion);)
 static inline count SampleExtra(This *t, cBounds *b)
 {
   number n = t->nextra;
-  t->peakfinder(&t->ndim, b, &n, t->xextra);
+  t->peakfinder(&t->ndim, b, &n, t->xextra, t->userdata);
   DoSample(t, n, t->xextra, t->fextra);
   return n;
 }
@@ -60,7 +75,7 @@ static inline count SampleExtra(This *t, cBounds *b)
 #include "common.c"
 
 #ifdef HAVE_FORK
-#include "Fork.c"
+#include "Parallel.c"
 #endif
 
 #include "Integrate.c"
diff --git a/src/external/libCuba/src/common/ChiSquare.c b/src/external/libCuba/src/common/ChiSquare.c
index baee4131..fc257887 100644
--- a/src/external/libCuba/src/common/ChiSquare.c
+++ b/src/external/libCuba/src/common/ChiSquare.c
@@ -3,7 +3,7 @@
 		the chi-square cdf
 		after W.J. Kennedy and J.E. Gentle,
 		Statistical computing, p. 116
-		last modified 9 Feb 05 th
+		last modified 12 Mar 15 th
 */
 
 #ifdef HAVE_ERF
@@ -31,7 +31,7 @@ static real ChiSquare(creal x, cint df)
   if( df > 1000 ) {
     if( x < 2 ) return 0;
     y = 2./(9*df);
-    y = (pow(x/df, 1/3.) - (1 - y))/sqrt(y);
+    y = (powx(x/df, 1/3.) - (1 - y))/sqrtx(y);
     if( y > 5 ) return 1;
     if( y < -18.8055 ) return 0;
     return Normal(y);
@@ -40,13 +40,13 @@ static real ChiSquare(creal x, cint df)
   y = .5*x;
 
   if( df & 1 ) {
-    creal sqrty = sqrt(y);
+    creal sqrty = sqrtx(y);
     real h = Erf(sqrty);
     count i;
 
     if( df == 1 ) return h;
 
-    y = sqrty*exp(-y)/.8862269254527579825931;
+    y = sqrty*expx(-y)/.8862269254527579825931;
     for( i = 3; i < df; i += 2 ) {
       h -= y;
       y *= x/i;
@@ -54,7 +54,7 @@ static real ChiSquare(creal x, cint df)
     y = h - y;
   }
   else {
-    real term = exp(-y), sum = term;
+    real term = expx(-y), sum = term;
     count i;
 
     for( i = 1; i < df/2; ++i )
diff --git a/src/external/libCuba/src/common/Erf.c b/src/external/libCuba/src/common/Erf.c
index 635b6e5c..c94d106d 100644
--- a/src/external/libCuba/src/common/Erf.c
+++ b/src/external/libCuba/src/common/Erf.c
@@ -4,7 +4,7 @@
 		= 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
+		last modified 12 Mar 15 th
 */
 
 
@@ -21,12 +21,12 @@ static real Erfc(creal x)
     1.66642447174307753e-07, 1.48455557345597957e+01,
     6.10399733098688199e+00, 1.26974899965115684e+01 };
   real y = x*x;
-  y = exp(-y)*x*(
+  y = expx(-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);
+  if( x < c[16] ) y += 2/(expx(c[17]*x) + 1);
   return y;
 }
 
@@ -40,7 +40,7 @@ static real Erf(creal x)
    -2.68661698447642378e-02,
     5.22387877685618101e-03,
    -8.49202435186918470e-04 };
-  real y = fabs(x);
+  real y = fabsx(x);
   if( y > .125 ) {
     y = 1 - Erfc(y);
     return (x > 0) ? y : -y;
diff --git a/src/external/libCuba/src/common/Fork.c b/src/external/libCuba/src/common/Fork.c
index 51b18884..99605d25 100644
--- a/src/external/libCuba/src/common/Fork.c
+++ b/src/external/libCuba/src/common/Fork.c
@@ -1,391 +1,98 @@
 /*
 	Fork.c
-		the parallel sampling routine
-		for the C versions of the Cuba routines
+		fork the cores for parallel sampling
+		(C version only)
 		by Thomas Hahn
-		last modified 25 Sep 13 th
+		last modified 23 Apr 15 th
 */
 
-#include <sys/select.h>
 
-#define MINSLICE 10
+#define ROUTINE "cubafork"
+#include "stddecl.h"
+
+#ifdef HAVE_FORK
+
+#include "sock.h"
+
 #define MINCORES 1
-/*#define MINCORES 2*/
 
-typedef struct {
-  number n, m, i;
-  VES_ONLY(count iter;)
-  DIV_ONLY(int phase SHM_ONLY(, shmid);)
-} Slice;
+coreinit cubafun_;
+extern int cubaverb_;
+extern corespec cubaworkers_;
 
-workerini cubaini;
+/*********************************************************************/
 
-#if defined HAVE_SHMGET && (defined SUAVE || defined DIVONNE)
-#define FRAMECOPY
-#endif
+static inline void Child(cint fd, cint core)
+{
+  dispatch d;
 
-#ifdef DEBUG
-#define TERM_RED "\e[31m"
-#define TERM_BLUE "\e[34m"
-#define TERM_RESET "\e[0m\n"
-#define MASTER(s, ...) \
-fprintf(stderr, TERM_RED ROUTINE " master %d(%d): " s TERM_RESET, core, getpid(), ##__VA_ARGS__)
-#define WORKER(s, ...) \
-fprintf(stderr, TERM_BLUE ROUTINE " worker %d(%d): " s TERM_RESET, core, getpid(), ##__VA_ARGS__)
+  while( readsock(fd, &d, sizeof d) == sizeof d ) {
+    if( d.thissize ) {
+      MemAlloc(d.thisptr, d.thissize);
+      WORKER("reading This (%lu)", d.thissize);
+      readsock(fd, d.thisptr, d.thissize);
+    }
+    WORKER("running %p on fd %d", d.thisptr, fd);
+    d.worker(d.thisptr, d.thissize, core, fd);
+    if( d.thissize ) free(d.thisptr);
+  }
+}
+
+/*********************************************************************/
+
+Extern void SUFFIX(cubafork)(Spin **pspin)
+{
+  char out[128];
+  int cores, core;
+  fdpid *pfp;
+  Spin *spin;
+
+  VerboseInit();
+
+  EnvInit(cubaworkers_.paccel, "CUBAACCELMAX", 1000);
+  EnvInit(cubaworkers_.pcores, "CUBACORESMAX", 10000);
+  EnvInit(cubaworkers_.naccel, "CUBAACCEL", 0);
+  EnvInit(cubaworkers_.ncores, "CUBACORES", -sysconf(_SC_NPROCESSORS_ONLN));
+
+#ifdef HAVE_GETLOADAVG
+  if( cubaworkers_.ncores < 0 ) {
+    static int load = uninitialized;
+    if( load == uninitialized ) {
+      double loadavg;
+      getloadavg(&loadavg, 1);
+      load = floor(loadavg);
+    }
+    cubaworkers_.ncores = IMax(-cubaworkers_.ncores - load, 0);
+  }
 #else
-#define MASTER(s, ...)
-#define WORKER(s, ...)
+  cubaworkers_.ncores = abs(cubaworkers_.ncores);
 #endif
 
-/*********************************************************************/
-
-#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 void DoSample(This *t, number n, creal *x, real *f
-  VES_ONLY(, creal *w, ccount iter))
-{
-  cint ncores = IMin(t->ncores, n/MINSLICE);
-
-  if( ncores < MINCORES ) DoSampleSerial(t, n, x, f VES_ONLY(, w, iter));
-  else {
-    Slice slice;
-    int core, abort;
-    number nx;
-    char s[128];
-
-    t->neval += n;
-
-    nx = n % ncores;
-    slice.m = slice.n = n/ncores + 1;
-    if( VERBOSE > 2 ) {
-      sprintf(s, "sampling " NUMBER " points each on %d cores",
-        slice.n, ncores);
-      Print(s);
-    }
-
-    slice.i = 0;
-    VES_ONLY(slice.iter = iter;)
-    DIV_ONLY(slice.phase = t->phase;)
-
-#ifdef DIVONNE
-    if( n > t->nframe ) {
-      FrameFree(t, ShmRm(t));
-      t->nframe = n;
-      FrameAlloc(t);
-    }
-    SHM_ONLY(slice.shmid = t->shmid;)
-#endif
-
-    SHM_ONLY(if( t->shmid != -1 ) {
-      slice.m = n;
-#ifdef FRAMECOPY
-      VES_ONLY(Copy(t->frame, w, n);)
-      Copy(t->frame + n*NW, x, n*t->ndim);
-#endif
-    })
-
-    for( core = 0; core < ncores; ++core ) {
-      cint fd = t->child[core];
-      slice.n -= (core == nx);
-      MASTER("sending samples (sli:%lu[+" VES_ONLY(NUMBER "w:%lu+")
-        NUMBER "x:%lu]) to fd %d",
-        sizeof slice, VES_ONLY(slice.n, sizeof *w,)
-        slice.n, t->ndim*sizeof *x, fd);
-      writesock(fd, &slice, sizeof slice);
-      SHM_ONLY(if( t->shmid == -1 )) {
-        VES_ONLY(writesock(fd, w, slice.n*sizeof *w);
-                 w += slice.n;)
-        writesock(fd, x, slice.n*t->ndim*sizeof *x);
-        x += slice.n*t->ndim;
-      }
-      slice.i += slice.n;
-    }
-
-    abort = 0;
-    for( core = 0; core < ncores; ++core ) {
-      cint fd = t->child[core];
-      readsock(fd, &slice, sizeof slice);
-      MASTER("reading samples (sli:%lu[+" NUMBER "f:%lu]) from fd %d",
-        sizeof slice, slice.n, t->ncomp*sizeof *f, fd);
-      if( slice.n == -1 ) abort = 1;
-      else SHM_ONLY(if( t->shmid == -1 )) readsock(fd,
-        f + slice.i*t->ncomp, slice.n*t->ncomp*sizeof *f);
-    }
-    if( abort ) longjmp(t->abort, -99);
-
-#ifdef FRAMECOPY
-    if( t->shmid != -1 )
-      Copy(f, t->frame + slice.m*(NW + t->ndim), slice.m*t->ncomp);
-#endif
-  }
-}
-
-/*********************************************************************/
-
-#ifdef DIVONNE
-
-static inline int ReadyCore(cThis *t)
-{
-  fd_set ready;
-  int core, n = 0;
-
-  FD_ZERO(&ready);
-  for( core = 0; core < t->ncores; ++core ) {
-    FD_SET(t->child[core], &ready);
-    n = IMax(n, t->child[core]);
-  }
-  select(n + 1, &ready, NULL, NULL, NULL);
-
-  for( core = 0; core < t->ncores; ++core )
-    if( FD_ISSET(t->child[core], &ready) ) break;
-
-  return core;
-}
-
-/*********************************************************************/
-
-typedef struct {
-  number neval, neval_opt, neval_cut;
-  count nregions, iregion, retval;
-} ExploreResult;
-
-static int Explore(This *t, cint iregion)
-{
-  csize_t regionsize = RegionSize;
-  Region *region;
-  int ireg = iregion, core = t->running;
-  Vector(Totals, totals, NCOMP);
-
-  if( t->ncores < MINCORES ) return ExploreSerial(t, iregion);
-
-  if( t->running >= ((iregion < 0) ? 1 : t->ncores) ) {
-    cint fd = t->child[core = ReadyCore(t)];
-    ExploreResult res;
-    count comp, succ;
-
-    --t->running;
-    MASTER("reading res + region (res:%lu+reg:%lu) from fd %d",
-      sizeof res, regionsize, fd);
-    readsock(fd, &res, sizeof res);
-    ireg = res.iregion;
-    region = RegionPtr(ireg);
-    succ = ireg + region->next;
-    readsock(fd, region, regionsize);
-    if( --res.nregions > 0 ) {
-      region->next = t->nregions - ireg;
-      EnlargeRegions(t, res.nregions);
-      MASTER("reading regions (%dreg:%lu) from fd %d",
-        res.nregions, regionsize, fd);
-      readsock(fd, RegionPtr(t->nregions), res.nregions*regionsize);
-      t->nregions += res.nregions;
-
-      RegionPtr(t->nregions-1)->next = succ - t->nregions + 1;
-    }
-
-    MASTER("reading totals (tot:%lu) from fd %d",
-      t->ncomp*sizeof(Totals), fd);
-    readsock(fd, totals, t->ncomp*sizeof(Totals));
-    for( comp = 0; comp < t->ncomp; ++comp )
-      t->totals[comp].secondspread =
-        Max(t->totals[comp].secondspread, totals[comp].secondspread);
-
-    t->neval += res.neval;
-    t->neval_opt += res.neval_opt;
-    t->neval_cut += res.neval_cut;
-
-    if( res.retval == -1 ) return -1;
+  cores = cubaworkers_.naccel + cubaworkers_.ncores;
+  if( cores < MINCORES ) {
+    *pspin = NULL;
+    return;
   }
 
-  if( iregion >= 0 ) {
-    Slice slice;
-    cint fd = t->child[core];
-    slice.n = 0;
-    slice.i = iregion;
-    slice.phase = t->phase;
-    region = RegionPtr(iregion);
-    MASTER("writing region (sli:%lu+sam:%lu+reg:%lu+tot:%lu) to fd %d",
-      sizeof slice, sizeof(Samples), regionsize,
-      t->ncomp*sizeof(Totals), fd);
-    writesock(fd, &slice, sizeof slice);
-    writesock(fd, &t->samples[region->isamples], sizeof(Samples));
-    writesock(fd, region, regionsize);
-    writesock(fd, t->totals, t->ncomp*sizeof(Totals));
-    region->depth = 0;
-    ++t->running;
-  }
-
-  return ireg;
-}
-#endif
-
-/*********************************************************************/
-
-static void DoChild(This *t, cint core, cint fd)
-{
-  Slice slice;
-
-#ifdef DIVONNE
-  csize_t regionsize = RegionSize;
-  Vector(Totals, totals, NCOMP);
-  ExploreResult res;
-
-  t->totals = totals;
-  t->ncores = 0;	/* no recursive forks */
-  t->size = 2*t->ndim + 2;
-  AllocRegions(t);
-#endif
-
-#ifdef SUAVE
-  SHM_ONLY(if( t->shmid == -1 ))
-    MemAlloc(t->frame, t->nframe*SAMPLESIZE);
-#endif
-
-  if( cubaini.initfun ) cubaini.initfun(cubaini.initarg);
-
-  while( readsock(fd, &slice, sizeof slice) ) {
-    number n = slice.n;
-    DIV_ONLY(t->phase = slice.phase;)
-    if( n > 0 ) {
-      real VES_ONLY(*w,) *x, *f;
-      WORKER("reading samples (sli:%lu[+" VES_ONLY(NUMBER "w:%lu+")
-        NUMBER "x:%lu]) from fd %d",
-        sizeof slice, VES_ONLY(n, sizeof *w,) n, t->ndim*sizeof *x, fd);
-
-#ifdef DIVONNE
-      if( slice.m > t->nframe ) {
-        FrameFree(t);
-        t->nframe = slice.m;
-        SHM_ONLY(t->shmid = slice.shmid; ShmMap(t) else)
-        MemAlloc(t->frame, t->nframe*SAMPLESIZE);
-      }
-#endif
-
-      VES_ONLY(w = t->frame;)
-      x = t->frame + slice.m*NW;
-      f = x + slice.m*t->ndim;
-
-      SHM_ONLY(if( t->shmid != -1 ) {
-        VES_ONLY(w += slice.i;)
-        x += slice.i*t->ndim;
-        f += slice.i*t->ncomp;
-      }
-      else) {
-        VES_ONLY(readsock(fd, w, n*sizeof *w);)
-        readsock(fd, x, n*t->ndim*sizeof *x);
-      }
-
-      slice.n |= SampleRaw(t, n, x, f VES_ONLY(, w, slice.iter));
-      WORKER("writing samples (sli:%lu[+" NUMBER "f:%lu]) to fd %d",
-        sizeof slice, slice.n, t->ncomp*sizeof *f, fd);
-      writesock(fd, &slice, sizeof slice);
-      if( SHM_ONLY(t->shmid == -1 &&) slice.n != -1 )
-        writesock(fd, f, slice.n*t->ncomp*sizeof *f);
-    }
-#ifdef DIVONNE
-    else {
-      Samples *samples, psamples;
-
-      WORKER("reading region (sli:%lu+sam:%lu+reg:%lu+tot:%lu) from fd %d",
-        sizeof slice, sizeof(Samples), regionsize,
-        t->ncomp*sizeof(Totals), fd);
-      readsock(fd, &psamples, sizeof(Samples));
-      readsock(fd, t->region, regionsize);
-      readsock(fd, totals, t->ncomp*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);
-      }
-
-      res.retval = ExploreSerial(t, 0);
-      res.neval = t->neval;
-      res.neval_opt = t->neval_opt;
-      res.neval_cut = t->neval_cut;
-      res.nregions = t->nregions;
-      res.iregion = slice.i;
-      WORKER("writing regions (res:%lu+%dreg:%lu+tot:%lu) to fd %d",
-        sizeof res, t->nregions, regionsize,
-        t->ncomp*sizeof(Totals), fd);
-      writesock(fd, &res, sizeof res);
-      writesock(fd, t->region, t->nregions*regionsize);
-      writesock(fd, totals, t->ncomp*sizeof(Totals));
-    }
-#endif
-  }
-
-  WORKER("wrapping up");
-  if( cubaini.exitfun ) cubaini.exitfun(cubaini.exitarg);
-
-  exit(0);
-}
-
-/*********************************************************************/
-
-#ifdef HAVE_GETLOADAVG
-double cubaloadavg_;
-#endif
-
-static inline void ForkCores(This *t)
-{
-  int core;
-  char s[128];
-  cchar *env = getenv("CUBACORES");
-
-  t->ncores = env ? atoi(env) : sysconf(_SC_NPROCESSORS_ONLN);
-#ifdef HAVE_GETLOADAVG
-  if( env == NULL || t->ncores < 0 ) {
-    if( cubaloadavg_ < 0 ) getloadavg(&cubaloadavg_, 1);
-    t->ncores = abs(t->ncores) - floor(cubaloadavg_);
-  }
-#endif
-
-  DIV_ONLY(t->running = 0;)
-
-  if( t->ncores < MINCORES ) return;
-  if( VERBOSE ) {
-    sprintf(s, "using %d cores via "
+  if( cubaverb_ ) {
+    sprintf(out, "using %d cores %d accelerators via "
 #ifdef HAVE_SHMGET
       "shared memory",
 #else
       "pipes",
 #endif
-      t->ncores);
-    Print(s);
+      cubaworkers_.ncores, cubaworkers_.naccel);
+    Print(out);
   }
 
   fflush(NULL);		/* make sure all buffers are flushed,
 			   or else buffered content will be written
 			   out multiply, at each child's exit(0) */
 
-  Alloc(t->child, t->ncores);
-  for( core = 0; core < t->ncores; ++core ) {
+  MemAlloc(spin, sizeof *spin + cores*sizeof *spin->fp);
+  spin->spec = cubaworkers_;
+  pfp = spin->fp;
+  for( core = -spin->spec.naccel; core < spin->spec.ncores; ++core ) {
     int fd[2];
     pid_t pid;
     assert(
@@ -393,32 +100,65 @@ static inline void ForkCores(This *t)
       (pid = fork()) != -1 );
     if( pid == 0 ) {
       close(fd[0]);
-      DoChild(t, core, fd[1]);
+      free(spin);
+      Child(fd[1], core);
+      exit(0);
     }
     MASTER("forked pid %d pipe %d(master) -> %d(worker)",
       pid, fd[0], fd[1]);
     close(fd[1]);
-    t->child[core] = fd[0];
+    pfp->fd = fd[0];
+    pfp->pid = pid;
+    ++pfp;
   }
+
+  *pspin = spin;
 }
 
 /*********************************************************************/
 
-static inline void WaitCores(cThis *t)
+Extern void SUFFIX(cubawait)(Spin **pspin)
 {
-  if( t->ncores >= MINCORES ) {
-    int core;
-    pid_t pid;
-    for( core = 0; core < t->ncores; ++core ) {
-      MASTER("closing fd %d", t->child[core]);
-      close(t->child[core]);
-    }
-    free(t->child);
-    for( core = 0; core < t->ncores; ++core ) {
-      MASTER("waiting for child");
-      wait(&pid);
-      MASTER("pid %d terminated", pid);
-    }
+  int cores, core, status;
+  Spin *spin;
+
+  MasterExit();
+
+  if( Invalid(pspin) || (spin = *pspin) == NULL ) return;
+
+  cores = spin->spec.naccel + spin->spec.ncores;
+
+  for( core = 0; core < cores; ++core ) {
+    MASTER("closing fd %d", spin->fp[core].fd);
+    close(spin->fp[core].fd);
   }
+
+#ifdef KILL_WORKERS
+  for( core = 0; core < cores; ++core ) {
+    MASTER("killing pid %d", spin->fp[core].pid);
+    kill(spin->fp[core].pid, SIGKILL);
+  }
+#endif
+
+  for( core = 0; core < cores; ++core ) {
+    DEB_ONLY(pid_t pid;)
+    MASTER("waiting for child");
+    DEB_ONLY(pid =) wait(&status);
+    MASTER("pid %d terminated with exit code %d", pid, status);
+  }
+
+  free(spin);
+  *pspin = NULL;
 }
 
+#else
+
+Extern void SUFFIX(cubafork)(Spin **pspin) {}
+
+Extern void SUFFIX(cubawait)(Spin **pspin)
+{
+  MasterExit();
+}
+
+#endif
+
diff --git a/src/external/libCuba/src/common/MSample.c b/src/external/libCuba/src/common/MSample.c
index 4c146346..e0c48142 100644
--- a/src/external/libCuba/src/common/MSample.c
+++ b/src/external/libCuba/src/common/MSample.c
@@ -3,7 +3,7 @@
 		the sampling routine for the
 		Mathematica versions of the Cuba routines
 		by Thomas Hahn
-		last modified 19 Mar 12 th
+		last modified 13 Mar 15 th
 */
 
 
@@ -11,20 +11,20 @@ static void DoSample(This *t, cnumber n, real *x, real *f
   VES_ONLY(, real *w, ccount iter))
 {
   real *mma_f;
-  long mma_n;
+  int mma_n;
 
   if( MLAbort ) longjmp(t->abort, -99);
 
   MLPutFunction(stdlink, "EvaluatePacket", 1);
   MLPutFunction(stdlink, "Cuba`" ROUTINE "`sample", 1 VES_ONLY(+2) DIV_ONLY(+1));
-  MLPutRealList(stdlink, x, n*t->ndim);
-  VES_ONLY(MLPutRealList(stdlink, w, n);
+  MLPutRealxList(stdlink, x, n*t->ndim);
+  VES_ONLY(MLPutRealxList(stdlink, w, n);
            MLPutInteger(stdlink, iter);)
   DIV_ONLY(MLPutInteger(stdlink, t->phase);)
   MLEndPacket(stdlink);
 
   MLNextPacket(stdlink);
-  if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) {
+  if( !MLGetRealxList(stdlink, &mma_f, &mma_n) ) {
     MLClearError(stdlink);
     MLNewPacket(stdlink);
     longjmp(t->abort, -99);
@@ -33,12 +33,12 @@ static void DoSample(This *t, cnumber n, real *x, real *f
   t->neval += mma_n;
 
   if( mma_n != n*t->ncomp ) {
-    MLDisownRealList(stdlink, mma_f, mma_n);
+    MLReleaseRealxList(stdlink, mma_f, mma_n);
     longjmp(t->abort, -3);
   }
  
   Copy(f, mma_f, n*t->ncomp);
-  MLDisownRealList(stdlink, mma_f, mma_n);
+  MLReleaseRealxList(stdlink, mma_f, mma_n);
 }
 
 /*********************************************************************/
@@ -50,16 +50,16 @@ static count SampleExtra(This *t, cBounds *b)
 {
   count n, nget;
   real *mma_f;
-  long mma_n;
+  int mma_n;
 
   MLPutFunction(stdlink, "EvaluatePacket", 1);
   MLPutFunction(stdlink, "Cuba`Divonne`findpeak", 2);
-  MLPutRealList(stdlink, (real *)b, 2*t->ndim);
+  MLPutRealxList(stdlink, (real *)b, 2*t->ndim);
   MLPutInteger(stdlink, t->phase);
   MLEndPacket(stdlink);
 
   MLNextPacket(stdlink);
-  if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) {
+  if( !MLGetRealxList(stdlink, &mma_f, &mma_n) ) {
     MLClearError(stdlink);
     MLNewPacket(stdlink);
     longjmp(t->abort, -99);
@@ -73,7 +73,7 @@ static count SampleExtra(This *t, cBounds *b)
     Copy(t->fextra, mma_f + nget*t->ndim, n*t->ncomp);
   }
 
-  MLDisownRealList(stdlink, mma_f, mma_n);
+  MLReleaseRealxList(stdlink, mma_f, mma_n);
 
   return n;
 }
diff --git a/src/external/libCuba/src/common/Makefile.am b/src/external/libCuba/src/common/Makefile.am
index cdf87141..d7f70508 100644
--- a/src/external/libCuba/src/common/Makefile.am
+++ b/src/external/libCuba/src/common/Makefile.am
@@ -1,12 +1,14 @@
 ## Process this file with automake to create Makefile.in
 
-c_sources = WorkerIni.c
+c_sources = \
+  Global.c \
+  Data.c
 
 AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE
 AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS)
 AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS)
 
-noinst_LTLIBRARIES = libworkerini.la
+noinst_LTLIBRARIES = libcommon.la
 
-libworkerini_la_SOURCES = $(c_sources)
-libworkerini_la_LDFLAGS = $(AM_LDFLAGS)
+libcommon_la_SOURCES = $(c_sources)
+libcommon_la_LDFLAGS = $(AM_LDFLAGS)
diff --git a/src/external/libCuba/src/common/Random.c b/src/external/libCuba/src/common/Random.c
index 4c2777c2..6d606157 100644
--- a/src/external/libCuba/src/common/Random.c
+++ b/src/external/libCuba/src/common/Random.c
@@ -1,7 +1,7 @@
 /*
 	Random.c
 		quasi- and pseudo-random-number generation
-		last modified 7 Aug 13 th
+		last modified 18 Mar 14 th
 */
 
 
@@ -88,14 +88,13 @@ static inline void SobolIni(This *t)
     299,   1,   3,   3,   9,   9,  25, 107,  39 };
 
   count dim, bit, nbits;
-  number max, *pini = ini;
-  cnumber nmax = 2*t->maxeval;
+  number *pini = ini, max;
 
-  for( nbits = 0, max = 1; max <= nmax; max <<= 1 ) ++nbits;
-  t->rng.sobol.norm = 1./max;
+  for( nbits = 0, max = t->maxeval; max; max >>= 1 ) ++nbits;
+  t->rng.sobol.norm = ldexp(.5, -nbits);
 
-  for( bit = 0; bit < nbits; ++bit )
-    t->rng.sobol.v[0][bit] = (max >>= 1);
+  for( bit = 0; bit <= nbits; ++bit )
+    t->rng.sobol.v[0][bit] = (number)1 << (nbits - bit);
 
   for( dim = 1; dim < t->ndim; ++dim ) {
     number *pv = t->rng.sobol.v[dim], *pvv = pv;
@@ -103,10 +102,10 @@ static inline void SobolIni(This *t)
     int inibits = -1, bit;
     for( j = powers; j; j >>= 1 ) ++inibits;
 
-    memcpy(pv, pini, inibits*sizeof(*pini));
+    memcpy(pv, pini, inibits*sizeof *pini);
     pini += 8;
 
-    for( bit = inibits; bit < nbits; ++bit ) {
+    for( bit = inibits; bit <= nbits; ++bit ) {
       number newv = *pvv, j = powers;
       int b;
       for( b = 0; b < inibits; ++b ) {
@@ -117,8 +116,8 @@ static inline void SobolIni(This *t)
       ++pvv;
     }
 
-    for( bit = 0; bit < nbits - 1; ++bit )
-      pv[bit] <<= nbits - bit - 1;
+    for( bit = 0; bit < nbits; ++bit )
+      pv[bit] <<= nbits - bit;
   }
 
   t->rng.sobol.seq = 0;
diff --git a/src/external/libCuba/src/common/WorkerIni.c b/src/external/libCuba/src/common/WorkerIni.c
deleted file mode 100644
index 52c9978a..00000000
--- a/src/external/libCuba/src/common/WorkerIni.c
+++ /dev/null
@@ -1,37 +0,0 @@
-/*
-	WorkerIni.c
-		set/run the init/exit functions for worker processes
-		by Thomas Hahn
-		last modified 6 Sep 12 th
-*/
-
-
-#include "stddecl.h"
-
-extern workerini cubaini;
-
-Extern void SUFFIX(cubasetinit)(subroutine f, void *arg)
-{
-  cubaini.initfun = f;
-  cubaini.initarg = arg;
-}
-
-
-Extern void SUFFIX(cubasetexit)(subroutine f, void *arg)
-{
-  cubaini.exitfun = f;
-  cubaini.exitarg = arg;
-}
-
-
-Extern void SUFFIX(cubaruninit)()
-{
-  if( cubaini.initfun ) cubaini.initfun(cubaini.initarg);
-}
-
-
-Extern void SUFFIX(cubarunexit)()
-{
-  if( cubaini.exitfun ) cubaini.exitfun(cubaini.exitarg);
-}
-
diff --git a/src/external/libCuba/src/common/stddecl.h b/src/external/libCuba/src/common/stddecl.h
index 6118aa0a..e5d046ff 100644
--- a/src/external/libCuba/src/common/stddecl.h
+++ b/src/external/libCuba/src/common/stddecl.h
@@ -1,7 +1,7 @@
 /*
 	stddecl.h
 		declarations common to all Cuba routines
-		last modified 17 Sep 13 th
+		last modified 23 Apr 15 th
 */
 
 
@@ -13,7 +13,6 @@
 #endif
 
 #define _DEFAULT_SOURCE
-#define _XOPEN_SOURCE
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -30,6 +29,7 @@
 #ifdef HAVE_FORK
 #include <sys/wait.h>
 #include <sys/socket.h>
+#include <signal.h>
 #ifdef HAVE_SHMGET
 #include <sys/ipc.h>
 #include <sys/shm.h>
@@ -83,10 +83,31 @@ void *alloca (size_t);
 
 #define SAMPLESIZE (NW + t->ndim + t->ncomp)*sizeof(real)
 
+
+enum { uninitialized = 0x61627563 };
+
+#define EnvInit(var, name, default) \
+  if( var == uninitialized ) { \
+    cchar *env = getenv(name); \
+    if( env == NULL ) var = default; \
+    else { \
+      var = atoi(env); \
+      if( cubaverb_ ) { \
+        char out[64]; \
+        sprintf(out, "env " name " = %d", (int)var); \
+        Print(out); \
+      } \
+    } \
+  }
+
+#define VerboseInit() EnvInit(cubaverb_, "CUBAVERBOSE", 0)
+#define MaxVerbose(flags) (flags + IDim(IMin(cubaverb_, 3) - ((flags) & 3)))
+
 #define VERBOSE (t->flags & 3)
 #define LAST (t->flags & 4)
 #define SHARPEDGES (t->flags & 8)
 #define KEEPFILE (t->flags & 16)
+#define ZAPSTATE (t->flags & 32)
 #define REGIONS (t->flags & 128)
 #define RNG (t->flags >> 8)
 
@@ -120,7 +141,7 @@ void *alloca (size_t);
 
 #define Zap(d) memset(d, 0, sizeof(d))
 
-#define MaxErr(avg) Max(t->epsrel*fabs(avg), t->epsabs)
+#define MaxErr(avg) Max(t->epsrel*fabsx(avg), t->epsabs)
 
 #ifdef __cplusplus
 #define mallocset(p, n) (*(void **)&p = malloc(n))
@@ -157,8 +178,21 @@ void *alloca (size_t);
 
 #ifdef MLVERSION
 #define ML_ONLY(...) __VA_ARGS__
+#define ML_NOT(...)
 #else
 #define ML_ONLY(...)
+#define ML_NOT(...) __VA_ARGS__
+
+#define CORE_MASTER (int []){32768}
+#define MasterInit() do if( !cubafun_.init ) { \
+  cubafun_.init = true; \
+  if( cubafun_.initfun ) cubafun_.initfun(cubafun_.initarg, CORE_MASTER); \
+} while( 0 )
+#define MasterExit() do if( cubafun_.init ) { \
+  cubafun_.init = false; \
+  if( cubafun_.exitfun ) cubafun_.exitfun(cubafun_.exitarg, CORE_MASTER); \
+} while( 0 )
+#define Invalid(s) ((s) == NULL || *(int *)(s) == -1)
 
 #ifdef HAVE_FORK
 #undef FORK_ONLY
@@ -168,37 +202,40 @@ void *alloca (size_t);
 #undef SHM_ONLY
 #define SHM_ONLY(...) __VA_ARGS__
 
-#define ShmMap(t, ...) if( t->shmid != -1 ) { \
-  t->frame = shmat(t->shmid, NULL, 0); \
-  if( t->frame == (void *)-1 ) Abort("shmat"); \
-  __VA_ARGS__ \
-}
-
-#define ShmRm(t) shmctl(t->shmid, IPC_RMID, NULL);
+#define MasterAlloc(t) \
+  t->shmid = shmget(IPC_PRIVATE, t->nframe*SAMPLESIZE, IPC_CREAT | 0600)
+#define MasterFree(t) shmctl(t->shmid, IPC_RMID, NULL)
+#define WorkerAlloc(t)
+#define WorkerFree(r)
 
 #undef ShmAlloc
-#define ShmAlloc(t, ...) \
-  t->shmid = shmget(IPC_PRIVATE, t->nframe*SAMPLESIZE, IPC_CREAT | 0600); \
-  ShmMap(t, __VA_ARGS__)
+#define ShmAlloc(t, who) \
+  who##Alloc(t); \
+  if( t->shmid != -1 ) { \
+    t->frame = shmat(t->shmid, NULL, 0); \
+    if( t->frame == (void *)-1 ) Abort("shmat"); \
+  }
 
 #undef ShmFree
-#define ShmFree(t, ...) if( t->shmid != -1 ) { \
-  shmdt(t->frame); \
-  __VA_ARGS__ \
-}
+#define ShmFree(t, who) \
+  if( t->shmid != -1 ) { \
+    shmdt(t->frame); \
+    who##Free(t); \
+  }
 
 #endif
 #endif
 #endif
   
-#define FrameAlloc(t, ...) \
-  SHM_ONLY(ShmAlloc(t, __VA_ARGS__) else) \
+#define FrameAlloc(t, who) \
+  SHM_ONLY(ShmAlloc(t, who) else) \
   MemAlloc(t->frame, t->nframe*SAMPLESIZE);
 
-#define FrameFree(t, ...) DIV_ONLY(if( t->nframe )) { \
-  SHM_ONLY(ShmFree(t, __VA_ARGS__) else) \
-  free(t->frame); \
-}
+#define FrameFree(t, who) \
+  DIV_ONLY(if( t->nframe )) { \
+    SHM_ONLY(ShmFree(t, who) else) \
+    free(t->frame); \
+  }
 
 
 #define StateDecl \
@@ -219,7 +256,9 @@ struct stat st
 
 typedef long long int signature_t;
 
-#define StateSignature(t, i) (0x41425543 + \
+enum { signature = 0x41425543 };
+
+#define StateSignature(t, i) (signature + \
   ((signature_t)(i) << 60) + \
   ((signature_t)(t)->ncomp << 48) + \
   ((signature_t)(t)->ndim << 32))
@@ -306,39 +345,100 @@ typedef const count ccount;
 #define PREFIX(s) ll##s
 #define NUMBER "%lld"
 #define NUMBER7 "%7lld"
+#define NUMBER_MAX LLONG_MAX
 typedef long long int number;
 #else
 #define PREFIX(s) s
 #define NUMBER "%d"
 #define NUMBER7 "%7d"
+#define NUMBER_MAX INT_MAX
 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. */
+#define SHOW(r) (double)(r)
+	/* floating-point numbers are printed with SHOW */
+
+#if REALSIZE == 16
+#include <quadmath.h>
+typedef __float128 real;
+#define RC(x) x##Q
+#define sqrtx sqrtq
+#define expx expq
+#define powx powq
+#define erfx erfq
+#define fabsx fabsq
+#define ldexpx ldexpq
+#define REAL_MAX_EXP FLT128_MAX_EXP
+#define REAL_MAX FLT128_MAX
+#elif REALSIZE == 10
+typedef long double real;
+#define RC(x) x##L
+#define sqrtx sqrtl
+#define expx expl
+#define powx powl
+#define erfx erfl
+#define fabsx fabsl
+#define ldexpx ldexpl
+#define REAL_MAX_EXP LDBL_MAX_EXP
+#define REAL_MAX LDBL_MAX
+#define MLPutRealxList MLPutReal128List
+#define MLGetRealxList MLGetReal128List
+#define MLReleaseRealxList MLReleaseReal128List
+#else
+typedef double real;
+#define RC(x) x
+#define sqrtx sqrt
+#define expx exp
+#define powx pow
+#define erfx erf
+#define fabsx fabs
+#define ldexpx ldexp
+#define REAL_MAX_EXP DBL_MAX_EXP
+#define REAL_MAX DBL_MAX
+#define MLPutRealxList MLPutReal64List
+#define MLGetRealxList MLGetReal64List
+#define MLReleaseRealxList MLReleaseReal64List
+#endif
 
 typedef const real creal;
 
-typedef void (*subroutine)();
+typedef void (*subroutine)(void *, cint *);
 
 typedef struct {
   subroutine initfun;
   void *initarg;
   subroutine exitfun;
   void *exitarg;
-} workerini;
+  bool init;
+} coreinit;
+
+typedef struct {
+  int ncores, naccel;
+  int pcores, paccel;
+} corespec;
+
+typedef struct {
+  int fd, pid;
+} fdpid;
+
+typedef struct {
+  corespec spec;
+  fdpid fp[];
+} Spin;
 
 
 struct _this;
 
+typedef struct {
+  void (*worker)(struct _this *, csize_t, cint, cint);
+  struct _this *thisptr;
+  size_t thissize;
+} dispatch;
+
+
 typedef unsigned int state_t;
 
 #define SOBOL_MINDIM 1
@@ -407,7 +507,7 @@ static inline real Max(creal a, creal b) {
 }
 
 static inline real Weight(creal sum, creal sqsum, cnumber n) {
-  creal w = sqrt(sqsum*n);
+  creal w = sqrtx(sqsum*n);
   return (n - 1)/Max((w + sum)*(w - sum), NOTZERO);
 }
 
@@ -452,7 +552,7 @@ static inline void Print(MLCONST char *s)
 
 #else
 
-#define Print(s) puts(s); fflush(stdout)
+#define Print(s) { puts(s); fflush(stdout); }
 
 #endif
 
diff --git a/src/external/libCuba/src/cuba.h b/src/external/libCuba/src/cuba.h
index 48b1c7ac..5e87ef1a 100644
--- a/src/external/libCuba/src/cuba.h
+++ b/src/external/libCuba/src/cuba.h
@@ -2,116 +2,122 @@
 	cuba.h
 		Prototypes for the Cuba library
 		this file is part of Cuba
-		last modified 30 Apr 13 th
+		last modified 13 Mar 15 th
 */
 
+typedef double cubareal;
+
+	/* integrand_t is intentionally a minimalistic integrand type.
+	   It includes neither the nvec and core arguments nor the
+	   extra arguments passed by Vegas/Suave (weight, iter) and
+	   Divonne (phase).
+	   In most cases, integrand_t is just what you want, otherwise
+	   simply use an explicit typecast to integrand_t in the Cuba
+	   invocation. */
+typedef int (*integrand_t)(const int *ndim, const cubareal x[],
+  const int *ncomp, cubareal f[], void *userdata);
+
+typedef void (*peakfinder_t)(const int *ndim, const cubareal b[],
+  int *n, cubareal x[], void *userdata);
+
 #ifdef __cplusplus
 extern "C" {
 #endif
 
-	/* NB: Divonne actually passes a fifth argument, a const int *
-	   which points to the integration phase.  This is used only
-	   rarely and most users are confused by the warnings the
-	   compiler emits if the `correct' prototype is used.  Thus,
-	   if you need to access this argument, use an explicit cast
-	   to integrand_t when invoking Divonne. */
-typedef int (*integrand_t)(const int *ndim, const double x[],
-  const int *ncomp, double f[], void *userdata);
-
-typedef void (*peakfinder_t)(const int *ndim, const double b[],
-  int *n, double x[]);
-
 void Vegas(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
-  const double epsrel, const double epsabs,
+  integrand_t integrand, void *userdata, const int nvec,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const int mineval, const int maxeval,
   const int nstart, const int nincrease, const int nbatch,
-  const int gridno, const char *statefile,
+  const int gridno, const char *statefile, void *spin,
   int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void llVegas(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
-  const double epsrel, const double epsabs,
+  integrand_t integrand, void *userdata, const long long int nvec,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const long long int mineval, const long long int maxeval,
   const long long int nstart, const long long int nincrease,
   const long long int nbatch,
-  const int gridno, const char *statefile,
+  const int gridno, const char *statefile, void *spin,
   long long int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void Suave(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
-  const double epsrel, const double epsabs,
+  integrand_t integrand, void *userdata, const int nvec,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const int mineval, const int maxeval,
-  const int nnew, const double flatness,
-  const char *statefile,
+  const int nnew, const int nmin,
+  const cubareal flatness, const char *statefile, void *spin,
   int *nregions, int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void llSuave(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
-  const double epsrel, const double epsabs,
+  integrand_t integrand, void *userdata, const long long int nvec,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const long long int mineval, const long long int maxeval,
-  const long long int nnew, const double flatness,
-  const char *statefile,
+  const long long int nnew, const long long int nmin,
+  const cubareal flatness, const char *statefile, void *spin,
   int *nregions, long long int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void Divonne(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
-  const double epsrel, const double epsabs,
+  integrand_t integrand, void *userdata, const int nvec,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const int mineval, const int maxeval,
   const int key1, const int key2, const int key3, const int maxpass,
-  const double border, const double maxchisq, const double mindeviation,
-  const int ngiven, const int ldxgiven, double xgiven[],
+  const cubareal border, const cubareal maxchisq, const cubareal mindeviation,
+  const int ngiven, const int ldxgiven, cubareal xgiven[],
   const int nextra, peakfinder_t peakfinder,
-  const char *statefile,
+  const char *statefile, void *spin,
   int *nregions, int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void llDivonne(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
-  const double epsrel, const double epsabs,
+  integrand_t integrand, void *userdata, const long long int nvec,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const long long int mineval, const long long int maxeval,
   const int key1, const int key2, const int key3, const int maxpass,
-  const double border, const double maxchisq, const double mindeviation,
-  const long long int ngiven, const int ldxgiven, double xgiven[],
-  const long long int nextra,
-  void (*peakfinder)(const int *, const double [], int *, double []),
-  const char *statefile,
+  const cubareal border, const cubareal maxchisq, const cubareal mindeviation,
+  const long long int ngiven, const int ldxgiven, cubareal xgiven[],
+  const long long int nextra, peakfinder_t peakfinder,
+  const char *statefile, void *spin,
   int *nregions, long long int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void Cuhre(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
-  const double epsrel, const double epsabs,
+  integrand_t integrand, void *userdata, const int nvec,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int mineval, const int maxeval,
   const int key,
-  const char *statefile,
+  const char *statefile, void *spin,
   int *nregions, int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void llCuhre(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
-  const double epsrel, const double epsabs,
+  integrand_t integrand, void *userdata, const long long int nvec,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags,
   const long long int mineval, const long long int maxeval,
   const int key,
-  const char *statefile,
+  const char *statefile, void *spin,
   int *nregions, long long int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
-void cubasetinit(void (*)(), void *);
-void cubasetexit(void (*)(), void *);
-void cubaruninit(void);
-void cubaruninit(void);
+void cubafork(void *pspin);
+void cubawait(void *pspin);
+
+void cubacores(const int n, const int p);
+void cubaaccel(const int n, const int p);
+
+void cubainit(void (*f)(), void *arg);
+void cubaexit(void (*f)(), void *arg);
 
 #ifdef __cplusplus
 }
diff --git a/src/external/libCuba/src/cuhre/Cuhre.c b/src/external/libCuba/src/cuhre/Cuhre.c
index 33e17f62..6dc0aaa1 100644
--- a/src/external/libCuba/src/cuhre/Cuhre.c
+++ b/src/external/libCuba/src/cuhre/Cuhre.c
@@ -2,7 +2,7 @@
 	Cuhre.c
 		Adaptive integration using cubature rules
 		by Thomas Hahn
-		last modified 17 Sep 13 th
+		last modified 22 Jul 14 th
 */
 
 
@@ -15,55 +15,70 @@
 /*********************************************************************/
 
 Extern void EXPORT(Cuhre)(ccount ndim, ccount ncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber nvec,
   creal epsrel, creal epsabs,
   cint flags, cnumber mineval, cnumber maxeval,
-  ccount key, cchar *statefile,
+  ccount key, cchar *statefile, Spin **pspin,
   count *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = ndim;
   t.ncomp = ncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = nvec;
   t.epsrel = epsrel;
   t.epsabs = epsabs;
-  t.flags = flags;
+  t.flags = MaxVerbose(flags);
   t.mineval = mineval;
   t.maxeval = maxeval;
   t.key = key;
   t.statefile = statefile;
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
  
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
 /*********************************************************************/
 
 Extern void EXPORT(cuhre)(ccount *pndim, ccount *pncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber *pnvec,
   creal *pepsrel, creal *pepsabs,
   cint *pflags, cnumber *pmineval, cnumber *pmaxeval,
-  ccount *pkey, cchar *statefile,
+  ccount *pkey, cchar *statefile, Spin **pspin,
   count *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob, cint statefilelen)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = *pndim;
   t.ncomp = *pncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = *pnvec;
   t.epsrel = *pepsrel;
   t.epsabs = *pepsabs;
-  t.flags = *pflags;
+  t.flags = MaxVerbose(*pflags);
   t.mineval = *pmineval;
   t.maxeval = *pmaxeval;
   t.key = *pkey;
   CString(t.statefile, statefile, statefilelen);
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
+
diff --git a/src/external/libCuba/src/cuhre/Integrate.c b/src/external/libCuba/src/cuhre/Integrate.c
index 25c40225..27fc7e2c 100644
--- a/src/external/libCuba/src/cuhre/Integrate.c
+++ b/src/external/libCuba/src/cuhre/Integrate.c
@@ -3,7 +3,7 @@
 		integrate over the unit hypercube
 		this file is part of Cuhre
 		checkpointing by B. Chokoufe
-		last modified 17 Sep 13 th
+		last modified 14 Mar 15 th
 */
 
 
@@ -11,6 +11,9 @@
 
 typedef struct pool {
   struct pool *next;
+#if REALSIZE > 8
+  void *dummy; /* for alignment */
+#endif
   char region[];
 } Pool;
 
@@ -42,12 +45,14 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   if( VERBOSE > 1 ) {
     sprintf(out, "Cuhre input parameters:\n"
       "  ndim " COUNT "\n  ncomp " COUNT "\n"
+      ML_NOT("  nvec " NUMBER "\n")
       "  epsrel " REAL "\n  epsabs " REAL "\n"
       "  flags %d\n  mineval " NUMBER "\n  maxeval " NUMBER "\n"
       "  key " COUNT "\n"
       "  statefile \"%s\"",
       t->ndim, t->ncomp,
-      t->epsrel, t->epsabs,
+      ML_NOT(t->nvec,)
+      SHOW(t->epsrel), SHOW(t->epsabs),
       t->flags, t->mineval, t->maxeval,
       t->key,
       t->statefile);
@@ -61,7 +66,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
 
   RuleAlloc(t);
   t->mineval = IMax(t->mineval, t->rule.n + 1);
-  FrameAlloc(t, ShmRm(t));
+  FrameAlloc(t, Master);
   ForkCores(t);
 
   if( (fail = setjmp(t->abort)) ) goto abort;
@@ -125,7 +130,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       for( tot = state->totals, comp = 0; tot < Tot; ++tot )
         oe += sprintf(oe, "\n[" COUNT "] "
           REAL " +- " REAL "  \tchisq " REAL " (" COUNT " df)",
-          ++comp, tot->avg, tot->err, tot->chisq, t->nregions - 1);
+          ++comp, SHOW(tot->avg), SHOW(tot->err),
+          SHOW(tot->chisq), t->nregions - 1);
       Print(out);
     }
 
@@ -188,7 +194,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
 
       tot->lastavg += diff = resL->avg + resR->avg - res->avg;
 
-      diff = fabs(.25*diff);
+      diff = fabsx(.25*diff);
       err = resL->err + resR->err;
       if( err > 0 ) {
         creal c = 1 + 2*diff/err;
@@ -213,7 +219,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       }
       else {
         tot->avg = avg;
-        tot->err = sqrt(sigsq);
+        tot->err = sqrtx(sigsq);
       }
     }
     ++t->nregions;
@@ -249,13 +255,13 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
         Result *Res;
 
         MLPutFunction(stdlink, "Cuba`Cuhre`region", 2);
-        MLPutRealList(stdlink, (real *)region->bounds, 2*t->ndim);
+        MLPutRealxList(stdlink, (real *)region->bounds, 2*t->ndim);
 
         MLPutFunction(stdlink, "List", t->ncomp);
         for( Res = (res = RegionResult(region)) + t->ncomp;
              res < Res; ++res ) {
           real r[] = {res->avg, res->err};
-          MLPutRealList(stdlink, r, Elements(r));
+          MLPutRealxList(stdlink, r, Elements(r));
         }
       }
   }
@@ -266,9 +272,7 @@ abort:
     cur = cur->next;
     free(pool);
   }
-
-  WaitCores(t);
-  FrameFree(t);
+  FrameFree(t, Master);
   RuleFree(t);
 
   StateRemove(t);
diff --git a/src/external/libCuba/src/cuhre/Rule.c b/src/external/libCuba/src/cuhre/Rule.c
index e3e931f3..129e9689 100644
--- a/src/external/libCuba/src/cuhre/Rule.c
+++ b/src/external/libCuba/src/cuhre/Rule.c
@@ -4,11 +4,12 @@
 		code lifted with minor modifications from DCUHRE
 		by J. Berntsen, T. Espelid, and A. Genz
 		this file is part of Cuhre
-		last modified 5 Aug 13 th
+		last modified 7 May 15 th
 */
 
 
 #define NextSet(p) p = (Set *)((char *)p + setsize)
+#define IndexSet(p, n) ((Set *)((char *)p + n*setsize))
 
 /*********************************************************************/
 
@@ -153,7 +154,7 @@ static void Rule13Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -298,7 +299,7 @@ static void Rule11Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -309,32 +310,55 @@ static void Rule11Alloc(This *t)
 static void Rule9Alloc(This *t)
 {
   static creal w[] = {
-    -.0023611709677855117884,   .11415390023857325268,
-    -.63833920076702389094,     .74849988504685208004,
-    -.0014324017033399125142,   .057471507864489725949,
-    -.14225104571434243234,    -.062875028738286979989,
-     .254591133248959089,     -1.207328566678236261,
-     .89567365764160676508,    -.36479356986049146661,
-     .0035417564516782676826,  -.072609367395893679605,
-     .10557491625218991012,     .0021486025550098687713,
-    -.032268563892953949998,    .010636783990231217481,
-     .014689102496143490175,    .51134708346467591431,
-     .45976448120806344646,     .18239678493024573331,
-    -.04508628929435784076,     .21415883524352793401,
-    -.027351546526545644722,    .054941067048711234101,
-     .11937596202570775297,     .65089519391920250593,
-     .14744939829434460168,     .057693384490973483573,
-     .034999626602143583822,  -1.3868627719278281436,
-    -.2386668732575008879,      .015532417276607053264,
-     .0035328099607090870236,   .09231719987444221619,
-     .02254314464717892038,     .013675773263272822361,
-    -.32544759695960125297,     .0017708782258391338413,
-     .0010743012775049343856,   .25150011495314791996 };
+    RC(-.002361170967785511788400941242259231309691),
+    RC(.1141539002385732526821323741697655347686),
+    RC(-.6383392007670238909386026193674701393074),
+    RC(.7484998850468520800423030047583803945205),
+    RC(-.001432401703339912514196154599769007103671),
+    RC(.05747150786448972594860897296200006759892),
+    RC(-.1422510457143424323449521620935950679394),
+    RC(-.06287502873828697998942424881040490136987),
+    RC(.2545911332489590890011611142429070613156),
+    RC(-1.207328566678236261002219995185143356737),
+    RC(.8956736576416067650809467826488567200939),
+    RC(-.3647935698604914666100134551377381205297),
+    RC(.003541756451678267682601411863388846964536),
+    RC(-.07260936739589367960492815865074633743652),
+    RC(.1055749162521899101218622863269817454540),
+    RC(.002148602555009868771294231899653510655506),
+    RC(-.03226856389295394999786630399875134318006),
+    RC(.01063678399023121748083624225818915724455),
+    RC(.01468910249614349017540783437728097691502),
+    RC(.5113470834646759143109387357149329909126),
+    RC(.4597644812080634464633352781605214342691),
+    RC(.1823967849302457333050067275688690602649),
+    RC(-.04508628929435784075980562738240804429658),
+    RC(.2141588352435279340097929526588394300172),
+    RC(-.02735154652654564472203690086290223507436),
+    RC(.05494106704871123410060080562462135546101),
+    RC(.1193759620257077529708962121565290178730),
+    RC(.6508951939192025059314756320878023215278),
+    RC(.1474493982943446016775696826942585013243),
+    RC(.05769338449097348357291272840392627722165),
+    RC(.03499962660214358382244159694487155861542),
+    RC(-1.386862771927828143599782668709014266770),
+    RC(-.2386668732575008878964134721962088068396),
+    RC(.01553241727660705326386197156586357005224),
+    RC(.003532809960709087023561817517751309380604),
+    RC(.09231719987444221619017126187763868745587),
+    RC(.02254314464717892037990281369120402214829),
+    RC(.01367577326327282236101845043145111753718),
+    RC(-.3254475969596012529657378160439011607639),
+    RC(.001770878225839133841300705931694423482268),
+    RC(.001074301277504934385647115949826755327753),
+    RC(.2515001149531479199576969952416196054795) };
 
   static creal g[] = {
-     .47795365790226950619,     .20302858736911986780,
-     .44762735462617812882,     .125,
-     .34303789878087814570 };
+    RC(.4779536579022695061928604197171830064732),
+    RC(.2030285873691198677998034402373279133258),
+    RC(.4476273546261781288207704806530998539285),
+    RC(.125),
+    RC(.3430378987808781457001426145164678603407) };
 
   enum { nsets = 9 };
 
@@ -440,7 +464,7 @@ static void Rule9Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -451,21 +475,33 @@ static void Rule9Alloc(This *t)
 static void Rule7Alloc(This *t)
 {
   static creal w[] = {
-     .019417866674748388428,   -.40385257701150182546,
-     .64485668767465982223,     .01177982690775806141,
-    -.18041318740733609012,    -.088785828081335044443,
-     .056328645808285941374,   -.0097089333373741942142,
-    -.99129176779582358138,    -.17757165616267008889,
-     .12359398032043233572,     .074978148702033690681,
-     .55489147051423559776,     .088041241522692771226,
-     .021118358455513385083,   -.0099302203239653333087,
-    -.064100053285010904179,    .030381729038221007659,
-     .0058899134538790307051,  -.0048544666686870971071,
-     .35514331232534017777 };
+    RC(.01941786667474838842844534313920462333850),
+    RC(-.4038525770115018254611834753723880293161),
+    RC(.6448566876746598222277360730193089551024),
+    RC(.01177982690775806141012214458820955067854),
+    RC(-.1804131874073360901182293138710989490609),
+    RC(-.08878582808133504444306598174517276122439),
+    RC(.05632864580828594137378124255408286479947),
+    RC(-.009708933337374194214222671569602311669249),
+    RC(-.9912917677958235813775106862002319060386),
+    RC(-.1775716561626700888861319634903455224488),
+    RC(.1235939803204323357183625846672135876752),
+    RC(.07497814870203369068087999555157339703666),
+    RC(.5548914705142355977605994477355651401434),
+    RC(.08804124152269277122645182458858273865209),
+    RC(.02111835845551338508329573367808085283304),
+    RC(-.009930220323965333308685820460105538586058),
+    RC(-.06410005328501090417895544042025034295870),
+    RC(.03038172903822100765927778829870429682489),
+    RC(.005889913453879030705061072294104775339268),
+    RC(-.004854466668687097107111335784801155834624),
+    RC(.3551433123253401777722639269806910448976) };
 
   static creal g[] = {
-     .47795365790226950619,     .20302858736911986780,
-     .375,                      .34303789878087814570 };
+    RC(.4779536579022695061928604197171830064732),
+    RC(.2030285873691198677998034402373279133258),
+    RC(.375),
+    RC(.3430378987808781457001426145164678603407) };
 
   enum { nsets = 6 };
 
@@ -541,7 +577,7 @@ static void Rule7Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -630,7 +666,8 @@ static void Sample(This *t, Region *region)
   Bounds *b, *B = region->bounds + t->ndim;
   Result *result = RegionResult(region), *res, *Res = result + t->ncomp;
   creal *errcoeff = t->rule.errcoeff;
-  creal ratio = Sq(first[2].gen[0]/first[1].gen[0]);
+  creal ratio = Sq(IndexSet(first,2)->gen[0]/
+                   IndexSet(first,1)->gen[0]);
 
   ccount offset = 2*t->ndim*t->ncomp;
   count dim, rul, n, maxdim = 0;
@@ -659,7 +696,7 @@ static void Sample(This *t, Region *region)
     for( dim = 0; dim < t->ndim; ++dim ) {
       creal *fp = f1 + t->ncomp;
       creal *fm = fp + t->ncomp;
-      creal fourthdiff = fabs(base +
+      creal fourthdiff = fabsx(base +
         ratio*(fp[0] + fm[0]) - (fp[offset] + fm[offset]));
       f1 = fm;
       if( fourthdiff > maxdiff ) {
@@ -688,7 +725,7 @@ static void Sample(This *t, Region *region)
       real maxerr = 0;
       for( s = first; s <= last; NextSet(s) )
         maxerr = Max(maxerr,
-          fabs(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
+          fabsx(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
       sum[rul] = maxerr;
     }
 
@@ -712,7 +749,7 @@ static void Sample(This *t, Region *region)
 
     for( res = result, comp = 0; res < Res; ++res )
       oe += sprintf(oe, "\n[" COUNT "] "
-        REAL " +- " REAL, ++comp, res->avg, res->err);
+        REAL " +- " REAL, ++comp, SHOW(res->avg), SHOW(res->err));
 
     Print(out);
   }
diff --git a/src/external/libCuba/src/cuhre/decl.h b/src/external/libCuba/src/cuhre/decl.h
index c846e9d2..bf924ba0 100644
--- a/src/external/libCuba/src/cuhre/decl.h
+++ b/src/external/libCuba/src/cuhre/decl.h
@@ -2,7 +2,7 @@
 	decl.h
 		Type declarations
 		this file is part of Cuhre
-		last modified 26 Jul 13 th
+		last modified 21 Jul 14 th
 */
 
 
@@ -47,16 +47,18 @@ typedef struct {
 
 typedef const Rule cRule;
 
-typedef int (*Integrand)(ccount *, creal *, ccount *, real *, void *);
+typedef int (*Integrand)(ccount *, creal *, ccount *, real *,
+  void *, cnumber *, cint *);
 
 typedef struct _this {
   count ndim, ncomp;
 #ifndef MLVERSION
   Integrand integrand;
   void *userdata;
+  number nvec;
 #ifdef HAVE_FORK
-  int ncores, *child;
   SHM_ONLY(int shmid;)
+  Spin *spin;
 #endif
 #endif
   real *frame;
diff --git a/src/external/libCuba/src/divonne/Divonne.c b/src/external/libCuba/src/divonne/Divonne.c
index a0577cc6..444dddde 100644
--- a/src/external/libCuba/src/divonne/Divonne.c
+++ b/src/external/libCuba/src/divonne/Divonne.c
@@ -4,7 +4,7 @@
 		originally by J.H. Friedman and M.H. Wright
 		(CERNLIB subroutine D151)
 		this version by Thomas Hahn
-		last modified 17 Sep 13 th
+		last modified 22 Jul 14 th
 */
 
 #define DIVONNE
@@ -16,7 +16,7 @@
 /*********************************************************************/
 
 Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber nvec,
   creal epsrel, creal epsabs,
   cint flags, cint seed,
   cnumber mineval, cnumber maxeval,
@@ -24,18 +24,22 @@ Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp,
   creal border, creal maxchisq, creal mindeviation,
   cnumber ngiven, ccount ldxgiven, real *xgiven,
   cnumber nextra, PeakFinder peakfinder,
-  cchar *statefile,
+  cchar *statefile, Spin **pspin,
   int *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = ndim;
   t.ncomp = ncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = nvec;
   t.epsrel = epsrel;
   t.epsabs = epsabs;
-  t.flags = flags;
+  t.flags = MaxVerbose(flags);
   t.seed = seed;
   t.mineval = mineval;
   t.maxeval = maxeval;
@@ -52,16 +56,19 @@ Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp,
   t.nextra = nextra;
   t.peakfinder = peakfinder;
   t.statefile = statefile;
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
 /*********************************************************************/
 
 Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber *pnvec,
   creal *pepsrel, creal *pepsabs,
   cint *pflags, cint *pseed,
   cnumber *pmineval, cnumber *pmaxeval,
@@ -69,18 +76,22 @@ Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp,
   creal *pborder, creal *pmaxchisq, creal *pmindeviation,
   cnumber *pngiven, ccount *pldxgiven, real *xgiven,
   cnumber *pnextra, PeakFinder peakfinder,
-  cchar *statefile,
+  cchar *statefile, Spin **pspin,
   int *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob, cint statefilelen)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = *pndim;
   t.ncomp = *pncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = *pnvec;
   t.epsrel = *pepsrel;
   t.epsabs = *pepsabs;
-  t.flags = *pflags;
+  t.flags = MaxVerbose(*pflags);
   t.seed = *pseed;
   t.mineval = *pmineval;
   t.maxeval = *pmaxeval;
@@ -97,9 +108,12 @@ Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp,
   t.nextra = *pnextra;
   t.peakfinder = peakfinder;
   CString(t.statefile, statefile, statefilelen);
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
diff --git a/src/external/libCuba/src/divonne/Explore.c b/src/external/libCuba/src/divonne/Explore.c
index 4f27ac25..e29520cb 100644
--- a/src/external/libCuba/src/divonne/Explore.c
+++ b/src/external/libCuba/src/divonne/Explore.c
@@ -2,7 +2,7 @@
 	Explore.c
 		sample region, determine min and max, split if necessary
 		this file is part of Divonne
-		last modified 2 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -19,10 +19,11 @@ static int ExploreSerial(This *t, ccount iregion)
   Region *region = RegionPtr(iregion);
   cBounds *bounds = region->bounds;
   Result *result = RegionResult(region);
+  real *minmax = RegionMinMax(region);
 
   Vector(Extrema, extrema, NCOMP);
   Vector(real, xtmp, NDIM);
-  Result *r, *r0;
+  Result *r;
   creal *x;
   real *f;
   real halfvol, maxerr;
@@ -102,7 +103,7 @@ skip:
       ftmp = FindMinimum(t, bounds, xtmp, e->fmin);
       if( ftmp < r->fmin ) {
         r->fmin = ftmp;
-        XCopy(&r->xminmax[0], xtmp);
+        XCopy(&minmax[2*comp*t->ndim], xtmp);
       }
 
       t->selectedcomp = Tag(comp);
@@ -110,12 +111,12 @@ skip:
       ftmp = -FindMinimum(t, bounds, xtmp, -e->fmax);
       if( ftmp > r->fmax ) {
         r->fmax = ftmp;
-        XCopy(&r->xminmax[t->ndim], xtmp);
+        XCopy(&minmax[(2*comp + 1)*t->ndim], xtmp);
       }
     }
 
     r->spread = halfvol*(r->fmax - r->fmin);
-    err = r->spread/Max(fabs(r->avg), NOTZERO);
+    err = r->spread/Max(fabsx(r->avg), NOTZERO);
     if( err > maxerr ) {
       maxerr = err;
       maxcomp = comp;
@@ -130,22 +131,22 @@ skip:
   }
 
   region->cutcomp = maxcomp;
-  r0 = RegionResult(region);
-  r = r0 + maxcomp;
+  r = RegionResult(region) + maxcomp;
   if( halfvol*(r->fmin + r->fmax) > r->avg ) {
     region->fminor = r->fmin;
     region->fmajor = r->fmax;
-    region->xmajor = &r->xminmax[t->ndim] - (real *)r0;
+    region->xmajor = (2*maxcomp + 1)*t->ndim;
   }
   else {
     region->fminor = r->fmax;
     region->fmajor = r->fmin;
-    region->xmajor = &r->xminmax[0] - (real *)r0;
+    region->xmajor = 2*maxcomp*t->ndim;
   }
 
   if( region->isamples == 0 ) {
     if( (region->depth < INIDEPTH && r->spread < samples->neff*r->err) ||
-        r->spread < t->totals[maxcomp].secondspread ) region->depth = 0;
+        r->spread < t->totals[maxcomp].secondspread )
+      region->depth = 0;
     if( region->depth == 0 )
       for( comp = 0; comp < t->ncomp; ++comp )
         t->totals[comp].secondspread =
diff --git a/src/external/libCuba/src/divonne/FindMinimum.c b/src/external/libCuba/src/divonne/FindMinimum.c
index c7e00b2d..ac7fd901 100644
--- a/src/external/libCuba/src/divonne/FindMinimum.c
+++ b/src/external/libCuba/src/divonne/FindMinimum.c
@@ -2,7 +2,7 @@
 	FindMinimum.c
 		find minimum (maximum) of hyperrectangular region
 		this file is part of Divonne
-		last modified 7 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -42,7 +42,7 @@ static inline real Dot(ccount n, creal *a, creal *b)
 
 static inline real Length(ccount n, creal *vec)
 {
-  return sqrt(Dot(n, vec, vec));
+  return sqrtx(Dot(n, vec, vec));
 }
 
 /*********************************************************************/
@@ -120,7 +120,7 @@ static void UpdateCholesky(cThis *t, ccount n, real *hessian,
     p[i] = dir;
     gamma += Sq(dir)/Hessian(i, i);
   }
-  gamma = Max(fabs(1 - gamma), EPS);
+  gamma = Max(fabsx(1 - gamma), EPS);
 
   while( --i >= 0 ) {
     creal dir = z[i] = p[i];
@@ -154,7 +154,7 @@ static inline void BFGS(cThis *t, ccount n, real *hessian,
 
   c = Dot(n, g, p);
   if( c >= 0 ) return;
-  c = 1/sqrt(-c);
+  c = 1/sqrtx(-c);
   for( i = 0; i < n; ++i )
     y[i] = c*g[i];
   UpdateCholesky(t, n, hessian, y, p);
@@ -199,7 +199,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
      c) the gradient is positive, i.e. we'd move uphill */
 
   if( step > 0 && range > tol2 && grad <= 0 ) {
-    creal eps = RTEPS*fabs(range) + ftol;
+    creal eps = RTEPS*fabsx(range) + ftol;
     creal mingrad = -1e-4*grad, maxgrad = -gtol*grad;
 
     real end = range + eps;
@@ -225,7 +225,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
         maxstep = maxstep*(1 + .75*RTEPS) + .75*tol;
       }
 
-      cur.dx = (fabs(step) >= tol) ? step : (step > 0) ? tol : -tol;
+      cur.dx = (fabsx(step) >= tol) ? step : (step > 0) ? tol : -tol;
       dist = distmin + cur.dx;
       for( i = 0; i < nfree; ++i ) {
         ccount dim = ifree[i];
@@ -248,7 +248,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
         if( cur.dx < 0 ) b = w;
         else a = w;
 
-        tol = RTEPS*fabs(distmin) + ftol;
+        tol = RTEPS*fabsx(distmin) + ftol;
         tol2 = tol + tol;
       }
       else {
@@ -260,14 +260,14 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
 
       if( distmin + b.dx <= xtol ) break;
       if( min.f < fini &&
-          a.f - min.f <= fabs(a.dx)*maxgrad &&
-          (fabs(distmin - range) > tol || maxstep < b.dx) ) break;
+          a.f - min.f <= fabsx(a.dx)*maxgrad &&
+          (fabsx(distmin - range) > tol || maxstep < b.dx) ) break;
 
       mid = .5*(a.dx + b.dx);
-      if( fabs(mid) <= tol2 - .5*(b.dx - a.dx) ) break;
+      if( fabsx(mid) <= tol2 - .5*(b.dx - a.dx) ) break;
 
       r = q = s = 0;
-      if( fabs(end) > tol ) {
+      if( fabsx(end) > tol ) {
         if( first ) {
           creal s1 = w.dx*grad;
           creal s2 = w.f - min.f;
@@ -281,7 +281,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
           q = 2*(s2 - s1);
         }
         if( q > 0 ) s = -s;
-        q = fabs(q);
+        q = fabsx(q);
         r = end;
         if( step != b1 || b.dx <= maxstep ) end = step;
       }
@@ -290,15 +290,15 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
       else if( b.dx > maxstep ) step = (step < b.dx) ? -4*a.dx : maxstep;
       else {
         real num = a.dx, den = b.dx;
-        if( fabs(b.dx) <= tol || (w.dx > 0 && fabs(a.dx) > tol) )
+        if( fabsx(b.dx) <= tol || (w.dx > 0 && fabsx(a.dx) > tol) )
           num = b.dx, den = a.dx;
         num /= -den;
-        step = (num < 1) ? .5*den*sqrt(num) : 5/11.*den*(.1 + 1/num);
+        step = (num < 1) ? .5*den*sqrtx(num) : 5/11.*den*(.1 + 1/num);
       }
 
       if( step > 0 ) a1 = a.dx, b1 = step;
       else a1 = step, b1 = b.dx;
-      if( fabs(s) < fabs(.5*q*r) && s > q*a1 && s < q*b1 ) {
+      if( fabsx(s) < fabsx(.5*q*r) && s > q*a1 && s < q*b1 ) {
         step = s/q;
         if( step - a.dx < tol2 || b.dx - step < tol2 )
           step = (mid > 0) ? tol : -tol;
@@ -307,7 +307,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
     }
 
     first = true;
-    if( fabs(distmin - range) < tol ) {
+    if( fabsx(distmin - range) < tol ) {
       distmin = range;
       if( maxstep > b.dx ) first = false;
     }
@@ -372,7 +372,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
      or we come close to a border. */
 
   XCopy(y, x);
-  ftest = SUFTOL*(1 + fabs(fx));
+  ftest = SUFTOL*(1 + fabsx(fx));
   delta = RTDELTA/5;
   do {
     delta = Min(5*delta, smax);
@@ -381,7 +381,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
       y[dim] = x[dim] + delta*p[i];
     }
     fy = Sample(t, y);
-    if( fabs(fy - fx) > ftest ) break;
+    if( fabsx(fy - fx) > ftest ) break;
   } while( delta != smax );
 
   /* Construct a second direction p' orthogonal to p, i.e. p.p' = 0.
@@ -425,7 +425,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
      or we come close to a border. */
 
   XCopy(z, y);
-  ftest = SUFTOL*(1 + fabs(fy));
+  ftest = SUFTOL*(1 + fabsx(fy));
   delta = RTDELTA/5;
   do {
     delta = Min(5*delta, smax);
@@ -434,7 +434,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
       z[dim] = y[dim] + delta*p[i];
     }
     fz = Sample(t, z);
-    if( fabs(fz - fy) > ftest ) break;
+    if( fabsx(fz - fy) > ftest ) break;
   } while( delta != smax );
 
   if( fy != fz ) {
@@ -541,12 +541,12 @@ static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin)
     bool resample = false;
     nfree = nfix = 0;
     for( dim = 0; dim < t->ndim; ++dim ) {
-      if( xmin[dim] < b[dim].lower + (1 + fabs(b[dim].lower))*QEPS ) {
+      if( xmin[dim] < b[dim].lower + (1 + fabsx(b[dim].lower))*QEPS ) {
         xmin[dim] = b[dim].lower;
         ifix[nfix++] = dim;
         resample = true;
       }
-      else if( xmin[dim] > b[dim].upper - (1 + fabs(b[dim].upper))*QEPS ) {
+      else if( xmin[dim] > b[dim].upper - (1 + fabsx(b[dim].upper))*QEPS ) {
         xmin[dim] = b[dim].upper;
         ifix[nfix++] = Tag(dim);
         resample = true;
@@ -562,7 +562,7 @@ static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin)
     if( local || Length(nfree, gfree) > GTOL ) break;
 
     ftmp = LocalSearch(t, nfree, ifree, b, xmin, fmin, tmp);
-    if( ftmp > fmin - (1 + fabs(fmin))*RTEPS )
+    if( ftmp > fmin - (1 + fabsx(fmin))*RTEPS )
       goto releasebounds;
     fmin = ftmp;
     XCopy(xmin, tmp);
@@ -586,7 +586,7 @@ static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin)
       minstep = INFTY;
       for( i = 0; i < nfree; ++i ) {
         count dim = Untag(ifree[i]);
-        if( fabs(p[i]) > EPS ) {
+        if( fabsx(p[i]) > EPS ) {
           real step;
           count fix;
           if( p[i] < 0 ) {
@@ -642,11 +642,11 @@ fixbound:
         BFGS(t, nfree, hessian, tmp, gfree, p, low.dx);
         XCopy(gfree, tmp);
 
-        if( fabs(low.dx - minstep) < QEPS*minstep ) goto fixbound;
+        if( fabsx(low.dx - minstep) < QEPS*minstep ) goto fixbound;
 
         fdiff = fini - fmin;
         fini = fmin;
-        if( fdiff > (1 + fabs(fmin))*FTOL ||
+        if( fdiff > (1 + fabsx(fmin))*FTOL ||
             low.dx*plen > (1 + Length(t->ndim, xmin))*FTOL ) continue;
       }
     }
diff --git a/src/external/libCuba/src/divonne/Integrate.c b/src/external/libCuba/src/divonne/Integrate.c
index c40e9502..335aa57c 100644
--- a/src/external/libCuba/src/divonne/Integrate.c
+++ b/src/external/libCuba/src/divonne/Integrate.c
@@ -5,7 +5,7 @@
 		then do a main integration over all regions
 		this file is part of Divonne
 		checkpointing by B. Chokoufe
-		last modified 5 Aug 13 th
+		last modified 13 Mar 15 th
 */
 
 
@@ -29,8 +29,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   Totals *tot, *Tot = state->totals + t->ncomp;
   Bounds *b, *B;
   Result *res;
-  count comp, err, iregion;
-  number nwant;
+  count comp, iregion;
+  number nwant, err;
   real nneed;
   ML_ONLY(number neff;)
   int fail;
@@ -38,6 +38,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   if( VERBOSE > 1 ) {
     sprintf(out, "Divonne input parameters:\n"
       "  ndim " COUNT "\n  ncomp " COUNT "\n"
+      ML_NOT("  nvec " NUMBER "\n")
       "  epsrel " REAL "\n  epsabs " REAL "\n"
       "  flags %d\n  seed %d\n"
       "  mineval " NUMBER "\n  maxeval " NUMBER "\n"
@@ -46,11 +47,12 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       "  ngiven " NUMBER "\n  nextra " NUMBER "\n"
       "  statefile \"%s\"",
       t->ndim, t->ncomp,
-      t->epsrel, t->epsabs,
+      ML_NOT(t->nvec,)
+      SHOW(t->epsrel), SHOW(t->epsabs),
       t->flags, t->seed,
       t->mineval, t->maxeval,
       t->key1, t->key2, t->key3, t->maxpass,
-      t->border.lower, t->maxchisq, t->mindeviation,
+      SHOW(t->border.lower), SHOW(t->maxchisq), SHOW(t->mindeviation),
       t->ngiven, t->nextra,
       t->statefile);
     Print(out);
@@ -132,7 +134,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   /* Step 1: partition the integration region */
 
   if( t->phase == 1 ) {
-    if( VERBOSE ) Print("Partitioning phase:");
+    if( VERBOSE ) Print("\nPartitioning phase:");
 
     if( ini ) Iterate(t, 0, INIDEPTH, 0, NULL);
 
@@ -164,7 +166,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
         integral[comp] = tot->avg;
         valid += tot->avg == tot->avg;
         if( tot->spreadsq > maxtot->spreadsq ) maxtot = tot;
-        tot->spread = sqrt(tot->spreadsq);
+        tot->spread = sqrtx(tot->spreadsq);
         error[comp] = tot->spread/t->samples[0].neff;
       }
 
@@ -197,7 +199,7 @@ if( StateWriteTest(t) ) { \
         for( comp = 0; comp < t->ncomp; ++comp )
           oe += sprintf(oe, "\n[" COUNT "] "
             REAL " +- " REAL,
-            comp + 1, integral[comp], error[comp]);
+            comp + 1, SHOW(integral[comp]), SHOW(error[comp]));
         Print(out);
       }
 
@@ -230,7 +232,7 @@ if( StateWriteTest(t) ) { \
     tot->maxerrsq = Sq(maxerr);
     tot->mindevsq = tot->maxerrsq*Sq(t->mindeviation);
   }
-  nwant = (number)Min(ceil(nneed), MARKMASK/40.);
+  nwant = (number)Min(ceil(nneed), NWANTMAX/40.);
 
   err = SamplesLookup(t, &t->samples[1], t->key2, nwant,
     (t->maxeval - t->neval)/t->nregions + 1, t->samples[0].n + 1);
@@ -308,11 +310,11 @@ refine:
                 Iterate(t, state->iregion, POSTDEPTH, 1, state->totals);
 
                 if( can_adjust ) {
-                  cnumber nnew = (tot->spreadsq/Sq(MARKMASK) > tot->maxerrsq) ?
-                    MARKMASK :
-                    (number)ceil(sqrt(tot->spreadsq/tot->maxerrsq));
+                  cnumber nnew = (tot->spreadsq/Sq(NWANTMAX) > tot->maxerrsq) ?
+                    NWANTMAX :
+                    (number)ceil(sqrtx(tot->spreadsq/tot->maxerrsq));
                   if( nnew > nwant + nwant/64 ) {
-                    ccount err = SamplesLookup(t, &t->samples[1], t->key2, nnew,
+                    cnumber err = SamplesLookup(t, &t->samples[1], t->key2, nnew,
                       (t->maxeval - t->neval)/t->nregions + 1, t->samples[1].n);
                     fail += Unmark(err)*t->nregions;
                     nwant = nnew;
@@ -404,14 +406,14 @@ refine:
         if( chisq > EPS ) chisq /= Max(chiden, NOTZERO);
 
         if( VERBOSE > 2 ) {
-#define Out2(f, r) (r)->avg, res->spread/t->samples[f].neff, (r)->err
+#define Out2(f, r) SHOW((r)->avg), SHOW(res->spread/t->samples[f].neff), SHOW((r)->err)
 #define Out(f) Out2(f, &tot->phase[f])
           oe += sprintf(oe, "\n[" COUNT "] "
             REAL " +- " REAL "(" REAL ")\n    "
             REAL " +- " REAL "(" REAL ")", ++comp, Out(0), Out(1));
           if( todo == 3 ) oe += sprintf(oe, "\n    "
             REAL " +- " REAL "(" REAL ")", Out2(2, res));
-          oe += sprintf(oe, "  \tchisq " REAL, chisq);
+          oe += sprintf(oe, "  \tchisq " REAL, SHOW(chisq));
         }
 
         tot->integral += avg;
@@ -419,7 +421,7 @@ refine:
         tot->chisq += chisq;
 
         res->avg = avg;
-        res->spread = sqrt(sigsq);
+        res->spread = sqrtx(sigsq);
         res->chisq = chisq;
       }
 
@@ -433,7 +435,7 @@ refine:
 
     for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) {
       integral[comp] = tot->integral;
-      error[comp] = sqrt(tot->sigsq);
+      error[comp] = sqrtx(tot->sigsq);
       prob[comp] = ChiSquare(tot->chisq, df);
     }
 
@@ -442,7 +444,8 @@ refine:
       for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp )
         oe += sprintf(oe, "\n[" COUNT "] "
           REAL " +- " REAL "  \tchisq " REAL " (" COUNT " df)",
-          comp + 1, integral[comp], error[comp], tot->chisq, df);
+          comp + 1, SHOW(integral[comp]), SHOW(error[comp]),
+          SHOW(tot->chisq), df);
       Print(out);
     }
 
@@ -468,12 +471,12 @@ refine:
 
       MLPutFunction(stdlink, "Cuba`Divonne`region", 4);
 
-      MLPutRealList(stdlink, bounds, 2*t->ndim);
+      MLPutRealxList(stdlink, bounds, 2*t->ndim);
 
       MLPutFunction(stdlink, "List", t->ncomp);
       for( Res = (res = RegionResult(region)) + t->ncomp; res < Res; ++res ) {
         real r[] = {res->avg, res->spread/neff, res->chisq};
-        MLPutRealList(stdlink, r, Elements(r));
+        MLPutRealxList(stdlink, r, Elements(r));
       }
 
       MLPutInteger(stdlink, region->depth + 1);  /* misused for df */
@@ -482,8 +485,7 @@ refine:
 #endif
 
 abort:
-  WaitCores(t);
-  FORK_ONLY(FrameFree(t, ShmRm(t));)
+  FORK_ONLY(FrameFree(t, Master);)
 
   RuleFree(t);
   SamplesFree(&t->samples[2]);
diff --git a/src/external/libCuba/src/divonne/Iterate.c b/src/external/libCuba/src/divonne/Iterate.c
index 3d9488de..ac4a4cda 100644
--- a/src/external/libCuba/src/divonne/Iterate.c
+++ b/src/external/libCuba/src/divonne/Iterate.c
@@ -2,7 +2,7 @@
 	Iterate.c
 		recursion over regions
 		this file is part of Divonne
-		last modified 2 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -85,10 +85,10 @@ FORK_ONLY(more:)
 
       norm = 1./nsplit--;
       for( res = RegionResult(parent), c = corr; c < C; ++res, ++c ) {
-        creal diff = fabs(res->avg - c->avg)*norm;
+        creal diff = fabsx(res->avg - c->avg)*norm;
         c->avg = diff*norm*nsplit;
         c->err = (c->err == 0) ? 1 : 1 + diff/c->err;
-        c->spread = (c->spread == 0) ? 1 : 1 + diff/sqrt(c->spread);
+        c->spread = (c->spread == 0) ? 1 : 1 + diff/sqrtx(c->spread);
       }
 
       for( td = todo; td < tdmax; ++td )
diff --git a/src/external/libCuba/src/divonne/Rule.c b/src/external/libCuba/src/divonne/Rule.c
index 5bb6c28b..813b9184 100644
--- a/src/external/libCuba/src/divonne/Rule.c
+++ b/src/external/libCuba/src/divonne/Rule.c
@@ -4,14 +4,10 @@
 		code lifted with minor modifications from DCUHRE
 		by J. Berntsen, T. Espelid, and A. Genz
 		this file is part of Divonne
-		last modified 26 Jul 13 th
+		last modified 7 May 15 th
 */
 
 
-#define NextSet(p) p = (Set *)((char *)p + setsize)
-
-/*********************************************************************/
-
 static void Rule13Alloc(This *t)
 {
   static creal w[][nrules] = {
@@ -153,7 +149,7 @@ static void Rule13Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -298,7 +294,7 @@ static void Rule11Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -309,32 +305,55 @@ static void Rule11Alloc(This *t)
 static void Rule9Alloc(This *t)
 {
   static creal w[] = {
-    -.0023611709677855117884,   .11415390023857325268,
-    -.63833920076702389094,     .74849988504685208004,
-    -.0014324017033399125142,   .057471507864489725949,
-    -.14225104571434243234,    -.062875028738286979989,
-     .254591133248959089,     -1.207328566678236261,
-     .89567365764160676508,    -.36479356986049146661,
-     .0035417564516782676826,  -.072609367395893679605,
-     .10557491625218991012,     .0021486025550098687713,
-    -.032268563892953949998,    .010636783990231217481,
-     .014689102496143490175,    .51134708346467591431,
-     .45976448120806344646,     .18239678493024573331,
-    -.04508628929435784076,     .21415883524352793401,
-    -.027351546526545644722,    .054941067048711234101,
-     .11937596202570775297,     .65089519391920250593,
-     .14744939829434460168,     .057693384490973483573,
-     .034999626602143583822,  -1.3868627719278281436,
-    -.2386668732575008879,      .015532417276607053264,
-     .0035328099607090870236,   .09231719987444221619,
-     .02254314464717892038,     .013675773263272822361,
-    -.32544759695960125297,     .0017708782258391338413,
-     .0010743012775049343856,   .25150011495314791996 };
+    RC(-.002361170967785511788400941242259231309691),
+    RC(.1141539002385732526821323741697655347686),
+    RC(-.6383392007670238909386026193674701393074),
+    RC(.7484998850468520800423030047583803945205),
+    RC(-.001432401703339912514196154599769007103671),
+    RC(.05747150786448972594860897296200006759892),
+    RC(-.1422510457143424323449521620935950679394),
+    RC(-.06287502873828697998942424881040490136987),
+    RC(.2545911332489590890011611142429070613156),
+    RC(-1.207328566678236261002219995185143356737),
+    RC(.8956736576416067650809467826488567200939),
+    RC(-.3647935698604914666100134551377381205297),
+    RC(.003541756451678267682601411863388846964536),
+    RC(-.07260936739589367960492815865074633743652),
+    RC(.1055749162521899101218622863269817454540),
+    RC(.002148602555009868771294231899653510655506),
+    RC(-.03226856389295394999786630399875134318006),
+    RC(.01063678399023121748083624225818915724455),
+    RC(.01468910249614349017540783437728097691502),
+    RC(.5113470834646759143109387357149329909126),
+    RC(.4597644812080634464633352781605214342691),
+    RC(.1823967849302457333050067275688690602649),
+    RC(-.04508628929435784075980562738240804429658),
+    RC(.2141588352435279340097929526588394300172),
+    RC(-.02735154652654564472203690086290223507436),
+    RC(.05494106704871123410060080562462135546101),
+    RC(.1193759620257077529708962121565290178730),
+    RC(.6508951939192025059314756320878023215278),
+    RC(.1474493982943446016775696826942585013243),
+    RC(.05769338449097348357291272840392627722165),
+    RC(.03499962660214358382244159694487155861542),
+    RC(-1.386862771927828143599782668709014266770),
+    RC(-.2386668732575008878964134721962088068396),
+    RC(.01553241727660705326386197156586357005224),
+    RC(.003532809960709087023561817517751309380604),
+    RC(.09231719987444221619017126187763868745587),
+    RC(.02254314464717892037990281369120402214829),
+    RC(.01367577326327282236101845043145111753718),
+    RC(-.3254475969596012529657378160439011607639),
+    RC(.001770878225839133841300705931694423482268),
+    RC(.001074301277504934385647115949826755327753),
+    RC(.2515001149531479199576969952416196054795) };
 
   static creal g[] = {
-     .47795365790226950619,     .20302858736911986780,
-     .44762735462617812882,     .125,
-     .34303789878087814570 };
+    RC(.4779536579022695061928604197171830064732),
+    RC(.2030285873691198677998034402373279133258),
+    RC(.4476273546261781288207704806530998539285),
+    RC(.125),
+    RC(.3430378987808781457001426145164678603407) };
 
   enum { nsets = 9 };
 
@@ -440,7 +459,7 @@ static void Rule9Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -451,21 +470,33 @@ static void Rule9Alloc(This *t)
 static void Rule7Alloc(This *t)
 {
   static creal w[] = {
-     .019417866674748388428,   -.40385257701150182546,
-     .64485668767465982223,     .01177982690775806141,
-    -.18041318740733609012,    -.088785828081335044443,
-     .056328645808285941374,   -.0097089333373741942142,
-    -.99129176779582358138,    -.17757165616267008889,
-     .12359398032043233572,     .074978148702033690681,
-     .55489147051423559776,     .088041241522692771226,
-     .021118358455513385083,   -.0099302203239653333087,
-    -.064100053285010904179,    .030381729038221007659,
-     .0058899134538790307051,  -.0048544666686870971071,
-     .35514331232534017777 };
+    RC(.01941786667474838842844534313920462333850),
+    RC(-.4038525770115018254611834753723880293161),
+    RC(.6448566876746598222277360730193089551024),
+    RC(.01177982690775806141012214458820955067854),
+    RC(-.1804131874073360901182293138710989490609),
+    RC(-.08878582808133504444306598174517276122439),
+    RC(.05632864580828594137378124255408286479947),
+    RC(-.009708933337374194214222671569602311669249),
+    RC(-.9912917677958235813775106862002319060386),
+    RC(-.1775716561626700888861319634903455224488),
+    RC(.1235939803204323357183625846672135876752),
+    RC(.07497814870203369068087999555157339703666),
+    RC(.5548914705142355977605994477355651401434),
+    RC(.08804124152269277122645182458858273865209),
+    RC(.02111835845551338508329573367808085283304),
+    RC(-.009930220323965333308685820460105538586058),
+    RC(-.06410005328501090417895544042025034295870),
+    RC(.03038172903822100765927778829870429682489),
+    RC(.005889913453879030705061072294104775339268),
+    RC(-.004854466668687097107111335784801155834624),
+    RC(.3551433123253401777722639269806910448976) };
 
   static creal g[] = {
-     .47795365790226950619,     .20302858736911986780,
-     .375,                      .34303789878087814570 };
+    RC(.4779536579022695061928604197171830064732),
+    RC(.2030285873691198677998034402373279133258),
+    RC(.375),
+    RC(.3430378987808781457001426145164678603407) };
 
   enum { nsets = 6 };
 
@@ -541,7 +572,7 @@ static void Rule7Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -666,7 +697,7 @@ static void SampleRule(This *t, ccount iregion)
       real maxerr = 0;
       for( s = first; s <= last; NextSet(s) )
         maxerr = Max(maxerr,
-          fabs(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
+          fabsx(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
       sum[rul] = maxerr;
     }
 
diff --git a/src/external/libCuba/src/divonne/Sample.c b/src/external/libCuba/src/divonne/Sample.c
index 7778f581..113b6aed 100644
--- a/src/external/libCuba/src/divonne/Sample.c
+++ b/src/external/libCuba/src/divonne/Sample.c
@@ -2,14 +2,16 @@
 	Sample.c
 		most of what is related to sampling
 		this file is part of Divonne
-		last modified 30 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
-#define MARKMASK 0xfffffff
+#define MARKMASK NUMBER_MAX
 #define Marked(x) ((x) & ~MARKMASK)
 #define Unmark(x) ((x) & MARKMASK)
 
+#define NWANTMAX NUMBER_MAX
+
 #define EXTRAPOLATE_EPS (.25*t->border.lower)
 /*#define EXTRAPOLATE_EPS 0x1p-26*/
 
@@ -96,7 +98,7 @@ static void SampleKorobov(This *t, ccount iregion)
   }
 
   if( dist > 0 ) {
-    dist = sqrt(dist)/EXTRAPOLATE_EPS;
+    dist = sqrtx(dist)/EXTRAPOLATE_EPS;
     for( dim = 0; dim < t->ndim; ++dim ) {
       real x2 = x[dim], dx = x2 - t->border.upper;
       if( dx > 0 ) {
@@ -148,7 +150,7 @@ static void SampleKorobov(This *t, ccount iregion)
          1..39 = multiplicator,        Korobov numbers,
        40..inf = absolute # of points, Korobov numbers.  */
 
-static count SamplesLookup(This *t, Samples *samples, cint key,
+static number SamplesLookup(This *t, Samples *samples, cint key,
   cnumber nwant, cnumber nmax, number nmin)
 {
   number n;
@@ -191,7 +193,8 @@ static count SamplesLookup(This *t, Samples *samples, cint key,
 static void SamplesAlloc(cThis *t, Samples *samples)
 {
 #define FIRST -INT_MAX
-#define MarkLast(x) (x | Marked(INT_MAX))
+#define MarkLast(x) ((x) | 0x40000000)
+#define UnmarkLast(x) ((x) & 0x3fffffff)
 
 #include "KorobovCoeff.c"
 
@@ -205,12 +208,12 @@ static void SamplesAlloc(cThis *t, Samples *samples)
 
     while( i = IMin(IDim(i), max),
            n > (p = prime[i + 1]) || n <= prime[i] ) {
-      cint d = (n - Unmark(p)) >> ++shift;
+      cint d = (n - UnmarkLast(p)) >> ++shift;
       i += Min1(d);
     }
 
     samples->coeff = coeff[i][t->ndim - KOROBOV_MINDIM];
-    samples->neff = p = Unmark(p);
+    samples->neff = p = UnmarkLast(p);
     samples->n = p/2 + 1;
   }
 
@@ -240,7 +243,7 @@ static real Sample(This *t, creal *x0)
   }
 
   if( dist > 0 ) {
-    dist = sqrt(dist)/EXTRAPOLATE_EPS;
+    dist = sqrtx(dist)/EXTRAPOLATE_EPS;
     for( dim = 0; dim < t->ndim; ++dim ) {
       real x2 = xtmp[dim], dx, b;
       if( (dx = x2 - (b = t->border.lower)) < 0 ||
diff --git a/src/external/libCuba/src/divonne/Split.c b/src/external/libCuba/src/divonne/Split.c
index 570258f8..abffcbf8 100644
--- a/src/external/libCuba/src/divonne/Split.c
+++ b/src/external/libCuba/src/divonne/Split.c
@@ -2,7 +2,7 @@
 	Split.c
 		determine optimal cuts for splitting a region
 		this file is part of Divonne
-		last modified 31 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -31,7 +31,7 @@ typedef struct {
 
 static inline real Div(creal a, creal b)
 {
-  return (b != 0 /*&& fabs(a) > SMALL*fabs(b)*/) ? a/b : a;
+  return (b != 0 /*&& fabsx(a) > SMALL*fabsx(b)*/) ? a/b : a;
 }
 
 /*********************************************************************/
@@ -58,7 +58,7 @@ static void SomeCut(This *t, Cut *cut, Bounds *b)
     yupper = Sample(t, xmid);
     xmid[dim] = x;
 
-    dev = fabs(ymid - .5*(ylower + yupper));
+    dev = fabsx(ymid - .5*(ylower + yupper));
     if( dev >= maxdev ) {
       maxdev = dev;
       maxdim = dim;
@@ -178,7 +178,7 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol,
 
     for( icut = 0; icut < ncuts; ++icut ) {
       Cut *c = &cut[icut];
-      creal diff = fabs(fmajor - c->f);
+      creal diff = fabsx(fmajor - c->f);
       if( diff <= mindiff ) {
         mindiff = diff;
         mincut = c;
@@ -230,18 +230,18 @@ repeat:
     if( lhssqnew <= lhssq ) {
       real fmax;
 
-      if( fabs(gammanew - gamma) < GAMMATOL*gamma ) break;
+      if( fabsx(gammanew - gamma) < GAMMATOL*gamma ) break;
       gamma = gammanew;
 
-      fmax = fabs(fgamma);
+      fmax = fabsx(fgamma);
       for( icut = 0; icut < ncuts; ++icut ) {
         Cut *c = &cut[icut];
         creal dfmin = SINGTOL*c->df;
         creal sol = c->sol/div;
         real df = c->f - c->fold;
-        df = (fabs(df) > SMALL*fabs(sol)) ? df/sol : 1;
-        c->df = (fabs(df) < fabs(dfmin)) ? dfmin : df;
-        fmax = Max(fmax, fabs(c->f));
+        df = (fabsx(df) > SMALL*fabsx(sol)) ? df/sol : 1;
+        c->df = (fabsx(df) < fabsx(dfmin)) ? dfmin : df;
+        fmax = Max(fmax, fabsx(c->f));
         c->fold = c->f;
       }
 
@@ -276,7 +276,7 @@ static void Split(This *t, ccount iregion)
   t->selectedcomp = region->cutcomp;
   t->neval_cut -= t->neval;
   ncuts = FindCuts(t, cut, region->bounds, region->vol,
-    (real *)RegionResult(region) + region->xmajor, region->fmajor,
+    RegionMinMax(region) + region->xmajor, region->fmajor,
     region->fmajor - region->fminor);
   t->neval_cut += t->neval;
 
diff --git a/src/external/libCuba/src/divonne/decl.h b/src/external/libCuba/src/divonne/decl.h
index cf127851..98023ec5 100644
--- a/src/external/libCuba/src/divonne/decl.h
+++ b/src/external/libCuba/src/divonne/decl.h
@@ -2,7 +2,7 @@
 	decl.h
 		Type declarations
 		this file is part of Divonne
-		last modified 26 Jul 13 th
+		last modified 9 Oct 14 th
 */
 
 
@@ -44,6 +44,8 @@ typedef struct {
 
 #define SetSize (sizeof(Set) + t->ndim*sizeof(real))
 
+#define NextSet(p) p = (Set *)((char *)p + setsize)
+
 typedef struct {
   Set *first, *last;
   real errcoeff[3];
@@ -71,12 +73,11 @@ typedef const Errors cErrors;
 typedef struct {
   real avg, err, spread, chisq;
   real fmin, fmax;
-  real xminmax[];
 } Result;
 
 typedef const Result cResult;
 
-#define ResultSize (sizeof(Result) + t->ndim*2*sizeof(real))
+#define MinMaxSize (t->ncomp*t->ndim*2*sizeof(real))
 
 typedef struct region {
   int depth, next;
@@ -85,29 +86,34 @@ typedef struct region {
   Bounds bounds[];
 } Region;
 
-#define RegionSize (sizeof(Region) + t->ndim*sizeof(Bounds) + t->ncomp*ResultSize)
+#define RegionSize (sizeof(Region) + t->ndim*sizeof(Bounds) + t->ncomp*sizeof(Result) + MinMaxSize)
 
 #define RegionResult(r) ((Result *)(r->bounds + t->ndim))
 
+#define RegionMinMax(r) ((real *)(RegionResult(r) + t->ncomp))
+
 #define RegionPtr(n) ((Region *)((char *)t->region + (n)*regionsize))
 
 
-typedef int (*Integrand)(ccount *, creal *, ccount *, real *, void *, cint *);
+typedef int (*Integrand)(ccount *, creal *, ccount *, real *,
+  void *, cnumber *, cint *, cint *);
 
-typedef void (*PeakFinder)(ccount *, cBounds *, number *, real *);
+typedef void (*PeakFinder)(ccount *, cBounds *, number *, real *, void *);
 
 typedef struct _this {
   count ndim, ncomp;
 #ifndef MLVERSION
   Integrand integrand;
   void *userdata;
-  PeakFinder peakfinder;
+  number nvec;
 #ifdef HAVE_FORK
-  int ncores, running, *child;
+  SHM_ONLY(int shmid;)
+  Spin *spin;
   real *frame;
   number nframe;
-  SHM_ONLY(int shmid;)
+  int running;
 #endif
+  PeakFinder peakfinder;
 #endif
   real epsrel, epsabs;
   int flags, seed;
diff --git a/src/external/libCuba/src/suave/Fluct.c b/src/external/libCuba/src/suave/Fluct.c
index ad1180ab..c226d0cc 100644
--- a/src/external/libCuba/src/suave/Fluct.c
+++ b/src/external/libCuba/src/suave/Fluct.c
@@ -2,35 +2,43 @@
 	Fluct.c
 		compute the fluctuation in the left and right half
 		this file is part of Suave
-		last modified 29 Jul 13 th
+		last modified 14 Mar 15 th
 */
 
 
-#if defined(HAVE_LONG_DOUBLE) && defined(HAVE_POWL)
+#if defined(HAVE_LONG_DOUBLE) && defined(HAVE_POWL) && REALSIZE <= 10
 
-typedef long double realx;
-#define XDBL_MAX_EXP LDBL_MAX_EXP
-#define XDBL_MAX LDBL_MAX
-#define powx powl
-#define ldexpx ldexpl
+typedef long double realL;
+#define REALL_MAX_EXP LDBL_MAX_EXP
+#define REALL_MAX LDBL_MAX
+#define powL powl
+#define ldexpL ldexpl
 
 #else
 
-typedef double realx;
-#define XDBL_MAX_EXP DBL_MAX_EXP
-#define XDBL_MAX DBL_MAX
-#define powx pow
-#define ldexpx ldexp
+typedef real realL;
+#define REALL_MAX_EXP REAL_MAX_EXP
+#define REALL_MAX REAL_MAX
+#define powL powx
+#define ldexpL ldexpx
 
 #endif
 
-typedef const realx crealx;
+typedef const realL crealL;
 
 typedef struct {
-  realx fluct;
+  realL fluct;
   number n;
 } Var;
 
+static inline realL MinL(crealL a, crealL b) {
+  return (a < b) ? a : b;
+}
+
+static inline realL MaxL(crealL a, crealL b) {
+  return (a > b) ? a : b;
+}
+
 /*********************************************************************/
 
 static void Fluct(cThis *t, Var *var,
@@ -41,27 +49,27 @@ static void Fluct(cThis *t, Var *var,
   count nvar = 2*t->ndim;
   creal norm = 1/(err*Max(fabs(avg), err));
   creal flat = 2/3./t->flatness;
-  crealx max = ldexpx(1., (int)((XDBL_MAX_EXP - 2)/t->flatness));
+  crealL max = ldexpL(1., (int)((REALL_MAX_EXP - 2)/t->flatness));
 
   Clear(var, nvar);
 
   while( n-- ) {
     count dim;
-    crealx arg = 1 + fabs(*w++)*Sq(*f - avg)*norm;
-    crealx ft = powx(arg < max ? arg : max, t->flatness);
+    crealL arg = 1 + fabs(*w++)*Sq(*f - avg)*norm;
+    crealL ft = powL(MinL(arg, max), t->flatness);
 
     f += t->ncomp;
 
     for( dim = 0; dim < t->ndim; ++dim ) {
       Var *v = &var[2*dim + (*x++ >= .5*(b[dim].lower + b[dim].upper))];
-      crealx f = v->fluct + ft;
-      v->fluct = (f > XDBL_MAX/2) ? XDBL_MAX/2 : f;
+      crealL f = v->fluct + ft;
+      v->fluct = MaxL(f, REALL_MAX/2);
       ++v->n;
     }
   }
 
   while( nvar-- ) {
-    var->fluct = powx(var->fluct, flat);
+    var->fluct = powL(var->fluct, flat);
     ++var;
   }
 }
diff --git a/src/external/libCuba/src/suave/Grid.c b/src/external/libCuba/src/suave/Grid.c
index 8e52fde8..5ba2c3b5 100644
--- a/src/external/libCuba/src/suave/Grid.c
+++ b/src/external/libCuba/src/suave/Grid.c
@@ -2,7 +2,7 @@
 	Grid.c
 		utility functions for the Vegas grid
 		this file is part of Suave
-		last modified 7 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -33,7 +33,7 @@ static void RefineGrid(cThis *t, Grid grid, Grid margsum)
     real impfun = 0;
     if( margsum[bin] > 0 ) {
       creal r = margsum[bin]*norm;
-      avgperbin += impfun = pow((r - 1)/log(r), 1.5);
+      avgperbin += impfun = powx((r - 1)/log(r), 1.5);
     }
     imp[bin] = impfun;
   }
diff --git a/src/external/libCuba/src/suave/Integrate.c b/src/external/libCuba/src/suave/Integrate.c
index 9c06d81f..4345f764 100644
--- a/src/external/libCuba/src/suave/Integrate.c
+++ b/src/external/libCuba/src/suave/Integrate.c
@@ -3,7 +3,7 @@
 		integrate over the unit hypercube
 		this file is part of Suave
 		checkpointing by B. Chokoufe
-		last modified 5 Aug 13 th
+		last modified 13 Mar 15 th
 */
 
 
@@ -26,22 +26,27 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   Result *tot, *Tot = state->totals + t->ncomp;
   Result *res, *resL, *resR;
   Bounds *b, *B;
+  cnumber minsamples = IMax(t->nmin, MINSAMPLES);
   count dim, comp;
   int fail;
 
   if( VERBOSE > 1 ) {
     sprintf(out, "Suave input parameters:\n"
       "  ndim " COUNT "\n  ncomp " COUNT "\n"
+      ML_NOT("  nvec " NUMBER "\n")
       "  epsrel " REAL "\n  epsabs " REAL "\n"
       "  flags %d\n  seed %d\n"
       "  mineval " NUMBER "\n  maxeval " NUMBER "\n"
-      "  nnew " NUMBER "\n  flatness " REAL "\n"
-      "  statefile \"%s\"\n",
+      "  nnew " NUMBER "\n  nmin " NUMBER "\n"
+      "  flatness " REAL "\n"
+      "  statefile \"%s\"",
       t->ndim, t->ncomp,
-      t->epsrel, t->epsabs,
+      ML_NOT(t->nvec,)
+      SHOW(t->epsrel), SHOW(t->epsabs),
       t->flags, t->seed,
       t->mineval, t->maxeval,
-      t->nnew, t->flatness,
+      t->nnew, t->nmin,
+      SHOW(t->flatness),
       t->statefile);
     Print(out);
   }
@@ -49,7 +54,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   if( BadComponent(t) ) return -2;
   if( BadDimension(t) ) return -1;
 
-  ShmAlloc(t, ShmRm(t));
+  ShmAlloc(t, Master);
   ForkCores(t);
 
   if( (fail = setjmp(t->abort)) ) goto abort;
@@ -132,7 +137,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       for( tot = state->totals, comp = 0; tot < Tot; ++tot )
         oe += sprintf(oe, "\n[" COUNT "] "
           REAL " +- " REAL "  \tchisq " REAL " (" COUNT " df)",
-          ++comp, tot->avg, tot->err, tot->chisq, state->df);
+          ++comp, SHOW(tot->avg), SHOW(tot->err),
+          SHOW(tot->chisq), state->df);
       Print(out);
     }
 
@@ -173,7 +179,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       region->result[maxcomp].avg, Max(maxerr, t->epsabs));
 
     bias = (t->epsrel < 1e-50) ? 2 :
-      Max(pow(2., -(real)region->div/t->ndim)/t->epsrel, 2.);
+      Max(powx(2., -(real)region->div/t->ndim)/t->epsrel, 2.);
     minfluct = INFTY;
     bisectdim = 0;
     for( dim = 0; dim < t->ndim; ++dim ) {
@@ -190,9 +196,9 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
     minfluct = vLR[0].fluct + vLR[1].fluct;
     nnewL = IMax(
       (minfluct == 0) ? t->nnew/2 : (count)(vLR[0].fluct/minfluct*t->nnew),
-      MINSAMPLES );
+      minsamples );
     nL = vLR[0].n + nnewL;
-    nnewR = IMax(t->nnew - nnewL, MINSAMPLES);
+    nnewR = IMax(t->nnew - nnewL, minsamples);
     nR = vLR[1].n + nnewR;
 
     regionL = RegionAlloc(t, nL, nnewL);
@@ -211,7 +217,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
     while( n-- ) {
       cbool final = (*w < 0);
       if( x[bisectdim] < mid ) {
-        if( final && wR > RegionW(regionR) ) wR[-1] = -fabs(wR[-1]);
+        if( final && wR > RegionW(regionR) ) wR[-1] = -fabsx(wR[-1]);
         *wL++ = *w++;
         XCopy(xL, x);
         xL += t->ndim;
@@ -219,7 +225,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
         fL += t->ncomp;
       }
       else {
-        if( final && wL > RegionW(regionL) ) wL[-1] = -fabs(wL[-1]);
+        if( final && wL > RegionW(regionL) ) wL[-1] = -fabsx(wL[-1]);
         *wR++ = *w++;
         XCopy(xR, x);
         xR += t->ndim;
@@ -259,15 +265,15 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       diff = Sq(.25*diff);
       sigsq = resL->sigsq + resR->sigsq;
       if( sigsq > 0 ) {
-        creal c = Sq(1 + sqrt(diff/sigsq));
+        creal c = Sq(1 + sqrtx(diff/sigsq));
         resL->sigsq *= c;
         resR->sigsq *= c;
       }
-      resL->err = sqrt(resL->sigsq += diff);
-      resR->err = sqrt(resR->sigsq += diff);
+      resL->err = sqrtx(resL->sigsq += diff);
+      resR->err = sqrtx(resR->sigsq += diff);
 
       tot->sigsq += resL->sigsq + resR->sigsq - res->sigsq;
-      tot->err = sqrt(tot->sigsq);
+      tot->err = sqrtx(tot->sigsq);
 
       tot->chisq += resL->chisq + resR->chisq - res->chisq;
     }
@@ -319,12 +325,12 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
 
       MLPutFunction(stdlink, "Cuba`Suave`region", 3);
 
-      MLPutRealList(stdlink, bounds, 2*t->ndim);
+      MLPutRealxList(stdlink, bounds, 2*t->ndim);
 
       MLPutFunction(stdlink, "List", t->ncomp);
       for( Res = (res = region->result) + t->ncomp; res < Res; ++res ) {
         real r[] = {res->avg, res->err, res->chisq};
-        MLPutRealList(stdlink, r, Elements(r));
+        MLPutRealxList(stdlink, r, Elements(r));
       }
 
       MLPutInteger(stdlink, region->df);
@@ -338,8 +344,7 @@ abort:
     anchor = anchor->next;
     free(region);
   }
-  WaitCores(t);
-  ShmFree(t);
+  ShmFree(t, Master);
 
   StateRemove(t);
 
diff --git a/src/external/libCuba/src/suave/Sample.c b/src/external/libCuba/src/suave/Sample.c
index 624c816c..a02249e2 100644
--- a/src/external/libCuba/src/suave/Sample.c
+++ b/src/external/libCuba/src/suave/Sample.c
@@ -2,7 +2,7 @@
 	Sample.c
 		the sampling step of Suave
 		this file is part of Suave
-		last modified 30 Jul 13 th
+		last modified 13 Mar 15 th
 */
 
 
@@ -71,7 +71,7 @@ static void Sample(This *t, cnumber nnew, Region *region,
 
   while( w < lastw ) {
     cbool final = (*w < 0);
-    creal weight = fabs(*w++);
+    creal weight = fabsx(*w++);
     ++n;
 
     for( c = cumul, comp = 0; c < C; ++c ) {
@@ -80,20 +80,20 @@ static void Sample(This *t, cnumber nnew, Region *region,
       c->sqsum += Sq(wfun);
 
       if( final ) {
-        if( n > 1 ) {
+        if( n >= t->nmin ) {
           real w = Weight(c->sum, c->sqsum, n);
           c->weightsum += c->weight = w;
           c->avgsum += c->avg = w*c->sum;
 
           if( VERBOSE > 2 ) {
-            creal sig = sqrt(1/w);
+            creal sig = sqrtx(1/w);
             ss[comp] += (df == 0) ?
               sprintf(ss[comp], "\n[" COUNT "] "
                 REAL " +- " REAL " (" NUMBER ")", comp + 1,
-                c->sum, sig, n) :
+                SHOW(c->sum), SHOW(sig), n) :
               sprintf(ss[comp], "\n    "
                 REAL " +- " REAL " (" NUMBER ")",
-                c->sum, sig, n);
+                SHOW(c->sum), SHOW(sig), n);
           }
 
           if( df == 0 ) c->guess = c->sum;
@@ -107,7 +107,7 @@ static void Sample(This *t, cnumber nnew, Region *region,
       }
     }
 
-    if( final ) ++df, n = 0;
+    if( final ) df -= NegQ(t->nmin - n - 1), n = 0;
   }
 
   region->df = --df;
@@ -124,7 +124,7 @@ static void Sample(This *t, cnumber nnew, Region *region,
       res->sigsq = sigsq;
       res->avg = avg;
     }
-    res->err = sqrt(res->sigsq);
+    res->err = sqrtx(res->sigsq);
 
     res->chisq = (sigsq < .9*NOTZERO) ? 0 : c->chisqsum - avg*c->chisum;
       /* This catches the special case where the integrand is constant
@@ -153,7 +153,7 @@ static void Sample(This *t, cnumber nnew, Region *region,
     for( comp = 0, res = region->result;
          comp < t->ncomp; ++comp, ++res ) {
       p += sprintf(p, "%s  \tchisq " REAL " (" COUNT " df)",
-        p0, res->chisq, df);
+        p0, SHOW(res->chisq), df);
       p0 += chars;
     }
 
diff --git a/src/external/libCuba/src/suave/Suave.c b/src/external/libCuba/src/suave/Suave.c
index 79323c14..3bd55a02 100644
--- a/src/external/libCuba/src/suave/Suave.c
+++ b/src/external/libCuba/src/suave/Suave.c
@@ -1,8 +1,8 @@
 /*
 	Suave.c
-		Subregion-adaptive Vegas Monte-Carlo integration
+		Subregion-adaptive Vegas Monte Carlo integration
 		by Thomas Hahn
-		last modified 17 Sep 13 th
+		last modified 28 Nov 14 th
 */
 
 
@@ -15,64 +15,80 @@
 /*********************************************************************/
 
 Extern void EXPORT(Suave)(ccount ndim, ccount ncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber nvec,
   creal epsrel, creal epsabs,
   cint flags, cint seed,
   cnumber mineval, cnumber maxeval,
-  cnumber nnew, creal flatness,
-  cchar *statefile,
+  cnumber nnew, cnumber nmin, creal flatness,
+  cchar *statefile, Spin **pspin,
   count *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = ndim;
   t.ncomp = ncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = nvec;
   t.epsrel = epsrel;
   t.epsabs = epsabs;
-  t.flags = flags;
+  t.flags = MaxVerbose(flags);
   t.seed = seed;
   t.mineval = mineval;
   t.maxeval = maxeval;
   t.nnew = nnew;
+  t.nmin = IMax(nmin, 2);
   t.flatness = flatness;
   t.statefile = statefile;
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
 /*********************************************************************/
 
 Extern void EXPORT(suave)(ccount *pndim, ccount *pncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber *pnvec,
   creal *pepsrel, creal *pepsabs,
   cint *pflags, cint *pseed,
   cnumber *pmineval, cnumber *pmaxeval,
-  cnumber *pnnew, creal *pflatness,
-  cchar *statefile,
+  cnumber *pnnew, cnumber *pnmin, creal *pflatness,
+  cchar *statefile, Spin **pspin,
   count *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob, cint statefilelen)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = *pndim;
   t.ncomp = *pncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = *pnvec;
   t.epsrel = *pepsrel;
   t.epsabs = *pepsabs;
-  t.flags = *pflags;
+  t.flags = MaxVerbose(*pflags);
   t.seed = *pseed;
   t.mineval = *pmineval;
   t.maxeval = *pmaxeval;
   t.nnew = *pnnew;
+  t.nmin = IMax(*pnmin, 2);
   t.flatness = *pflatness;
   CString(t.statefile, statefile, statefilelen);
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
diff --git a/src/external/libCuba/src/suave/decl.h b/src/external/libCuba/src/suave/decl.h
index 8944e60f..4d738b37 100644
--- a/src/external/libCuba/src/suave/decl.h
+++ b/src/external/libCuba/src/suave/decl.h
@@ -2,7 +2,7 @@
 	decl.h
 		Type declarations
 		this file is part of Suave
-		last modified 29 Jul 13 th
+		last modified 25 Nov 14 th
 */
 
 
@@ -35,23 +35,24 @@ typedef struct {
 typedef const Bounds cBounds;
 
 typedef int (*Integrand)(ccount *, creal *, ccount *, real *,
-  void *, creal *, cint *);
+  void *, cnumber *, cint *, creal *, cint *);
 
 typedef struct _this {
   count ndim, ncomp;
 #ifndef MLVERSION
   Integrand integrand;
   void *userdata;
+  number nvec;
 #ifdef HAVE_FORK
-  int ncores, *child;
-  real *frame;
   SHM_ONLY(int shmid;)
+  Spin *spin;
+  real *frame;
 #endif
 #endif
   real epsrel, epsabs;
   int flags, seed;
   number mineval, maxeval;
-  number nnew;
+  number nnew, nmin;
   real flatness;
   cchar *statefile;
   count nregions;
diff --git a/src/external/libCuba/src/vegas/Grid.c b/src/external/libCuba/src/vegas/Grid.c
index 0ad8a79e..7fff0aa9 100644
--- a/src/external/libCuba/src/vegas/Grid.c
+++ b/src/external/libCuba/src/vegas/Grid.c
@@ -2,7 +2,7 @@
 	Grid.c
 		utility functions for the Vegas grid
 		this file is part of Vegas
-		last modified 13 Dec 11 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -70,7 +70,7 @@ static void RefineGrid(cThis *t, Grid grid, Grid margsum)
     real impfun = 0;
     if( margsum[bin] > 0 ) {
       creal r = margsum[bin]*norm;
-      avgperbin += impfun = pow((r - 1)/log(r), 1.5);
+      avgperbin += impfun = powx((r - 1)/log(r), 1.5);
     }
     imp[bin] = impfun;
   }
diff --git a/src/external/libCuba/src/vegas/Integrate.c b/src/external/libCuba/src/vegas/Integrate.c
index 14557871..bb1bae79 100644
--- a/src/external/libCuba/src/vegas/Integrate.c
+++ b/src/external/libCuba/src/vegas/Integrate.c
@@ -2,7 +2,7 @@
 	Integrate.c
 		integrate over the unit hypercube
 		this file is part of Vegas
-		last modified 8 Aug 13 th
+		last modified 13 Mar 15 th
 */
 
 
@@ -31,6 +31,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   if( VERBOSE > 1 ) {
     sprintf(out, "Vegas input parameters:\n"
       "  ndim " COUNT "\n  ncomp " COUNT "\n"
+      ML_NOT("  nvec " NUMBER "\n")
       "  epsrel " REAL "\n  epsabs " REAL "\n"
       "  flags %d\n  seed %d\n"
       "  mineval " NUMBER "\n  maxeval " NUMBER "\n"
@@ -38,7 +39,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       "  nbatch " NUMBER "\n  gridno %d\n"
       "  statefile \"%s\"",
       t->ndim, t->ncomp,
-      t->epsrel, t->epsabs,
+      ML_NOT(t->nvec,)
+      SHOW(t->epsrel), SHOW(t->epsabs),
       t->flags, t->seed,
       t->mineval, t->maxeval,
       t->nstart, t->nincrease, t->nbatch,
@@ -49,7 +51,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   if( BadComponent(t) ) return -2;
   if( BadDimension(t) ) return -1;
 
-  FrameAlloc(t, ShmRm(t));
+  FrameAlloc(t, Master);
   ForkCores(t);
   Alloc(bins, t->nbatch*t->ndim);
 
@@ -68,12 +70,12 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
     t->rng.skiprandom(t, t->neval);
   }
 
-  if( ini ) {
+  if( ini | ZAPSTATE ) {
+    t->neval = 0;
     state->niter = 0;
     state->nsamples = t->nstart;
     FClear(state->cumul);
-    GetGrid(t, state_grid);
-    t->neval = 0;
+    if( ini ) GetGrid(t, state_grid);
   }
 
   /* main iteration loop */
@@ -143,7 +145,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       real avg = sigsq*(c->avgsum += w*c->sum);
 
       c->avg = LAST ? (sigsq = 1/w, c->sum) : avg;
-      c->err = sqrt(sigsq);
+      c->err = sqrtx(sigsq);
       fail |= (c->err > MaxErr(c->avg));
 
       if( state->niter == 0 ) c->guess = c->sum;
@@ -163,7 +165,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       for( c = state->cumul, comp = 0; c < C; ++c )
         oe += sprintf(oe, "\n[" COUNT "] "
           REAL " +- " REAL "  \tchisq " REAL " (" COUNT " df)",
-          ++comp, c->avg, c->err, c->chisq, state->niter);
+          ++comp, SHOW(c->avg), SHOW(c->err),
+          SHOW(c->chisq), state->niter);
       Print(out);
     }
 
@@ -215,8 +218,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
 abort:
   PutGrid(t, state_grid);
   free(bins);
-  WaitCores(t);
-  FrameFree(t);
+  FrameFree(t, Master);
 
   StateRemove(t);
 
diff --git a/src/external/libCuba/src/vegas/Vegas.c b/src/external/libCuba/src/vegas/Vegas.c
index 43a608f7..46a3fd85 100644
--- a/src/external/libCuba/src/vegas/Vegas.c
+++ b/src/external/libCuba/src/vegas/Vegas.c
@@ -1,8 +1,8 @@
 /*
 	Vegas.c
-		Vegas Monte-Carlo integration
+		Vegas Monte Carlo integration
 		by Thomas Hahn
-		last modified 17 Sep 13 th
+		last modified 25 Nov 14 th
 */
 
 
@@ -15,22 +15,27 @@
 /*********************************************************************/
 
 Extern void EXPORT(Vegas)(ccount ndim, ccount ncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber nvec,
   creal epsrel, creal epsabs, cint flags, cint seed,
   cnumber mineval, cnumber maxeval,
-  cnumber nstart, cnumber nincrease, cnumber nbatch,
-  cint gridno, cchar *statefile,
+  cnumber nstart, cnumber nincrease,
+  cnumber nbatch, cint gridno,
+  cchar *statefile, Spin **pspin,
   number *pneval, int *pfail,
   real *integral, real *error, real *prob)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = ndim;
   t.ncomp = ncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = nvec;
   t.epsrel = epsrel;
   t.epsabs = epsabs;
-  t.flags = flags;
+  t.flags = MaxVerbose(flags);
   t.seed = seed;
   t.mineval = mineval;
   t.maxeval = maxeval;
@@ -39,30 +44,38 @@ Extern void EXPORT(Vegas)(ccount ndim, ccount ncomp,
   t.nbatch = nbatch;
   t.gridno = gridno;
   t.statefile = statefile;
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
 /*********************************************************************/
 
 Extern void EXPORT(vegas)(ccount *pndim, ccount *pncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber *pnvec,
   creal *pepsrel, creal *pepsabs, cint *pflags, cint *pseed,
   cnumber *pmineval, cnumber *pmaxeval,
   cnumber *pnstart, cnumber *pnincrease, 
-  cnumber *pnbatch, cint *pgridno, cchar *statefile,
+  cnumber *pnbatch, cint *pgridno,
+  cchar *statefile, Spin **pspin,
   number *pneval, int *pfail,
   real *integral, real *error, real *prob, cint statefilelen)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = *pndim;
   t.ncomp = *pncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = *pnvec;
   t.epsrel = *pepsrel;
   t.epsabs = *pepsabs;
-  t.flags = *pflags;
+  t.flags = MaxVerbose(*pflags);
   t.seed = *pseed;
   t.mineval = *pmineval;
   t.maxeval = *pmaxeval;
@@ -71,8 +84,11 @@ Extern void EXPORT(vegas)(ccount *pndim, ccount *pncomp,
   t.nbatch = *pnbatch;
   t.gridno = *pgridno;
   CString(t.statefile, statefile, statefilelen);
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
diff --git a/src/external/libCuba/src/vegas/decl.h b/src/external/libCuba/src/vegas/decl.h
index 0bd2614d..d528c2b8 100644
--- a/src/external/libCuba/src/vegas/decl.h
+++ b/src/external/libCuba/src/vegas/decl.h
@@ -2,7 +2,7 @@
 	decl.h
 		Type declarations
 		this file is part of Vegas
-		last modified 21 Dec 11 th
+		last modified 21 Jul 14 th
 */
 
 
@@ -29,16 +29,17 @@ typedef struct {
 typedef const Cumulants cCumulants;
 
 typedef int (*Integrand)(ccount *, creal *, ccount *, real *,
-  void *, creal *, cint *);
+  void *, cnumber *, cint *, creal *, cint *);
 
 typedef struct _this {
   count ndim, ncomp;
 #ifndef MLVERSION
   Integrand integrand;
   void *userdata;
+  number nvec;
 #ifdef HAVE_FORK
-  int ncores, *child;
   SHM_ONLY(int shmid;)
+  Spin *spin;
 #endif
 #endif
   real *frame;