add more GNU R sources from stats package

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/HoltWinters.c

+/*  R : A Computer Language for Statistical Data Analysis
+ *
+ *  Copyright (C) 2003-7  The R Core Team
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, a copy is available at
+ *  https://www.R-project.org/Licenses/.
+ */
+
+/* Originally contributed by David Meyer */
+
+#include <stdlib.h>
+#include <string.h>  // memcpy
+
+#include <R.h>
+#include "ts.h"
+
+void HoltWinters (double *x,
+		  int    *xl,
+		  double *alpha,
+		  double *beta,
+		  double *gamma,
+		  int    *start_time,
+		  int    *seasonal,
+		  int    *period,
+		  int    *dotrend,
+		  int    *doseasonal,
+
+		  double *a,
+		  double *b,
+		  double *s,
+
+		  /* return values */
+		  double *SSE,
+		  double *level,
+		  double *trend,
+		  double *season
+    )
+
+{
+    double res = 0, xhat = 0, stmp = 0;
+    int i, i0, s0;
+
+    /* copy start values to the beginning of the vectors */
+    level[0] = *a;
+    if (*dotrend == 1) trend[0] = *b;
+    if (*doseasonal == 1) memcpy(season, s, *period * sizeof(double));
+
+    for (i = *start_time - 1; i < *xl; i++) {
+	/* indices for period i */
+	i0 = i - *start_time + 2;
+	s0 = i0 + *period - 1;
+
+	/* forecast *for* period i */
+	xhat = level[i0 - 1] + (*dotrend == 1 ? trend[i0 - 1] : 0);
+	stmp = *doseasonal == 1 ? season[s0 - *period] : (*seasonal != 1);
+	if (*seasonal == 1)
+	    xhat += stmp;
+	else
+	    xhat *= stmp;
+
+	/* Sum of Squared Errors */
+	res   = x[i] - xhat;
+	*SSE += res * res;
+
+	/* estimate of level *in* period i */
+	if (*seasonal == 1)
+	    level[i0] = *alpha       * (x[i] - stmp)
+		      + (1 - *alpha) * (level[i0 - 1] + trend[i0 - 1]);
+	else
+	    level[i0] = *alpha       * (x[i] / stmp)
+		      + (1 - *alpha) * (level[i0 - 1] + trend[i0 - 1]);
+
+	/* estimate of trend *in* period i */
+	if (*dotrend == 1)
+	    trend[i0] = *beta        * (level[i0] - level[i0 - 1])
+		      + (1 - *beta)  * trend[i0 - 1];
+
+	/* estimate of seasonal component *in* period i */
+	if (*doseasonal == 1) {
+	    if (*seasonal == 1)
+		season[s0] = *gamma       * (x[i] - level[i0])
+			   + (1 - *gamma) * stmp;
+	    else
+		season[s0] = *gamma       * (x[i] / level[i0])
+			   + (1 - *gamma) * stmp;
+	}
+    }
+}

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/Srunmed.c

+/*
+ *  R : A Computer Language for Statistical Data Analysis
+ *  Copyright (C) 1995--2002 Martin Maechler <maechler@stat.math.ethz.ch>
+ *  Copyright (C) 2003       The R Foundation
+ *  Copyright (C) 2012-2016  The R Core Team
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, a copy is available at
+ *  https://www.R-project.org/Licenses/
+ */
+
+#include "modreg.h"
+
+#include "Trunmed.c"
+
+static void Srunmed(double* y, double* smo, R_xlen_t n, int bw,
+		    int end_rule, int debug)
+{
+/*
+ *  Computes "Running Median" smoother with medians of 'band'
+
+ *  Input:
+ *	y(n)	- responses in order of increasing predictor values
+ *	n	- number of observations
+ *	bw	- span of running medians (MUST be ODD !!)
+ *	end_rule -- 0: Keep original data at ends {j; j < b2 | j > n-b2}
+ *		 -- 1: Constant ends = median(y[1],..,y[bw]) "robust"
+ *  Output:
+ *	smo(n)	- smoothed responses
+
+ * NOTE:  The 'end' values are just copied !! this is fast but not too nice !
+ */
+
+/* Local variables */
+    double rmed, rmin, temp, rnew, yout, yi;
+    double rbe, rtb, rse, yin, rts;
+    int imin, ismo, i, j, first, last, band2, kminus, kplus;
+
+
+    double *scrat = (double *) R_alloc(bw, sizeof(double));
+    /*was  malloc( (unsigned) bw * sizeof(double));*/
+
+    if(bw > n)
+	error(_("bandwidth/span of running medians is larger than n"));
+
+/* 1. Compute  'rmed' := Median of the first 'band' values
+   ======================================================== */
+
+    for (int i = 0; i < bw; ++i)
+	scrat[i] = y[i];
+
+    /* find minimal value  rmin = scrat[imin] <= scrat[j] */
+    rmin = scrat[0]; imin = 0;
+    for (int i = 1; i < bw; ++i)
+	if (scrat[i] < rmin) {
+	    rmin = scrat[i]; imin = i;
+	}
+    /* swap scrat[0] <-> scrat[imin] */
+    temp = scrat[0]; scrat[0] = rmin; scrat[imin] = temp;
+    /* sort the rest of of scrat[] by bubble (?) sort -- */
+    for (int i = 2; i < bw; ++i) {
+	if (scrat[i] < scrat[i - 1]) {/* find the proper place for scrat[i] */
+	    temp = scrat[i];
+	    j = i;
+	    do {
+		scrat[j] = scrat[j - 1];
+		--j;
+	    } while (scrat[j - 1] > temp); /* now:  scrat[j-1] <= temp */
+	    scrat[j] = temp;
+	}
+    }
+
+    band2 = bw / 2;
+    rmed = scrat[band2];/* == Median( y[(1:band2)-1] ) */
+    /* "malloc" had  free( (char*) scrat);*/ /*-- release scratch memory --*/
+
+    if(end_rule == 0) { /*-- keep DATA at end values */
+	for (i = 0; i < band2; ++i)
+	    smo[i] = y[i];
+    }
+    else { /* if(end_rule == 1)  copy median to CONSTANT end values */
+	for (i = 0; i < band2; ++i)
+	    smo[i] = rmed;
+    }
+    smo[band2] = rmed;
+    band2++; /* = bw / 2 + 1*/;
+
+    if(debug) REprintf("(bw,b2)= (%d,%d)\n", bw, band2);
+
+/* Big	FOR Loop: RUNNING median, update the median 'rmed'
+   ======================================================= */
+    for (first = 1, last = bw, ismo = band2;
+	 last < n;
+	 ++first, ++last, ++ismo) {
+
+	yin = y[last];
+	yout = y[first - 1];
+
+	if(debug) REprintf(" is=%d, y(in/out)= %10g, %10g", ismo ,yin, yout);
+
+	rnew = rmed; /* New median = old one   in all the simple cases --*/
+
+	if (yin < rmed) {
+	    if (yout >= rmed) {
+		kminus = 0;
+		if (yout > rmed) {/*	--- yin < rmed < yout --- */
+		    if(debug) REprintf(": yin < rmed < yout ");
+		    rnew = yin;/* was -rinf */
+		    for (i = first; i <= last; ++i) {
+			yi = y[i];
+			if (yi < rmed) {
+			    ++kminus;
+			    if (yi > rnew)	rnew = yi;
+			}
+		    }
+		    if (kminus < band2)		rnew = rmed;
+		}
+		else {/*		--- yin < rmed = yout --- */
+		    if(debug) REprintf(": yin < rmed == yout ");
+		    rse = rts = yin;/* was -rinf */
+		    for (i = first; i <= last; ++i) {
+			yi = y[i];
+			if (yi <= rmed) {
+			    if (yi < rmed) {
+				++kminus;
+				if (yi > rts)	rts = yi;
+				if (yi > rse)	rse = yi;
+			    } else		rse = yi;
+
+			}
+		    }
+		    rnew = (kminus == band2) ? rts : rse ;
+		    if(debug) REprintf("k- : %d,", kminus);
+		}
+	    } /* else: both  yin, yout < rmed -- nothing to do .... */
+	}
+	else if (yin != rmed) { /* yin > rmed */
+	    if (yout <= rmed) {
+		kplus = 0;
+		if (yout < rmed) {/* -- yout < rmed < yin --- */
+		    if(debug) REprintf(": yout < rmed < yin ");
+		    rnew = yin; /* was rinf */
+		    for (i = first; i <= last; ++i) {
+			yi = y[i];
+			if (yi > rmed) {
+			    ++kplus;
+			    if (yi < rnew)	rnew = yi;
+			}
+		    }
+		    if (kplus < band2)	rnew = rmed;
+
+		} else { /* -- yout = rmed < yin --- */
+		    if(debug) REprintf(": yout == rmed < yin ");
+		    rbe = rtb = yin; /* was rinf */
+		    for (i = first; i <= last; ++i) {
+			yi = y[i];
+			if (yi >= rmed) {
+			    if (yi > rmed) {
+				++kplus;
+				if (yi < rtb)	rtb = yi;
+				if (yi < rbe)	rbe = yi;
+			    } else		rbe = yi;
+			}
+		    }
+		    rnew = (kplus == band2) ? rtb : rbe;
+		    if(debug) REprintf("k+ : %d,", kplus);
+		}
+	    } /* else: both  yin, yout > rmed --> nothing to do */
+	} /* else: yin == rmed -- nothing to do .... */
+	if(debug) REprintf("=> %12g, %12g\n", rmed, rnew);
+	rmed = rnew;
+	smo[ismo] = rmed;
+    } /*     END FOR ------------ big Loop -------------------- */
+
+    if(end_rule == 0) { /*-- keep DATA at end values */
+	for (i = ismo; i < n; ++i)
+	    smo[i] = y[i];
+    }
+    else { /* if(end_rule == 1)  copy median to CONSTANT end values */
+	for (i = ismo; i < n; ++i)
+	    smo[i] = rmed;
+    }
+} /* Srunmed */
+
+SEXP runmed(SEXP x, SEXP stype, SEXP sk, SEXP end, SEXP print_level)
+{
+    if (TYPEOF(x) != REALSXP) error("numeric 'x' required");
+    R_xlen_t n = XLENGTH(x);
+    int type = asInteger(stype), k = asInteger(sk),
+	iend = asInteger(end), pl = asInteger(print_level);
+    SEXP ans = PROTECT(allocVector(REALSXP, n));
+    if (type == 1) {
+	if (IS_LONG_VEC(x)) 
+	    error("long vectors are not supported for algorithm = \"Turlach\"");
+	int *i1 = (int *) R_alloc(k + 1, sizeof(int)), 
+	    *i2 = (int *) R_alloc(2*k + 1, sizeof(int));
+	double *d1 = (double *) R_alloc(2*k + 1, sizeof(double));
+	Trunmed(n, k, REAL(x), REAL(ans), i1, i2, d1, iend, pl);
+    } else {
+	Srunmed(REAL(x), REAL(ans), n, k, iend, pl > 0);
+    }
+    UNPROTECT(1);
+    return ans;
+}

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/Trunmed.c

