//extern "C"
SEXP log_marg_A0k(SEXP WpostR, SEXP A0R, SEXP N2R, SEXP consttermR,
			     SEXP bfR, SEXP UTR, SEXP TinvR, SEXP dfR, SEXP n0R)
{
  int i, j, *dTi, db, m, N2, df, len;
  N2=INTEGER(N2R)[0]; df=INTEGER(dfR)[0]; m=INTEGER(coerceVector(listElt(WpostR,"m"),INTSXP))[0];

  double *pbfi, *lpa0, *cterm, lN2, tol, maxvlog, lqlog;
  lN2=log((double)N2); tol=1E-12; cterm=REAL(consttermR); //Rprintf("m: %d\nN2: %d\ndf: %d\n",m,N2,df);

  // Initialize Tinv/b.free/vlog variables
  SEXP Ti, bfi; Matrix Tinv; ColumnVector bfree, vlog(N2), qlog;

  // Initialize Wlist/W/Wmat objects and populate Wlist from WpostR
  Wlist Wall(WpostR,N2); Wobj W; Matrix Wmat;

  // Initialize SEXP/ptr to store/access log marginal A0k values
  SEXP lpa0yao; PROTECT(lpa0yao=allocVector(REALSXP,m)); lpa0=REAL(lpa0yao);
  for(i=0;i<m-1;i++){
    PROTECT(Ti=VECTOR_ELT(TinvR,i));
    dTi=getdims(Ti); Tinv=R2Cmat(Ti,dTi[0],dTi[1]);
    UNPROTECT(1); //Rprintf("Tinv[[%d]](%dx%d) initialized\n",i,dTi[0],dTi[1]);

    PROTECT(bfi=VECTOR_ELT(bfR,i));
    db=length(bfi); bfree.ReSize(db); pbfi=REAL(bfi); bfree<<pbfi;
    UNPROTECT(1); //Rprintf("bfree[[%d]](%d) initialized\n",i,db);

    for(j=1;j<=N2;j++){
      Wall.getWobj(W,j); Wmat=W.getWelt(i+1); W.clear();
      vlog(j)=getvlog(Wmat,Tinv,bfree,cterm[i],df,tol); //Rprintf("vlog(%d): %f\n",j,vlog(j));
    }

    // Modified harmonic mean of the max
    maxvlog=vlog.Maximum(); qlog=vlog-maxvlog; len=qlog.Storage();
    lqlog=0; for(j=1;j<=len;j++) lqlog+=exp(qlog(j)); lqlog=log(lqlog); // log(sum(exp(qlog)))
    lpa0[i]=maxvlog-lN2+lqlog; //Rprintf("lpa0[%d] = %f\n", i, lpa0[i]);
  }

  // Computations for last column
  PROTECT(Ti=VECTOR_ELT(TinvR,m-1));
  dTi=getdims(Ti); Tinv=R2Cmat(Ti,dTi[0],dTi[1]);
  UNPROTECT(1); //Rprintf("Tinv[[%d]](%dx%d) initialized\n",i,dTi[0],dTi[1]);

  PROTECT(bfi=VECTOR_ELT(bfR,m-1));
  pbfi=REAL(bfi); bfree.ReSize(length(bfi)); bfree<<pbfi;
  UNPROTECT(1); //Rprintf("bfree[[%d]](%d) initialized\n",i,db);

  UTobj UT(UTR); Matrix A0=R2Cmat(A0R,m,m);
  A0=drawA0cpp(A0,UT,df,INTEGER(n0R),W); Wmat=W.getWelt(m);
  lpa0[m-1]=getvlog(Wmat,Tinv,bfree,cterm[m-1],df,tol); //Rprintf("lpa0[%d] = %f\n",m-1,lpa0[m-1]);

  // Return R object lpa0yao
  UNPROTECT(1); return lpa0yao;
}

