void massdist3d(double *x1, double *x2,	double *x3, int *n, 
		double *a1, double *a2, double *a3, 
		double *b1, double *b2, double *b3,
		int *M1, int *M2, int *M3, double *weight, double *est)
{
  double fx1, fx2, fx3, xdelta1, xdelta2, xdelta3, xpos1, xpos2, xpos3, wi;   
  int i, ix1, ix2, ix3, ixmax1, ixmin1, ixmax2, ixmax3, ixmin2, ixmin3, MM1, MM2, MM3;
  
  MM1 = M1[0];
  MM2 = M2[0];
  MM3 = M3[0];
  ixmin1 = 0;
  ixmax1 = MM1 - 2;
  ixmin2 = 0;
  ixmax2 = MM2 - 2;
  ixmin3 = 0;
  ixmax3 = MM3 - 2;
  xdelta1 = (b1[0] - a1[0]) / (MM1 - 1);
  xdelta2 = (b2[0] - a2[0]) / (MM2 - 1);
  xdelta3 = (b3[0] - a3[0]) / (MM3 - 1);
 
  // set all est = 0 
  for (i=0; i < MM1*MM2*MM3; i++)  
    est[i] = 0.0;

  // assign linear binning weights
  for(i=0; i < n[0]; i++) {
    if(R_FINITE(x1[i]) && R_FINITE(x2[i]) && R_FINITE(x3[i])) {
      xpos1 = (x1[i] - a1[0]) / xdelta1;
      xpos2 = (x2[i] - a2[0]) / xdelta2;
      xpos3 = (x3[i] - a3[0]) / xdelta3;
      ix1 = floor(xpos1);
      ix2 = floor(xpos2);
      ix3 = floor(xpos3);
      fx1 = xpos1 - ix1;
      fx2 = xpos2 - ix2;
      fx3 = xpos3 - ix3;
      wi = weight[i];   
      
      if(ixmin1 <= ix1 && ix1 <= ixmax1 && ixmin2 <= ix2 && ix2 <= ixmax2 && ixmin3 <= ix3 && ix3 <= ixmax3) {
	est[ix3*MM1*MM2 + ix2*MM1 + ix1]             += wi*(1-fx1)*(1-fx2)*(1-fx3);   
	est[ix3*MM1*MM2 + ix2*MM1 + ix1 + 1]         += wi*fx1*(1-fx2)*(1-fx3);
	est[ix3*MM1*MM2 + (ix2+1)*MM1 + ix1]         += wi*(1-fx1)*fx2*(1-fx3);   
	est[ix3*MM1*MM2 + (ix2+1)*MM1 + ix1 + 1]     += wi*fx1*fx2*(1-fx3);
	est[(ix3+1)*MM1*MM2 + ix2*MM1 + ix1]         += wi*(1-fx1)*(1-fx2)*fx3;   
	est[(ix3+1)*MM1*MM2 + ix2*MM1 + ix1 + 1]     += wi*fx1*(1-fx2)*fx3;
	est[(ix3+1)*MM1*MM2 + (ix2+1)*MM1 + ix1]     += wi*(1-fx1)*fx2*fx3;   
	est[(ix3+1)*MM1*MM2 + (ix2+1)*MM1 + ix1 + 1] += wi*fx1*fx2*fx3;
      }
    }
  }
}

double fprec(double x, double digits)
{
    double l10, pow10, sgn, p10, P10;
    int e10, e2, do_round, dig;
    /* Max.expon. of 10 (=308.2547) */
    const double max10e = numeric_limits<double>::max_exponent * M_LOG10_2;

#ifdef IEEE_754
    if (ISNAN(x) || ISNAN(digits))
	return x + digits;
    if (!R_FINITE(x)) return x;
    if (!R_FINITE(digits)) {
	if(digits > 0) return x;
	else return 0;
    }
#endif
    if(x == 0) return x;
    dig = (int)FLOOR(digits+0.5);
    if (dig > MAX_DIGITS) {
	return x;
    } else if (dig < 1)
	dig = 1;

    sgn = 1.0;
    if(x < 0.0) {
	sgn = -sgn;
	x = -x;
    }
    l10 = log10(x);
    e10 = (int)(dig-1-FLOOR(l10));
    if(fabs(l10) < max10e - 2) {
	p10 = 1.0;
	if(e10 > max10e) {
	    p10 =  std::pow(10., e10-max10e);
	    e10 = static_cast<int>(max10e);
	} else if(e10 < - max10e) {
	    p10 =  std::pow(10., e10+max10e);
	    e10 = static_cast<int>(-max10e);	
	}
	pow10 = std::pow(10., e10);
	return(sgn*(FLOOR((x*pow10)*p10+0.5)/pow10)/p10);
    } else { /* -- LARGE or small -- */
	do_round = max10e - l10	 >= std::pow(10., -dig);
	e2 = dig + ((e10>0)? 1 : -1) * MAX_DIGITS;
	p10 = std::pow(10., e2);	x *= p10;
	P10 = std::pow(10., e10-e2);	x *= P10;
	/*-- p10 * P10 = 10 ^ e10 */
	if(do_round) x += 0.5;
	x = FLOOR(x) / p10;
	return(sgn*x/P10);
    }
}

