void R_approx(double *x, double *y, int *nxy, double *xout, int *nout,
	      int *method, double *yleft, double *yright, double *f)
{
    int i;
    appr_meth M = {0.0, 0.0, 0.0, 0.0, 0}; /* -Wall */

    /* check interpolation method */

    switch(*method) {
    case 1: /* linear */
	    break;
    case 2: /* constant */
      if(!R_FINITE(*f) || *f < 0 || *f > 1)
        error(_("approx(): invalid f value"));
      M.f2 = *f;
      M.f1 = 1 - *f;
      break;
    default:
    	error(_("approx(): invalid interpolation method"));
	    break;
    }

    for(i = 0 ; i < *nxy ; i++)
	    if(ISNA(x[i]) || ISNA(y[i]))
	      error(_("approx(): attempted to interpolate NA values"));

    M.kind = *method;
    M.ylow = *yleft;
    M.yhigh = *yright;

    for(i = 0 ; i < *nout; i++)
	    if(!ISNA(xout[i]))
	      xout[i] = approx1(xout[i], x, y, *nxy, &M);
}

SEXP geoddist_alongpath(SEXP lat, SEXP lon, SEXP a, SEXP f)
{
  if (!isReal(lat))
    error("latitude must be a numeric (floating-point) vector");
  if (!isReal(lon))
    error("longitude must be a numeric (floating-point) vector");
  SEXP res;
  //int n = INTEGER(GET_LENGTH(lat));
  //int nlon = INTEGER(GET_LENGTH(lon));
  int n = GET_LENGTH(lat);
  int nlon = GET_LENGTH(lon);
  if (n != nlon)
    error("lengths of latitude and longitude vectors must match, but they are %d and %d, respectively", n, nlon);
  double *latp = REAL(lat);
  double *lonp = REAL(lon);
  double *ap = REAL(a);
  double *fp = REAL(f);
  PROTECT(res = allocVector(REALSXP, n));
  double *resp = REAL(res);
  double last = 0.0;
  resp[0] = ISNA(lonp[0]) ? NA_REAL : 0.0;
  for (int i = 0; i < n-1; i++) {
    double faz, baz, s;
    if (ISNA(latp[i]) || ISNA(lonp[i]) || ISNA(latp[i+1]) || ISNA(lonp[i+1])) {
      resp[i+1] = NA_REAL;
      last = 0.0; // reset
    } else {
      geoddist_core(latp+i, lonp+i, latp+i+1, lonp+i+1, ap, fp, &faz, &baz, &s);
      resp[i+1] = last + s;
      last = resp[i+1];
    }
  }
  UNPROTECT(1);
  return(res);
}

// TODO: NEED TO DEAL WITH Inf/NA/NaN
void rtexpon_rate(double *x, double *left, double *right, double *rate, int *num)
{
  RNG r;

  #ifdef USE_R
  GetRNGstate();
  #endif

  for(int i=0; i < *num; ++i){
    #ifdef USE_R
    if (i%SAMPCHECK==0) R_CheckUserInterrupt();

    if (ISNAN(left[i]) || ISNAN(right[i]) || ISNAN(rate[i]) || !R_FINITE(left[i]) ||
	ISNA(left[i])  || ISNA(right[i])  || ISNA(rate[i]) ) 
      {
	x[i] = R_NaN; 
	fprintf(stderr, "rtexpon_rate: caught non finite left value: %g; x[i] = %g.\n", left[i], x[i]);
      }

    #endif

    // TODO: Really, I need to check if it is +/- Inf.
    if (MYFINITE(right[i]))
      x[i] = r.texpon_rate(left[i], right[i], rate[i]);
    else 
      x[i] = r.texpon_rate(left[i], rate[i]);

  }

  #ifdef USE_R
  PutRNGstate();
  #endif
}