SEXP icav_Matrix_sort(SEXP _M,SEXP Param)
{
    // sorts a cloned version
    //copies values in (Col_index and Row_index) of _M to a new matrix mimat. size of mimat < size of _M

    // param[0] rows of _M
    // param[1] cols of _M
    // param[2] chunk size
    SEXP mimat;
    int i,j,_COL;
    int MAX_T,Tnum;
    int one = 1;

    double *M,*mm;
    int *param;

    PROTECT(_M = coerceVector(_M, REALSXP));
    PROTECT(Param = coerceVector(Param, INTSXP));

    M = REAL(_M);
    param = INTEGER(Param);
    PROTECT(mimat = allocMatrix(REALSXP,param[0],param[1]));
    mm=REAL(mimat);
    MAX_T = omp_get_max_threads();

    //we order the cols
    for(i=0; i<param[1]; i++) {
        _COL = i*param[0];
        for(j=0; j<param[0]; j++) {
            mm[j+_COL]=M[j+_COL];
        }
    }
    for(i=0; i<param[1]; i++) {
        _COL = i*param[0];
        qsort (&mm[_COL], param[0], sizeof(double), icav_compare_doubles);
    }

    UNPROTECT(3);
    return(mimat);
}

SEXP binomial_dev_resids(SEXP y, SEXP mu, SEXP wt)
{
    int i, n = LENGTH(y), lmu = LENGTH(mu), lwt = LENGTH(wt), nprot = 1;
    SEXP ans;
    double mui, yi, *rmu, *ry, *rwt, *rans;

    if (!isReal(y)) {y = PROTECT(coerceVector(y, REALSXP)); nprot++;}
    ry = REAL(y);
    ans = PROTECT(duplicate(y));
    rans = REAL(ans);
    if (!isReal(mu)) {mu = PROTECT(coerceVector(mu, REALSXP)); nprot++;}
    if (!isReal(wt)) {wt = PROTECT(coerceVector(wt, REALSXP)); nprot++;}
    rmu = REAL(mu);
    rwt = REAL(wt);
    if (lmu != n && lmu != 1)
	error(_("argument %s must be a numeric vector of length 1 or length %d"),
	      "mu", n);
    if (lwt != n && lwt != 1)
	error(_("argument %s must be a numeric vector of length 1 or length %d"),
	      "wt", n);
    /* Written separately to avoid an optimization bug on Solaris cc */
    if(lmu > 1) {
	for (i = 0; i < n; i++) {
	    mui = rmu[i];
	    yi = ry[i];
	    rans[i] = 2 * rwt[lwt > 1 ? i : 0] *
		(y_log_y(yi, mui) + y_log_y(1 - yi, 1 - mui));
	}
    } else {
	mui = rmu[0];
	for (i = 0; i < n; i++) {
	    yi = ry[i];
	    rans[i] = 2 * rwt[lwt > 1 ? i : 0] *
		(y_log_y(yi, mui) + y_log_y(1 - yi, 1 - mui));
	}
    }

    UNPROTECT(nprot);
    return ans;
}

/* all, any */
SEXP attribute_hidden do_logic3(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP ans, s, t, call2;
    int narm, has_na = 0;
    /* initialize for behavior on empty vector
       all(logical(0)) -> TRUE
       any(logical(0)) -> FALSE
     */
    Rboolean val = PRIMVAL(op) == _OP_ALL ? TRUE : FALSE;

    PROTECT(args = fixup_NaRm(args));
    PROTECT(call2 = duplicate(call));
    SETCDR(call2, args);

    if (DispatchGroup("Summary", call2, op, args, env, &ans)) {
	UNPROTECT(2);
	return(ans);
    }

    ans = matchArgExact(R_NaRmSymbol, &args);
    narm = asLogical(ans);

    for (s = args; s != R_NilValue; s = CDR(s)) {
	t = CAR(s);
	/* Avoid memory waste from coercing empty inputs, and also
	   avoid warnings with empty lists coming from sapply */
	if(xlength(t) == 0) continue;
	/* coerceVector protects its argument so this actually works
	   just fine */
	if (TYPEOF(t) != LGLSXP) {
	    /* Coercion of integers seems reasonably safe, but for
	       other types it is more often than not an error.
	       One exception is perhaps the result of lapply, but
	       then sapply was often what was intended. */
	    if(TYPEOF(t) != INTSXP)
		warningcall(call,
			    _("coercing argument of type '%s' to logical"),
			    type2char(TYPEOF(t)));
	    t = coerceVector(t, LGLSXP);
	}
	val = checkValues(PRIMVAL(op), narm, LOGICAL(t), XLENGTH(t));
        if (val != NA_LOGICAL) {
            if ((PRIMVAL(op) == _OP_ANY && val)
                || (PRIMVAL(op) == _OP_ALL && !val)) {
                has_na = 0;
                break;
            }
        } else has_na = 1;
    }
    UNPROTECT(2);
    return has_na ? ScalarLogical(NA_LOGICAL) : ScalarLogical(val);
}

