Implement nbeta distribution functions in stats

com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNBeta.java

+/*
+ * 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);
+    }
+}

com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNBeta.java

+/*
+ * 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;
+    }
+}

com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNBeta.java

+/*
+ * 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);
+    }
+}

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);
         }
     }
 

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":

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
 

mx.fastr/copyrights/gnu_r_ihaka_core.copyright.star.regex

-/\*\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.*

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