From a71931b73e333ec4145de433b280dcc3b59e6d70 Mon Sep 17 00:00:00 2001
From: stepan <stepan.sindelar@oracle.com>
Date: Wed, 28 Dec 2016 16:10:51 +0100
Subject: [PATCH] Implement nbeta distribution functions in stats

---
 .../truffle/r/library/stats/DNBeta.java       |  94 ++++
 .../truffle/r/library/stats/PNBeta.java       | 126 ++++++
 .../truffle/r/library/stats/QNBeta.java       |  76 ++++
 .../r/library/stats/StatsFunctions.java       | 108 ++++-
 .../foreign/CallAndExternalFunctions.java     |  10 +
 .../truffle/r/test/ExpectedTestOutput.test    | 428 ++++++++++++++++++
 .../test/library/stats/TestDistributions.java |   9 +-
 .../gnu_r_ihaka_core.copyright.star.regex     |   2 +-
 mx.fastr/copyrights/overrides                 |   4 +
 9 files changed, 843 insertions(+), 14 deletions(-)
 create mode 100644 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNBeta.java
 create mode 100644 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNBeta.java
 create mode 100644 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNBeta.java

diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNBeta.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNBeta.java
new file mode 100644
index 0000000000..2a6fe146eb
--- /dev/null
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNBeta.java
@@ -0,0 +1,94 @@
+/*
+ * 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-12, The R Core Team
+ * Copyright (c) 2016, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+package com.oracle.truffle.r.library.stats;
+
+import static com.oracle.truffle.r.library.stats.GammaFunctions.dpoisRaw;
+
+import com.oracle.truffle.r.library.stats.StatsFunctions.Function4_1;
+
+public class DNBeta implements Function4_1 {
+    private static final double eps = 1.e-15;
+    private final DBeta dbeta = new DBeta();
+
+    @Override
+    public double evaluate(double x, double a, double b, double ncp, boolean giveLog) {
+        if (Double.isNaN(x) || Double.isNaN(a) || Double.isNaN(b) || Double.isNaN(ncp)) {
+            return x + a + b + ncp;
+        }
+        if (ncp < 0 || a <= 0 || b <= 0 || !Double.isFinite(a) || !Double.isFinite(b) || !Double.isFinite(ncp)) {
+            return RMathError.defaultError();
+        }
+
+        if (x < 0 || x > 1) {
+            return DPQ.rd0(giveLog);
+        }
+        if (ncp == 0) {
+            return dbeta.evaluate(x, a, b, giveLog);
+        }
+
+        /* New algorithm, starting with *largest* term : */
+        double ncp2 = 0.5 * ncp;
+        double dx2 = ncp2 * x;
+        double d = (dx2 - a - 1) / 2;
+        double capD = d * d + dx2 * (a + b) - a;
+        int kMax;
+        if (capD <= 0) {
+            kMax = 0;
+        } else {
+            capD = Math.ceil(d + Math.sqrt(capD));
+            kMax = (capD > 0) ? (int) capD : 0;
+        }
+
+        /* The starting "middle term" --- first look at it's log scale: */
+        double term = dbeta.evaluate(x, a + kMax, b, /* log = */ true);
+        /* LDOUBLE */double pK = dpoisRaw(kMax, ncp2, true);
+        if (x == 0. || !Double.isFinite(term) || !Double.isFinite(pK)) {
+            /* if term = +Inf */
+            return DPQ.rdexp(pK + term, giveLog);
+        }
+
+        /*
+         * Now if s_k := pK * t_k {here = Math.exp(pK + term)} would underflow, we should rather
+         * scale everything and re-scale at the end:
+         */
+
+        pK += term; /*
+                     * = Math.log(pK) + Math.log(t_k) == Math.log(s_k) -- used at end to rescale
+                     */
+        /* mid = 1 = the rescaled value, instead of mid = Math.exp(pK); */
+
+        /* Now sum from the inside out */
+        /* LDOUBLE */double sum = term = 1. /* = mid term */;
+        /* middle to the left */
+        int k = kMax;
+        while (k > 0 && term > sum * eps) {
+            k--;
+            /* LDOUBLE */double q = /* 1 / r_k = */ (k + 1) * (k + a) / (k + a + b) / dx2;
+            term *= q;
+            sum += term;
+        }
+        /* middle to the right */
+        term = 1.;
+        k = kMax;
+        do {
+            /* LDOUBLE */double q = /* r_{old k} = */ dx2 * (k + a + b) / (k + a) / (k + 1);
+            k++;
+            term *= q;
+            sum += term;
+        } while (term > sum * eps);
+
+        // #ifdef HAVE_LONG_DOUBLE
+        // return DPQ.rdMath.exp((double)(pK + logl(sum)));
+        // #else
+        return DPQ.rdexp(pK + Math.log(sum), giveLog);
+    }
+}
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNBeta.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNBeta.java
new file mode 100644
index 0000000000..58795411a1
--- /dev/null
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNBeta.java
@@ -0,0 +1,126 @@
+/*
+ * 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) 2000-2013, The R Core Team
+ * Copyright (c) 2016, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+package com.oracle.truffle.r.library.stats;
+
+import static com.oracle.truffle.r.library.stats.GammaFunctions.lgammafn;
+
+import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
+import com.oracle.truffle.r.library.stats.RMathError.MLError;
+import com.oracle.truffle.r.library.stats.StatsFunctions.Function4_2;
+import com.oracle.truffle.r.library.stats.TOMS708.Bratio;
+
+public class PNBeta implements Function4_2 {
+    @Override
+    public double evaluate(double x, double a, double b, double ncp, boolean lowerTail, boolean logP) {
+        if (Double.isNaN(x) || Double.isNaN(a) || Double.isNaN(b) || Double.isNaN(ncp)) {
+            return x + a + b + ncp;
+        }
+        try {
+            DPQ.rpbounds01(x, 0., 1., lowerTail, logP);
+        } catch (EarlyReturn e) {
+            return e.result;
+        }
+        return pnbeta2(x, 1 - x, a, b, ncp, lowerTail, logP);
+    }
+
+    private double pnbeta2(double x, double oX, double a, double b, double ncp, boolean lowerTail, boolean logP) {
+        /* LDOUBLE */
+        double ans = pnbetaRaw(x, oX, a, b, ncp);
+
+        /* return DPQ.rdtval(ans), but we want to warn about cancellation here */
+        if (lowerTail) {
+            // #ifdef HAVE_LONG_DOUBLE
+            // return (double) (logP ? logl(ans) : ans);
+            // #else
+            return logP ? Math.log(ans) : ans;
+        } else {
+            if (ans > 1. - 1e-10) {
+                RMathError.error(MLError.PRECISION, "pnbeta");
+            }
+            if (ans > 1.0) {
+                ans = 1.0;
+            } /* Precaution */
+            // #if defined(HAVE_LONG_DOUBLE) && defined(HAVE_LOG1PL)
+            // return (double) (logP ? log1pl(-ans) : (1. - ans));
+            // #else
+            /* include standalone case */
+            return (logP ? RMath.log1p(-ans) : (1. - ans));
+        }
+
+    }
+
+    /*
+     * GnuR: change errmax and itrmax if desired; original (AS 226, R84) had (errmax; itrmax) =
+     * (1e-6; 100)
+     */
+    private static final double errmax = 1.0e-9;
+    private static final int itrmax = 10000; /*
+                                              * GnuR: 100 is not enough for pf(ncp=200) see PR#11277
+                                              */
+
+    double pnbetaRaw(double x, double oX, double a, double b, double ncp) {
+        /* oX == 1 - x but maybe more accurate */
+        if (ncp < 0. || a <= 0. || b <= 0.) {
+            return RMathError.defaultError();
+        }
+
+        if (x < 0. || oX > 1. || (x == 0. && oX == 1.)) {
+            return 0.;
+        }
+        if (x > 1. || oX < 0. || (x == 1. && oX == 0.)) {
+            return 1.;
+        }
+
+        double c = ncp / 2.;
+
+        /* initialize the series */
+        double x0 = Math.floor(RMath.fmax2(c - 7. * Math.sqrt(c), 0.));
+        double a0 = a + x0;
+        double lbeta = lgammafn(a0) + lgammafn(b) - lgammafn(a0 + b);
+
+        /* temp = pbeta_raw(x, a0, b, true, false), but using (x, oX): */
+        double temp = Bratio.bratio(a0, b, x, oX, false).w;
+
+        /* LDOUBLE */double gx = Math.exp(a0 * Math.log(x) + b * (x < .5 ? RMath.log1p(-x) : Math.log(oX)) - lbeta - Math.log(a0));
+        /* LDOUBLE */double q;
+        if (a0 > a) {
+            q = Math.exp(-c + x0 * Math.log(c) - lgammafn(x0 + 1.));
+        } else {
+            q = Math.exp(-c);
+        }
+
+        /* LDOUBLE */double sumq = 1. - q;
+        /* LDOUBLE */double ans = q * temp;
+        /* LDOUBLE */double ax = ans;
+
+        /* recurse over subsequent terms until convergence is achieved */
+        double j = Math.floor(x0); // x0 could be billions, and is in package EnvStats
+        double errbd;
+        do {
+            j++;
+            temp -= gx;
+            gx *= x * (a + b + j - 1.) / (a + j);
+            q *= c / j;
+            sumq -= q;
+            ax = temp * q;
+            ans += ax;
+            errbd = ((temp - gx) * sumq);
+        } while (errbd > errmax && j < itrmax + x0);
+
+        if (errbd > errmax) {
+            RMathError.error(MLError.PRECISION, "pnbeta");
+        }
+        if (j >= itrmax + x0) {
+            RMathError.error(MLError.NOCONV, "pnbeta");
+        }
+        return ans;
+    }
+}
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNBeta.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNBeta.java
new file mode 100644
index 0000000000..ea00bc27cf
--- /dev/null
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNBeta.java
@@ -0,0 +1,76 @@
+/*
+ * 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) 2006, The R Core Team
+ * Copyright (c) 2016, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+package com.oracle.truffle.r.library.stats;
+
+import static com.oracle.truffle.r.library.stats.MathConstants.DBL_EPSILON;
+import static com.oracle.truffle.r.library.stats.MathConstants.DBL_MIN;
+
+import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
+import com.oracle.truffle.r.library.stats.StatsFunctions.Function4_2;
+
+public class QNBeta implements Function4_2 {
+    private static final double accu = 1e-15;
+    private static final double Eps = 1e-14; /* must be > accu */
+
+    private final PNBeta pnbeta = new PNBeta();
+
+    @Override
+    public double evaluate(double p, double a, double b, double ncp, boolean lowerTail, boolean logP) {
+        if (Double.isNaN(p) || Double.isNaN(a) || Double.isNaN(b) || Double.isNaN(ncp)) {
+            return p + a + b + ncp;
+        }
+        if (!Double.isFinite(a)) {
+            return RMathError.defaultError();
+        }
+
+        if (ncp < 0. || a <= 0. || b <= 0.) {
+            return RMathError.defaultError();
+        }
+
+        try {
+            DPQ.rqp01boundaries(p, 0, 1, lowerTail, logP);
+        } catch (EarlyReturn e) {
+            return e.result;
+        }
+
+        p = DPQ.rdtqiv(p, lowerTail, logP);
+
+        /*
+         * Invert pnbeta(.) : 1. finding an upper and lower bound
+         */
+        if (p > 1 - DBL_EPSILON) {
+            return 1.0;
+        }
+        double pp = RMath.fmin2(1 - DBL_EPSILON, p * (1 + Eps));
+        double ux = 0.5;
+        while (ux < 1 - DBL_EPSILON && pnbeta.evaluate(ux, a, b, ncp, true, false) < pp) {
+            ux = 0.5 * (1 + ux);
+        }
+        pp = p * (1 - Eps);
+        double lx = 0.5;
+        while (lx > DBL_MIN && pnbeta.evaluate(lx, a, b, ncp, true, false) > pp) {
+            lx *= 0.5;
+        }
+
+        /* 2. interval (lx,ux) halving : */
+        double nx;
+        do {
+            nx = 0.5 * (lx + ux);
+            if (pnbeta.evaluate(nx, a, b, ncp, true, false) > p) {
+                ux = nx;
+            } else {
+                lx = nx;
+            }
+        } while ((ux - lx) / nx > accu);
+
+        return 0.5 * (ux + lx);
+    }
+}
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/StatsFunctions.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/StatsFunctions.java
index 48c80044ba..9d5484b763 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/StatsFunctions.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/StatsFunctions.java
@@ -21,6 +21,8 @@ import com.oracle.truffle.api.nodes.Node;
 import com.oracle.truffle.api.profiles.BranchProfile;
 import com.oracle.truffle.api.profiles.ConditionProfile;
 import com.oracle.truffle.api.profiles.LoopConditionProfile;