SEXP convolve2(SEXP a, SEXP b) {
    int na, nb, nab;
    double *xa, *xb, *xab;
    SEXP ab;

    a = PROTECT(coerceVector(a, REALSXP));
    b = PROTECT(coerceVector(b, REALSXP));
    na = length(a);
    nb = length(b);
    nab = na + nb - 1;
    ab = PROTECT(allocVector(REALSXP, nab));
    xa = REAL(a);
    xb = REAL(b);
    xab = REAL(ab);
    for(int i = 0; i < nab; i++)
        xab[i] = 0.0;
    for(int i = 0; i < na; i++)
        for(int j = 0; j < nb; j++)
            xab[i + j] += xa[i] * xb[j];
    UNPROTECT(3);
    return ab;
}

SEXP createPopulation(SEXP popSize_ext, SEXP funcSet, SEXP inSet, SEXP maxDepth_ext, SEXP constProb_ext, SEXP subtreeProb_ext, SEXP constScaling_ext) {
  SEXP pop;
  PROTECT(popSize_ext = coerceVector(popSize_ext, INTSXP));
  int popSize= INTEGER(popSize_ext)[0];

  PROTECT(pop= allocVector(VECSXP, popSize));
  for(int i=0; i < popSize; i++)
  {
    SET_VECTOR_ELT(pop, i, randFuncGrow(funcSet, inSet, maxDepth_ext, constProb_ext, subtreeProb_ext, constScaling_ext));
  }
  UNPROTECT(2);
  return pop;
}

double getScalarReal(SEXP foo)
{
    if (! isNumeric(foo))
        error("argument must be numeric");
    if (LENGTH(foo) != 1)
        error("argument must be scalar");
    if (isReal(foo)) {
        return REAL(foo)[0];
    } else {
        SEXP bar = coerceVector(foo, REALSXP);
        return REAL(bar)[0];
    }
}

SEXP pacf1(SEXP acf, SEXP lmax)
{
    int lagmax = asInteger(lmax);
    acf = PROTECT(coerceVector(acf, REALSXP));
    SEXP ans = PROTECT(allocVector(REALSXP, lagmax));
    uni_pacf(REAL(acf), REAL(ans), lagmax);
    SEXP d = PROTECT(allocVector(INTSXP, 3));
    INTEGER(d)[0] = lagmax;
    INTEGER(d)[1] = INTEGER(d)[2] = 1;
    setAttrib(ans, R_DimSymbol, d);
    UNPROTECT(3);
    return ans;
}

// Returns the proper integer chrom target (as is or through factor levels)
// Negative values should skip further treatments (level not found in factor)
int chromTarget(SEXP chrom, SEXP targetChrom) {
	int output;
	if(isFactor(chrom)) {
		// Convert 'targetChrom' to a character vector
		if(isFactor(targetChrom)) { targetChrom = PROTECT(asCharacterFactor(targetChrom));
		} else                    { targetChrom = PROTECT(coerceVector(targetChrom, STRSXP));
		}
		
		if(LENGTH(targetChrom) != 1 || STRING_ELT(targetChrom, 0) == NA_STRING) {
			error("As 'chrom' is factor, target 'chrom' must be a single non-NA character-coercible value");
		}
		
		// From character to integer position in the index (0+)
		SEXP levels = PROTECT(getAttrib(chrom, R_LevelsSymbol));
		for(int i = 0; i < LENGTH(levels); i++) {
			if(strcmp(CHAR(STRING_ELT(levels, i)), CHAR(STRING_ELT(targetChrom, 0))) == 0) {
				// Early exit when found
				output = i;
				UNPROTECT(2);
				return output;
			}
		}
		
		// Was not found in levels
		output = -1;
		UNPROTECT(2);
	} else {
		// Integer chrom
		targetChrom = PROTECT(coerceVector(targetChrom, INTSXP));
		if(LENGTH(targetChrom) != 1 || INTEGER(targetChrom)[0] == NA_INTEGER || INTEGER(targetChrom)[0] < 0) {
			error("As 'chrom' is integer, target 'chrom' must be a single non-NA integer-coercible strictly positive value");
		}
		
		// Position in the index (0+)
		output = INTEGER(targetChrom)[0] - 1;
		UNPROTECT(1);
	}
	return output;
}