double qsignrank(double x, double n, int lower_tail, int log_p)
{
    double f, p;

#ifdef IEEE_754
    if (ISNAN(x) || ISNAN(n))
	return(x + n);
#endif
    if (!R_FINITE(x) || !R_FINITE(n))
	ML_ERR_return_NAN;
    R_Q_P01_check(x);

    n = floor(n + 0.5);
    if (n <= 0)
	ML_ERR_return_NAN;

    if (x == R_DT_0)
	return(0);
    if (x == R_DT_1)
	return(n * (n + 1) / 2);

    if(log_p || !lower_tail)
	x = R_DT_qIv(x); /* lower_tail,non-log "p" */

    int nn = (int) n;
    w_init_maybe(nn);
    f = exp(- n * M_LN2);
    p = 0;
    int q = 0;
    if (x <= 0.5) {
	x = x - 10 * DBL_EPSILON;
	for (;;) {
	    p += csignrank(q, nn) * f;
	    if (p >= x)
		break;
	    q++;
	}
    }
    else {
	x = 1 - x + 10 * DBL_EPSILON;
	for (;;) {
	    p += csignrank(q, nn) * f;
	    if (p > x) {
		q = (int)(n * (n + 1) / 2 - q);
		break;
	    }
	    q++;
	}
    }

    return(q);
}

double qunif(double p, double a, double b, int lower_tail, int log_p)
{
#ifdef IEEE_754
    if (ISNAN(p) || ISNAN(a) || ISNAN(b))
	return p + a + b;
#endif
    R_Q_P01_check(p);
    if (!R_FINITE(a) || !R_FINITE(b)) ML_ERR_return_NAN;
    if (b < a) ML_ERR_return_NAN;
    if (b == a) return a;

    return a + R_DT_qIv(p) * (b - a);
}

double pnchisq(double x, double f, double theta, int lower_tail, int log_p)
{
#ifdef IEEE_754
    if (ISNAN(x) || ISNAN(f) || ISNAN(theta))
	return x + f + theta;
    if (!R_FINITE(f) || !R_FINITE(theta))
	ML_ERR_return_NAN;
#endif

    if (f < 0. || theta < 0.) ML_ERR_return_NAN;

    return (R_DT_val(pnchisq_raw(x, f, theta, 1e-12, 8*DBL_EPSILON, 1000000)));
}

double beta(double a, double b)
{
#ifdef NOMORE_FOR_THREADS
    static double xmin, xmax = 0;/*-> typically = 171.61447887 for IEEE */
    static double lnsml = 0;/*-> typically = -708.3964185 */

    if (xmax == 0) {
	    gammalims(&xmin, &xmax);
	    lnsml = log(d1mach(1));
    }
#else
/* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
 *   xmin, xmax : see ./gammalims.c
 *   lnsml = log(DBL_MIN) = log(2 ^ -1022) = -1022 * log(2)
*/
# define xmin  -170.5674972726612
# define xmax   171.61447887182298
# define lnsml -708.39641853226412
#endif


#ifdef IEEE_754
    /* NaNs propagated correctly */
    if(ISNAN(a) || ISNAN(b)) return a + b;
#endif

    if (a < 0 || b < 0)
	ML_ERR_return_NAN
    else if (a == 0 || b == 0)
	return ML_POSINF;
    else if (!R_FINITE(a) || !R_FINITE(b))
	return 0;

    if (a + b < xmax) {/* ~= 171.61 for IEEE */
//	return gammafn(a) * gammafn(b) / gammafn(a+b);
	/* All the terms are positive, and all can be large for large
	   or small arguments.  They are never much less than one.
	   gammafn(x) can still overflow for x ~ 1e-308, 
	   but the result would too. 
	*/
	return (1 / gammafn(a+b)) * gammafn(a) * gammafn(b);
    } else {
	double val = lbeta(a, b);
	if (val < lnsml) {
	    /* a and/or b so big that beta underflows */
	    ML_ERROR(ME_UNDERFLOW, "beta");
	    /* return ML_UNDERFLOW; pointless giving incorrect value */
	}
	return exp(val);
    }
}