+import com.oracle.truffle.r.library.stats.StatsFunctionsFactory.Function4_1NodeGen;
+import com.oracle.truffle.r.library.stats.StatsFunctionsFactory.Function4_2NodeGen;
 import com.oracle.truffle.r.nodes.attributes.UnaryCopyAttributesNode;
 import com.oracle.truffle.r.nodes.builtin.CastBuilder;
 import com.oracle.truffle.r.nodes.builtin.RExternalBuiltinNode;
@@ -44,7 +46,25 @@ public final class StatsFunctions {
         // private
     }
 
-    public interface Function3_2 {
+    public interface Function4_2 {
+        double evaluate(double a, double b, double c, double d, boolean x, boolean y);
+    }
+
+    public interface Function4_1 extends Function4_2 {
+        @Override
+        default double evaluate(double a, double b, double c, double d, boolean x, boolean y) {
+            return evaluate(a, b, c, d, x);
+        }
+
+        double evaluate(double a, double b, double c, double d, boolean x);
+    }
+
+    public interface Function3_2 extends Function4_2 {
+        @Override
+        default double evaluate(double a, double b, double c, double d, boolean x, boolean y) {
+            return evaluate(a, b, c, x, y);
+        }
+
         double evaluate(double a, double b, double c, boolean x, boolean y);
     }
 
@@ -80,9 +100,11 @@ public final class StatsFunctions {
         final NACheck aCheck = NACheck.create();
         final NACheck bCheck = NACheck.create();
         final NACheck cCheck = NACheck.create();
+        final NACheck dCheck = NACheck.create();
         final ConditionProfile copyAttrsFromA = ConditionProfile.createBinaryProfile();
         final ConditionProfile copyAttrsFromB = ConditionProfile.createBinaryProfile();
         final ConditionProfile copyAttrsFromC = ConditionProfile.createBinaryProfile();
+        final ConditionProfile copyAttrsFromD = ConditionProfile.createBinaryProfile();
         final VectorLengthProfile resultVectorLengthProfile = VectorLengthProfile.create();
         final LoopConditionProfile loopConditionProfile = LoopConditionProfile.createCountingProfile();
 
@@ -91,15 +113,16 @@ public final class StatsFunctions {
         }
     }
 
-    private static RAbstractDoubleVector evaluate3(Node node, Function3_2 function, RAbstractDoubleVector a, RAbstractDoubleVector b, RAbstractDoubleVector c, boolean x, boolean y,
-                    StatFunctionProfiles profiles, UnaryCopyAttributesNode copyAttributesNode) {
+    private static RAbstractDoubleVector evaluate4(Node node, Function4_2 function, RAbstractDoubleVector a, RAbstractDoubleVector b, RAbstractDoubleVector c, RAbstractDoubleVector d, boolean x,
+                    boolean y, StatFunctionProfiles profiles, UnaryCopyAttributesNode copyAttributesNode) {
         int aLength = a.getLength();
         int bLength = b.getLength();
         int cLength = c.getLength();
-        if (aLength == 0 || bLength == 0 || cLength == 0) {
+        int dLength = d.getLength();
+        if (aLength == 0 || bLength == 0 || cLength == 0 || dLength == 0) {
             return RDataFactory.createEmptyDoubleVector();
         }
-        int length = profiles.resultVectorLengthProfile.profile(Math.max(aLength, Math.max(bLength, cLength)));
+        int length = profiles.resultVectorLengthProfile.profile(Math.max(aLength, Math.max(bLength, Math.max(cLength, dLength))));
         RNode.reportWork(node, length);
         double[] result = new double[length];
 
@@ -108,22 +131,24 @@ public final class StatsFunctions {
         profiles.aCheck.enable(a);
         profiles.bCheck.enable(b);
         profiles.cCheck.enable(c);
+        profiles.dCheck.enable(d);
         profiles.loopConditionProfile.profileCounted(length);
         for (int i = 0; profiles.loopConditionProfile.inject(i < length); i++) {
             double aValue = a.getDataAt(i % aLength);
             double bValue = b.getDataAt(i % bLength);
             double cValue = c.getDataAt(i % cLength);
+            double dValue = d.getDataAt(i % dLength);
             double value;
-            if (Double.isNaN(aValue) || Double.isNaN(bValue) || Double.isNaN(cValue)) {
+            if (Double.isNaN(aValue) || Double.isNaN(bValue) || Double.isNaN(cValue) || Double.isNaN(dValue)) {
                 profiles.nan.enter();
-                if (profiles.aCheck.check(aValue) || profiles.bCheck.check(bValue) || profiles.cCheck.check(cValue)) {
+                if (profiles.aCheck.check(aValue) || profiles.bCheck.check(bValue) || profiles.cCheck.check(cValue) || profiles.cCheck.check(dValue)) {
                     value = RRuntime.DOUBLE_NA;
                     complete = false;
                 } else {
                     value = Double.NaN;
                 }
             } else {
-                value = function.evaluate(aValue, bValue, cValue, x, y);
+                value = function.evaluate(aValue, bValue, cValue, dValue, x, y);
                 if (Double.isNaN(value)) {
                     profiles.nan.enter();
                     nans = true;
@@ -143,6 +168,8 @@ public final class StatsFunctions {
             copyAttributesNode.execute(resultVec, b);
         } else if (profiles.copyAttrsFromC.profile(cLength == length)) {
             copyAttributesNode.execute(resultVec, c);
+        } else if (profiles.copyAttrsFromD.profile((dLength == length))) {
+            copyAttributesNode.execute(resultVec, d);
         }
 
         return resultVec;
@@ -168,7 +195,64 @@ public final class StatsFunctions {
         protected RAbstractDoubleVector evaluate(RAbstractDoubleVector a, RAbstractDoubleVector b, RAbstractDoubleVector c, boolean x, boolean y,
                         @Cached("create()") StatFunctionProfiles profiles,
                         @Cached("create()") UnaryCopyAttributesNode copyAttributesNode) {
-            return evaluate3(this, function, a, b, c, x, y, profiles, copyAttributesNode);
+            return evaluate4(this, function, a, b, c, DUMMY_VECTOR, x, y, profiles, copyAttributesNode);
+        }
+    }
+
+    public abstract static class Function4_1Node extends RExternalBuiltinNode.Arg5 {
+        private final Function4_1 function;
+
+        public Function4_1Node(Function4_1 function) {
+            this.function = function;
+        }
+
+        public static Function4_1Node create(Function4_1 function) {
+            return Function4_1NodeGen.create(function);
+        }
+
+        @Override
+        protected void createCasts(CastBuilder casts) {
+            casts.arg(0).asDoubleVector();
+            casts.arg(1).asDoubleVector();
+            casts.arg(2).asDoubleVector();
+            casts.arg(3).asDoubleVector();
+            casts.arg(4).asLogicalVector().findFirst().map(toBoolean());
+        }
+
+        @Specialization
+        protected RAbstractDoubleVector evaluate(RAbstractDoubleVector a, RAbstractDoubleVector b, RAbstractDoubleVector c, RAbstractDoubleVector d, boolean x,
+                        @Cached("create()") StatFunctionProfiles profiles,
+                        @Cached("create()") UnaryCopyAttributesNode copyAttributesNode) {
+            return evaluate4(this, function, a, b, c, d, x, false /* dummy */, profiles, copyAttributesNode);
+        }
+    }
+
+    public abstract static class Function4_2Node extends RExternalBuiltinNode.Arg6 {
+        private final Function4_2 function;
+
+        public Function4_2Node(Function4_2 function) {
+            this.function = function;
+        }
+
+        public static Function4_2Node create(Function4_2 function) {
+            return Function4_2NodeGen.create(function);
+        }
+
+        @Override
+        protected void createCasts(CastBuilder casts) {
+            casts.arg(0).asDoubleVector();
+            casts.arg(1).asDoubleVector();
+            casts.arg(2).asDoubleVector();
+            casts.arg(3).asDoubleVector();
+            casts.arg(4).asLogicalVector().findFirst().map(toBoolean());
+            casts.arg(5).asLogicalVector().findFirst().map(toBoolean());
+        }
+
+        @Specialization
+        protected RAbstractDoubleVector evaluate(RAbstractDoubleVector a, RAbstractDoubleVector b, RAbstractDoubleVector c, RAbstractDoubleVector d, boolean x, boolean y,
+                        @Cached("create()") StatFunctionProfiles profiles,
+                        @Cached("create()") UnaryCopyAttributesNode copyAttributesNode) {
+            return evaluate4(this, function, a, b, c, d, x, y, profiles, copyAttributesNode);
         }
     }
 
@@ -191,7 +275,7 @@ public final class StatsFunctions {
         protected RAbstractDoubleVector evaluate(RAbstractDoubleVector a, RAbstractDoubleVector b, RAbstractDoubleVector c, boolean x,
                         @Cached("create()") StatFunctionProfiles profiles,
                         @Cached("create()") UnaryCopyAttributesNode copyAttributesNode) {
-            return evaluate3(this, function, a, b, c, x, false /* dummy */, profiles, copyAttributesNode);
+            return evaluate4(this, function, a, b, c, DUMMY_VECTOR, x, false /* dummy */, profiles, copyAttributesNode);
         }
     }
 
@@ -213,7 +297,7 @@ public final class StatsFunctions {
         protected RAbstractDoubleVector evaluate(RAbstractDoubleVector a, RAbstractDoubleVector b, boolean x,
                         @Cached("create()") StatFunctionProfiles profiles,
                         @Cached("create()") UnaryCopyAttributesNode copyAttributesNode) {
-            return evaluate3(this, function, a, b, DUMMY_VECTOR, x, false /* dummy */, profiles, copyAttributesNode);
+            return evaluate4(this, function, a, b, DUMMY_VECTOR, DUMMY_VECTOR, x, false /* dummy */, profiles, copyAttributesNode);
         }
     }
 
@@ -236,7 +320,7 @@ public final class StatsFunctions {
         protected RAbstractDoubleVector evaluate(RAbstractDoubleVector a, RAbstractDoubleVector b, boolean x, boolean y,
                         @Cached("create()") StatFunctionProfiles profiles,
                         @Cached("create()") UnaryCopyAttributesNode copyAttributesNode) {
-            return evaluate3(this, function, a, b, DUMMY_VECTOR, x, y, profiles, copyAttributesNode);
+            return evaluate4(this, function, a, b, DUMMY_VECTOR, DUMMY_VECTOR, x, y, profiles, copyAttributesNode);
         }
     }
 
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 2186960df1..78b960ddcc 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
@@ -51,6 +51,7 @@ import com.oracle.truffle.r.library.stats.CompleteCases;
 import com.oracle.truffle.r.library.stats.CovcorNodeGen;
 import com.oracle.truffle.r.library.stats.CutreeNodeGen;
 import com.oracle.truffle.r.library.stats.DBeta;
+import com.oracle.truffle.r.library.stats.DNBeta;
 import com.oracle.truffle.r.library.stats.DNChisq;
 import com.oracle.truffle.r.library.stats.DNorm;
 import com.oracle.truffle.r.library.stats.DPois;
@@ -75,6 +76,7 @@ import com.oracle.truffle.r.library.stats.Logis;
 import com.oracle.truffle.r.library.stats.Logis.DLogis;
 import com.oracle.truffle.r.library.stats.Logis.RLogis;
 import com.oracle.truffle.r.library.stats.PGamma;
+import com.oracle.truffle.r.library.stats.PNBeta;
 import com.oracle.truffle.r.library.stats.PNChisq;
 import com.oracle.truffle.r.library.stats.PPois;
 import com.oracle.truffle.r.library.stats.Pbeta;
@@ -83,6 +85,7 @@ import com.oracle.truffle.r.library.stats.Pf;
 import com.oracle.truffle.r.library.stats.Pnorm;
 import com.oracle.truffle.r.library.stats.Pt;
 import com.oracle.truffle.r.library.stats.QBeta;
+import com.oracle.truffle.r.library.stats.QNBeta;
 import com.oracle.truffle.r.library.stats.QNChisq;
 import com.oracle.truffle.r.library.stats.QPois;
 import com.oracle.truffle.r.library.stats.Qbinom;
@@ -108,6 +111,7 @@ import com.oracle.truffle.r.library.stats.Rt;
 import com.oracle.truffle.r.library.stats.Signrank.RSignrank;
 import com.oracle.truffle.r.library.stats.SplineFunctionsFactory.SplineCoefNodeGen;
 import com.oracle.truffle.r.library.stats.SplineFunctionsFactory.SplineEvalNodeGen;
+import com.oracle.truffle.r.library.stats.StatsFunctions;
 import com.oracle.truffle.r.library.stats.StatsFunctionsFactory;
 import com.oracle.truffle.r.library.stats.Unif.DUnif;
 import com.oracle.truffle.r.library.stats.Unif.PUnif;
@@ -372,6 +376,12 @@ public class CallAndExternalFunctions {
                     return StatsFunctionsFactory.Function2_1NodeGen.create(new DPois());
                 case "dbeta":
                     return StatsFunctionsFactory.Function3_1NodeGen.create(new DBeta());
+                case "dnbeta":
+                    return StatsFunctions.Function4_1Node.create(new DNBeta());
+                case "qnbeta":
+                    return StatsFunctions.Function4_2Node.create(new QNBeta());
+                case "pnbeta":
+                    return StatsFunctions.Function4_2Node.create(new PNBeta());
                 case "dt":
                     return StatsFunctionsFactory.Function2_1NodeGen.create(new Dt());
                 case "rlnorm":
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 194f3d6540..275f7e542b 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
@@ -111561,16 +111561,34 @@ attr(,"is.truffle.object")
 Warning message:
 In dbeta(0, -1, 0.5) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, -4,  15,  0)
+[1] NaN
+Warning message:
+In dbeta(0, -4, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
 #dbeta(0, -Inf,  0.5)
 [1] NaN
 Warning message:
 In dbeta(0, -Inf, 0.5) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, -Inf,  15,  0)
+[1] NaN
+Warning message:
+In dbeta(0, -Inf, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
 #dbeta(0, 0,  0.5)
 [1] Inf
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 0,  15,  0)
+[1] NaN
+Warning message:
+In dbeta(0, 0, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
 #dbeta(0, 0.5, -1)
 [1] NaN
@@ -111595,14 +111613,78 @@ In dbeta(0, 0.5, -Inf) : NaNs produced
 #dbeta(0, 0.5, NaN)
 [1] NaN
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 10,  15, -4)
+[1] NaN
+Warning message:
+In dbeta(0, 10, 15, -4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 10,  15, -Inf)
+[1] NaN
+Warning message:
+In dbeta(0, 10, 15, -Inf) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 10,  15, 0)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 10,  15, Inf)
+[1] NaN
+Warning message:
+In dbeta(0, 10, 15, Inf) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 10,  15, NaN)
+[1] NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 10, -4,  0)
+[1] NaN
+Warning message:
+In dbeta(0, 10, -4, 0) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 10, -Inf,  0)
+[1] NaN
+Warning message:
+In dbeta(0, 10, -Inf, 0) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 10, 0,  0)
+[1] NaN
+Warning message:
+In dbeta(0, 10, 0, 0) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 10, Inf,  0)
+[1] NaN
+Warning message:
+In dbeta(0, 10, Inf, 0) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, 10, NaN,  0)
+[1] NaN
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
 #dbeta(0, Inf,  0.5)
 [1] 0
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, Inf,  15,  0)
+[1] NaN
+Warning message:
+In dbeta(0, Inf, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
 #dbeta(0, NaN,  0.5)
 [1] NaN
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(0, NaN,  15,  0)
+[1] NaN
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
 #dbeta(c(-Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 0.5, 0.5, log=F)
 [1] 0.000000e+00 0.000000e+00          Inf 4.911628e+14 0.000000e+00
@@ -111662,6 +111744,42 @@ In dbeta(0, 0.5, -Inf) : NaNs produced
 #dbeta(c(0.6, 0.1, 42e-33), 6, 3, log=T)
 [1]    0.7372544   -6.5996825 -356.1142283
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(c(10, 15, 100), 7, 11, 0.37e-10, log=F)
+[1] 0 0 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(c(10, 15, 100), 7, 11, 0.37e-10, log=T)
+[1] -Inf -Inf -Inf
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(c(10, 15, 100), 7, 113e11, 1, log=F)
+[1] 0 0 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(c(10, 15, 100), 7, 113e11, 1, log=T)
+[1] -Inf -Inf -Inf
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 10, 15, 0, log=F)
+[1]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
+[6]  0.000000e+00 7.975867e-267  0.000000e+00           NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 10, 15, 0, log=T)
+[1]      -Inf      -Inf      -Inf      -Inf      -Inf      -Inf -612.7138
+[8]      -Inf       NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 7, 13, 3, log=F)
+[1]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
+[6]  0.000000e+00 4.319955e-178  0.000000e+00           NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
+#dbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 7, 13, 3, log=T)
+[1]      -Inf      -Inf      -Inf      -Inf      -Inf      -Inf -408.3969
+[8]      -Inf       NaN
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext#
 #dcauchy(0, -Inf,  -1)
 [1] NaN
