From 5569b0514fbed053fa762d9542a10175d9f299d2 Mon Sep 17 00:00:00 2001
From: stepan <stepan.sindelar@oracle.com>
Date: Fri, 6 Jan 2017 12:06:39 +0100
Subject: [PATCH] Stats: port {q/p}tukey functions

---
 .../foreign/CallAndExternalFunctions.java     |   6 +
 .../truffle/r/runtime/nmath/distr/PTukey.java | 310 ++++++++++++++++++
 .../truffle/r/runtime/nmath/distr/QTukey.java | 155 +++++++++
 .../truffle/r/test/ExpectedTestOutput.test    | 254 ++++++++++++++
 .../test/library/stats/TestDistributions.java |  16 +-
 mx.fastr/copyrights/overrides                 |   2 +
 6 files changed, 742 insertions(+), 1 deletion(-)
 create mode 100644 com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/PTukey.java
 create mode 100644 com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/QTukey.java

diff --git a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java
index b21920e7d2..1c9a963a3b 100644
--- a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java
+++ b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java
@@ -127,6 +127,7 @@ import com.oracle.truffle.r.runtime.nmath.distr.PNBinom.PNBinomFunc;
 import com.oracle.truffle.r.runtime.nmath.distr.PNBinom.PNBinomMu;
 import com.oracle.truffle.r.runtime.nmath.distr.PNChisq;
 import com.oracle.truffle.r.runtime.nmath.distr.PPois;
+import com.oracle.truffle.r.runtime.nmath.distr.PTukey;
 import com.oracle.truffle.r.runtime.nmath.distr.Pbeta;
 import com.oracle.truffle.r.runtime.nmath.distr.Pbinom;
 import com.oracle.truffle.r.runtime.nmath.distr.Pf;
@@ -142,6 +143,7 @@ import com.oracle.truffle.r.runtime.nmath.distr.QNBinom.QNBinomFunc;
 import com.oracle.truffle.r.runtime.nmath.distr.QNBinom.QNBinomMu;
 import com.oracle.truffle.r.runtime.nmath.distr.QNChisq;
 import com.oracle.truffle.r.runtime.nmath.distr.QPois;
+import com.oracle.truffle.r.runtime.nmath.distr.QTukey;
 import com.oracle.truffle.r.runtime.nmath.distr.Qbinom;
 import com.oracle.truffle.r.runtime.nmath.distr.Qf;
 import com.oracle.truffle.r.runtime.nmath.distr.Qnf;