// variant: means by /nrow rather than /npairs when dealing with missing data
// deviates from R's pairwise.complete.obs, but highly similar.
// good results in context of particular biological validations,
// but formal correctness undetermined
void
pearson_distances_pairwise_complete_obs_variant(
	double * const d,
	const double * const matrix,
	int const nrow,
	int const ncol
){
	std::ptrdiff_t p(0);
	t_float EX(0), EY(0), EXX(0), EYY(0), EXY(0), x(0), y(0);
	for(int col1(0), end(ncol); col1<(end-1); ++col1){
		for(int col2(col1+1); col2<end; ++col2){
			// Pearson correlation distance
			EX=0, EY=0, EXX=0, EYY=0, EXY=0, x=0, y=0;
			unsigned npairs(0);
			for(int row(0); row<nrow; ++row){
				// R indexes its arrays BY COLUMN
				x = matrix[col1*nrow+row];
				y = matrix[col2*nrow+row];
				if(ISNA(x) || ISNA(y)) continue;
				++npairs;
				EX += x;
				EY += y;
				EXX += x*x;
				EYY += y*y;
				EXY += x*y;
			}
			if(npairs<1) d[p++] = 2.0;
			else d[p++] = 1.0 - (EXY - EX*EY/nrow) / sqrt( (EXX - EX*EX/nrow)*(EYY - EY*EY/nrow) );
		}
	}
}

const char
*EncodeComplex(Rcomplex x, int wr, int dr, int er, int wi, int di, int ei,
	       const char *dec)
{
    static char buff[NB];

    /* IEEE allows signed zeros; strip these here */
    if (x.r == 0.0) x.r = 0.0;
    if (x.i == 0.0) x.i = 0.0;

    if (ISNA(x.r) || ISNA(x.i)) {
	snprintf(buff, NB,
		 "%*s", /* was "%*s%*s", R_print.gap, "", */
		 min(wr+wi+2, (NB-1)), CHAR(R_print.na_string));
    } else {
	char Re[NB];
	const char *Im, *tmp;
	int flagNegIm = 0;
	Rcomplex y;
	/* formatComplex rounded, but this does not, and we need to
	   keep it that way so we don't get strange trailing zeros.
	   But we do want to avoid printing small exponentials that
	   are probably garbage.
	 */
	z_prec_r(&y, &x, R_print.digits);
	/* EncodeReal has static buffer, so copy */
	tmp = EncodeReal0(y.r == 0. ? y.r : x.r, wr, dr, er, dec);
	strcpy(Re, tmp);
	if ( (flagNegIm = (x.i < 0)) ) x.i = -x.i;
	Im = EncodeReal0(y.i == 0. ? y.i : x.i, wi, di, ei, dec);
	snprintf(buff, NB, "%s%s%si", Re, flagNegIm ? "-" : "+", Im);
    }
    buff[NB-1] = '\0';
    return buff;
}

void rtnorm(double *x, double *left, double* right, double *mu, double *sig, int *num)
{
  // TODO: NEED TO DEAL WITH Inf/NA/NaN

  RNG r;
  
  #ifdef USE_R
  GetRNGstate();
  #endif
  
  for(int i=0; i < *num; ++i){
    #ifdef USE_R
    if (i%SAMPCHECK==0) R_CheckUserInterrupt(); 

    if (ISNAN(left[i]) || ISNAN(right[i]) || ISNAN(mu[i]) || ISNAN(sig[i]))
      x[i] = R_NaN;

    if (ISNA(left[i]) || ISNA(right[i]) || ISNA(mu[i]) || ISNA(sig[i]))
      x[i] = NA_REAL;
    
    #endif
    
    #ifdef USE_R

    if (left[i] != R_NegInf && right[i] != R_PosInf) {
      x[i] = r.tnorm(left[i], right[i], mu[i], sig[i]);
    } else if (left[i] != R_NegInf && right[i] == R_PosInf) {
      x[i] = r.tnorm(left[i], mu[i], sig[i]);
    } else if (left[i] == R_NegInf && right[i] != R_PosInf) {
      x[i] = -1.0 * r.tnorm(-1.0 * right[i], -1.0 * mu[i], sig[i]);
    } else if (left[i] == R_NegInf && right[i] == R_PosInf) {
      x[i] = r.norm(mu[i], sig[i]);
    } else {
      x[i] = R_NaN;
    }

    #else

    // TODO: Need to adjust so that we deal with +/- inf.

    if (MYFINITE(left[i]) && MYFINITE(right[i]))
      x[i] = r.tnorm(left[i], right[i], mu[i], sig[i]);
    else if (MYFINITE(left[i]))
      x[i] = r.tnorm(left[i], mu[i], sig[i]);
    else if (MYFINITE(right[i]))
      x[i] = -1.0 * r.tnorm(-1.0 * right[i], -1.0 * mu[i], sig[i]);
    else 
      x[i] = r.norm(mu[i], sig[i]);

    #endif
	
  }
  
  #ifdef USE_R
  PutRNGstate();
  #endif
}