+/* Copyright (C) 1995   Berwin A. Turlach <berwin@alphasun.anu.edu.au>
+ * Copyright (C) 2000-2 Martin Maechler <maechler@stat.math.ethz.ch>
+ * Copyright (C) 2003   The R Foundation
+ * Copyright (C) 2012-2016   The R Core Team
+
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the
+ * GNU General Public License for more details.
+
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, a copy is available at
+ * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ */
+
+/* These routines implement a running median smoother according to the
+ * algorithm described in Haerdle und Steiger (1995).
+ *
+ * A tech-report of that paper is available under
+ * ftp://amadeus.wiwi.hu-berlin.de/pub/papers/sfb/sfb1994/dpsfb940015.ps.Z
+ *
+ * The current implementation does not use any global variables!
+ */
+
+/* Changes for R port by Martin Maechler ((C) above):
+ *
+ *  s/long/int/			R uses int, not long (as S does)
+ *  s/void/static void/		most routines are internal
+ *
+ * Added  print_level  and end_rule  arguments
+ */
+
+/* Variable	name	descri- | Identities from paper
+ * name here	paper	ption   | (1-indexing)
+ * ---------    -----	-----------------------------------
+ * window[]      H      the array containing the double heap
+ * data[]        X      the data (left intact)
+ * ...		 i	1st permuter  H[i[m]]    == X[i + m]
+ * ...		 j	2nd permuter  X[i +j[m]] == H[m]
+ */
+
+#include <math.h>
+
+static void
+swap(int l, int r, double *window, int *outlist, int *nrlist, int print_level)
+{
+    /* swap positions `l' and `r' in window[] and nrlist[]
+     *
+     * ---- Used in R_heapsort() and many other routines
+     */
+    int nl, nr;
+    double tmp;
+
+    if(print_level >= 3) Rprintf("SW(%d,%d) ", l,r);
+    tmp = window[l]; window[l]	    = window[r];  window[r] = tmp;
+    nl	= nrlist[l]; nrlist[l] = (nr= nrlist[r]); nrlist[r] = nl;
+
+    outlist[nl/* = nrlist[r] */] = r;
+    outlist[nr/* = nrlist[l] */] = l;
+}
+
+static void
+siftup(int l, int r, double *window, int *outlist, int *nrlist, int print_level)
+{
+/* Used only in R_heapsort() */
+    int i, j, nrold;
+    double x;
+
+    if(print_level >= 2) Rprintf("siftup(%d,%d) ", l,r);
+    i	  = l;
+    j	  = 2*i;
+    x	  = window[i];
+    nrold = nrlist[i];
+    while (j <= r) {
+	if (j < r)
+	    if (window[j] < window[j+1])
+		j++;
+	if (x >= window[j])
+	    break;
+
+	window[i]	   = window[j];
+	outlist[nrlist[j]] = i;
+	nrlist[i]	   = nrlist[j];
+	i = j;
+	j = 2*i;
+    }
+    window[i]	   = x;
+    outlist[nrold] = i;
+    nrlist[i]	   = nrold;
+}
+
+static void
+R_heapsort(int low, int up, double *window, int *outlist, int *nrlist,
+	 int print_level)
+{
+    int l, u;
+
+    l = (up/2) + 1;
+    u = up;
+    while(l > low) {
+	l--;
+	siftup(l, u, window, outlist, nrlist, print_level);
+    }
+    while(u > low) {
+	swap(l, u, window, outlist, nrlist, print_level);
+	u--;
+	siftup(l, u, window, outlist, nrlist, print_level);
+    }
+}
+
+static void
+inittree(R_xlen_t n, int k, int k2, const double *data, double *window,
+	 int *outlist, int *nrlist, int print_level)
+{
+    int i, k2p1;
+    double big;
+
+    for(i=1; i <= k; i++) { /* use 1-indexing for our three arrays !*/
+	window[i] = data[i-1];
+	nrlist[i] = outlist[i] = i;
+    }
+
+    /* sort the window[] -- sort *only* called here */
+    R_heapsort(1, k, window, outlist, nrlist, print_level);
+
+    big = fabs(window[k]);
+    if (big < fabs(window[1]))
+	big = fabs(window[1]);
+    /* big := max | X[1..k] | */
+    for(i=k; i < n; i++)
+	if (big < fabs(data[i]))
+	    big = fabs(data[i]);
+    /* big == max(|data_i|,  i = 1,..,n) */
+    big = 1 + 2. * big;/* such that -big < data[] < +big (strictly !) */
+
+    for(i=k; i > 0; i--) {
+	window[i+k2] = window[i];
+	nrlist[i+k2] = nrlist[i]-1;
+    }
+
+    for(i=0; i<k; i++)
+	outlist[i]=outlist[i+1]+k2;
+
+    k2p1 = k2+1;
+    for(i=0; i<k2p1; i++) {
+	window[i]	 = -big;
+	window[k+k2p1+i] =  big;
+    }
+} /* inittree*/
+
+static void
+toroot(int outvirt, int k, R_xlen_t nrnew, int outnext,
+       const double *data, double *window, int *outlist, int *nrlist,
+       int print_level)
+{
+    int father;
+
+    if(print_level >= 2) Rprintf("toroot(%d, %d,%d) ", k, (int) nrnew, outnext);
+
+    do {
+	father			  = outvirt/2;
+	window[outvirt+k]	  = window[father+k];
+	outlist[nrlist[father+k]] = outvirt+k;
+	nrlist[outvirt+k]	  = nrlist[father+k];
+	outvirt			  = father;
+    } while (father != 0);
+    window[k]	     = data[nrnew];
+    outlist[outnext] = k;
+    nrlist[k]	     = outnext;
+}
+
+static void
+downtoleave(int outvirt, int k,
+	    double *window, int *outlist, int *nrlist, int print_level)
+{
+    int childl, childr;
+
+    if(print_level >= 2) Rprintf("\n downtoleave(%d, %d)\n   ", outvirt,k);
+    for(;;) {
+	childl = outvirt*2;
+	childr = childl-1;
+	if (window[childl+k] < window[childr+k])
+	    childl = childr;
+	if (window[outvirt+k] >= window[childl+k])
+	    break;
+	/* seg.fault happens here: invalid outvirt/childl ? */
+	swap(outvirt+k, childl+k, window, outlist, nrlist, print_level);
+	outvirt = childl;
+    }
+}
+
+static void
+uptoleave(int outvirt, int k,
+	  double *window, int *outlist, int *nrlist, int print_level)
+{
+    int childl, childr;
+
+    if(print_level >= 2) Rprintf("\n uptoleave(%d, %d)\n   ", outvirt,k);
+    for(;;) {
+	childl = outvirt*2;
+	childr = childl+1;
+	if (window[childl+k] > window[childr+k])
+	    childl = childr;
+	if (window[outvirt+k] <= window[childl+k])
+	    break;
+	swap(outvirt+k, childl+k, window, outlist, nrlist, print_level);
+	outvirt = childl;
+    }
+}
+
+static void
+upperoutupperin(int outvirt, int k,
+		double *window, int *outlist, int *nrlist, int print_level)
+{
+    int father;
+
+    if(print_level >= 2) Rprintf("\nUpperoutUPPERin(%d, %d)\n  ", outvirt,k);
+    uptoleave(outvirt, k, window, outlist, nrlist, print_level);
+    father = outvirt/2;
+    while (window[outvirt+k] < window[father+k]) {
+	swap(outvirt+k, father+k, window, outlist, nrlist, print_level);
+	outvirt = father;
+	father	= outvirt/2;
+    }
+    if(print_level >= 2) Rprintf("\n");
+}
+
+static void
+upperoutdownin(int outvirt, int k, R_xlen_t nrnew, int outnext,
+	       const double *data, double *window, int *outlist, int *nrlist,
+	       int print_level)
+{
+    if(print_level >= 2) Rprintf("\n__upperoutDOWNin(%d, %d)\n  ", outvirt,k);
+    toroot(outvirt, k, nrnew, outnext, data, window, outlist, nrlist, print_level);
+    if(window[k] < window[k-1]) {
+	swap(k, k-1, window, outlist, nrlist, print_level);
+	downtoleave(/*outvirt = */ -1, k, window, outlist, nrlist, print_level);
+    }
+}
+
+static void
+downoutdownin(int outvirt, int k,
+	      double *window, int *outlist, int *nrlist, int print_level)
+{
+    int father;
+
+    if(print_level >= 2) Rprintf("\nDownoutDOWNin(%d, %d)\n  ", outvirt,k);
+    downtoleave(outvirt, k, window, outlist, nrlist, print_level);
+    father = outvirt/2;
+    while (window[outvirt+k] > window[father+k]) {
+	swap(outvirt+k, father+k, window, outlist, nrlist, print_level);
+	outvirt = father;
+	father	= outvirt/2;
+    }
+    if(print_level >= 2) Rprintf("\n");
+}
+
+static void
+downoutupperin(int outvirt, int k, R_xlen_t nrnew, int outnext,
+	       const double *data, double *window, int *outlist, int *nrlist,
+	       int print_level)
+{
+    if(print_level >= 2) Rprintf("\n__downoutUPPERin(%d, %d)\n  ", outvirt,k);
+    toroot(outvirt, k, nrnew, outnext, data, window, outlist, nrlist, print_level);
+    if(window[k] > window[k+1]) {
+	swap(k, k+1, window, outlist, nrlist, print_level);
+	uptoleave(/*outvirt = */ +1, k, window, outlist, nrlist, print_level);
+    }
+}
+
+static void
+wentoutone(int k, double *window, int *outlist, int *nrlist, int print_level)
+{
+    if(print_level >= 2) Rprintf("\nwentOUT_1(%d)\n  ", k);
+    swap(k, k+1, window, outlist, nrlist, print_level);
+    uptoleave(/*outvirt = */ +1, k, window, outlist, nrlist, print_level);
+}
+
+static void
+wentouttwo(int k, double *window, int *outlist, int *nrlist, int print_level)
+{
+    if(print_level >= 2) Rprintf("\nwentOUT_2(%d)\n  ", k);
+    swap(k, k-1, window, outlist, nrlist, print_level);
+    downtoleave(/*outvirt = */ -1, k, window, outlist, nrlist, print_level);
+}
+
+/* For Printing `diagnostics' : */
+#define Rm_PR(_F_,_A_)		for(j = 0; j <= 2*k; j++) Rprintf(_F_, _A_)
+#define RgPRINT_j(A_J)		Rm_PR("%6g", A_J); Rprintf("\n")
+#define RdPRINT_j(A_J)		Rm_PR("%6d", A_J); Rprintf("\n")
+
+#define R_PRINT_4vec()							\
+	Rprintf(" %9s: ","j");        RdPRINT_j(j);			\
+	Rprintf(" %9s: ","window []");RgPRINT_j(window[j]);		\
+	Rprintf(" %9s: "," nrlist[]");RdPRINT_j(nrlist[j]);		\
+	Rprintf(" %9s: ","outlist[]");RdPRINT_j((j <= k2|| j > k+k2)? -9\
+						: outlist[j - k2])
+
+static void
+runmedint(R_xlen_t n, int k, int k2, const double *data, double *median,
+	  double *window, int *outlist, int *nrlist,
+	  int end_rule, int print_level)
+{
+    /* Running Median of `k' ,  k == 2*k2 + 1 *
+     * end_rule == 0: leave values at the end,
+     *          otherwise: "constant" end values
+     */
+    int outnext, out, outvirt;
+
+    if(end_rule)
+	for(int i = 0; i <= k2; median[i++] = window[k]);
+    else {
+	for(int i = 0; i <  k2; median[i] = data[i], i++);
+	median[k2] = window[k];
+    }
+    outnext = 0;
+    for(R_xlen_t i = k2+1; i < n-k2; i++) {/* compute (0-index) median[i] == X*_{i+1} */
+	out  = outlist[outnext];
+	R_xlen_t nrnew = i+k2;
+	window[out] = data[nrnew];
+	outvirt	= out-k;
+	if (out > k)
+	    if(data[nrnew] >= window[k])
+		upperoutupperin(outvirt, k, window, outlist, nrlist, print_level);
+	    else
+		upperoutdownin(outvirt, k, nrnew, outnext,
+			       data, window, outlist, nrlist, print_level);
+	else if(out < k)
+	    if(data[nrnew] < window[k])
+		downoutdownin(outvirt, k, window, outlist, nrlist, print_level);
+	    else
+		downoutupperin(outvirt, k, nrnew, outnext,
+			       data, window, outlist, nrlist, print_level);
+	else if(window[k] > window[k+1])
+	    wentoutone(k, window, outlist, nrlist, print_level);
+	else if(window[k] < window[k-1])
+	    wentouttwo(k, window, outlist, nrlist, print_level);
+	median[i] = window[k];
+	outnext	  = (outnext+1)%k;
+    }
+    if(end_rule)
+	for(R_xlen_t i = n-k2; i < n; median[i++] = window[k]);
+    else
+	for(R_xlen_t i = n-k2; i < n; median[i] = data[i], i++);
+}/* runmedint() */
+
+/* This is the function called from R or S: */
+static void Trunmed(R_xlen_t n,/* = length(data) */
+		    int k,/* is odd <= n */
+		    const double *data,
+		    double *median, /* (n) */
+		    int   *outlist,/* (k+1) */
+		    int   *nrlist,/* (2k+1) */
+		    double *window,/* (2k+1) */
+		    int   end_rule,
+		    int   print_level)
+{
+    int k2 = (k - 1)/2, /* k == *kk == 2 * k2 + 1 */
+	j;
+
+    inittree (n, k, k2, data,
+	      /* initialize these: */
+	      window, (int *)outlist, (int *)nrlist, (int) print_level);
+
+    /* window[], outlist[], and nrlist[] are all 1-based (indices) */
+
+    if(print_level) {
+	Rprintf("After inittree():\n");
+	R_PRINT_4vec();
+    }
+    runmedint(n, k, k2, data, median, window, outlist, nrlist,
+	      end_rule, print_level);
+}

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/ansari.c