@@ -112027,16 +112145,28 @@ In dunif(0, Inf, 3.3) : NaNs produced
 Warning message:
 In pbeta(0, -1, 0.5) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, -4,  15,  0)
+[1] 0
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
 #pbeta(0, -Inf,  0.5)
 [1] NaN
 Warning message:
 In pbeta(0, -Inf, 0.5) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, -Inf,  15,  0)
+[1] 0
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
 #pbeta(0, 0,  0.5)
 [1] 0
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 0,  15,  0)
+[1] 0
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
 #pbeta(0, 0.5, -1)
 [1] NaN
@@ -112061,14 +112191,62 @@ In pbeta(0, 0.5, -Inf) : NaNs produced
 #pbeta(0, 0.5, NaN)
 [1] NaN
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 10,  15, -4)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 10,  15, -Inf)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 10,  15, 0)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 10,  15, Inf)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 10,  15, NaN)
+[1] NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 10, -4,  0)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 10, -Inf,  0)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 10, 0,  0)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 10, Inf,  0)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, 10, NaN,  0)
+[1] NaN
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
 #pbeta(0, Inf,  0.5)
 [1] 0
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, Inf,  15,  0)
+[1] 0
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
 #pbeta(0, NaN,  0.5)
 [1] NaN
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(0, NaN,  15,  0)
+[1] NaN
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
 #pbeta(c(-Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 0.5, 0.5, lower.tail=F, log.p=F)
 [1]   1   1   1   1   0 NaN