/* zeroin2(f, ax, bx, f.ax, f.bx, tol, maxiter) */
SEXP zeroin2(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    double f_ax, f_bx;
    double xmin, xmax, tol;
    int iter;
    SEXP v, res;
    struct callinfo info;

    args = CDR(args);
    PrintDefaults();

    /* the function to be minimized */
    v = CAR(args);
    if (!isFunction(v)) error(_("attempt to minimize non-function"));
    args = CDR(args);

    /* xmin */
    xmin = asReal(CAR(args));
    if (!R_FINITE(xmin)) error(_("invalid '%s' value"), "xmin");
    args = CDR(args);

    /* xmax */
    xmax = asReal(CAR(args));
    if (!R_FINITE(xmax)) error(_("invalid '%s' value"), "xmax");
    if (xmin >= xmax) error(_("'xmin' not less than 'xmax'"));
    args = CDR(args);

    /* f(ax) = f(xmin) */
    f_ax = asReal(CAR(args));
    if (ISNA(f_ax)) error(_("NA value for '%s' is not allowed"), "f.lower");
    args = CDR(args);

    /* f(bx) = f(xmax) */
    f_bx = asReal(CAR(args));
    if (ISNA(f_bx)) error(_("NA value for '%s' is not allowed"), "f.upper");
    args = CDR(args);

    /* tol */
    tol = asReal(CAR(args));
    if (!R_FINITE(tol) || tol <= 0.0) error(_("invalid '%s' value"), "tol");
    args = CDR(args);

    /* maxiter */
    iter = asInteger(CAR(args));
    if (iter <= 0) error(_("'maxiter' must be positive"));

    info.R_env = rho;
    PROTECT(info.R_fcall = lang2(v, R_NilValue)); /* the info used in fcn2() */
    PROTECT(res = allocVector(REALSXP, 3));
    REAL(res)[0] =
	R_zeroin2(xmin, xmax, f_ax, f_bx, (double (*)(double, void*)) fcn2,
		 (void *) &info, &tol, &iter);
    REAL(res)[1] = (double)iter;
    REAL(res)[2] = tol;
    UNPROTECT(2);
    return res;
}

GBMRESULT CDataset::SetData
(
    double *adX,
    int *aiXOrder,
    double *adY,
    double *adOffset,
    double *adWeight,
    double *adMisc,
    int cRows,
    int cCols,
    int *acVarClasses,
    int *alMonotoneVar
)
{
    GBMRESULT hr = GBM_OK;

    if((adX == NULL) || (adY == NULL))
    {
        hr = GBM_INVALIDARG;
        goto Error;
    }

    this->cRows = cRows;
    this->cCols = cCols;

    this->adX = adX;
    this->aiXOrder = aiXOrder;
    this->adY = adY;
    this->adOffset = adOffset;
    this->adWeight = adWeight;
    this->acVarClasses = acVarClasses;
    this->alMonotoneVar = alMonotoneVar;

    if((adOffset != NULL) && !ISNA(*adOffset))
    {
        this->adOffset = adOffset;
        fHasOffset = true;
    }
    else
    {
        this->adOffset = NULL;
        fHasOffset = false;
    }
    if((adMisc != NULL) && !ISNA(*adMisc))
    {
        this->adMisc = adMisc;
    }
    else
    {
        this->adMisc = NULL;
    }

Cleanup:
   return hr;
Error:
    goto Cleanup;
}

// bivariate normal measurement error density
void ou2_dmeasure (double *lik, double *y, double *x, double *p, int give_log, 
		   int *obsindex, int *stateindex, int *parindex, int *covindex,
		   int covdim, double *covar, double t) 
{
  double sd = fabs(TAU);
  double f = 0.0;
  f += (ISNA(Y1)) ? 0.0 : dnorm(Y1,x[X1],sd,1);
  f += (ISNA(Y2)) ? 0.0 : dnorm(Y2,x[X2],sd,1);
  *lik = (give_log) ? f : exp(f);
}