SEXP lsq_dense_QR(SEXP X, SEXP y)
{
    SEXP ans;
    int info, n, p, k, *Xdims, *ydims, lwork;
    double *work, tmp, *xvals;

    if (!(isReal(X) & isMatrix(X)))
	error(_("X must be a numeric (double precision) matrix"));
    Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
    n = Xdims[0];
    p = Xdims[1];
    if (!(isReal(y) & isMatrix(y)))
	error(_("y must be a numeric (double precision) matrix"));
    ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));
    if (ydims[0] != n)
	error(_(
	    "number of rows in y (%d) does not match number of rows in X (%d)"),
	    ydims[0], n);
    k = ydims[1];
    if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);
    xvals = (double *) R_alloc(n * p, sizeof(double));
    Memcpy(xvals, REAL(X), n * p);
    ans = PROTECT(duplicate(y));
    lwork = -1;
    F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
		    &tmp, &lwork, &info);
    if (info)
	error(_("First call to Lapack routine dgels returned error code %d"),
	      info);
    lwork = (int) tmp;
    work = (double *) R_alloc(lwork, sizeof(double));
    F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
		    work, &lwork, &info);
    if (info)
	error(_("Second call to Lapack routine dgels returned error code %d"),
	      info);
    UNPROTECT(1);
    return ans;
}

SEXP SWilk(SEXP x)
{
    int n, ifault = 0;
    double W = 0, pw;  /* original version tested W on entry */
    x = PROTECT(coerceVector(x, REALSXP));
    n = LENGTH(x);
    swilk(REAL(x), n, &W, &pw, &ifault);
    if (ifault > 0 && ifault != 7)
	error("ifault=%d. This should not happen", ifault);
    SEXP ans = PROTECT(allocVector(REALSXP, 2));
    REAL(ans)[0] = W, REAL(ans)[1] = pw;
    UNPROTECT(2);
    return ans;
}

SEXP acf(SEXP x, SEXP lmax, SEXP sCor)
{
    int nx = nrows(x), ns = ncols(x), lagmax = asInteger(lmax),
        cor = asLogical(sCor);
    x = PROTECT(coerceVector(x, REALSXP));
    SEXP ans = PROTECT(allocVector(REALSXP, (lagmax + 1)*ns*ns));
    acf0(REAL(x), nx, ns, lagmax, cor, REAL(ans));
    SEXP d = PROTECT(allocVector(INTSXP, 3));
    INTEGER(d)[0] = lagmax + 1;
    INTEGER(d)[1] = INTEGER(d)[2] = ns;
    setAttrib(ans, R_DimSymbol, d);
    UNPROTECT(3);
    return ans;
}

size_t dplRlength(SEXP x) {
    size_t xlength;
    SEXP sn, tmp, ncall;
    PROTECT_INDEX ipx;
    PROTECT(tmp = ncall = allocList(2));
    SET_TYPEOF(ncall, LANGSXP);
    SETCAR(tmp, install("length")); tmp = CDR(tmp);
    SETCAR(tmp, x);
    PROTECT_WITH_INDEX(sn = eval(ncall, R_BaseEnv), &ipx);
    REPROTECT(sn = coerceVector(sn, REALSXP), ipx);
    xlength = (size_t) *REAL(sn);
    UNPROTECT(2);
    return xlength;
}

/*
 * Calculate piece-wise linear partial value function
 * @param x: vector of measurements (length N)
 * @param y: vector of cut-offs (length n) defining the PVF
 * @param v: vector of values corresponding to the cut-offs (length n)
 * @return vector of values (length N)
 */
SEXP smaa_pvf(SEXP _x, SEXP _y, SEXP _v) {
	int const N = length(_x);
	int const n = length(_y);
	
	_x = PROTECT(coerceVector(_x, REALSXP));
	_y = PROTECT(coerceVector(_y, REALSXP));
	_v = PROTECT(coerceVector(_v, REALSXP));
	double const *x = REAL(_x);
	double const *y = REAL(_y);
	double const *v = REAL(_v);

	SEXP _out = PROTECT(allocVector(REALSXP, N));
	double *out = REAL(_out);

	for (unsigned i = 0; i < N; ++i) {
		unsigned j;
		for (j = 1; j < (n - 1) && y[j] < x[i]; ++j) ;
		out[i] = v[j - 1] + (x[i] - y[j - 1]) * ((v[j] - v[j - 1]) / (y[j] - y[j - 1]));
	}

	UNPROTECT(4);
	return _out;
}