+/*
+ *  R : A Computer Language for Statistical Data Analysis
+ *  Copyright (C) 1999-2016  The R Core Team
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, a copy is available at
+ *  https://www.R-project.org/Licenses/
+ *
+ */
+
+/* ansari.c
+   Compute the exact distribution of the Ansari-Bradley test statistic.
+   */
+
+#include <string.h>
+#include <R.h>
+#include <math.h>  // for floor
+#include <Rmath.h>		/* uses choose() */
+#include "stats.h"
+
+static double ***
+w_init(int m, int n)
+{
+    int i;
+    double ***w;
+
+    w = (double ***) R_alloc(m + 1, sizeof(double **));
+    memset(w, '\0', (m+1) * sizeof(double**));
+    for (i = 0; i <= m; i++) {
+	w[i] = (double**) R_alloc(n + 1, sizeof(double *));
+	memset(w[i], '\0', (n+1) * sizeof(double*));
+    }
+    return(w);
+}
+
+
+static double
+cansari(int k, int m, int n, double ***w)
+{
+    int i, l, u;
+
+    l = (m + 1) * (m + 1) / 4;
+    u = l + m * n / 2;
+
+    if ((k < l) || (k > u))
+	return(0);
+
+    if (w[m][n] == 0) {
+	w[m][n] = (double *) R_alloc(u + 1, sizeof(double));
+	memset(w[m][n], '\0', (u + 1) * sizeof(double));
+	for (i = 0; i <= u; i++)
+	    w[m][n][i] = -1;
+    }
+
+    if (w[m][n][k] < 0) {
+	if (m == 0)
+	    w[m][n][k] = (k == 0);
+	else if (n == 0)
+	    w[m][n][k] = (k == l);
+	else
+	    w[m][n][k] = cansari(k, m, n - 1, w)
+		+ cansari(k - (m + n + 1) / 2, m - 1, n, w);
+    }
+
+    return(w[m][n][k]);
+}
+
+
+static void
+pansari(int len, double *Q, double *P, int m, int n)
+{
+    int i, j, l, u;
+    double c, p, q;
+    double ***w;
+
+    w = w_init(m, n);
+    l = (m + 1) * (m + 1) / 4;
+    u = l + m * n / 2;
+    c = choose(m + n, m);
+    for (i = 0; i < len; i++) {
+	q = floor(Q[i] + 1e-7);
+	if (q < l)
+	    P[i] = 0;
+	else if (q > u)
+	    P[i] = 1;
+	else {
+	    p = 0;
+	    for (j = l; j <= q; j++) p += cansari(j, m, n, w);
+	    P[i] = p / c;
+	}
+    }
+}
+
+static void
+qansari(int len, double *P, double *Q, int m, int n)
+{
+    int i, l, u;
+    double c, p, xi;
+    double ***w;
+
+    w = w_init(m, n);
+    l = (m + 1) * (m + 1) / 4;
+    u = l + m * n / 2;
+    c = choose(m + n, m);
+    for (i = 0; i < len; i++) {
+	xi = P[i];
+	if(xi < 0 || xi > 1)
+	    error(_("probabilities outside [0,1] in qansari()"));
+	if(xi == 0)
+	    Q[i] = l;
+	else if(xi == 1)
+	    Q[i] = u;
+	else {
+	    p = 0.;
+	    int q = 0;
+	    for(;;) {
+		p += cansari(q, m, n, w) / c;
+		if (p >= xi) break;
+		q++;
+	    }
+	    Q[i] = q;
+	}
+    }
+}
+
+#include <Rinternals.h>
+SEXP pAnsari(SEXP q, SEXP sm, SEXP sn)
+{
+    int m = asInteger(sm), n = asInteger(sn);
+    q = PROTECT(coerceVector(q, REALSXP));
+    int len = LENGTH(q);
+    SEXP p = PROTECT(allocVector(REALSXP, len));
+    pansari(len, REAL(q), REAL(p), m, n);
+    UNPROTECT(2);
+    return p;
+}
+
+SEXP qAnsari(SEXP p, SEXP sm, SEXP sn)
+{
+    int m = asInteger(sm), n = asInteger(sn);
+    p = PROTECT(coerceVector(p, REALSXP));
+    int len = LENGTH(p);
+    SEXP q = PROTECT(allocVector(REALSXP, len));
+    qansari(len, REAL(p), REAL(q), m, n);
+    UNPROTECT(2);
+    return q;
+}

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/bandwidths.c