SEXP attribute_hidden complex_math2(SEXP call, SEXP op, SEXP args, SEXP env)
{
    R_xlen_t i, n, na, nb;
    Rcomplex ai, bi, *a, *b, *y;
    SEXP sa, sb, sy;
    Rboolean naflag = FALSE;
    cm2_fun f;

    switch (PRIMVAL(op)) {
    case 0: /* atan2 */
	f = z_atan2; break;
    case 10001: /* round */
	f = z_rround; break;
    case 2: /* passed from do_log1arg */
    case 10:
    case 10003: /* passed from do_log */
	f = z_logbase; break;
    case 10004: /* signif */
	f = z_prec; break;
    default:
	errorcall_return(call, _("unimplemented complex function"));
    }

    PROTECT(sa = coerceVector(CAR(args), CPLXSXP));
    PROTECT(sb = coerceVector(CADR(args), CPLXSXP));
    na = XLENGTH(sa); nb = XLENGTH(sb);
    if ((na == 0) || (nb == 0)) {
        UNPROTECT(2);
        return(allocVector(CPLXSXP, 0));
    }
    n = (na < nb) ? nb : na;
    PROTECT(sy = allocVector(CPLXSXP, n));
    a = COMPLEX(sa); b = COMPLEX(sb); y = COMPLEX(sy);
    for (i = 0; i < n; i++) {
	ai = a[i % na]; bi = b[i % nb];
	if(ISNA(ai.r) && ISNA(ai.i) &&
	   ISNA(bi.r) && ISNA(bi.i)) {
	    y[i].r = NA_REAL; y[i].i = NA_REAL;
	} else {
	    f(&y[i], &ai, &bi);
	    if ( (ISNAN(y[i].r) || ISNAN(y[i].i)) &&
		 !(ISNAN(ai.r) || ISNAN(ai.i) || ISNAN(bi.r) || ISNAN(bi.i)) )
		naflag = TRUE;
	}
    }
    if (naflag)
	warningcall(call, "NaNs produced in function \"%s\"", PRIMNAME(op));
    if(n == na) {
	DUPLICATE_ATTRIB(sy, sa);
    } else if(n == nb) {
	DUPLICATE_ATTRIB(sy, sb);
    }
    UNPROTECT(3);
    return sy;
}

static void printComplexMatrix(SEXP sx, int offset, int r_pr, int r, int c,
                               SEXP rl, SEXP cl, const char *rn, const char *cn)
{
    _PRINT_INIT_rl_rn;
    Rcomplex *x = COMPLEX(sx) + offset;

    int *dr = (int *) R_alloc(c, sizeof(int)),
         *er = (int *) R_alloc(c, sizeof(int)),
          *wr = (int *) R_alloc(c, sizeof(int)),
           *di = (int *) R_alloc(c, sizeof(int)),
            *ei = (int *) R_alloc(c, sizeof(int)),
             *wi = (int *) R_alloc(c, sizeof(int));

    /* Determine the column widths */

    for (j = 0; j < c; j++) {
        formatComplex(&x[j * r], (R_xlen_t) r,
                      &wr[j], &dr[j], &er[j],
                      &wi[j], &di[j], &ei[j], 0);
        _PRINT_SET_clabw;
        w[j] = wr[j] + wi[j] + 2;
        if (w[j] < clabw)
            w[j] = clabw;
        w[j] += R_print.gap;
    }

    _PRINT_DEAL_c_eq_0;
    while (jmin < c) {
        width = rlabw;
        do {
            width += w[jmax];
            jmax++;
        }
        while (jmax < c && width + w[jmax] < R_print.width);

        _PRINT_ROW_LAB;

        for (j = jmin; j < jmax ; j++)
            MatrixColumnLabel(cl, j, w[j]);
        for (i = 0; i < r_pr; i++) {
            MatrixRowLabel(rl, i, rlabw, lbloff);
            for (j = jmin; j < jmax; j++) {
                if (ISNA(x[i + j * r].r) || ISNA(x[i + j * r].i))
                    Rprintf("%s", EncodeReal(NA_REAL, w[j], 0, 0, OutDec));
                else
                    Rprintf("%s",
                            EncodeComplex(x[i + j * r],
                                          wr[j] + R_print.gap, dr[j], er[j],
                                          wi[j], di[j], ei[j], OutDec));
            }
        }
        Rprintf("\n");
        jmin = jmax;
    }
}