@@ -112189,6 +112367,76 @@ In pbeta(0, 0.5, -Inf) : NaNs produced
 #pbeta(c(0.6, 0.1, 42e-33), 6, 3, lower.tail=T, log.p=T)
 [1]   -1.153931  -10.662347 -430.153626
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100), 7, 11, 0.37e-10, lower.tail=F, log.p=F)
+[1] 0 0 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100), 7, 11, 0.37e-10, lower.tail=F, log.p=T)
+[1] -Inf -Inf -Inf
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100), 7, 11, 0.37e-10, lower.tail=T, log.p=F)
+[1] 1 1 1
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100), 7, 11, 0.37e-10, lower.tail=T, log.p=T)
+[1] 0 0 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100), 7, 113e11, 1, lower.tail=F, log.p=F)
+[1] 0 0 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100), 7, 113e11, 1, lower.tail=F, log.p=T)
+[1] -Inf -Inf -Inf
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100), 7, 113e11, 1, lower.tail=T, log.p=F)
+[1] 1 1 1
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100), 7, 113e11, 1, lower.tail=T, log.p=T)
+[1] 0 0 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 10, 15, 0, lower.tail=F, log.p=F)
+[1]   0   0   0   1   1   1   1   0 NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 10, 15, 0, lower.tail=F, log.p=T)
+[1]           -Inf           -Inf           -Inf   0.000000e+00   0.000000e+00
+[6]   0.000000e+00 -3.349864e-298           -Inf            NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 10, 15, 0, lower.tail=T, log.p=F)
+[1]  1.000000e+00  1.000000e+00  1.000000e+00  0.000000e+00  0.000000e+00
+[6]  0.000000e+00 3.349864e-298  1.000000e+00           NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 10, 15, 0, lower.tail=T, log.p=T)
+[1]    0.0000    0.0000    0.0000      -Inf      -Inf      -Inf -684.9614
+[8]    0.0000       NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 7, 13, 3, lower.tail=F, log.p=F)
+[1]   0   0   0   1   1   1   1   0 NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 7, 13, 3, lower.tail=F, log.p=T)
+[1]           -Inf           -Inf           -Inf   0.000000e+00   0.000000e+00
+[6]   0.000000e+00 -2.591973e-209           -Inf            NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 7, 13, 3, lower.tail=T, log.p=F)
+[1]  1.000000e+00  1.000000e+00  1.000000e+00  0.000000e+00  0.000000e+00
+[6]  0.000000e+00 2.591973e-209  1.000000e+00           NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
+#pbeta(c(10, 15, 100, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 7, 13, 3, lower.tail=T, log.p=T)
+[1]    0.0000    0.0000    0.0000      -Inf      -Inf      -Inf -480.2879
+[8]    0.0000       NaN
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext#
 #pcauchy(0, -Inf,  -1)
 [1] NaN