+/*
+ *  R : A Computer Language for Statistical Data Analysis
+ *  bandwidth.c by W. N. Venables and B. D. Ripley  Copyright (C) 1994-2001
+ *  Copyright (C) 2012-2017  The R Core Team
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, a copy is available at
+ *  https://www.R-project.org/Licenses/
+ */
+
+#include <stdlib.h> //abs
+#include <math.h>
+#include <Rmath.h> // M_* constants
+#include <Rinternals.h>
+
+// or include "stats.h"
+#ifdef ENABLE_NLS
+#include <libintl.h>
+#define _(String) dgettext ("stats", String)
+#else
+#define _(String) (String)
+#endif
+
+#define DELMAX 1000
+/* Avoid slow and possibly error-producing underflows by cutting off at
+   plus/minus sqrt(DELMAX) std deviations */
+/* Formulae (6.67) and (6.69) of Scott (1992), the latter corrected. */
+
+SEXP bw_ucv(SEXP sn, SEXP sd, SEXP cnt, SEXP sh)
+{
+    double h = asReal(sh), d = asReal(sd), sum = 0.0, term, u;
+    int n = asInteger(sn), nbin = LENGTH(cnt);
+    double *x = REAL(cnt);
+    for (int i = 0; i < nbin; i++) {
+	double delta = i * d / h;
+	delta *= delta;
+	if (delta >= DELMAX) break;
+	term = exp(-delta / 4.) - sqrt(8.0) * exp(-delta / 2.);
+	sum += term * x[i];
+    }
+    u = (0.5 + sum/n) / (n * h * M_SQRT_PI);
+    // = 1 / (2 * n * h * sqrt(PI)) + sum / (n * n * h * sqrt(PI));
+    return ScalarReal(u);
+}
+
+SEXP bw_bcv(SEXP sn, SEXP sd, SEXP cnt, SEXP sh)
+{
+    double h = asReal(sh), d = asReal(sd), sum = 0.0, term, u;
+    int n = asInteger(sn), nbin = LENGTH(cnt);
+    double *x = REAL(cnt);
+
+    sum = 0.0;
+    for (int i = 0; i < nbin; i++) {
+	double delta = i * d / h; delta *= delta;
+	if (delta >= DELMAX) break;
+	term = exp(-delta / 4) * (delta * delta - 12 * delta + 12);
+	sum += term * x[i];
+    }
+    u = (1 + sum/(32.0*n)) / (2.0 * n * h * M_SQRT_PI);
+    // = 1 / (2 * n * h * sqrt(PI)) + sum / (64 * n * n * h * sqrt(PI));
+    return ScalarReal(u);
+}
+
+SEXP bw_phi4(SEXP sn, SEXP sd, SEXP cnt, SEXP sh)
+{
+    double h = asReal(sh), d = asReal(sd), sum = 0.0, term, u;
+    int n = asInteger(sn), nbin = LENGTH(cnt);
+    double *x = REAL(cnt);
+
+    for (int i = 0; i < nbin; i++) {
+	double delta = i * d / h; delta *= delta;
+	if (delta >= DELMAX) break;
+	term = exp(-delta / 2.) * (delta * delta - 6. * delta + 3.);
+	sum += term * x[i];
+    }
+    sum = 2. * sum + n * 3.;	/* add in diagonal */
+    u = sum / ((double)n * (n - 1) * pow(h, 5.0)) * M_1_SQRT_2PI;
+    // = sum / (n * (n - 1) * pow(h, 5.0) * sqrt(2 * PI));
+    return ScalarReal(u);
+}
+
+SEXP bw_phi6(SEXP sn, SEXP sd, SEXP cnt, SEXP sh)
+{
+    double h = asReal(sh), d = asReal(sd), sum = 0.0, term, u;
+    int n = asInteger(sn), nbin = LENGTH(cnt);
+    double *x = REAL(cnt);
+
+    for (int i = 0; i < nbin; i++) {
+	double delta = i * d / h; delta *= delta;
+	if (delta >= DELMAX) break;
+	term = exp(-delta / 2) *
+	    (delta * delta * delta - 15 * delta * delta + 45 * delta - 15);
+	sum += term * x[i];
+    }
+    sum = 2. * sum - 15. * n;	/* add in diagonal */
+    u = sum / ((double)n * (n - 1) * pow(h, 7.0)) * M_1_SQRT_2PI;
+    // = sum / (n * (n - 1) * pow(h, 7.0) * sqrt(2 * PI));
+    return ScalarReal(u);
+}
+
+/*
+   Use double cnt as from R 3.4.0, as counts can exceed INT_MAX for
+   large n (65537 in the worse case but typically not at n = 1 million
+   for a smooth distribution -- and this is by default no longer used
+   for n > 500).
+*/
+
+SEXP bw_den(SEXP nbin, SEXP sx)
+{
+    int nb = asInteger(nbin), n = LENGTH(sx);
+    double xmin, xmax, rang, dd, *x = REAL(sx);
+
+    xmin = R_PosInf; xmax = R_NegInf;
+    for (int i = 0; i < n; i++) {
+	if(!R_FINITE(x[i]))
+	    error(_("non-finite x[%d] in bandwidth calculation"), i+1);
+	if(x[i] < xmin) xmin = x[i];
+	if(x[i] > xmax) xmax = x[i];
+    }
+    rang = (xmax - xmin) * 1.01;
+    dd = rang / nb;
+
+    SEXP ans = PROTECT(allocVector(VECSXP, 2)),
+	sc = SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, nb));
+    SET_VECTOR_ELT(ans, 0, ScalarReal(dd));
+    double *cnt = REAL(sc);
+    for (int i = 0; i < nb; i++) cnt[i] = 0.0;
+
+    for (int i = 1; i < n; i++) {
+	int ii = (int)(x[i] / dd);
+	for (int j = 0; j < i; j++) {
+	    int jj = (int)(x[j] / dd);
+	    cnt[abs(ii - jj)] += 1.0;
+	}
+    }
+
+    UNPROTECT(1);
+    return ans;
+}
+
+/* Input: counts for nb bins */
+SEXP bw_den_binned(SEXP sx)
+{
+    int nb = LENGTH(sx);
+    int *x = INTEGER(sx);
+
+    SEXP ans = PROTECT(allocVector(REALSXP, nb));
+    double *cnt = REAL(ans);
+    for (int ib = 0; ib < nb; ib++) cnt[ib] = 0.0;
+
+    for (int ii = 0; ii < nb; ii++) {
+	int w = x[ii];
+	cnt[0] += w*(w-1.); // don't count distances to self
+	for (int jj = 0; jj < ii; jj++)
+	    cnt[ii - jj] += w * x[jj];
+    }
+    cnt[0] *= 0.5; // counts in the same bin got double-counted
+
+    UNPROTECT(1);
+    return ans;
+}

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/burg.c