/**Cross Mean. Function takes in nxd matrix x and generates a dxd matrix with means of each column in x*/
SEXP crossMean(SEXP x){	
	int n, d;
	int *xdims;
	double sum;
	double *meanptr, *xptr;
	SEXP mean; 
	
	xdims = getDims(x);
	n = xdims[0];
	d = xdims[1];
	
	PROTECT(x = coerceVector(x, REALSXP)); //PROTECT 1
	xptr = REAL(x);	
	PROTECT(mean = allocMatrix(REALSXP,1,d)); //PROTECT 2
	meanptr = REAL(mean);
	
	int k =0; //index for mean	
	for(int i=0;i<=(d-1)*n;i=i+n){
		sum = 0;		
		for(int j=i;j<=i+n-1;j++){			
			sum = sum+xptr[j];			
			}//end inner for		
		meanptr[k] = sum/n;	//calculate mean and store value
		k++;	
	}//end outer for
			
	/**Mean vector completed. Now generate crossMean vector.*/
	
	SEXP crossM;
	double *crossMptr;
	PROTECT(crossM = allocMatrix(REALSXP,d,d)); //PROTECT 3
	crossMptr = REAL(crossM);
	double prod;
	
	int c = 0; //index for crossMean	
	/**Calculates the mean. First fills the diagonal, then fills the upper and lower triangular portions using symmetry.*/
	for(int i=0; i<d;i++){		
		c = i*(d+1);
		crossMptr[c] = meanptr[i] * meanptr[i];		
		for(int j=i+1;j<d;j++){			
			prod = meanptr[i] * meanptr[j];					
			c = i + j*d;			
			crossMptr[c] = prod;			
			c = j + i*d;
			crossMptr[c] = prod; 
		}//end inner for
	}//end outer for		
	UNPROTECT(3);
	return(crossM);		
}//end crossMean

// used in gap.inter.R
SEXP icav_wss_matrix(SEXP _M, SEXP Param)
{
    //copies values in (Col_index and Row_index) of _M to a new matrix mimat. size of mimat < size of _M

    // param[0] rows of _M
    // param[1] cols of _M
//param[0] == param[1] this matrix must be symmetric
    SEXP mimat;
    int i,j,_COL;
    int MAX_T,Tnum;

    double *M,*mm,tmp_max;
    double tmp_acc;
    int *param;

    PROTECT(_M = coerceVector(_M, REALSXP));
    PROTECT(Param = coerceVector(Param, INTSXP));

    M = REAL(_M);
    param = INTEGER(Param);

    PROTECT(mimat = allocVector(REALSXP, 1));

    mm = REAL(mimat);
    tmp_acc = 0.;
    //we order the cols
    for(i=0; i<param[1]; i++) {
        _COL = i*param[0];
        for(j=0; j<param[0]; j++) {
            tmp_acc += M[_COL+j]*M[_COL+j];
        }
    }

    *mm = tmp_acc/(2*param[0]);
    UNPROTECT(3);
    return(mimat);
}

SEXP R_colwiseProd(SEXP V, SEXP X)
{

  int    *xdims, nrx, ncx, ii, jj, kk, len_V;
  double *xptr, *vptr;
  xdims = getDims(X);
  nrx = xdims[0];
  ncx = xdims[1];
  len_V = length(V);

  PROTECT(X = coerceVector(X, REALSXP));
  xptr  = REAL(X);

  PROTECT(V = coerceVector(V, REALSXP));
  vptr = REAL(V);

  double *ansptr;
  SEXP ans;
  PROTECT(ans = allocMatrix(REALSXP, nrx, ncx));
  ansptr = REAL(ans);
  
  kk = 0;
  for (jj = 0; jj < ncx; jj++){
    //Rprintf("kk=%i len_V=%i\n", kk, len_V);
    for (ii=0; ii<nrx; ii++){
      ansptr[ii + nrx * jj] = vptr[kk] * xptr[ii + nrx * jj];
    }
    kk++;
    if (kk == len_V){
      //Rprintf("HERE: kk=%i len_V=%i\n", kk, len_V);
      kk=0;
    }
  }
  
  UNPROTECT(3);
  return(ans);
}