@@ -112698,34 +112946,70 @@ In punif(0, Inf, 3.3) : NaNs produced
 Warning message:
 In qbeta(-4.2e-39, 0.5, 0.5) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(-0.42e-38, 10, 15, 0)
+[1] NaN
+Warning message:
+In qbeta(-4.2e-39, 10, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(-42, 0.5, 0.5)
 [1] NaN
 Warning message:
 In qbeta(-42, 0.5, 0.5) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(-42, 10, 15, 0)
+[1] NaN
+Warning message:
+In qbeta(-42, 10, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(-Inf, 0.5, 0.5)
 [1] NaN
 Warning message:
 In qbeta(-Inf, 0.5, 0.5) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(-Inf, 10, 15, 0)
+[1] NaN
+Warning message:
+In qbeta(-Inf, 10, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(0, -1,  0.5)
 [1] NaN
 Warning message:
 In qbeta(0, -1, 0.5) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, -4,  15,  0)
+[1] NaN
+Warning message:
+In qbeta(0, -4, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(0, -Inf,  0.5)
 [1] NaN
 Warning message:
 In qbeta(0, -Inf, 0.5) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, -Inf,  15,  0)
+[1] NaN
+Warning message:
+In qbeta(0, -Inf, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(0, 0,  0.5)
 [1] 0
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 0,  15,  0)
+[1] NaN
+Warning message:
+In qbeta(0, 0, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(0, 0.5, -1)
 [1] NaN
@@ -112750,24 +113034,94 @@ In qbeta(0, 0.5, -Inf) : NaNs produced
 #qbeta(0, 0.5, NaN)
 [1] NaN
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 10,  15, -4)
+[1] NaN
+Warning message:
+In qbeta(0, 10, 15, -4) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 10,  15, -Inf)
+[1] NaN
+Warning message:
+In qbeta(0, 10, 15, -Inf) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 10,  15, 0)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 10,  15, Inf)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 10,  15, NaN)
+[1] NaN
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 10, -4,  0)
+[1] NaN
+Warning message:
+In qbeta(0, 10, -4, 0) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 10, -Inf,  0)
+[1] NaN
+Warning message:
+In qbeta(0, 10, -Inf, 0) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 10, 0,  0)
+[1] NaN
+Warning message:
+In qbeta(0, 10, 0, 0) : NaNs produced
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 10, Inf,  0)
+[1] 0
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, 10, NaN,  0)
+[1] NaN
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(0, Inf,  0.5)
 [1] 0
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, Inf,  15,  0)
+[1] NaN
+Warning message:
+In qbeta(0, Inf, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(0, NaN,  0.5)
 [1] NaN
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(0, NaN,  15,  0)
+[1] NaN
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(Inf, 0.5, 0.5)
 [1] NaN
 Warning message:
 In qbeta(Inf, 0.5, 0.5) : NaNs produced
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(Inf, 10, 15, 0)
+[1] NaN
+Warning message:
+In qbeta(Inf, 10, 15, 0) : NaNs produced
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(NaN, 0.5, 0.5)
 [1] NaN
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(NaN, 10, 15, 0)
+[1] NaN
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 0, 0, lower.tail=F, log.p=F)
 [1] 1.0 0.0 0.0 0.5 1.0 0.0 0.0