+/*
+ *  R : A Computer Language for Statistical Data Analysis
+
+ *  Copyright (C) 1999-2016  The R Core Team
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, a copy is available at
+ *  https://www.R-project.org/Licenses/.
+ */
+
+#include <R.h>
+
+static void
+burg(int n, double*x, int pmax, double *coefs, double *var1, double *var2)
+{
+    double d, phii, *u, *v, *u0, sum;
+
+    u = (double *) R_alloc(n, sizeof(double));
+    v = (double *) R_alloc(n, sizeof(double));
+    u0 = (double *) R_alloc(n, sizeof(double));
+
+    for(int i = 0; i < pmax*pmax; i++) coefs[i] = 0.0;
+    sum = 0.0;
+    for(int t = 0; t < n; t++) {
+	u[t] = v[t] = x[n - 1 - t];
+	sum += x[t] * x[t];
+    }
+    var1[0] = var2[0] = sum/n;
+    for(int p = 1; p <= pmax; p++) { /* do AR(p) */
+	sum = 0.0;
+	d = 0;
+	for(int t = p; t < n; t++) {
+	    sum += v[t]*u[t-1];
+	    d += v[t]*v[t] + u[t-1]*u[t-1];
+	}
+	phii = 2*sum/d;
+	coefs[pmax*(p-1) + (p-1)] = phii;
+	if(p > 1)
+	    for(int j = 1; j < p; j++)
+		coefs[p-1 + pmax*(j-1)] =
+		    coefs[p-2 + pmax*(j-1)] - phii* coefs[p-2 + pmax*(p-j-1)];
+	/* update u and v */
+	for(int t = 0; t < n; t++)
+	    u0[t] = u[t];
+	for(int t = p; t < n; t++) {
+	    u[t] = u0[t-1] - phii * v[t];
+	    v[t] = v[t] - phii * u0[t-1];
+	}
+	var1[p] = var1[p-1] * (1 - phii * phii);
+	d = 0.0;
+	for(int t = p; t < n; t++) d += v[t]*v[t] + u[t]*u[t];
+	var2[p] = d/(2.0*(n-p));
+    }
+}
+
+#include <Rinternals.h>
+
+SEXP Burg(SEXP x, SEXP order)
+{
+    x = PROTECT(coerceVector(x, REALSXP));
+    int n = LENGTH(x), pmax = asInteger(order);
+    SEXP coefs = PROTECT(allocVector(REALSXP, pmax * pmax)),
+	var1 = PROTECT(allocVector(REALSXP, pmax + 1)),
+	var2 = PROTECT(allocVector(REALSXP, pmax + 1));
+    burg(n, REAL(x), pmax, REAL(coefs), REAL(var1), REAL(var2));
+    SEXP ans = PROTECT(allocVector(VECSXP, 3));
+    SET_VECTOR_ELT(ans, 0, coefs);
+    SET_VECTOR_ELT(ans, 1, var1);
+    SET_VECTOR_ELT(ans, 2, var2);
+    UNPROTECT(5);
+    return ans;
+}

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/d2x2xk.c