SEXP Expect_matrix(SEXP S1, SEXP S, SEXP lambda, SEXP time, SEXP theta0, SEXP theta1, SEXP matdiag){
    
    int nvar, npoints, nprotect=0;
    
    nvar = GET_LENGTH(lambda);
    npoints = GET_LENGTH(time);
    
    PROTECT(time = coerceVector(time,REALSXP)); nprotect++;
    PROTECT(theta0 = coerceVector(theta0,REALSXP)); nprotect++;
    PROTECT(theta1 = coerceVector(theta1,REALSXP)); nprotect++;
    // results
    SEXP expectation = PROTECT(allocVector(REALSXP,nvar*npoints)); nprotect++;
    
    if(!isComplex(lambda)){
    // eigenvectors
    PROTECT(S1 = coerceVector(S1,REALSXP)); nprotect++;
    PROTECT(S = coerceVector(S,REALSXP)); nprotect++;
    // matrix exponential
    SEXP matexp = PROTECT(allocVector(REALSXP,nvar*nvar*npoints)); nprotect++;

    // Compute the exponential matrix
    multi_exp_matrix (&nvar, &npoints, REAL(time), REAL(lambda), REAL(S), REAL(S1), REAL(matexp));
    
    // Compute the expectations
    optimum (&nvar, &npoints, REAL(time), REAL(theta0), REAL(theta1), REAL(expectation), REAL(matexp), REAL(matdiag));
    // Done.
        
    }else{
    
    double complex *matexp;
    // complex eigenvalues & eigenvectors
    PROTECT(S1 = coerceVector(S1,CPLXSXP)); nprotect++;
    PROTECT(S = coerceVector(S,CPLXSXP)); nprotect++;
        
    // alloc a complex vector in C rather than R structure...
    matexp = Calloc(nvar*nvar*npoints,double complex);
        
    // Compute the exponential matrix
    multi_exp_matrix_complex (&nvar, &npoints, REAL(time), COMPLEX(lambda), COMPLEX(S), COMPLEX(S1), matexp);
        
    // Compute the expectations
    optimum_complex(&nvar, &npoints, REAL(time), REAL(theta0), REAL(theta1), REAL(expectation), matexp, REAL(matdiag));
    // Done.
    // Free the memory
    Free(matexp);
    }


    UNPROTECT(nprotect);
    return expectation;
    
}

SEXP do_backtrack(SEXP n, SEXP rs, SEXP cs)
{
    int i, fill, nr = length(rs), nc = length(cs), nmat = asInteger(n);
    int N = nr * nc;

    /* check & cast */
    if(TYPEOF(rs) != INTSXP)
	rs = coerceVector(rs, INTSXP);
    PROTECT(rs);
    if(TYPEOF(cs) != INTSXP)
	cs = coerceVector(cs, INTSXP);
    PROTECT(cs);
    int *rowsum = INTEGER(rs);
    int *colsum = INTEGER(cs);

    /* initialize work arrays for backtrack()*/
    int *ind = (int *) R_alloc(nr * nc, sizeof(int));
    int *rfill = (int *) R_alloc(nr, sizeof(int));
    int * cfill = (int *) R_alloc(nc, sizeof(int));
    for (i = 0, fill = 0; i < nr; i++)
	fill += rowsum[i];
    int *x = (int *) R_alloc(nr * nc, sizeof(int));

    SEXP out =  PROTECT(alloc3DArray(INTSXP, nr, nc, nmat));
    int *iout = INTEGER(out);

    GetRNGstate();
    /* Call static C function */
    for(i = 0; i < nmat; i++) {
	backtrack(x, rowsum, colsum, fill, nr, nc, rfill, cfill, ind);
	memcpy(iout + i * N, x, N * sizeof(int));
    }
    PutRNGstate();

    UNPROTECT(3);
    return out;
}

SEXP clamp(SEXP d, SEXP rr, SEXP usevals) {
					
	R_len_t i;
	SEXP val;
	double *xd, *r, *xval;
	PROTECT(d = coerceVector(d, REALSXP));
	PROTECT(rr = coerceVector(rr, REALSXP));
  	int uv = INTEGER(usevals)[0];
	r = REAL(rr);
	xd = REAL(d);

	int n = length(d);
	PROTECT( val = allocVector(REALSXP, n) );
	xval = REAL(val);
	
	if (uv) {
		for (i=0; i<n; i++) {
			if ( xd[i] < r[0] ) { 
				xval[i] = r[0];
			} else if ( xd[i] > r[1] ) { 
				xval[i] = r[1];
			} else {
				xval[i] = xd[i];
			}
		}		
	} else {
		for (i=0; i<n; i++) {
			if ( (xd[i] < r[0]) | (xd[i] > r[1])) {
				xval[i] = R_NaReal;
			} else {
				xval[i] = xd[i];
			}
		}	
	}
	UNPROTECT(3);
	return(val);
}