@@ -112803,6 +113157,15 @@ In qbeta(Inf, 0.5, 0.5) : NaNs produced
 [1]  0.000000e+00 4.352496e-157  2.447174e-02  5.000000e-01  7.938926e-01
 [6]  1.000000e+00  1.000000e+00
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 10, 15, 0, lower.tail=F, log.p=F)
+[1] 1.0000000 1.0000000 0.5264118 0.3972924 0.3463781 0.0000000 0.0000000
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 10, 15, 0, lower.tail=T, log.p=F)
+[1] 0.000000e+00 3.412491e-09 2.772130e-01 3.972924e-01 4.497463e-01
+[6] 1.000000e+00 1.000000e+00
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 2, 5, lower.tail=F, log.p=F)
 [1] 1.0000000 1.0000000 0.5103163 0.2644500 0.1818035 0.0000000 0.0000000
@@ -112821,6 +113184,34 @@ In qbeta(Inf, 0.5, 0.5) : NaNs produced
 [1] 0.000000e+00 4.966097e-14 4.617846e-01 6.794810e-01 7.586657e-01
 [6] 1.000000e+00 1.000000e+00
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 7, 11, 0.37e-10, lower.tail=F, log.p=F)
+[1] 1.0000000 1.0000000 0.5373549 0.3846872 0.3251782 0.0000000 0.0000000
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 7, 11, 0.37e-10, lower.tail=T, log.p=F)
+[1] 0.000000e+00 1.551047e-12 2.461368e-01 3.846872e-01 4.466173e-01
+[6] 1.000000e+00 1.000000e+00
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 7, 113e11, 1, lower.tail=F, log.p=F)
+[1] 1.000000e+00 1.000000e+00 9.979557e-13 6.326635e-13 5.134397e-13
+[6] 0.000000e+00 0.000000e+00
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 7, 113e11, 1, lower.tail=T, log.p=F)
+[1] 0.000000e+00 2.042740e-24 3.697562e-13 6.326635e-13 7.690718e-13
+[6] 1.000000e+00 1.000000e+00
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 7, 13, 3, lower.tail=F, log.p=F)
+[1] 1.0000000 1.0000000 0.5361327 0.3904742 0.3326837 0.0000000 0.0000000
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 7, 13, 3, lower.tail=T, log.p=F)
+[1] 0.000000e+00 1.677349e-12 2.547294e-01 3.904742e-01 4.499427e-01
+[6] 1.000000e+00 1.000000e+00
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), Inf, 0, lower.tail=F, log.p=F)
 [1] 1 1 1 1 1 0 0