SEXP adjRatios (SEXP split, SEXP div, SEXP close) {

    /* Initialize REAL pointers to function arguments */
    double *real_close = REAL(close);
    double *real_split = REAL(split);
    double *real_div   = REAL(div);
    
    /* Initalize loop and PROTECT counters */
    int i, P = 0;
    /* Initalize object length (NOTE: all arguments are the same length) */
    int N = length(close);

    /* Initalize result R objects */
    SEXP result;    PROTECT(result  = allocVector(VECSXP, 2)); P++;
    SEXP s_ratio;   PROTECT(s_ratio = allocVector(REALSXP,N)); P++;
    SEXP d_ratio;   PROTECT(d_ratio = allocVector(REALSXP,N)); P++;
    
    /* Initialize REAL pointers to R objects and set their last value to '1' */
    double *rs_ratio = REAL(s_ratio);
    double *rd_ratio = REAL(d_ratio);
    rs_ratio[N-1] = 1;
    rd_ratio[N-1] = 1;

    /* Loop over split/div vectors from newest period to oldest */
    for(i = N-1; i > 0; i--) {
        /* Carry newer ratio value backward */
        if(ISNA(real_split[i])) {
            rs_ratio[i-1] = rs_ratio[i];
        /* Update split ratio */
        } else {
            rs_ratio[i-1] = rs_ratio[i] * real_split[i];
        }
        /* Carry newer ratio value backward */
        if(ISNA(real_div[i])) {
            rd_ratio[i-1] = rd_ratio[i];
        } else {
        /* Update dividend ratio */
            rd_ratio[i-1] = rd_ratio[i] *
                (1.0 - real_div[i] / real_close[i-1]);
        }
    }
    
    /* Assign results to list */
    SET_VECTOR_ELT(result, 0, s_ratio);
    SET_VECTOR_ELT(result, 1, d_ratio);

    /* UNPROTECT R objects and return result */
    UNPROTECT(P);
    return(result);
}

int Engine::gradient()
{
    double repsL, repsR;
    std::vector<double> x1(xBuffer_.size());
    std::vector<double> x2(xBuffer_.size());

    for (unsigned int i = 0; i < xBuffer_.size(); ++i)
    {
        std::copy(xBuffer_.begin(), xBuffer_.end(), x1.begin());
        std::copy(xBuffer_.begin(), xBuffer_.end(), x2.begin());
        repsL = reps_;
        repsR = reps_;

        x1[i] = xBuffer_[i] + repsR;
        if (x1[i] > upper_[i])
        {
            x1[i] = upper_[i];
            repsR = x1[i] - xBuffer_[i];
        }

        x2[i] = xBuffer_[i] - repsL;
        if (x2[i] < lower_[i])
        {
            x2[i] = lower_[i];
            repsL = xBuffer_[i] - x2[i];
        }
        g_[i] = (fObjective(x1) - fObjective(x2)) / (repsL + repsR);

        if (ISNA(g_[i]) || !R_FINITE(g_[i]))
        {
            g_[i] = 101.;
        }
    }
    return 0;
}

attribute_hidden
const char *EncodeReal2(double x, int w, int d, int e)
{
    static char buff[NB];
    char fmt[20];

    /* IEEE allows signed zeros (yuck!) */
    if (x == 0.0) x = 0.0;
    if (!R_FINITE(x)) {
	if(ISNA(x)) snprintf(buff, NB, "%*s", w, CHAR(R_print.na_string));
	else if(ISNAN(x)) snprintf(buff, NB, "%*s", w, "NaN");
	else if(x > 0) snprintf(buff, NB, "%*s", w, "Inf");
	else snprintf(buff, NB, "%*s", w, "-Inf");
    }
    else if (e) {
	if(d) {
	    sprintf(fmt,"%%#%d.%de", min(w, (NB-1)), d);
	    snprintf(buff, NB, fmt, x);
	}
	else {
	    sprintf(fmt,"%%%d.%de", min(w, (NB-1)), d);
	    snprintf(buff, NB, fmt, x);
	}
    }
    else { /* e = 0 */
	sprintf(fmt,"%%#%d.%df", min(w, (NB-1)), d);
	snprintf(buff, NB, fmt, x);
    }
    buff[NB-1] = '\0';
    return buff;
}

SEXP attribute_hidden StringFromReal(double x, int *warn)
{
    int w, d, e;
    formatReal(&x, 1, &w, &d, &e, 0);
    if (ISNA(x)) return NA_STRING;
    else return mkChar(EncodeRealDrop0(x, w, d, e, OutDec));
}