double mgamma(double order, double shape, double scale, int give_log)
{
    if (!R_FINITE(shape) ||
        !R_FINITE(scale) ||
        !R_FINITE(order) ||
        shape <= 0.0 ||
        scale <= 0.0)
        return R_NaN;

    if (order <= -shape)
	return R_PosInf;

    return R_pow(scale, order) * gammafn(order + shape) / gammafn(shape);
}

double mgfgamma(double x, double shape, double scale, int give_log)
{
    if (!R_FINITE(shape) ||
        !R_FINITE(scale) ||
        shape <= 0.0 ||
        scale <= 0.0 ||
        scale * x > 1.)
        return R_NaN;

    if (x == 0.0)
        return ACT_D_exp(0.0);

    return ACT_D_exp(-shape * log1p(-scale * x));
}

double qcauchy(double p, double location, double scale,
               int lower_tail, int log_p)
{
#ifdef IEEE_754
    if (ISNAN(p) || ISNAN(location) || ISNAN(scale))
        return p + location + scale;
#endif
    if(!R_FINITE(p) || !R_FINITE(location) || !R_FINITE(scale))
        ML_ERR_return_NAN;
    R_Q_P01_check(p);
    if (scale <= 0) ML_ERR_return_NAN;

    return location + scale * tan(M_PI * (R_DT_qIv(p) - 0.5));
}

double rinvparalogis(double shape, double scale)
{
    double tmp;

    if (!R_FINITE(shape) ||
        !R_FINITE(scale) ||
        shape <= 0.0 ||
        scale <= 0.0)
        return R_NaN;;

    tmp = -1.0 / shape;

    return scale * R_pow(R_pow(unif_rand(), tmp) - 1.0, tmp);
}

double runif(rng_t unif_rand, double a, double b)
{
    if (!R_FINITE(a) || !R_FINITE(b) || b < a)	ML_ERR_return_NAN;

    if (a == b)
	return a;
    else {
	double u;
	/* This is true of all builtin generators, but protect against
	   user-supplied ones */
	do {u = unif_rand();} while (u <= 0 || u >= 1);
	return a + (b - a) * u;
    }
}

double qinvpareto(double p, double shape, double scale, int lower_tail,
                  int log_p)
{
    if (!R_FINITE(shape) ||
        !R_FINITE(scale) ||
        shape <= 0.0 ||
        scale <= 0.0)
        return R_NaN;;

    ACT_Q_P01_boundaries(p, 0, R_PosInf);
    p = ACT_D_qIv(p);

    return scale / (R_pow(ACT_D_Lval(p), -1.0 / shape) - 1.0);
}

SEXP mom_calc_int2(SEXP is, SEXP m, SEXP nb, SEXP weights, SEXP card) {
    SEXP Omega;
    int hm = INTEGER_POINTER(m)[0];
    int n = length(card);
    double *eta, *zeta, *omega, sum, res;
    int i, ii, j, k1, k2, k3;
    int iis = length(is);

    omega = (double *) R_alloc((size_t) hm, sizeof(double));
    eta = (double *) R_alloc((size_t) n, sizeof(double));
    zeta = (double *) R_alloc((size_t) n, sizeof(double));
    for (j=0; j<hm; j++) omega[j] = 0.0;

    for (ii=0; ii<iis; ii++) {
        R_CheckUserInterrupt();
        i = INTEGER_POINTER(is)[ii]-ROFFSET;
        for (j=0; j<n; j++) eta[j] = 0.0;
        eta[i] = 1.0;
        for (j=1; j<hm; j=j+2) {
            for (k1=0; k1<n; k1++) {
                k3 = INTEGER_POINTER(card)[k1];
                if (k3 == 0) {
                    zeta[k1] = 0.0;
                } else {
                    sum = 0.0;
                    for (k2=0; k2<k3; k2++) {
                        sum += eta[INTEGER_POINTER(VECTOR_ELT(nb, k1))[k2]
                            - ROFFSET] * NUMERIC_POINTER(VECTOR_ELT(weights,
                            k1))[k2];
                    }
                    zeta[k1] = sum;
                }
            }
            res = F77_CALL(ddot)(&n, zeta, &c__1, eta, &c__1);
            if (R_FINITE(res)) omega[(j-1)] += res;
            else error("non-finite dot product %d, %d", i, j);
            res = F77_CALL(ddot)(&n, zeta, &c__1, zeta, &c__1);
            if (R_FINITE(res)) omega[j] += res;
            else error("non-finite dot product %d, %d", i, j);
            for (k1=0; k1<n; k1++) eta[k1] = zeta[k1];
        }
    }

    PROTECT(Omega = NEW_NUMERIC(hm));
    for (j=0; j<hm; j++) NUMERIC_POINTER(Omega)[j] = omega[j];

    UNPROTECT(1);
    return(Omega);
}