+/*
+ *  R : A Computer Language for Statistical Data Analysis
+ *  Copyright (C) 2000-2016  The R Core Team.
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, a copy is available at
+ *  https://www.R-project.org/Licenses/
+ */
+
+/* for mantelhaen.test */
+
+#include <R.h>
+#include <Rmath.h>
+
+static void
+int_d2x2xk(int K, double *m, double *n, double *t, double *d)
+{
+    int i, j, l, w, y, z;
+    double u, **c;
+
+    c = (double **) R_alloc(K + 1, sizeof(double *));
+    l = y = z = 0;
+    c[0] = (double *) R_alloc(1, sizeof(double));
+    c[0][0] = 1;
+    for(i = 0; i < K; i++) {
+	y = imax2(0,  (int)(*t - *n));
+	z = imin2((int)*m, (int)*t);
+	c[i + 1] = (double *) R_alloc(l + z - y + 1, sizeof(double));
+	for(j = 0; j <= l + z - y; j++) c[i + 1][j] = 0;
+	for(j = 0; j <= z - y; j++) {
+	    u = dhyper(j + y, *m, *n, *t, FALSE);
+	    for(w = 0; w <= l; w++) c[i + 1][w + j] += c[i][w] * u;
+	}
+	l = l + z - y;
+	m++; n++; t++;
+    }
+
+    u = 0;
+    for(j = 0; j <= l; j++) u += c[K][j];
+    for(j = 0; j <= l; j++) d[j] = c[K][j] / u;
+}
+
+#include <Rinternals.h>
+
+SEXP d2x2xk(SEXP sK, SEXP m, SEXP n, SEXP t, SEXP srn)
+{
+    int K = asInteger(sK), rn = asInteger(srn);
+    m = PROTECT(coerceVector(m, REALSXP));
+    n = PROTECT(coerceVector(n, REALSXP));
+    t = PROTECT(coerceVector(t, REALSXP));
+    SEXP ans = PROTECT(allocVector(REALSXP, rn));
+    int_d2x2xk(K, REAL(m), REAL(n), REAL(t), REAL(ans));
+    UNPROTECT(4);
+    return ans;
+}

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/family.c