Rboolean tmin(double *x, index_type n, double *value, Rboolean narm, 
              double NA_VALUE)
{
  double s = 0.0; /* -Wall */
  Rboolean updated = (Rboolean)FALSE;

  for (index_type i = 0; i < n; i++) {
    if (ISNAN(x[i])) {/* Na(N) */
      if (!narm) {
        if(!ISNA(s)) s = x[i]; /* so any NA trumps all NaNs */
        if(!updated) updated = (Rboolean)TRUE;
      } // narm = TRUE then nothing.
    } // When for all i: (ISNAN && narm) = TRUE then updated = FALSE
    else if (!updated || x[i] < s) {/* Never true if s is NA/NaN */
      s = x[i];
      if(!updated) updated = (Rboolean)TRUE;
    }
  }
  
  if (!updated) {
    // *value = NA_REAL;
    if (narm) {  // To Make consistent w/ R-min
      *value = R_PosInf;
    } else {
      *value = NA_REAL;
    }
  } else {
    *value = s;
  }
  return((Rboolean)TRUE);
}

// [[register]]
SEXP any_na( SEXP x ) {
	SEXP out;
	PROTECT(out = allocVector(LGLSXP, 1));
	int len = length(x);
	switch( TYPEOF(x) ) {
	case REALSXP: {
		double* ptr = REAL(x);
		for( int i=0; i < len; ++i ) {
			if( ISNA( ptr[i] ) || ISNAN( ptr[i] ) ) {
				LOGICAL(out)[0] = TRUE;
				UNPROTECT(1);
				return out;
			}
		}
		LOGICAL(out)[0] = FALSE;
		UNPROTECT(1);
		return out;
	}
	case INTSXP: {
		int* ptr = INTEGER(x);
		for( int i=0; i < len; ++i ) {
			if( ptr[i] == NA_INTEGER ) {
				LOGICAL(out)[0] = TRUE;
				UNPROTECT(1);
				return out;
			}
		}
		LOGICAL(out)[0] = FALSE;
		UNPROTECT(1);
		return out;
	}
	case LGLSXP: {
		int* ptr = LOGICAL(x);
		for( int i=0; i < len; ++i ) {
			if( ptr[i] == NA_LOGICAL ) {
				LOGICAL(out)[0] = TRUE;
				UNPROTECT(1);
				return out;
			}
		}
		LOGICAL(out)[0] = FALSE;
		UNPROTECT(1);
		return out;
	}
	case STRSXP: {
		for( int i=0; i < len; ++i ) {
			if( STRING_ELT(x, i) == NA_STRING ) {
				LOGICAL(out)[0] = TRUE;
				UNPROTECT(1);
				return out;
			}
		}
		LOGICAL(out)[0] = FALSE;
		UNPROTECT(1);
		return out;
	}
	}
	error("argument is of incompatible type '%s'", type2char( TYPEOF(x) ) );
	return x;
}

static Rboolean cmath1(double complex (*f)(double complex),
		       Rcomplex *x, Rcomplex *y, R_xlen_t n)
{
    R_xlen_t i;
    Rboolean naflag = FALSE;
    for (i = 0 ; i < n ; i++) {
	if (ISNA(x[i].r) || ISNA(x[i].i)) {
	    y[i].r = NA_REAL; y[i].i = NA_REAL;
	} else {
	    SET_C99_COMPLEX(y, i, f(toC99(x + i)));
	    if ( (ISNAN(y[i].r) || ISNAN(y[i].i)) &&
		!(ISNAN(x[i].r) || ISNAN(x[i].i)) ) naflag = TRUE;
	}
    }
    return naflag;
}

signed char CNodeContinuous::WhichNode
(
    double *adX,
    unsigned long cRow,
    unsigned long cCol,
    unsigned long iRow
)
{
    signed char ReturnValue = 0;
    double dX = adX[iSplitVar*cRow + iRow];

    if(!ISNA(dX))
    {
        if(dX < dSplitValue)
        {
            ReturnValue = -1;
        }
        else
        {
            ReturnValue = 1;
        }
    }
    // if missing value returns 0

    return ReturnValue;
}

extern SEXP gini_matrix(SEXP value,
		 SEXP nrowR,
		 SEXP ncolR) {
  double *pmat = REAL(value);
  int nrow = INTEGER(nrowR)[0];
  int ncol = INTEGER(ncolR)[0];
  double rowvec[ncol];
  double curr;
  int i, j, k;
  
  SEXP res;
  PROTECT(res = allocVector(REALSXP,
			    nrow));
  for (i=0; i<nrow; i++) {
    k=0;
    for(j=0; j<ncol; j++) {
      curr=pmat[i+j*nrow];
      if(!ISNA(curr)) {
	rowvec[k++]=curr;
      }
    }
    REAL(res)[i] = stat_gini(rowvec, k);
  }
  UNPROTECT(1);
  return(res);
}