//compute the observed hands
// r[i] is the number of hands with i+1 (different) value(s)
void pokerTest(int *hands, int nbh, int d, int *res)
{
        int i, j; //loop indexes
        int nbzero; //zero counter
        
        int * temp = (int *) R_alloc(d, sizeof(int) );
	
        if (!R_FINITE(nbh) || !R_FINITE(d))
                error(_("non finite argument"));
            
        //init
        for(j = 0; j < d; j++)
                res[j] = 0;
                   
        for(i = 0; i < nbh; i++) 
        {
                //erase previous line
                for(j = 0; j < d; j++)
                        temp[j] = 0;
            
                //browse the i+1th hand
                for(j = 0; j < d; j++)
                { 
                        //if(hands[i + j * nb] > -1 && hands[i + j * nb] <d)
                            temp[ hands[i + j * nbh] ] ++;
                        //else
                            //error(_("internal error in pokertest"));
                }
            /*
            Rprintf("temp : ");
            for(j = 0; j < d; j++)
                Rprintf(" %d\t", temp[j]);
            Rprintf("\n");
            */
                
                nbzero = 0;
            
                //find the i+1 th hand
                for(j = 0; j < d; j++)
                {
                        if(temp[j] == 0)
                                nbzero++;
                }
            //nb of different value is d-nbzero
            res[d - nbzero - 1] ++;
                    
        }
      
}

double rnchisq(double df, double lambda)
{
    if (!R_FINITE(df) || !R_FINITE(lambda) || df < 0. || lambda < 0.)
	ML_ERR_return_NAN;

    if(lambda == 0.) {
	return (df == 0.) ? 0. : rgamma(df / 2., 2.);
    }
    else {
	double r = rpois( lambda / 2.);
	if (r > 0.)  r = rchisq(2. * r);
	if (df > 0.) r += rgamma(df / 2., 2.);
	return r;
    }
}

/* winProgressBar(width, title, label, min, max, initial) */
SEXP winProgressBar(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP tmp, ptr;
    int width, iv;
    double d;
    const char *title, *label;
    winprogressbar *pbar;
    Rboolean haveLabel;

    args = CDR(args);
    pbar = Calloc(1, winprogressbar);
    width = asInteger(CAR(args)); args = CDR(args);
    if(width == NA_INTEGER || width < 0) width = 200;
    tmp = CAR(args); args = CDR(args);
    if(!isString(tmp) || length(tmp) < 1 || STRING_ELT(tmp, 0) == NA_STRING)
	errorcall(call, "invalid '%s' argument", "title");
    title = translateChar(STRING_ELT(tmp, 0));
    tmp = CAR(args); args = CDR(args);
    if(!isString(tmp) || length(tmp) < 1 || STRING_ELT(tmp, 0) == NA_STRING)
	errorcall(call, "invalid '%s' argument", "Label");
    label = translateChar(STRING_ELT(tmp, 0));
    haveLabel = strlen(label) > 0;
    d = asReal(CAR(args)); args = CDR(args);
    if (!R_FINITE(d)) errorcall(call, "invalid '%s' argument", "min");
    pbar->min = d;
    d = asReal(CAR(args)); args = CDR(args);
    if (!R_FINITE(d)) errorcall(call, "invalid '%s' argument", "max");
    pbar->max = d;
    d = asReal(CAR(args)); args = CDR(args);
    if (!R_FINITE(d)) errorcall(call, "invalid '%s' argument", "initial");
    pbar->val = d;

    pbar->width = width;
    pbar->wprog = newwindow(title, rect(0, 0, width+40, haveLabel ? 100: 80),
			    Titlebar | Centered);
    setbackground(pbar->wprog, dialog_bg());
    if(haveLabel)
	pbar->lab = newlabel(label, rect(10, 15, width+20, 25), AlignCenter);
    pbar->pb = newprogressbar(rect(20, haveLabel ? 50 : 30, width, 20),
			      0, width, 1, 1);
    iv = pbar->width * (pbar->val - pbar->min)/(pbar->max - pbar->min);
    setprogressbar(pbar->pb, iv);
    show(pbar->wprog);
    ptr = R_MakeExternalPtr(pbar, install("winProgressBar"), R_NilValue);
    R_RegisterCFinalizerEx(ptr, pbarFinalizer, TRUE);

    return ptr;
}