+/*
+ *  R : A Computer Language for Statistical Data Analysis
+ *  Copyright (C) 2005-2016  The R Core Team
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, a copy is available at
+ *  https://www.R-project.org/Licenses/
+ *
+ *  Quartz Quartz device module header file
+ *
+ */
+
+#include <Rinternals.h>
+#include <Rconfig.h>
+#include <R_ext/Constants.h>
+#include <float.h>
+#include <math.h>
+#include "stats.h"
+#include "statsR.h"
+
+static const double THRESH = 30.;
+static const double MTHRESH = -30.;
+static const double INVEPS = 1/DOUBLE_EPS;
+
+/**
+ * Evaluate x/(1 - x). An inline function is used so that x is
+ * evaluated once only.
+ *
+ * @param x input in the range (0, 1)
+ *
+ * @return x/(1 - x)
+ */
+static R_INLINE double x_d_omx(double x) {
+    if (x < 0 || x > 1)
+	error(_("Value %g out of range (0, 1)"), x);
+    return x/(1 - x);
+}
+
+/**
+ * Evaluate x/(1 + x). An inline function is used so that x is
+ * evaluated once only. [but inlining is optional!]
+ *
+ * @param x input
+ *
+ * @return x/(1 + x)
+ */
+static R_INLINE double x_d_opx(double x) {return x/(1 + x);}
+
+SEXP logit_link(SEXP mu)
+{
+    int i, n = LENGTH(mu);
+    SEXP ans = PROTECT(shallow_duplicate(mu));
+    double *rans = REAL(ans), *rmu=REAL(mu);
+
+    if (!n || !isReal(mu))
+	error(_("Argument %s must be a nonempty numeric vector"), "mu");
+    for (i = 0; i < n; i++)
+	rans[i] = log(x_d_omx(rmu[i]));
+    UNPROTECT(1);
+    return ans;
+}
+
+SEXP logit_linkinv(SEXP eta)
+{
+    SEXP ans = PROTECT(shallow_duplicate(eta));
+    int i, n = LENGTH(eta);
+    double *rans = REAL(ans), *reta = REAL(eta);
+
+    if (!n || !isReal(eta))
+	error(_("Argument %s must be a nonempty numeric vector"), "eta");
+    for (i = 0; i < n; i++) {
+	double etai = reta[i], tmp;
+	tmp = (etai < MTHRESH) ? DOUBLE_EPS :
+	    ((etai > THRESH) ? INVEPS : exp(etai));
+	rans[i] = x_d_opx(tmp);
+    }
+    UNPROTECT(1);
+    return ans;
+}
+
+SEXP logit_mu_eta(SEXP eta)
+{
+    SEXP ans = PROTECT(shallow_duplicate(eta));
+    int i, n = LENGTH(eta);
+    double *rans = REAL(ans), *reta = REAL(eta);
+
+    if (!n || !isReal(eta))
+	error(_("Argument %s must be a nonempty numeric vector"), "eta");
+    for (i = 0; i < n; i++) {
+	double etai = reta[i];
+	double opexp = 1 + exp(etai);
+
+	rans[i] = (etai > THRESH || etai < MTHRESH) ? DOUBLE_EPS :
+	    exp(etai)/(opexp * opexp);
+    }
+    UNPROTECT(1);
+    return ans;
+}
+
+static R_INLINE
+double y_log_y(double y, double mu)
+{
+    return (y != 0.) ? (y * log(y/mu)) : 0;
+}
+
+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(shallow_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;
+}

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/filter.c

+/*
+ *  R : A Computer Language for Statistical Data Analysis
+
+ *  Copyright (C) 1999-2016   The R Core Team
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, a copy is available at
+ *  https://www.R-project.org/Licenses/.
+ */
+
+#ifdef HAVE_CONFIG_H
+# include <config.h>
+#endif
+
+#include <R.h>
+#include "ts.h"
+
+#ifndef min
+#define min(a, b) ((a < b)?(a):(b))
+#define max(a, b) ((a < b)?(b):(a))
+#endif
+
+// currently ISNAN includes NAs
+#define my_isok(x) (!ISNA(x) & !ISNAN(x))
+
+SEXP cfilter(SEXP sx, SEXP sfilter, SEXP ssides, SEXP scircular)
+{
+   if (TYPEOF(sx) != REALSXP || TYPEOF(sfilter) != REALSXP)
+       error("invalid input");
+    R_xlen_t nx = XLENGTH(sx), nf = XLENGTH(sfilter);
+    int sides = asInteger(ssides), circular = asLogical(scircular);
+    if(sides == NA_INTEGER || circular == NA_LOGICAL)  error("invalid input");
+
+    SEXP ans = allocVector(REALSXP, nx);
+
+    R_xlen_t i, j, nshift;
+    double z, tmp, *x = REAL(sx), *filter = REAL(sfilter), *out = REAL(ans);
+
+    if(sides == 2) nshift = nf /2; else nshift = 0;
+    if(!circular) {
+	for(i = 0; i < nx; i++) {
+	    z = 0;
+	    if(i + nshift - (nf - 1) < 0 || i + nshift >= nx) {
+		out[i] = NA_REAL;
+		continue;
+	    }
+	    for(j = max(0, nshift + i - nx); j < min(nf, i + nshift + 1) ; j++) {
+		tmp = x[i + nshift - j];
+		if(my_isok(tmp)) z += filter[j] * tmp;
+		else { out[i] = NA_REAL; goto bad; }
+	    }
+	    out[i] = z;
+	bad:
+	    continue;
+	}
+    } else { /* circular */
+	for(i = 0; i < nx; i++)
+	{
+	    z = 0;
+	    for(j = 0; j < nf; j++) {
+		R_xlen_t ii = i + nshift - j;
+		if(ii < 0) ii += nx;
+		if(ii >= nx) ii -= nx;
+		tmp = x[ii];
+		if(my_isok(tmp)) z += filter[j] * tmp;
+		else { out[i] = NA_REAL; goto bad2; }
+	    }
+	    out[i] = z;
+	bad2:
+	    continue;
+	}
+    }
+    return ans;
+}
+
+/* recursive filtering */
+SEXP rfilter(SEXP x, SEXP filter, SEXP out)
+{
+   if (TYPEOF(x) != REALSXP || TYPEOF(filter) != REALSXP
+       || TYPEOF(out) != REALSXP) error("invalid input");
+    R_xlen_t nx = XLENGTH(x), nf = XLENGTH(filter);
+    double sum, tmp, *r = REAL(out), *rx = REAL(x), *rf = REAL(filter);
+
+    for(R_xlen_t i = 0; i < nx; i++) {
+	sum = rx[i];
+	for (R_xlen_t j = 0; j < nf; j++) {
+	    tmp = r[nf + i - j - 1];
+	    if(my_isok(tmp)) sum += tmp * rf[j];
+	    else { r[nf + i] = NA_REAL; goto bad3; }
+	}
+	r[nf + i] = sum;
+    bad3:
+	continue;
+    }
+    return out;
+}
+
+/* now allows missing values */
+static void
+acf0(double *x, int n, int ns, int nl, Rboolean correlation, double *acf)
+{
+    int d1 = nl+1, d2 = ns*d1;
+
+    for(int u = 0; u < ns; u++)
+	for(int v = 0; v < ns; v++)
+	    for(int lag = 0; lag <= nl; lag++) {
+		double sum = 0.0; int nu = 0;
+		for(int i = 0; i < n-lag; i++)
+		    if(!ISNAN(x[i + lag + n*u]) && !ISNAN(x[i + n*v])) {
+			nu++;
+			sum += x[i + lag + n*u] * x[i + n*v];
+		    }
+		acf[lag + d1*u + d2*v] = (nu > 0) ? sum/(nu + lag) : NA_REAL;
+	    }
+    if(correlation) {
+	if(n == 1) {
+	    for(int u = 0; u < ns; u++)
+		acf[0 + d1*u + d2*u] = 1.0;
+	} else {
+	    double *se = (double *) R_alloc(ns, sizeof(double));
+	    for(int u = 0; u < ns; u++)
+		se[u] = sqrt(acf[0 + d1*u + d2*u]);
+	    for(int u = 0; u < ns; u++)
+		for(int v = 0; v < ns; v++)
+		    for(int lag = 0; lag <= nl; lag++) { // ensure correlations remain in  [-1,1] :
+			double a = acf[lag + d1*u + d2*v] / (se[u]*se[v]);
+			acf[lag + d1*u + d2*v] = (a > 1.) ? 1. : ((a < -1.) ? -1. : a);
+		    }
+	}
+    }
+}
+
+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;
+}

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/kendall.c