@@ -476,6 +478,10 @@ public class CallAndExternalFunctions {
                     return StatsFunctionsNodes.Function3_2Node.create(new Qnt());
                 case "qsignrank":
                     return StatsFunctionsNodes.Function2_2Node.create(new QSignrank());
+                case "qtukey":
+                    return StatsFunctionsNodes.Function4_2Node.create(new QTukey());
+                case "ptukey":
+                    return StatsFunctionsNodes.Function4_2Node.create(new PTukey());
                 case "rmultinom":
                     return RMultinomNode.create();
                 case "Approx":
diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/PTukey.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/PTukey.java
new file mode 100644
index 0000000000..b238d42c77
--- /dev/null
+++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/PTukey.java
@@ -0,0 +1,310 @@
+/*
+ * This material is distributed under the GNU General Public License
+ * Version 2. You may review the terms of this license at
+ * http://www.gnu.org/licenses/gpl-2.0.html
+ *
+ * Copyright (C) 1998 Ross Ihaka
+ * Copyright (c) 2000--2007, The R Core Team
+ * Copyright (c) 2017, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+package com.oracle.truffle.r.runtime.nmath.distr;
+
+import static com.oracle.truffle.r.runtime.nmath.GammaFunctions.lgammafn;
+import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_1_SQRT_2PI;
+import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_LN2;
+
+import com.oracle.truffle.r.runtime.nmath.DPQ;
+import com.oracle.truffle.r.runtime.nmath.MathFunctions.Function4_2;
+import com.oracle.truffle.r.runtime.nmath.RMathError;
+import com.oracle.truffle.r.runtime.nmath.RMathError.MLError;
+
+public class PTukey implements Function4_2 {
+    private static final int nlegq = 16;
+    private static final int ihalfq = 8;
+    private static final int nleg = 12;
+    private static final int ihalf = 6;
+    private static final double eps1 = -30.0;
+    private static final double eps2 = 1.0e-14;
+    private static final double dhaf = 100.0;
+    private static final double dquar = 800.0;
+    private static final double deigh = 5000.0;
+    private static final double dlarg = 25000.0;
+    private static final double ulen1 = 1.0;
+    private static final double ulen2 = 0.5;
+    private static final double ulen3 = 0.25;
+    private static final double ulen4 = 0.125;
+    private static final double[] xlegq = {
+                    0.989400934991649932596154173450,
+                    0.944575023073232576077988415535,
+                    0.865631202387831743880467897712,
+                    0.755404408355003033895101194847,
+                    0.617876244402643748446671764049,
+                    0.458016777657227386342419442984,
+                    0.281603550779258913230460501460,
+                    0.950125098376374401853193354250e-1
+    };
+    private static final double[] alegq = {
+                    0.271524594117540948517805724560e-1,
+                    0.622535239386478928628438369944e-1,
+                    0.951585116824927848099251076022e-1,
+                    0.124628971255533872052476282192,
+                    0.149595988816576732081501730547,
+                    0.169156519395002538189312079030,
+                    0.182603415044923588866763667969,
+                    0.189450610455068496285396723208
+    };
+
+    private final Pnorm pnorm = new Pnorm();
+
+    @Override
+    public double evaluate(double q, double rr, double cc, double df, boolean lowerTail, boolean logP) {
+        if (Double.isNaN(q) || Double.isNaN(rr) || Double.isNaN(cc) || Double.isNaN(df)) {
+            return RMathError.defaultError();
+        }
+
+        if (q <= 0) {
+            return DPQ.rdt0(lowerTail, logP);
+        }
+
+        /* df must be > 1 */
+        /* there must be at least two values */
+
+        if (df < 2 || rr < 1 || cc < 2) {
+            return RMathError.defaultError();
+        }
+
+        if (!Double.isFinite(q)) {
+            return DPQ.rdt1(lowerTail, logP);
+        }
+
+        if (df > dlarg) {
+            return DPQ.rdtval(wprob(q, rr, cc), lowerTail, logP);
+        }
+
+        /* calculate leading constant */
+
+        double f2 = df * 0.5;
+        /* lgammafn(u) = Math.log(gamma(u)) */
+        double f2lf = ((f2 * Math.log(df)) - (df * M_LN2)) - lgammafn(f2);
+        double f21 = f2 - 1.0;
+
+        /* integral is divided into unit, half-unit, quarter-unit, or */
+        /* eighth-unit length intervals depending on the value of the */
+        /* degrees of freedom. */
+
+        double ff4 = df * 0.25;
+        double ulen = getULen(df);
+        f2lf += Math.log(ulen);
+
+        /* integrate over each subinterval */
+        double ans = 0.0;
+
+        double otsum = 0.0;
+        for (int i = 1; i <= 50; i++) {
+            otsum = 0.0;
+
+            /* legendre quadrature with order = nlegq */
+            /* nodes (stored in xlegq) are symmetric around zero. */
+
+            double twa1 = (2 * i - 1) * ulen;
+
+            for (int jj = 1; jj <= nlegq; jj++) {
+                int j;
+                double t1;
+                if (ihalfq < jj) {
+                    j = jj - ihalfq - 1;
+                    t1 = (f2lf + (f21 * Math.log(twa1 + (xlegq[j] * ulen)))) - (((xlegq[j] * ulen) + twa1) * ff4);
+                } else {
+                    j = jj - 1;
+                    t1 = (f2lf + (f21 * Math.log(twa1 - (xlegq[j] * ulen)))) + (((xlegq[j] * ulen) - twa1) * ff4);
+
+                }
+
+                /* if Math.exp(t1) < 9e-14, then doesn't contribute to integral */
+                if (t1 >= eps1) {
+                    double qsqz;
+                    if (ihalfq < jj) {
+                        qsqz = q * Math.sqrt(((xlegq[j] * ulen) + twa1) * 0.5);
+                    } else {
+                        qsqz = q * Math.sqrt(((-(xlegq[j] * ulen)) + twa1) * 0.5);
+                    }
+
+                    /* call wprob to find integral of range portion */
+
+                    double wprb = wprob(qsqz, rr, cc);
+                    double rotsum = (wprb * alegq[j]) * Math.exp(t1);
+                    otsum += rotsum;
+                }
+                /* end legendre integral for interval i */
+                /* L200: */
+            }
+
+            /*
+             * if integral for interval i < 1e-14, then stop. However, in order to avoid small area
+             * under left tail, at least 1 / ulen intervals are calculated.
+             */
+            if (i * ulen >= 1.0 && otsum <= eps2) {
+                break;
+            }
+
+            /* end of interval i */
+            /* L330: */
+
+            ans += otsum;
+        }
+
+        if (otsum > eps2) { /* not converged */
+            RMathError.error(MLError.PRECISION, "ptukey");
+        }
+        if (ans > 1.) {
+            ans = 1.;
+        }
+        return DPQ.rdtval(ans, lowerTail, logP);
+    }
+
+    private double getULen(double df) {
+        if (df <= dhaf) {
+            return ulen1;
+        } else if (df <= dquar) {
+            return ulen2;
+        } else if (df <= deigh) {
+            return ulen3;
+        }
+        return ulen4;
+    }
+
+    private static final double C1 = -30.;
+    private static final double C2 = -50.;
+    private static final double C3 = 60.;
+    private static final double bb = 8.;
+    private static final double wlar = 3.;
+    private static final double wincr1 = 2.;
+    private static final double wincr2 = 3.;
+    private static final double[] xleg = {
+                    0.981560634246719250690549090149,
+                    0.904117256370474856678465866119,
+                    0.769902674194304687036893833213,
+                    0.587317954286617447296702418941,
+                    0.367831498998180193752691536644,
+                    0.125233408511468915472441369464
+    };
+    private static final double[] aleg = {
+                    0.047175336386511827194615961485,
+                    0.106939325995318430960254718194,
+                    0.160078328543346226334652529543,
+                    0.203167426723065921749064455810,
+                    0.233492536538354808760849898925,
+                    0.249147045813402785000562436043
+    };
+
+    private double wprob(double w, double rr, double cc) {
+        // double a, ac, prW, b, binc, c, cc1,
+        // pminus, pplus, qexpo, qsqz, rinsum, wi, wincr, xx;
+        // LDOUBLE blb, bub, einsum, elsum;
+        // int j, jj;
+
+        double qsqz = w * 0.5;
+
+        /* if w >= 16 then the integral lower bound (occurs for c=20) */
+        /* is 0.99999999999995 so return a value of 1. */
+
+        if (qsqz >= bb) {
+            return 1.0;
+        }
+
+        /* find (f(w/2) - 1) ^ cc */
+        /* (first term in integral of hartley's form). */
+
+        double prW = 2 * pnorm.evaluate(qsqz, 0., 1., true, false) - 1.; /* erf(qsqz / M_SQRT2) */
+        /* if prW ^ cc < 2e-22 then set prW = 0 */
+        if (prW >= Math.exp(C2 / cc)) {
+            prW = Math.pow(prW, cc);
+        } else {
+            prW = 0.0;
+        }
+
+        /* if w is large then the second component of the */
+        /* integral is small, so fewer intervals are needed. */
+
+        double wincr = w > wlar ? wincr1 : wincr2;
+
+        /* find the integral of second term of hartley's form */
+        /* for the integral of the range for equal-length */
+        /* intervals using legendre quadrature. limits of */
+        /* integration are from (w/2, 8). two or three */
+        /* equal-length intervals are used. */
+
+        /* blb and bub are lower and upper limits of integration. */
+
+        /* LDOUBLE */double blb = qsqz;
+        double binc = (bb - qsqz) / wincr;
+        /* LDOUBLE */double bub = blb + binc;
+        /* LDOUBLE */double einsum = 0.0;
+
+        /* integrate over each interval */
+
+        double cc1 = cc - 1.0;
+        for (double wi = 1; wi <= wincr; wi++) {
+            /* LDOUBLE */double elsum = 0.0;
+            double a = 0.5 * (bub + blb);
+
+            /* legendre quadrature with order = nleg */
+
+            double b = 0.5 * (bub - blb);
+
+            for (int jj = 1; jj <= nleg; jj++) {
+                double xx;
+                int j;
+                if (ihalf < jj) {
+                    j = (nleg - jj) + 1;
+                    xx = xleg[j - 1];
+                } else {
+                    j = jj;
+                    xx = -xleg[j - 1];
+                }
+                double c = b * xx;
+                double ac = a + c;
+
+                /* if Math.exp(-qexpo/2) < 9e-14, */
+                /* then doesn't contribute to integral */
+
+                double qexpo = ac * ac;
+                if (qexpo > C3) {
+                    break;
+                }
+
+                double pplus = 2 * pnorm.evaluate(ac, 0., 1., true, false);
+                double pminus = 2 * pnorm.evaluate(ac, w, 1., true, false);
+
+                /* if rinsum ^ (cc-1) < 9e-14, */
+                /* then doesn't contribute to integral */
+
+                double rinsum = (pplus * 0.5) - (pminus * 0.5);
+                if (rinsum >= Math.exp(C1 / cc1)) {
+                    rinsum = (aleg[j - 1] * Math.exp(-(0.5 * qexpo))) * Math.pow(rinsum, cc1);
+                    elsum += rinsum;
+                }
+            }
+            elsum *= (((2.0 * b) * cc) * M_1_SQRT_2PI);
+            einsum += elsum;
+            blb = bub;
+            bub += binc;
+        }
+
+        /* if prW ^ rr < 9e-14, then return 0 */
+        prW += einsum;
+        if (prW <= Math.exp(C1 / rr)) {
+            return 0.;
+        }
+
+        prW = Math.pow(prW, rr);
+        if (prW >= 1.) {
+            /* 1 was iMax was eps */
+            return 1.;
+        }
+
+        return prW;
+    }
+}
diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/QTukey.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/QTukey.java
new file mode 100644
index 0000000000..dbd69d4df5
--- /dev/null
+++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/QTukey.java
@@ -0,0 +1,155 @@
+/*
+ * This material is distributed under the GNU General Public License
+ * Version 2. You may review the terms of this license at
+ * http://www.gnu.org/licenses/gpl-2.0.html
+ *
+ * Copyright (C) 1998 Ross Ihaka
+ * Copyright (c) 2000--2005, The R Core Team
+ * Copyright (c) 2017, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+/*
+ *  based in part on AS70 (C) 1974 Royal Statistical Society
+ */
+package com.oracle.truffle.r.runtime.nmath.distr;
+
+import com.oracle.truffle.r.runtime.nmath.DPQ;
+import com.oracle.truffle.r.runtime.nmath.DPQ.EarlyReturn;
+import com.oracle.truffle.r.runtime.nmath.MathFunctions.Function4_2;
+import com.oracle.truffle.r.runtime.nmath.RMath;
+import com.oracle.truffle.r.runtime.nmath.RMathError;
+import com.oracle.truffle.r.runtime.nmath.RMathError.MLError;
+import com.oracle.truffle.r.runtime.nmath.TOMS708;
+
+/*
+ *  Copenhaver, Margaret Diponzio & Holland, Burt S.
+ *  Multiple comparisons of simple effects in
+ *  the two-way analysis of variance with fixed effects.
+ *  Journal of Statistical Computation and Simulation,
+ *  Vol.30, pp.1-15, 1988.
+ *
+ *  Uses the secant method to find critical values.
+ *
+ *  p = confidence level (1 - alpha)
+ *  rr = no. of rows or groups
+ *  cc = no. of columns or treatments
+ *  df = degrees of freedom of error term
+ *
+ *  ir(1) = error flag = 1 if wprob probability > 1
+ *  ir(2) = error flag = 1 if ptukey probability > 1
+ *  ir(3) = error flag = 1 if convergence not reached in 50 iterations
+ *       = 2 if df < 2
+ *
+ *  qtukey = returned critical value
+ *
+ *  If the difference between successive iterates is less than eps,
+ *  the search is terminated
+ */
+
+public final class QTukey implements Function4_2 {
+    private static final double eps = 0.0001;
+    private static final int maxiter = 50;
+
+    private final PTukey ptukey = new PTukey();
+
+    @Override
+    public double evaluate(double pIn, double rr, double cc, double df, boolean lowerTail, boolean logP) {
+        if (Double.isNaN(pIn) || Double.isNaN(rr) || Double.isNaN(cc) || Double.isNaN(df)) {
+            RMathError.error(MLError.DOMAIN, "qtukey");
+            return pIn + rr + cc + df;
+        }
+
+        /* df must be > 1 ; there must be at least two values */
+        if (df < 2 || rr < 1 || cc < 2) {
+            return RMathError.defaultError();
+        }
+
+        try {
+            DPQ.rqp01boundaries(pIn, 0, Double.POSITIVE_INFINITY, lowerTail, logP);
+        } catch (EarlyReturn e) {
+            return e.result;
+        }
+
+        double p = DPQ.rdtqiv(pIn, lowerTail, logP); /* lowerTail,non-log "p" */
+
+        /* Initial value */
+
+        double x0 = qinv(p, cc, df);
+
+        /* Find prob(value < x0) */
+
+        double valx0 = ptukey.evaluate(x0, rr, cc, df, /* LOWER */true, /* LOG_P */false) - p;
+
+        /* Find the second iterate and prob(value < x1). */
+        /* If the first iterate has probability value */
+        /* exceeding p then second iterate is 1 less than */
+        /* first iterate; otherwise it is 1 greater. */
+
+        double x1 = valx0 > 0.0 ? RMath.fmax2(0.0, x0 - 1.0) : x0 + 1.0;
+        double valx1 = ptukey.evaluate(x1, rr, cc, df, /* LOWER */true, /* LOG_P */false) - p;
+
+        /* Find new iterate */
+
+        double ans = 0.;
+        for (int iter = 1; iter < maxiter; iter++) {
+            ans = x1 - ((valx1 * (x1 - x0)) / (valx1 - valx0));
+            valx0 = valx1;
+
+            /* New iterate must be >= 0 */
+
+            x0 = x1;
+            if (ans < 0.0) {
+                ans = 0.0;
+                valx1 = -p;
+            }
+            /* Find prob(value < new iterate) */
+
+            valx1 = ptukey.evaluate(ans, rr, cc, df, /* LOWER */true, /* LOG_P */false) - p;
+            x1 = ans;
+
+            /* If the difference between two successive */
+            /* iterates is less than eps, stop */
+            double xabs = TOMS708.fabs(x1 - x0);
+            if (xabs < eps) {
+                return ans;
+            }
+        }
+
+        /* The process did not converge in 'maxiter' iterations */
+        RMathError.error(MLError.NOCONV, "qtukey");
+        return ans;
+
+    }
+
+    private static final double p0 = 0.322232421088;
+    private static final double q0 = 0.993484626060e-01;
+    private static final double p1 = -1.0;
+    private static final double q1 = 0.588581570495;
+    private static final double p2 = -0.342242088547;
+    private static final double q2 = 0.531103462366;
+    private static final double p3 = -0.204231210125;
+    private static final double q3 = 0.103537752850;
+    private static final double p4 = -0.453642210148e-04;
+    private static final double q4 = 0.38560700634e-02;
+    private static final double c1 = 0.8832;
+    private static final double c2 = 0.2368;
+    private static final double c3 = 1.214;
+    private static final double c4 = 1.208;
+    private static final double c5 = 1.4142;
+    private static final double vmax = 120.0;
+
+    static double qinv(double p, double c, double v) {
+        double ps = 0.5 - 0.5 * p;
+        double yi = Math.sqrt(Math.log(1.0 / (ps * ps)));
+        double t = yi + ((((yi * p4 + p3) * yi + p2) * yi + p1) * yi + p0) / ((((yi * q4 + q3) * yi + q2) * yi + q1) * yi + q0);
+        if (v < vmax) {
+            t += (t * t * t + t) / v / 4.0;
+        }
+        double q = c1 - c2 * t;
+        if (v < vmax) {
+            q += -c3 / v + c4 * t / v;
+        }
+        return t * (q * Math.log(c - 1.0) + c5);
+    }
+}
diff --git a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test
index 1b894cb4b7..10eed45227 100644
--- a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test
+++ b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test
@@ -114624,6 +114624,102 @@ In pt(c(1), -10, ncp = 2, lower.tail = T, log.p = F) : NaNs produced
 Warning message:
 In pt(c(1), -10, ncp = 2, lower.tail = T, log.p = T) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, -10,  5,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, -Inf,  5,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 0,  5,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 1,  5,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10,  5, -10)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10,  5, -Inf)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10,  5, 0)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10,  5, 1)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10,  5, Inf)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10,  5, NaN)