@@ -112864,6 +113255,15 @@ In qbeta(Inf, 0.5, 0.5) : NaNs produced
 [1]  0.000000e+00 4.352496e-157  2.447174e-02  5.000000e-01  7.938926e-01
 [6]  1.000000e+00  1.000000e+00
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 10, 15, 0, lower.tail=F, log.p=T)
+[1] 1.0000000 1.0000000 0.5264118 0.3972924 0.3463781 0.0000000 0.0000000
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 10, 15, 0, lower.tail=T, log.p=T)
+[1] 0.000000e+00 3.412491e-09 2.772130e-01 3.972924e-01 4.497463e-01
+[6] 1.000000e+00 1.000000e+00
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 2, 5, lower.tail=F, log.p=T)
 [1] 1.0000000 1.0000000 0.5103163 0.2644500 0.1818035 0.0000000 0.0000000
@@ -112882,6 +113282,34 @@ In qbeta(Inf, 0.5, 0.5) : NaNs produced
 [1] 0.000000e+00 4.966097e-14 4.617846e-01 6.794810e-01 7.586657e-01
 [6] 1.000000e+00 1.000000e+00
 
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 7, 11, 0.37e-10, lower.tail=F, log.p=T)
+[1] 1.0000000 1.0000000 0.5373549 0.3846872 0.3251782 0.0000000 0.0000000
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 7, 11, 0.37e-10, lower.tail=T, log.p=T)
+[1] 0.000000e+00 1.551047e-12 2.461368e-01 3.846872e-01 4.466173e-01
+[6] 1.000000e+00 1.000000e+00
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 7, 113e11, 1, lower.tail=F, log.p=T)
+[1] 1.000000e+00 1.000000e+00 9.979557e-13 6.326635e-13 5.134397e-13
+[6] 0.000000e+00 0.000000e+00
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 7, 113e11, 1, lower.tail=T, log.p=T)
+[1] 0.000000e+00 2.042740e-24 3.697562e-13 6.326635e-13 7.690718e-13
+[6] 1.000000e+00 1.000000e+00
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 7, 13, 3, lower.tail=F, log.p=T)
+[1] 1.0000000 1.0000000 0.5361327 0.3904742 0.3326837 0.0000000 0.0000000
+
+##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
+#qbeta(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 7, 13, 3, lower.tail=T, log.p=T)
+[1] 0.000000e+00 1.677349e-12 2.547294e-01 3.904742e-01 4.499427e-01
+[6] 1.000000e+00 1.000000e+00
+
 ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext#
 #qbeta(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), Inf, 0, lower.tail=F, log.p=T)
 [1] 1 1 1 1 1 0 0
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 d7ef66b59b..f745736a92 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
@@ -89,7 +89,14 @@ public class TestDistributions extends TestBase {
                     test("1, 1", withDefaultQ("0.5", "2")).
                     test("420, 4", withQuantiles("0.42e-10", "100", "13e10", "11e111")).
                     test("0.13e-8, 1", withQuantiles("0.42e-10", "100", "13e10", "11e111")).
-                    test("1, 0.13e-8", withQuantiles("0.42e-10", "100", "13e10", "11e111"))
+                    test("1, 0.13e-8", withQuantiles("0.42e-10", "100", "13e10", "11e111")),
+            // tests of nbeta, which is called in beta when third param is not missing
+            distr("beta").
+                    addErrorParamValues("-4", "0").
+                    test("10, 15, 0", withDefaultQ("10", "15", "100")).
+                    test("7, 13, 3", withDefaultQ("10", "15", "100")).
+                    test("7, 11, 0.37e-10", withQuantiles("10", "15", "100")).
+                    test("7, 113e11, 1", withQuantiles("10", "15", "100"))
     };
     // @formatter:on
 