JNIEXPORT jboolean JNICALL Java_r_ext_SystemLibs_fmod
  (JNIEnv *jenv, jclass jcls, jdoubleArray xArg, jdoubleArray yArg, jdoubleArray resArg, jint size) {
  
  double *x = (*jenv)->GetPrimitiveArrayCritical(jenv, xArg, 0);
  double *y = (*jenv)->GetPrimitiveArrayCritical(jenv, yArg, 0);
  double *res = (*jenv)->GetPrimitiveArrayCritical(jenv, resArg, 0);
  int i;
  int lostAccuracy = 0;

  for (i = 0; i < size; i++) {
    double a = x[i];
    double b = y[i];

    if (b == 0) {  // LICENSE: transcribed from GNU-R, which is licensed under GPL
      res[i] = R_NaN;
    } else {
      double q = a / b;
      double tmp = a - floor(q) * b;
      if ((fabs(q) > 1/DBL_EPSILON ) && R_FINITE(q)) {   // R_AccuracyInfo.eps
        lostAccuracy = 1;
      }
      res[i] = tmp - floor(tmp/b) * b;  
    }    
  }
   
  (*jenv)->ReleasePrimitiveArrayCritical(jenv, xArg, x, 0);
  (*jenv)->ReleasePrimitiveArrayCritical(jenv, yArg, y, 0);
  (*jenv)->ReleasePrimitiveArrayCritical(jenv, resArg, res, 0);  
  
  return lostAccuracy ? JNI_TRUE : JNI_FALSE;
}

double fround(double x, double digits) {
#define MAX_DIGITS DBL_MAX_10_EXP
    /* = 308 (IEEE); was till R 0.99: (DBL_DIG - 1) */
    /* Note that large digits make sense for very small numbers */
    LDOUBLE pow10, sgn, intx;
    int dig;

#ifdef IEEE_754
    if (ISNAN(x) || ISNAN(digits))
	return x + digits;
    if(!R_FINITE(x)) return x;
#endif

    if (digits > MAX_DIGITS)
	digits = MAX_DIGITS;
    dig = (int)floor(digits + 0.5);
    if(x < 0.) {
	sgn = -1.;
	x = -x;
    } else
	sgn = 1.;
    if (dig == 0) {
	return sgn * R_rint(x);
    } else if (dig > 0) {
        pow10 = R_pow_di(10., dig);
	intx = floor(x);
	return sgn * (intx + R_rint((x-intx) * pow10) / pow10);
    } else {
        pow10 = R_pow_di(10., -dig);
        return sgn * R_rint(x/pow10) * pow10;
    }
}

double logspaceAdd(const double loga, const double logb) {
    if (!R_FINITE(loga))
        return logb;
    if (loga > logb)
        return logspaceAdd(logb, loga);
    return logb + log1p(exp(loga - logb));
}

double psignrank(double x, double n, int lower_tail, int log_p)
{
    int i;
    double f, p;

#ifdef IEEE_754
    if (ISNAN(x) || ISNAN(n))
    return(x + n);
#endif
    if (!R_FINITE(n)) ML_ERR_return_NAN;
    n = floor(n + 0.5);
    if (n <= 0) ML_ERR_return_NAN;

    x = floor(x + 1e-7);
    if (x < 0.0)
	return(R_DT_0);
    if (x >= n * (n + 1) / 2)
	return(R_DT_1);

    w_init_maybe(n);
    f = exp(- n * M_LN2);
    p = 0;
    if (x <= (n * (n + 1) / 4)) {
	for (i = 0; i <= x; i++)
	    p += csignrank(i, n) * f;
    }
    else {
	x = n * (n + 1) / 2 - x;
	for (i = 0; i < x; i++)
	    p += csignrank(i, n) * f;
	lower_tail = !lower_tail; /* p = 1 - p; */
    }

    return(R_DT_val(p));
} /* psignrank() */