+[1] NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10, -10,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10, -Inf,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10, 0,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10, 1,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10, Inf,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, 10, NaN,  4)
+[1] NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, Inf,  5,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(0, NaN,  5,  4)
+[1] NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(c(-1, 1, 1.9, 3, 5, 10, 15, 20, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 10, 5, 4, lower.tail=F, log.p=F)
+ [1] 1.000000e+00 1.000000e+00 9.934805e-01 8.305223e-01 3.211940e-01
+ [6] 2.527372e-02 4.050344e-03 1.054660e-03 5.155335e-07 1.000000e+00
+[11] 1.000000e+00 1.000000e+00 1.000000e+00 0.000000e+00          NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(c(-1, 1, 1.9, 3, 5, 10, 15, 20, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 10, 5, 4, lower.tail=F, log.p=T)
+ [1]   0.000000000   0.000000000  -0.006540828  -0.185700543  -1.135710025
+ [6]  -3.677990327  -5.508953426  -6.854537201 -14.478063579   0.000000000
+[11]   0.000000000   0.000000000   0.000000000          -Inf           NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(c(-1, 1, 1.9, 3, 5, 10, 15, 20, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 10, 5, 4, lower.tail=T, log.p=F)
+ [1] 0.000000000 0.000000000 0.006519484 0.169477737 0.678806015 0.974726284
+ [7] 0.995949656 0.998945340 0.999999484 0.000000000 0.000000000 0.000000000
+[13] 0.000000000 1.000000000         NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#ptukey(c(-1, 1, 1.9, 3, 5, 10, 15, 20, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 10, 5, 4, lower.tail=T, log.p=T)
+ [1]          -Inf          -Inf -5.032960e+00 -1.775034e+00 -3.874199e-01
+ [6] -2.559858e-02 -4.058569e-03 -1.055216e-03 -5.155336e-07          -Inf
+[11]          -Inf          -Inf          -Inf  0.000000e+00           NaN
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
 #punif(0, -3, -Inf)
 [1] NaN
@@ -116951,6 +117047,164 @@ In qt(log(c(0, 4.2e-79, 0.1, 0.5, 0.7, 1 - 4.2e-79, 1)), -10, ncp = 2,  :
 [1]        -Inf -19.7943519  -2.2815516  -1.0000000  -0.4755995         Inf
 [7]         Inf
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(-0.42e-38, 10, 5, 4)
+[1] NaN
+Warning message:
+In qtukey(-4.2e-39, 10, 5, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(-42, 10, 5, 4)
+[1] NaN
+Warning message:
+In qtukey(-42, 10, 5, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(-Inf, 10, 5, 4)
+[1] NaN
+Warning message:
+In qtukey(-Inf, 10, 5, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, -10,  5,  4)
+[1] NaN
+Warning message:
+In qtukey(0, -10, 5, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, -Inf,  5,  4)
+[1] NaN
+Warning message:
+In qtukey(0, -Inf, 5, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 0,  5,  4)
+[1] NaN
+Warning message:
+In qtukey(0, 0, 5, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 1,  5,  4)
+[1] NaN
+Warning message:
+In qtukey(0, 1, 5, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10,  5, -10)
+[1] NaN
+Warning message:
+In qtukey(0, 10, 5, -10) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10,  5, -Inf)
+[1] NaN
+Warning message:
+In qtukey(0, 10, 5, -Inf) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10,  5, 0)
+[1] NaN
+Warning message:
+In qtukey(0, 10, 5, 0) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10,  5, 1)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10,  5, Inf)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10,  5, NaN)
+[1] NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10, -10,  4)
+[1] NaN
+Warning message:
+In qtukey(0, 10, -10, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10, -Inf,  4)
+[1] NaN
+Warning message:
+In qtukey(0, 10, -Inf, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10, 0,  4)
+[1] NaN
+Warning message:
+In qtukey(0, 10, 0, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10, 1,  4)
+[1] NaN
+Warning message:
+In qtukey(0, 10, 1, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10, Inf,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, 10, NaN,  4)
+[1] NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, Inf,  5,  4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(0, NaN,  5,  4)
+[1] NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(Inf, 10, 5, 4)
+[1] NaN
+Warning message:
+In qtukey(Inf, 10, 5, 4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(NaN, 10, 5, 4)
+[1] NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 10, 5, 4, lower.tail=F, log.p=F)
+[1]      Inf      NaN 7.124296 4.179388 3.465450 0.000000 0.000000
+Warning messages:
+1: In qtukey(c(0, 4.2e-79, 0.1, 0.5, 0.7, 1 - 4.2e-79, 1), 10, 5, 4,  :
+  convergence failed in 'qtukey'
+2: In qtukey(c(0, 4.2e-79, 0.1, 0.5, 0.7, 1 - 4.2e-79, 1), 10, 5, 4,  :
+  NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 10, 5, 4, lower.tail=T, log.p=F)
+[1] 0.0000000 0.6000925       NaN 4.1793878 5.1226189       Inf       Inf
+Warning messages:
+1: In qtukey(c(0, 4.2e-79, 0.1, 0.5, 0.7, 1 - 4.2e-79, 1), 10, 5, 4,  :
+  convergence failed in 'qtukey'
+2: In qtukey(c(0, 4.2e-79, 0.1, 0.5, 0.7, 1 - 4.2e-79, 1), 10, 5, 4,  :
+  NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 10, 5, 4, lower.tail=F, log.p=T)
+[1]      Inf      NaN 7.124296 4.179388 3.465450 0.000000 0.000000
+Warning messages:
+1: In qtukey(log(c(0, 4.2e-79, 0.1, 0.5, 0.7, 1 - 4.2e-79, 1)), 10,  :
+  convergence failed in 'qtukey'
+2: In qtukey(log(c(0, 4.2e-79, 0.1, 0.5, 0.7, 1 - 4.2e-79, 1)), 10,  :
+  NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qtukey(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 10, 5, 4, lower.tail=T, log.p=T)
+[1] 0.0000000 0.6000925       NaN 4.1793878 5.1226189       Inf       Inf
+Warning messages:
+1: In qtukey(log(c(0, 4.2e-79, 0.1, 0.5, 0.7, 1 - 4.2e-79, 1)), 10,  :
+  convergence failed in 'qtukey'
+2: In qtukey(log(c(0, 4.2e-79, 0.1, 0.5, 0.7, 1 - 4.2e-79, 1)), 10,  :
+  NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qunif(-0.42e-38, -3, 3.3)
 [1] NaN
diff --git a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestDistributions.java b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestDistributions.java
index c142d36cf8..ce30d364d5 100644
--- a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestDistributions.java
+++ b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestDistributions.java
@@ -160,13 +160,21 @@ public class TestDistributions extends TestBase {
                     test("3e100, ncp=0.5", withDefaultQ("-30", "-20", "-4", "0.5", "1.3", "2", "3", "4", "10", "100")).
                     test("Inf, ncp=-1", withQuantiles("-10", "-5", "-4", "-3", "-2", "-1", "0", "1.1", "2", "3", "4", "10", "100")).
                     // negative first parameter => error
-                    test("-10, ncp=2", withQuantiles("1"))
+                    test("-10, ncp=2", withQuantiles("1")),
+            distr("tukey").
+                    hasNoDensityFunction().
+                    addErrorParamValues("-10", "0", "1").
+                    test("10, 5, 4", withDefaultQ("-1", "1", "1.9", "3", "5", "10", "15", "20", "100"))
     };
     // @formatter:on
 
     @Test
     public void testDensityFunctions() {
         for (DistrTest testCase : testCases) {
+            if (!testCase.hasDensityFunction) {
+                continue;
+            }
+
             for (ParamsAndQuantiles paramsAndQ : testCase.paramsAndQuantiles) {
                 testDensityFunction("d" + testCase.name, paramsAndQ.params, paramsAndQ.quantiles);
             }
@@ -250,6 +258,7 @@ public class TestDistributions extends TestBase {
         public final String name;
         public final ArrayList<ParamsAndQuantiles> paramsAndQuantiles = new ArrayList<>();
         private int paramsCount = -1;
+        private boolean hasDensityFunction = true;
         /**
          * Set of single R values that are supposed to produce error when used as any of the
          * parameters.
@@ -282,6 +291,11 @@ public class TestDistributions extends TestBase {
             errorParamValues.clear();
             return this;
         }
+
+        public DistrTest hasNoDensityFunction() {
+            hasDensityFunction = false;
+            return this;
+        }
     }
 
     /**
diff --git a/mx.fastr/copyrights/overrides b/mx.fastr/copyrights/overrides
index 5b8f4bcc37..c4d5072ca3 100644
--- a/mx.fastr/copyrights/overrides
+++ b/mx.fastr/copyrights/overrides
@@ -33,6 +33,8 @@ com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/DNorm.
 com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/RBeta.java,gnu_r_gentleman_ihaka.copyright
 com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/RMultinom.java,gnu_r_gentleman_ihaka.copyright
 com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Cauchy.java,gnu_r_ihaka_core.copyright
+com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/QTukey.java,gnu_r_ihaka_core.copyright
+com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/PTukey.java,gnu_r_ihaka_core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Cdist.java,gnu_r_gentleman_ihaka.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/CompleteCases.java,gnu_r_gentleman_ihaka2.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Covcor.java,gnu_r.copyright
-- 
GitLab