diff --git a/mx.fastr/copyrights/gnu_r_ihaka_core.copyright.star.regex b/mx.fastr/copyrights/gnu_r_ihaka_core.copyright.star.regex
index 1bc53ecfd3..8ede89dd35 100644
--- a/mx.fastr/copyrights/gnu_r_ihaka_core.copyright.star.regex
+++ b/mx.fastr/copyrights/gnu_r_ihaka_core.copyright.star.regex
@@ -1 +1 @@
-/\*\n \* This material is distributed under the GNU General Public License\n \* Version 2. You may review the terms of this license at\n \* http://www.gnu.org/licenses/gpl-2.0.html\n \*\n \* Copyright \(C\) 1998 Ross Ihaka\n \* Copyright \(c\) (?:[1-2][09][0-9][0-9]--?)?[1-2][09][0-9][0-9], The R Core Team\n \* Copyright \(c\) (?:(20[0-9][0-9]), )?(20[0-9][0-9]), Oracle and/or its affiliates\n \*\n \* All rights reserved.\n \*/\n.*
+/\*\n \* This material is distributed under the GNU General Public License\n \* Version 2. You may review the terms of this license at\n \* http://www.gnu.org/licenses/gpl-2.0.html\n \*\n \* Copyright \(C\) 1998 Ross Ihaka\n \* Copyright \(c\) (?:[1-2][09][0-9][0-9]--?)?(?:[1-2][09])?[0-9][0-9], The R Core Team\n \* Copyright \(c\) (?:(20[0-9][0-9]), )?(20[0-9][0-9]), Oracle and/or its affiliates\n \*\n \* All rights reserved.\n \*/\n.*
diff --git a/mx.fastr/copyrights/overrides b/mx.fastr/copyrights/overrides
index 6b838fa31c..a4760a69f7 100644
--- a/mx.fastr/copyrights/overrides
+++ b/mx.fastr/copyrights/overrides
@@ -49,6 +49,9 @@ com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pt.java,gnu_
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RHyper.java,gnu_r_ihaka.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathInit.java,gnu_r.core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbeta.java,gnu_r.core.copyright
+com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNBeta.java,gnu_r.core.copyright
+com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNBeta.java,gnu_r.core.copyright
+com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNBeta.java,gnu_r_ihaka_core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/LBeta.java,gnu_r_ihaka.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbinom.java,gnu_r_ihaka.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pf.java,gnu_r_ihaka_core.copyright
@@ -82,6 +85,7 @@ com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RGamma.java,
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PGamma.java,gnu_r_welinder.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RNbinomMu.java,gnu_r_ihaka_core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Logis.java,gnu_r_ihaka_core.copyright
+com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pnf.java,gnu_r_ihaka_core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rf.java,gnu_r_ihaka_core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Exp.java,gnu_r_ihaka_core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Geom.java,gnu_r_ihaka_core.copyright
-- 
GitLab