+/*
+ *  R : A Computer Language for Statistical Data Analysis
+ *  Copyright (C) 1999-2016   The R Core Team.
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, a copy is available at
+ *  https://www.R-project.org/Licenses/
+ */
+
+/* Kendall's rank correlation tau and its exact distribution in case of no ties
+*/
+
+#include <string.h>
+#include <R.h>
+#include <math.h> // for floor
+#include <Rmath.h>
+
+/*
+   and the exact distribution of  T = (n * (n - 1) * tau + 1) / 4,
+   which is -- if there are no ties -- the number of concordant ordered pairs
+*/
+
+static double
+ckendall(int k, int n, double **w) 
+{
+    int i, u;
+    double s;
+
+    u =  (n * (n - 1) / 2);
+    if ((k < 0) || (k > u)) return(0);
+    if (w[n] == 0) {
+	w[n] = (double *) R_alloc(u + 1, sizeof(double));
+	memset(w[n], '\0', sizeof(double) * (u+1));
+	for (i = 0; i <= u; i++) w[n][i] = -1;
+    }
+    if (w[n][k] < 0) {
+	if (n == 1)
+	    w[n][k] = (k == 0);
+	else {
+	    s = 0;
+	    for (i = 0; i < n; i++)
+		s += ckendall(k - i, n - 1, w);
+	    w[n][k] = s;
+	}
+    }
+    return(w[n][k]);
+}
+
+#if 0
+void
+dkendall(int *len, double *x, int *n) 
+{
+    int i;
+    double **w;
+
+    w = (double **) R_alloc(*n + 1, sizeof(double *));
+
+    for (i = 0; i < *len; i++)
+	if (fabs(x[i] - floor(x[i] + 0.5)) > 1e-7) {
+	    x[i] = 0;
+	} else {
+	    x[i] = ckendall((int)x[i], *n) / gammafn(*n + 1, w);
+	}
+}
+#endif
+
+static void
+pkendall(int len, double *Q, double *P, int n) 
+{
+    int i, j;
+    double p, q;
+    double **w;
+
+    w = (double **) R_alloc(n + 1, sizeof(double *));
+    memset(w, '\0', sizeof(double*) * (n+1));
+
+    for (i = 0; i < len; i++) {
+	q = floor(Q[i] + 1e-7);
+	if (q < 0)
+	    P[i] = 0;
+	else if (q > (n * (n - 1) / 2))
+	    P[i] = 1;
+	else {
+	    p = 0;
+	    for (j = 0; j <= q; j++) p += ckendall(j, n, w);
+	    P[i] = p / gammafn(n + 1);
+	}
+    }
+}
+
+#include <Rinternals.h>
+SEXP pKendall(SEXP q, SEXP sn)
+{
+    q = PROTECT(coerceVector(q, REALSXP));
+    int len = LENGTH(q), n = asInteger(sn);
+    SEXP p = PROTECT(allocVector(REALSXP, len));
+    pkendall(len, REAL(q), REAL(p), n);
+    UNPROTECT(2);
+    return p;
+}

com.oracle.truffle.r.native/gnur/patch/src/library/stats/src/ksmooth.c

+/*
+ *  R : A Computer Language for Statistical Data Analysis
+ *  Copyright (C) 1998-2016	The R Foundation
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, a copy is available at
+ *  https://www.R-project.org/Licenses/
+ */
+
+#include <math.h>
+#include <R.h>			/* for NA_REAL, includes math.h */
+#include <Rinternals.h>
+
+#ifdef ENABLE_NLS
+#include <libintl.h>
+#define _(String) dgettext ("stats", String)
+#else
+#define _(String) (String)
+#endif
+
+static double dokern(double x, int kern)
+{
+    if(kern == 1) return(1.0);
+    if(kern == 2) return(exp(-0.5*x*x));
+    return(0.0); /* -Wall */
+}
+
+static void BDRksmooth(double *x, double *y, R_xlen_t n,
+		       double *xp, double *yp, R_xlen_t np,
+		       int kern, double bw)
+{
+    R_xlen_t imin = 0;
+    double cutoff = 0.0, num, den, x0, w;
+
+    /* bandwidth is in units of half inter-quartile range. */
+    if(kern == 1) {bw *= 0.5; cutoff = bw;}
+    if(kern == 2) {bw *= 0.3706506; cutoff = 4*bw;}
+    while(x[imin] < xp[0] - cutoff && imin < n) imin++;
+    for(R_xlen_t j = 0; j < np; j++) {
+	num = den = 0.0;
+	x0 = xp[j];
+	for(R_xlen_t i = imin; i < n; i++) {
+	    if(x[i] < x0 - cutoff) imin = i;
+	    else {
+		if(x[i] > x0 + cutoff) break;
+		w = dokern(fabs(x[i] - x0)/bw, kern);
+		num += w*y[i];
+		den += w;
+	    }
+	}
+	if(den > 0) yp[j] = num/den; else yp[j] = NA_REAL;
+    }
+}
+
+
+// called only from  spline()  in ./ppr.f
+void NORET F77_SUB(bdrsplerr)(void)
+{
+    error(_("only 2500 rows are allowed for sm.method=\"spline\""));
+}
+
+void F77_SUB(splineprt)(double* df, double* gcvpen, int* ismethod,
+			      double* lambda, double *edf)
+{
+    Rprintf("spline(df=%5.3g, g.pen=%11.6g, ismeth.=%+2d) -> (lambda, edf) = (%.7g, %5.2f)\n",
+	    *df, *gcvpen, *ismethod, *lambda, *edf);
+    return;
+}
+
+// called only from smooth(..., trace=TRUE)  in ./ppr.f :
+void F77_SUB(smoothprt)(double* span, int* iper, double* var, double* cvar)
+{
+    Rprintf("smooth(span=%4g, iper=%+2d) -> (var, cvar) = (%g, %g)\n",
+	    *span, *iper, *var, *cvar);
+    return;
+}
+
+
+SEXP ksmooth(SEXP x, SEXP y, SEXP xp, SEXP skrn, SEXP sbw)
+{
+    int krn = asInteger(skrn);
+    double bw = asReal(sbw);
+    x = PROTECT(coerceVector(x, REALSXP));
+    y = PROTECT(coerceVector(y, REALSXP));
+    xp = PROTECT(coerceVector(xp, REALSXP));
+    R_xlen_t nx = XLENGTH(x), np = XLENGTH(xp);
+    SEXP yp = PROTECT(allocVector(REALSXP, np));
+
+    BDRksmooth(REAL(x), REAL(y), nx, REAL(xp), REAL(yp), np, krn, bw);
+    SEXP ans = PROTECT(allocVector(VECSXP, 2));
+    SET_VECTOR_ELT(ans, 0, xp);
+    SET_VECTOR_ELT(ans, 1, yp);
+    SEXP nm = allocVector(STRSXP, 2);
+    setAttrib(ans, R_NamesSymbol, nm);
+    SET_STRING_ELT(nm, 0, mkChar("x"));
+    SET_STRING_ELT(nm, 1, mkChar("y"));
+    UNPROTECT(5);
+    return ans;
+}