diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java
index 57479bf6173ead3e861672edda459a8b8618e7af..f3286dd935c434ed1e9e09711e3d5b0350fce863 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java
@@ -16,137 +16,108 @@ import com.oracle.truffle.api.nodes.ControlFlowException;
 import com.oracle.truffle.r.runtime.RError;
 import com.oracle.truffle.r.runtime.RError.Message;
 
-// transcribed from dpq.h
-
+/**
+ * Contains macros transcribed from dpq.h. Naming convention is all lowercase and remove underscores
+ * from GnuR name. Some macros change control flow by explicitly returning, this is handled by
+ * throwing {@link EarlyReturn} exception that encapsulates the desired return value.
+ */
 public final class DPQ {
     private DPQ() {
-        // private
+        // only static methods
+    }
+
+    public static final class EarlyReturn extends ControlFlowException {
+        private static final long serialVersionUID = 1182697355931636213L;
+        public final double result;
+
+        private EarlyReturn(double result) {
+            this.result = result;
+        }
     }
 
     // R >= 3.1.0: # define R_nonint(x) (fabs((x) - R_forceint(x)) > 1e-7)
-    // Note: if true should be followed by "return d0(logP)"
-    // Consider using dNointCheck instead
+    // Note: if true should be followed by "return d0(logP)", consider using nointCheck instead
     public static boolean nonint(double x) {
         return Math.abs(x - Math.round(x)) > 1e-7 * Math.max(1., Math.abs(x));
     }
 
     // R_D__0
-    public static double d0(boolean logP) {
+    public static double rd0(boolean logP) {
         return logP ? Double.NEGATIVE_INFINITY : 0.;
     }
 
     // R_D__1
-    public static double d1(boolean logP) {
+    public static double rd1(boolean logP) {
         return logP ? 0. : 1.;
     }
 
     // R_DT_0
-    public static double dt0(boolean logP, boolean lowerTail) {
-        return lowerTail ? d0(logP) : d1(logP);
+    public static double rdt0(boolean lowerTail, boolean logP) {
+        return lowerTail ? rd0(logP) : rd1(logP);
     }
 
     // R_D_log
-    public static double dLog(double p, boolean logP) {
+    public static double rdlog(double p, boolean logP) {
         return logP ? p : Math.log(p);
     }
 
-    // R_DT_1
-    public static double dt1(boolean logP, boolean lowerTail) {
-        return lowerTail ? d1(logP) : d0(logP);
+    public static double rdtlog(double p, boolean lowerTail, boolean logp) {
+        return lowerTail ? rdlog(p, logp) : rdlexp(p, logp);
     }
 
-    /* Use 0.5 - p + 0.5 to perhaps gain 1 bit of accuracy */
+    // R_DT_1
+    public static double rdt1(boolean lowerTail, boolean logP) {
+        return lowerTail ? rd1(logP) : rd0(logP);
+    }
 
     // R_D_Lval
-    public static double dLval(boolean lowerTail, double p) {
+    public static double rdlval(double p, boolean lowerTail) {
         return lowerTail ? p : 0.5 - p + 0.5;
     }
 
-    public static double dCval(double p, boolean lowerTail) {
+    public static double rdcval(double p, boolean lowerTail) {
         return lowerTail ? 0.5 - p + 0.5 : p; /* 1 - p */
     }
 
-    //
-    public static double dVal(double x, boolean logP) {
+    public static double dval(double x, boolean logP) {
         return logP ? Math.log(x) : x; /* x in pF(x,..) */
     }
 
-    public static double dExp(double x, boolean logP) {
+    public static double rdexp(double x, boolean logP) {
         return logP ? x : Math.exp(x); /* exp(x) */
     }
 
-    /* log(1-exp(x)): R_D_LExp(x) == (log1p(- R_D_qIv(x))) but even more stable: */
-    // #define R_D_LExp(x) (log_p ? R_Log1_Exp(x) : log1p(-x))
-    public static double dLExp(double x, boolean logP) {
-        return (logP ? log1Exp(x, logP) : StatsUtil.log1p(-x));
+    // R_D_LExp
+    public static double rdlexp(double x, boolean logP) {
+        return (logP ? rlog1exp(x) : RMath.log1p(-x));
     }
 
-    // #define R_Log1_Exp(x) ((x) > -M_LN2 ? log(-expm1(x)) : log1p(-exp(x)))
-    public static double log1Exp(double x, boolean logP) {
-        return ((x) > -M_LN2 ? Math.log(-StatsUtil.expm1(x)) : StatsUtil.log1p(-Math.exp(x)));
+    public static double rdfexp(double f, double x, boolean giveLog) {
+        return giveLog ? -0.5 * Math.log(f) + x : Math.exp(x) / Math.sqrt(f);
     }
 
-    // #define R_D_log(p) (log_p ? (p) : log(p)) /* log(p) */
-    public static double dClog(double p, boolean logP) {
-        return logP ? StatsUtil.log1p(-p) : 0.5 - p + 0.5; /* [log](1-p) */
-    }
-
-    // #define R_DT_Clog(p) (lower_tail? R_D_LExp(p): R_D_log(p))/* log(1-p) in qF*/
-    public static double dtCLog(double p, boolean lowerTail, boolean logP) {
-        return lowerTail ? dLExp(p, logP) : dLog(p, logP);
+    // R_Log1_Exp
+    public static double rlog1exp(double x) {
+        return ((x) > -M_LN2 ? Math.log(-RMath.expm1(x)) : RMath.log1p(-Math.exp(x)));
     }
 
-    //
-    // // log(1 - exp(x)) in more stable form than log1p(- R_D_qIv(x)) :
-    // #define R_Log1_Exp(x) ((x) > -M_LN2 ? log(-expm1(x)) : log1p(-exp(x)))
-    //
-    // #define R_D_LExp(x) (log_p ? R_Log1_Exp(x) : log1p(-x))
-    //
-    // #define R_DT_val(x) (lower_tail ? R_D_val(x) : R_D_Clog(x))
-    public static double dtCval(double x, boolean lowerTail, boolean logP) {
-        return lowerTail ? dClog(x, logP) : dVal(x, logP);
+    // R_DT_Clog
+    public static double rdtclog(double p, boolean lowerTail, boolean logP) {
+        return lowerTail ? rdlexp(p, logP) : rdlog(p, logP);
     }
 
-    //
     // R_DT_qIv
-    public static double dtQIv(double p, boolean lowerTail, boolean logP) {
-        return logP ? lowerTail ? Math.exp(p) : -Math.expm1(p) : dLval(lowerTail, p);
-    }
-
-    // /*#define R_DT_CIv(p) R_D_Cval(R_D_qIv(p)) * 1 - p in qF */
-
-    public static double dtCIv(double p, boolean lowerTail, boolean logP) {
-        return logP ? lowerTail ? -Math.expm1(p) : Math.exp(p) : dCval(p, lowerTail);
+    public static double rdtqiv(double p, boolean lowerTail, boolean logP) {
+        return logP ? lowerTail ? Math.exp(p) : -Math.expm1(p) : rdlval(p, lowerTail);
     }
 
-    //
-    // #define R_DT_exp(x) R_D_exp(R_D_Lval(x)) /* exp(x) */
-    // #define R_DT_Cexp(x) R_D_exp(R_D_Cval(x)) /* exp(1 - x) */
-    //
-    // #define R_DT_log(p) (lower_tail? R_D_log(p) : R_D_LExp(p))/* log(p) in qF */
-    // #define R_DT_Clog(p) (lower_tail? R_D_LExp(p): R_D_log(p))/* log(1-p) in qF*/
-    // #define R_DT_Log(p) (lower_tail? (p) : R_Log1_Exp(p))
-
-    public static final class EarlyReturn extends ControlFlowException {
-        private static final long serialVersionUID = 1182697355931636213L;
-        public final double result;
-
-        private EarlyReturn(double result) {
-            this.result = result;
-        }
+    // R_DT_CIv
+    public static double rdtciv(double p, boolean lowerTail, boolean logP) {
+        return logP ? lowerTail ? -Math.expm1(p) : Math.exp(p) : rdcval(p, lowerTail);
     }
 
-    /*
-     * Do the boundaries exactly for q*() functions : Often _LEFT_ = ML_NEGINF , and very often
-     * _RIGHT_ = ML_POSINF;
-     *
-     * R_Q_P01_boundaries(p, _LEFT_, _RIGHT_) :<==>
-     *
-     * R_Q_P01_check(p); if (p == R_DT_0) return _LEFT_ ; if (p == R_DT_1) return _RIGHT_;
-     *
-     * the following implementation should be more efficient (less tests):
-     */
-    public static void qP01Boundaries(double p, double left, double right, boolean lowerTail, boolean logP) throws EarlyReturn {
+    // R_Q_P01_boundaries
+    public static void rqp01boundaries(double p, double left, double right, boolean lowerTail, boolean logP) throws EarlyReturn {
         if (logP) {
             if (p > 0) {
                 throw new EarlyReturn(Double.NaN);
@@ -171,20 +142,21 @@ public final class DPQ {
         }
     }
 
-    // #define R_Q_P01_check(p) \
-    // if ((log_p && p > 0) || \
-    // (!log_p && (p < 0 || p > 1)) ) \
-    // ML_ERR_return_NAN
-    public static void qQP01Check(double p, boolean logP) throws EarlyReturn {
+    // R_Q_P01_check
+    public static void rqp01check(double p, boolean logP) throws EarlyReturn {
         if ((logP && p > 0) || (!logP && (p < 0 || p > 1))) {
-            throw new EarlyReturn(StatsUtil.mlError());
+            throw new EarlyReturn(RMath.mlError());
         }
     }
 
-    /* [neg]ative or [non int]eger : */
-    public static boolean dNegInonint(double x) {
-        return x < 0 || nonint(x);
-    }
+    // Unimplemented macros:
+    //
+    // #define R_DT_exp(x) R_D_exp(R_D_Lval(x)) /* exp(x) */
+    // #define R_DT_Cexp(x) R_D_exp(R_D_Cval(x)) /* exp(1 - x) */
+    //
+    // #define R_DT_log(p) (lower_tail? R_D_log(p) : R_D_LExp(p))/* log(p) in qF */
+    // #define R_DT_Clog(p) (lower_tail? R_D_LExp(p): R_D_log(p))/* log(1-p) in qF*/
+    // #define R_DT_Log(p) (lower_tail? (p) : R_Log1_Exp(p))
 
     // FastR helpers:
 
@@ -192,10 +164,10 @@ public final class DPQ {
         RError.warning(RError.SHOW_CALLER, Message.NON_INTEGER_N, varName, x);
     }
 
-    public static void dNonintCheck(double x, boolean giveLog) throws EarlyReturn {
+    public static void nonintCheck(double x, boolean giveLog) throws EarlyReturn {
         if (nonint(x)) {
             nointCheckWarning(x, "x");
-            throw new EarlyReturn(d0(giveLog));
+            throw new EarlyReturn(rd0(giveLog));
         }
     }
 }
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPois.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPois.java
index 7c5086d02224a59aaea4d6895e77dc52a14740a7..b01d2a4d75603a4b862849e498ef24a41ce0e167 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPois.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPois.java
@@ -25,16 +25,16 @@ public final class DPois implements Function2_1 {
             return x + lambda;
         }
         if (lambda < 0) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
 
         try {
-            DPQ.dNonintCheck(x, giveLog);
+            DPQ.nonintCheck(x, giveLog);
         } catch (EarlyReturn e) {
             return e.result;
         }
         if (x < 0 || !Double.isFinite(x)) {
-            return DPQ.d0(giveLog);
+            return DPQ.rd0(giveLog);
         }
 
         return dpoisRaw(forceint(x), lambda, giveLog);
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Dbinom.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Dbinom.java
index d59f1cde171977849191627f6f5377d6b48d7cca..509f543c81aeb163e1eb86ccf6725ea40a88f95d 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Dbinom.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Dbinom.java
@@ -28,25 +28,25 @@ public final class Dbinom implements StatsFunctions.Function3_1 {
     public static double dbinomRaw(double x, double n, double p, double q, boolean giveLog) {
 
         if (p == 0) {
-            return ((x == 0) ? DPQ.d1(giveLog) : DPQ.d0(giveLog));
+            return ((x == 0) ? DPQ.rd1(giveLog) : DPQ.rd0(giveLog));
         }
         if (q == 0) {
-            return ((x == n) ? DPQ.d1(giveLog) : DPQ.d0(giveLog));
+            return ((x == n) ? DPQ.rd1(giveLog) : DPQ.rd0(giveLog));
         }
 
         if (x == 0) {
             if (n == 0) {
-                return DPQ.d1(giveLog);
+                return DPQ.rd1(giveLog);
             }
             double lc = (p < 0.1) ? -GammaFunctions.bd0(n, n * q) - n * p : n * Math.log(q);
-            return DPQ.dExp(lc, giveLog);
+            return DPQ.rdexp(lc, giveLog);
         }
         if (x == n) {
             double lc = (q < 0.1) ? -GammaFunctions.bd0(n, n * p) - n * q : n * Math.log(p);
-            return DPQ.dExp(lc, giveLog);
+            return DPQ.rdexp(lc, giveLog);
         }
         if (x < 0 || x > n) {
-            return DPQ.d0(giveLog);
+            return DPQ.rd0(giveLog);
         }
 
         /*
@@ -62,7 +62,7 @@ public final class Dbinom implements StatsFunctions.Function3_1 {
          */
         double lf = MathConstants.M_LN_2PI + Math.log(x) + Math.log1p(-x / n);
 
-        return DPQ.dExp(lc - 0.5 * lf, giveLog);
+        return DPQ.rdexp(lc - 0.5 * lf, giveLog);
     }
 
     @Override
@@ -73,19 +73,19 @@ public final class Dbinom implements StatsFunctions.Function3_1 {
             return x + n + p;
         }
 
-        if (p < 0 || p > 1 || DPQ.dNegInonint(n)) {
+        if (p < 0 || p > 1 || n < 0 || DPQ.nonint(n)) {
             nanProfile.enter();
             return Double.NaN;
         }
 
         try {
-            DPQ.dNonintCheck(x, giveLog);
+            DPQ.nonintCheck(x, giveLog);
         } catch (EarlyReturn e) {
             return e.result;
         }
 
         if (x < 0 || !Double.isFinite(x)) {
-            return DPQ.d0(giveLog);
+            return DPQ.rd0(giveLog);
         }
 
         return dbinomRaw(MathConstants.forceint(x), MathConstants.forceint(n), p, 1 - p, giveLog);
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Dt.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Dt.java
index 71c723c5535b3b77a3f2300419e3b318e2027d39..62245897573ddf9e784b89fd9b97acd1d21fa249 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Dt.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Dt.java
@@ -29,11 +29,11 @@ public final class Dt implements Function2_1 {
         }
 
         if (n <= 0) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
 
         if (!Double.isFinite(x)) {
-            return DPQ.d0(giveLog);
+            return DPQ.rd0(giveLog);
         }
         if (!Double.isFinite(n)) {
             return dnorm(x, 0., 1., giveLog);
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Exp.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Exp.java
index 95c5b68a0042bace4ebc2f772bcd03310853dbdd..05c81882607860c2964fcd497ca817555579a8b3 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Exp.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Exp.java
@@ -27,11 +27,11 @@ public class Exp {
             }
 
             if (scale <= 0.0) {
-                return StatsUtil.mlError();
+                return RMath.mlError();
             }
 
             if (x < 0.) {
-                return DPQ.d0(giveLog);
+                return DPQ.rd0(giveLog);
             }
             return (giveLog ? (-x / scale) - Math.log(scale) : Math.exp(-x / scale) / scale);
         }
@@ -41,7 +41,7 @@ public class Exp {
         @Override
         public double evaluate(double scale, RandomNumberProvider rand) {
             if (!Double.isFinite(scale) || scale <= 0.0) {
-                return scale == 0. ? 0. : StatsUtil.mlError();
+                return scale == 0. ? 0. : RMath.mlError();
             }
             return scale * rand.expRand();
         }
@@ -54,16 +54,16 @@ public class Exp {
                 return x + scale;
             }
             if (scale < 0) {
-                return StatsUtil.mlError();
+                return RMath.mlError();
             }
 
             if (x <= 0.) {
-                return DPQ.dt0(logP, lowerTail);
+                return DPQ.rdt0(lowerTail, logP);
             }
 
             /* same as weibull( shape = 1): */
             x = -(x / scale);
-            return lowerTail ? (logP ? DPQ.log1Exp(x, logP) : -StatsUtil.expm1(x)) : DPQ.dExp(x, logP);
+            return lowerTail ? (logP ? DPQ.rlog1exp(x) : -RMath.expm1(x)) : DPQ.rdexp(x, logP);
         }
     }
 
@@ -75,20 +75,20 @@ public class Exp {
             }
 
             if (scale < 0) {
-                return StatsUtil.mlError();
+                return RMath.mlError();
             }
 
             try {
-                DPQ.qQP01Check(p, logP);
+                DPQ.rqp01check(p, logP);
             } catch (EarlyReturn e) {
                 return e.result;
             }
 
-            if (p == DPQ.dt0(logP, lowerTail)) {
+            if (p == DPQ.rdt0(lowerTail, logP)) {
                 return 0;
             }
 
-            return -scale * DPQ.dtCLog(p, lowerTail, logP);
+            return -scale * DPQ.rdtclog(p, lowerTail, logP);
 
         }
     }
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java
index e04341a0ed35cb7847ac00ebe5c48d2e1144cecf..a00d9241c07bd37df1dc9914a9ffa83f45c1878e 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java
@@ -18,29 +18,25 @@
  */
 package com.oracle.truffle.r.library.stats;
 
-import static com.oracle.truffle.r.library.stats.StatsUtil.DBLEPSILON;
-import static com.oracle.truffle.r.library.stats.StatsUtil.M_1_SQRT_2PI;
-import static com.oracle.truffle.r.library.stats.StatsUtil.M_2PI;
-import static com.oracle.truffle.r.library.stats.StatsUtil.M_LN2;
-import static com.oracle.truffle.r.library.stats.StatsUtil.M_SQRT_32;
-import static com.oracle.truffle.r.library.stats.StatsUtil.chebyshevEval;
-import static com.oracle.truffle.r.library.stats.StatsUtil.expm1;
-import static com.oracle.truffle.r.library.stats.StatsUtil.fmax2;
-import static com.oracle.truffle.r.library.stats.StatsUtil.log1p;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rd0;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rd1;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rdexp;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rdfexp;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rdt0;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rdt1;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rdtclog;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rdtlog;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rdtqiv;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rlog1exp;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rqp01check;
+import static com.oracle.truffle.r.library.stats.DPQ.rd0;
+import static com.oracle.truffle.r.library.stats.DPQ.rd1;
+import static com.oracle.truffle.r.library.stats.DPQ.rdexp;
+import static com.oracle.truffle.r.library.stats.DPQ.rdfexp;
+import static com.oracle.truffle.r.library.stats.DPQ.rdt0;
+import static com.oracle.truffle.r.library.stats.DPQ.rdt1;
+import static com.oracle.truffle.r.library.stats.DPQ.rdtclog;
+import static com.oracle.truffle.r.library.stats.DPQ.rdtqiv;
+import static com.oracle.truffle.r.library.stats.DPQ.rlog1exp;
+import static com.oracle.truffle.r.library.stats.MathConstants.DBL_EPSILON;
+import static com.oracle.truffle.r.library.stats.MathConstants.M_1_SQRT_2PI;
+import static com.oracle.truffle.r.library.stats.MathConstants.M_2PI;
+import static com.oracle.truffle.r.library.stats.MathConstants.M_LN2;
+import static com.oracle.truffle.r.library.stats.MathConstants.M_SQRT_32;
+import static com.oracle.truffle.r.library.stats.RMath.fmax2;
 
 import com.oracle.truffle.api.CompilerDirectives.CompilationFinal;
 import com.oracle.truffle.api.CompilerDirectives.TruffleBoundary;
+import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
 import com.oracle.truffle.r.library.stats.StatsFunctions.Function3_1;
 import com.oracle.truffle.r.library.stats.StatsFunctions.Function3_2;
 import com.oracle.truffle.r.runtime.RError;
@@ -186,7 +182,7 @@ public abstract class GammaFunctions {
             /* allow to underflow below */
         } else if (x < xbig) {
             tmp = 10 / x;
-            return chebyshevEval(tmp * tmp * 2 - 1, ALGMCS, nalgm) / x;
+            return RMath.chebyshevEval(tmp * tmp * 2 - 1, ALGMCS, nalgm) / x;
         }
         return 1 / (x * 12);
     }
@@ -257,7 +253,7 @@ public abstract class GammaFunctions {
             }
             y = x - n; /* n = floor(x) ==> y in [ 0, 1 ) */
             --n;
-            value = chebyshevEval(y * 2 - 1, GAMCS, ngam) + .9375;
+            value = RMath.chebyshevEval(y * 2 - 1, GAMCS, ngam) + .9375;
             if (n == 0) {
                 return value; /* x = 1.dddd = 1+y */
             }
@@ -461,9 +457,10 @@ public abstract class GammaFunctions {
             return p + nu;
         }
 
-        if (rqp01check(p, logp)) {
-            // TODO ML_ERR_return_NAN
-            return Double.NaN;
+        try {
+            DPQ.rqp01check(p, logp);
+        } catch (EarlyReturn e) {
+            return e.result;
         }
 
         if (nu <= 0) {
@@ -474,7 +471,7 @@ public abstract class GammaFunctions {
         alpha = 0.5 * nu; /* = [pq]gamma() shape */
         c = alpha - 1;
 
-        if (nu < (-1.24) * (p1 = rdtlog(p, lowerTail, logp))) { /* for small chi-squared */
+        if (nu < (-1.24) * (p1 = DPQ.rdtlog(p, lowerTail, logp))) { /* for small chi-squared */
             /*
              * log(alpha) + g = log(alpha) + log(gamma(alpha)) = = log(alpha*gamma(alpha)) =
              * lgamma(alpha+1) suffers from catastrophic cancellation when alpha << 1
@@ -777,7 +774,7 @@ public abstract class GammaFunctions {
     /* Accurate calculation of log(1+x)-x, particularly for small x. */
     private static double log1pmx(double x) {
         if (x > 1 || x < minLog1Value) {
-            return log1p(x) - x;
+            return RMath.log1p(x) - x;
         } else { /*
                   * -.791 <= x <= 1 -- expand in [x/(2+x)]^2 =: y : log(1+x) - x = x/(2+x) * [ 2 * y
                   * * S(y) - x], with --------------------------------------------- S(y) = 1/3 + y/5
@@ -839,7 +836,7 @@ public abstract class GammaFunctions {
     } /* lgamma1p */
 
     /* If |x| > |k| * M_cutoff, then log[ exp(-x) * k^x ] =~= -x */
-    private static final double M_cutoff = M_LN2 * Double.MAX_EXPONENT / DBLEPSILON;
+    private static final double M_cutoff = M_LN2 * Double.MAX_EXPONENT / DBL_EPSILON;
 
     /*
      * dpois_wrap (x_P_1, lambda, g_log) == dpois (x_P_1 - 1, lambda, g_log) := exp(-L) L^k /
@@ -879,10 +876,10 @@ public abstract class GammaFunctions {
             c *= -x / n;
             term = c / (alph + n);
             sum += term;
-        } while (Math.abs(term) > DBLEPSILON * Math.abs(sum));
+        } while (Math.abs(term) > DBL_EPSILON * Math.abs(sum));
 
         if (lowerTail) {
-            double f1 = logp ? log1p(sum) : 1 + sum;
+            double f1 = logp ? RMath.log1p(sum) : 1 + sum;
             double f2;
             if (alph > 1) {
                 f2 = dpoisRaw(alph, x, logp);
@@ -896,10 +893,10 @@ public abstract class GammaFunctions {
         } else {
             double lf2 = alph * Math.log(x) - lgamma1p(alph);
             if (logp) {
-                return rlog1exp(log1p(sum) + lf2);
+                return rlog1exp(RMath.log1p(sum) + lf2);
             } else {
                 double f1m1 = sum;
-                double f2m1 = expm1(lf2);
+                double f2m1 = RMath.expm1(lf2);
                 return -(f1m1 + f2m1 + f1m1 * f2m1);
             }
         }
@@ -914,7 +911,7 @@ public abstract class GammaFunctions {
             localY++;
             term *= x / localY;
             sum += term;
-        } while (term > sum * DBLEPSILON);
+        } while (term > sum * DBL_EPSILON);
 
         /*
          * sum = \sum_{n=1}^ oo x^n / (y*(y+1)*...*(y+n-1)) = \sum_{n=0}^ oo x^(n+1) /
@@ -956,7 +953,7 @@ public abstract class GammaFunctions {
 
         f0 = y / d;
         /* Needed, e.g. for pgamma(10^c(100,295), shape= 1.1, log=TRUE): */
-        if (Math.abs(y - 1) < Math.abs(d) * DBLEPSILON) { /* includes y < d = Inf */
+        if (Math.abs(y - 1) < Math.abs(d) * DBL_EPSILON) { /* includes y < d = Inf */
             return f0;
         }
 
@@ -1007,7 +1004,7 @@ public abstract class GammaFunctions {
             if (b2 != 0) {
                 f = a2 / b2;
                 /* convergence check: relative; "absolute" for very small f : */
-                if (Math.abs(f - of) <= DBLEPSILON * fmax2(f0, Math.abs(f))) {
+                if (Math.abs(f - of) <= DBL_EPSILON * fmax2(f0, Math.abs(f))) {
                     return f;
                 }
                 of = f;
@@ -1023,7 +1020,7 @@ public abstract class GammaFunctions {
         double term = 1;
         double sum = 0;
 
-        while (localY >= 1 && term > sum * DBLEPSILON) {
+        while (localY >= 1 && term > sum * DBL_EPSILON) {
             term *= localY / lambda;
             sum += term;
             localY--;
@@ -1084,7 +1081,7 @@ public abstract class GammaFunctions {
                 term *= -i / x2;
                 sum += term;
                 i += 2;
-            } while (Math.abs(term) > DBLEPSILON * sum);
+            } while (Math.abs(term) > DBL_EPSILON * sum);
 
             return 1 / sum;
         } else {
@@ -1162,7 +1159,7 @@ public abstract class GammaFunctions {
 
         if (logp) {
             double ndOverP = dpnorm(s2pt, !lowerTail, np);
-            return np + log1p(f * ndOverP);
+            return np + RMath.log1p(f * ndOverP);
         } else {
             double nd = dnorm(s2pt, 0., 1., logp);
             return np + f * nd;
@@ -1198,7 +1195,7 @@ public abstract class GammaFunctions {
             double sum;
             double d = dpoisWrap(alph, x, logp);
             if (alph < 1) {
-                if (x * DBLEPSILON > 1 - alph) {
+                if (x * DBL_EPSILON > 1 - alph) {
                     sum = rd1(logp);
                 } else {
                     double f = pdLowerCf(alph, x - (alph - 1)) * x / alph;
@@ -1207,7 +1204,7 @@ public abstract class GammaFunctions {
                 }
             } else {
                 sum = pdLowerSeries(x, alph - 1); /* = (alph-1)/x + o((alph-1)/x) */
-                sum = logp ? log1p(sum) : 1 + sum;
+                sum = logp ? RMath.log1p(sum) : 1 + sum;
             }
             if (!lowerTail) {
                 res = logp ? sum + d : sum * d;
@@ -1222,7 +1219,7 @@ public abstract class GammaFunctions {
          * We lose a fair amount of accuracy to underflow in the cases where the final result is
          * very close to DBL_MIN. In those cases, simply redo via log space.
          */
-        if (!logp && res < Double.MIN_VALUE / DBLEPSILON) {
+        if (!logp && res < Double.MIN_VALUE / DBL_EPSILON) {
             /* with(.Machine, double.xmin / double.eps) #|-> 1.002084e-292 */
             return Math.exp(pgammaRaw(x, alph, lowerTail, true));
         } else {
@@ -1439,7 +1436,7 @@ public abstract class GammaFunctions {
         if (logp) {
             cum[0] = (-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp);
             if ((lower && x > 0.) || (upper && x <= 0.)) {
-                ccum[0] = log1p(-Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp);
+                ccum[0] = RMath.log1p(-Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp);
             }
         } else {
             cum[0] = Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp;
@@ -1480,7 +1477,7 @@ public abstract class GammaFunctions {
         }
 
         /* Consider changing these : */
-        eps = DBLEPSILON * 0.5;
+        eps = DBL_EPSILON * 0.5;
 
         /* i_tail in {0,1,2} =^= {lower, upper, both} */
         lower = iTail != 1;
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Geom.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Geom.java
index 41954c08979efe0ff8b0958be98e4643d9827508..ca8331b6c4f278d88f6fca1b4bce1bb4ded1a96e 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Geom.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Geom.java
@@ -39,11 +39,11 @@ public class Geom {
     public static final class QGeom implements Function2_2 {
         public double evaluate(double p, double prob, boolean lowerTail, boolean logP) {
             if (prob <= 0 || prob > 1) {
-                return StatsUtil.mlError();
+                return RMath.mlError();
             }
 
             try {
-                DPQ.qP01Boundaries(p, 0, Double.POSITIVE_INFINITY, lowerTail, logP);
+                DPQ.rqp01boundaries(p, 0, Double.POSITIVE_INFINITY, lowerTail, logP);
             } catch (EarlyReturn e) {
                 return e.result;
             }
@@ -56,7 +56,7 @@ public class Geom {
                 return 0;
             }
             /* add a fuzz to ensure left continuity, but value must be >= 0 */
-            return StatsUtil.fmax2(0, Math.ceil(DPQ.dtCLog(p, lowerTail, logP) / StatsUtil.log1p(-prob) - 1 - 1e-12));
+            return RMath.fmax2(0, Math.ceil(DPQ.rdtclog(p, lowerTail, logP) / RMath.log1p(-prob) - 1 - 1e-12));
         }
 
     }
@@ -68,17 +68,17 @@ public class Geom {
                 return x + p;
             }
             if (p <= 0 || p > 1) {
-                return StatsUtil.mlError();
+                return RMath.mlError();
             }
 
             try {
-                DPQ.dNonintCheck(x, giveLog);
+                DPQ.nonintCheck(x, giveLog);
             } catch (EarlyReturn e) {
                 return e.result;
             }
 
             if (x < 0 || !Double.isFinite(x) || p == 0) {
-                return DPQ.d0(giveLog);
+                return DPQ.rd0(giveLog);
             }
             /* prob = (1-p)^x, stable for small p */
             double prob = Dbinom.dbinomRaw(0., forceint(x), p, 1 - p, giveLog);
@@ -90,7 +90,7 @@ public class Geom {
         @Override
         public double evaluate(double p, RandomNumberProvider rand) {
             if (!Double.isFinite(p) || p <= 0 || p > 1) {
-                return StatsUtil.mlError();
+                return RMath.mlError();
             }
             return RPois.rpois(rand.expRand() * ((1 - p) / p), rand);
         }
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/LBeta.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/LBeta.java
index 6cc0de5a33d977f8dec7d56bf66456d975a15ee3..8f997c3f586187da864b6cb32346d7ea03826c01 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/LBeta.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/LBeta.java
@@ -38,7 +38,7 @@ public final class LBeta {
 
         /* both arguments must be >= 0 */
         if (p < 0) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         } else if (p == 0) {
             return Double.POSITIVE_INFINITY;
         } else if (!Double.isFinite(q)) { /* q == +Inf */
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathConstants.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathConstants.java
index 1fbefb2a995f10460e01aef039b2ad94135b163c..0885fbeecf503bb056f1c289b3c2096f5ef3e5a3 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathConstants.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathConstants.java
@@ -3,9 +3,10 @@
  * Version 2. You may review the terms of this license at
  * http://www.gnu.org/licenses/gpl-2.0.html
  *
- * Copyright (c) 1995-2012, The R Core Team
- * Copyright (c) 2003, The R Foundation
- * Copyright (c) 2016, 2016, Oracle and/or its affiliates
+ * Copyright (C) 1998 Ross Ihaka
+ * Copyright (c) 1998--2012, The R Core Team
+ * Copyright (c) 2004, The R Foundation
+ * Copyright (c) 2013, 2016, Oracle and/or its affiliates
  *
  * All rights reserved.
  */
@@ -55,7 +56,7 @@ public final class MathConstants {
     // 1/sqrt(2)
     public static final double M_SQRT1_2 = 0.707106781186547524400844362105;
 
-    /* R-Specific Constants */
+    /* R-Specific Constants from dpq.h and Rmath.h and others */
     // sqrt(3)
     public static final double M_SQRT_3 = 1.732050807568877293527446341506;
     // sqrt(32)
@@ -77,6 +78,12 @@ public final class MathConstants {
     // log(sqrt(pi/2)) == log(pi/2)/2
     public static final double M_LN_SQRT_PId2 = 0.225791352644727432363097614947;
 
+    public static final double DBL_MANT_DIG = 53;
+
+    public static final int DBL_MAX_EXP = 1024;
+
+    public static final int DBL_MIN_EXP = -1021;
+
     public static final double DBL_EPSILON = Math.ulp(1.0);
 
     /**
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbeta.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbeta.java
index 4c950e6d35727e6fc3511907fc7f621500f399af..d7bfc2f1a7c7e07017967aa46dfef307478a4375 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbeta.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbeta.java
@@ -32,19 +32,19 @@ public abstract class Pbeta extends RExternalBuiltinNode.Arg5 {
             }
             if (a == 0 || a / b == 0) {
                 // point mass 1 at 0 ==> P(X <= x) = 1, all x > 0
-                return DPQ.dt1(logProb, lowerTail);
+                return DPQ.rdt1(lowerTail, logProb);
             }
             if (b == 0 || b / a == 0) {
                 // point mass 1 at 1 ==> P(X <= x) = 0, all x < 1
-                return DPQ.dt0(logProb, lowerTail);
+                return DPQ.rdt0(lowerTail, logProb);
             }
 
             // else, remaining case: a = b = Inf : point mass 1 at 1/2
             if (x < 0.5) {
-                return DPQ.dt0(logProb, lowerTail);
+                return DPQ.rdt0(lowerTail, logProb);
             }
             // else, x >= 0.5 :
-            return DPQ.dt1(logProb, lowerTail);
+            return DPQ.rdt1(lowerTail, logProb);
         }
         // Now: 0 < a < Inf; 0 < b < Inf
 
@@ -77,10 +77,10 @@ public abstract class Pbeta extends RExternalBuiltinNode.Arg5 {
         // allowing a==0 and b==0 <==> treat as one- or two-point mass
 
         if (x <= 0) {
-            return DPQ.dt0(logP, lowerTail);
+            return DPQ.rdt0(lowerTail, logP);
         }
         if (x >= 1) {
-            return DPQ.dt1(logP, lowerTail);
+            return DPQ.rdt1(lowerTail, logP);
         }
 
         return pbetaRaw(x, a, b, lowerTail, logP);
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbinom.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbinom.java
index 1891d8f9107a65b8ed0d79278f88373d722b6eb6..70fcce282db80d30395effba2881a35c155aafcb 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbinom.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbinom.java
@@ -34,7 +34,7 @@ public final class Pbinom implements StatsFunctions.Function3_2 {
         if (DPQ.nonint(size)) {
             nanProfile.enter();
             DPQ.nointCheckWarning(size, "n");
-            return DPQ.d0(logP);
+            return DPQ.rd0(logP);
         }
         size = Math.round(size);
         /* PR#8560: n=0 is a valid value */
@@ -44,11 +44,11 @@ public final class Pbinom implements StatsFunctions.Function3_2 {
         }
 
         if (q < 0) {
-            return DPQ.dt0(logP, lowerTail);
+            return DPQ.rdt0(lowerTail, logP);
         }
         q = Math.floor(q + 1e-7);
         if (size <= q) {
-            return DPQ.dt1(logP, lowerTail);
+            return DPQ.rdt1(lowerTail, logP);
         }
         return Pbeta.pbeta(prob, q + 1, size - q, !lowerTail, logP, nanProfile);
     }
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pf.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pf.java
index 51cf4f7e06d1f2833ce943b59dc2296041260add..f595cf1cc339e26326e79c2d9ebe5d3984e182cf 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pf.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pf.java
@@ -11,7 +11,8 @@
  */
 package com.oracle.truffle.r.library.stats;
 
-import static com.oracle.truffle.r.library.stats.StatsUtil.*;
+import static com.oracle.truffle.r.library.stats.DPQ.rdt0;
+import static com.oracle.truffle.r.library.stats.DPQ.rdt1;
 
 import com.oracle.truffle.api.profiles.BranchProfile;
 
@@ -45,7 +46,7 @@ public final class Pf implements StatsFunctions.Function3_2 {
                     return rdt0(lowerTail, logP);
                 }
                 if (x == 1) {
-                    return logP ? -M_LN2 : 0.5;
+                    return logP ? -MathConstants.M_LN2 : 0.5;
                 }
                 if (x > 1) {
                     return rdt1(lowerTail, logP);
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pnorm.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pnorm.java
index bae5a546448f141bb8a75e4f012792e19b7fe32b..a8c29ddde26e7c17a1aba7d83c9d29887136e321 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pnorm.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pnorm.java
@@ -41,11 +41,11 @@ public final class Pnorm implements StatsFunctions.Function3_2 {
                 return Double.NaN;
             }
             /* sigma = 0 : */
-            return (x < mu) ? DPQ.d0(logP) : DPQ.d1(logP);
+            return (x < mu) ? DPQ.rd0(logP) : DPQ.rd1(logP);
         }
         double p = (x - mu) / sigma;
         if (!Double.isFinite(p)) {
-            return (x < mu) ? DPQ.d0(logP) : DPQ.d1(logP);
+            return (x < mu) ? DPQ.rd0(logP) : DPQ.rd1(logP);
         }
 
         PnormBoth pnormBoth = new PnormBoth(p);
@@ -204,11 +204,11 @@ public final class Pnorm implements StatsFunctions.Function3_2 {
                 swapTail(x, lower);
             } else { /* large x such that probs are 0 or 1 */
                 if (x > 0) {
-                    cum = DPQ.d1(logP);
-                    ccum = DPQ.d0(logP);
+                    cum = DPQ.rd1(logP);
+                    ccum = DPQ.rd0(logP);
                 } else {
-                    cum = DPQ.d0(logP);
-                    ccum = DPQ.d1(logP);
+                    cum = DPQ.rd0(logP);
+                    ccum = DPQ.rd1(logP);
                 }
             }
 
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QHyper.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QHyper.java
index d3e23b0e73276ac48c384b17a182ec639278afcb..c3a43d3259eba623afcabfc64774ea901a893bbf 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QHyper.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QHyper.java
@@ -13,9 +13,9 @@ 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.forceint;
-import static com.oracle.truffle.r.library.stats.StatsUtil.fmax2;
-import static com.oracle.truffle.r.library.stats.StatsUtil.fmin2;
-import static com.oracle.truffle.r.library.stats.StatsUtil.lfastchoose;
+import static com.oracle.truffle.r.library.stats.RMath.fmax2;
+import static com.oracle.truffle.r.library.stats.RMath.fmin2;
+import static com.oracle.truffle.r.library.stats.RMath.lfastchoose;
 
 import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
 
@@ -34,7 +34,7 @@ public final class QHyper {
             return p + nr + nb + n;
         }
         if (!Double.isFinite(p) || !Double.isFinite(nr) || !Double.isFinite(nb) || !Double.isFinite(n)) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
 
         nr = forceint(nr);
@@ -42,7 +42,7 @@ public final class QHyper {
         capN = nr + nb;
         n = forceint(n);
         if (nr < 0 || nb < 0 || n < 0 || n > capN) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
 
         /*
@@ -54,7 +54,7 @@ public final class QHyper {
         xend = fmin2(n, nr);
 
         try {
-            DPQ.qP01Boundaries(p, xstart, xend, lowerTail, logP);
+            DPQ.rqp01boundaries(p, xstart, xend, lowerTail, logP);
         } catch (EarlyReturn ex) {
             return ex.result;
         }
@@ -75,7 +75,7 @@ public final class QHyper {
         nb -= xb;
 
         if (!lowerTail || logP) {
-            p = DPQ.dtQIv(p, lowerTail, logP);
+            p = DPQ.rdtqiv(p, lowerTail, logP);
         }
         p *= 1 - 1000 * DBL_EPSILON; /* was 64, but failed on FreeBSD sometimes */
         sum = smallN ? term : Math.exp(term);
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Qbinom.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Qbinom.java
index 6bab507fdc59415010679410aaf6d0cdabfcb6c9..32e4f1df57df5f8f87b110d4cdef6dff2ae41402 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Qbinom.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Qbinom.java
@@ -126,7 +126,7 @@ public final class Qbinom implements StatsFunctions.Function3_2 {
          * [cancellation for p ~= 1, etc]:
          */
         if (!lowerTail || logProb) {
-            p = DPQ.dtQIv(p, lowerTail, logProb); /* need check again (cancellation!): */
+            p = DPQ.rdtqiv(p, lowerTail, logProb); /* need check again (cancellation!): */
             if (p == 0) {
                 return 0;
             }
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Qnorm.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Qnorm.java
index 73dc865c65ef8d64be9ea18e57f867475d5e3a9d..42efb456233e808e5f0119665e833637675143b1 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Qnorm.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Qnorm.java
@@ -31,7 +31,7 @@ public final class Qnorm implements StatsFunctions.Function3_2 {
             return p + mu + sigma;
         }
         try {
-            DPQ.qP01Boundaries(p, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, lowerTail, logP);
+            DPQ.rqp01boundaries(p, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, lowerTail, logP);
         } catch (EarlyReturn early) {
             return early.result;
         }
@@ -44,7 +44,7 @@ public final class Qnorm implements StatsFunctions.Function3_2 {
             return mu;
         }
 
-        double p2 = DPQ.dtQIv(p, lowerTail, logP); /* real lower_tail prob. p */
+        double p2 = DPQ.rdtqiv(p, lowerTail, logP); /* real lower_tail prob. p */
         double q = p2 - 0.5;
 
         debugPrintf("qnorm(p=%10.7g, m=%g, s=%g, l.t.= %d, log= %d): q = %g\n", p, mu, sigma, lowerTail, logP, q);
@@ -76,7 +76,7 @@ public final class Qnorm implements StatsFunctions.Function3_2 {
             /* r = min(p, 1-p) < 0.075 */
             double r;
             if (q > 0) {
-                r = DPQ.dtCIv(p, lowerTail, logP); /* 1-p */
+                r = DPQ.rdtciv(p, lowerTail, logP); /* 1-p */
             } else {
                 r = p2; /* = R_DT_Iv(p) ^= p */
             }
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RBeta.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RBeta.java
index 9ba8e64a12ac14fa2575067e12d4dc4b878b8eec..d62ab89e7e74c7d8962352525c10588f022e4c2d 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RBeta.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RBeta.java
@@ -12,10 +12,10 @@
  */
 package com.oracle.truffle.r.library.stats;
 
+import static com.oracle.truffle.r.library.stats.MathConstants.DBL_MAX_EXP;
 import static com.oracle.truffle.r.library.stats.MathConstants.M_LN2;
-import static com.oracle.truffle.r.library.stats.StatsUtil.DBL_MAX_EXP;
-import static com.oracle.truffle.r.library.stats.StatsUtil.fmax2;
-import static com.oracle.truffle.r.library.stats.StatsUtil.fmin2;
+import static com.oracle.truffle.r.library.stats.RMath.fmax2;
+import static com.oracle.truffle.r.library.stats.RMath.fmin2;
 
 import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandFunction2_Double;
 import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandomNumberProvider;
@@ -27,7 +27,7 @@ public final class RBeta implements RandFunction2_Double {
     @Override
     public double evaluate(double aa, double bb, RandomNumberProvider rand) {
         if (Double.isNaN(aa) || Double.isNaN(bb) || aa < 0. || bb < 0.) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
         if (!Double.isFinite(aa) && !Double.isFinite(bb)) { // a = b = Inf : all mass at 1/2
             return 0.5;
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RCauchy.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RCauchy.java
index accf46e7c246e8914ab9075568984edbe7194507..6792699dccf791bc5636ed24cf401d6882ba782e 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RCauchy.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RCauchy.java
@@ -20,7 +20,7 @@ public final class RCauchy implements RandFunction2_Double {
     @Override
     public double evaluate(double location, double scale, RandomNumberProvider rand) {
         if (Double.isNaN(location) || !Double.isFinite(scale) || scale < 0) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
         if (scale == 0. || !Double.isFinite(location)) {
             return location;
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RChisq.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RChisq.java
index 3b4c21ac102e287b7f0c21c0fdd01a6dd7eff473..67ba496753caf017b672e0003654a11f1eb4904d 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RChisq.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RChisq.java
@@ -17,7 +17,7 @@ import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandomNumberPr
 public final class RChisq implements RandFunction1_Double {
     public static double rchisq(double df, RandomNumberProvider rand) {
         if (!Double.isFinite(df) || df < 0.0) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
         return new RGamma().evaluate(df / 2.0, 2.0, rand);
     }
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RGamma.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RGamma.java
index 0e8581b5035469b7af5bc26d512b224447bf5dfd..eea9b5aea8027a913db009a0f96260fa22327fed 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RGamma.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RGamma.java
@@ -66,13 +66,13 @@ public class RGamma implements RandFunction2_Double {
         double retVal;
 
         if (Double.isNaN(a) || Double.isNaN(scale)) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
         if (a <= 0.0 || scale <= 0.0) {
             if (scale == 0. || a == 0.) {
                 return 0.;
             }
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
         if (!Double.isFinite(a) || !Double.isFinite(scale)) {
             return Double.POSITIVE_INFINITY;
@@ -192,7 +192,7 @@ public class RGamma implements RandFunction2_Double {
                 /* Step 11: hat acceptance (h) */
                 /* (if q not positive go to step 8) */
                 if (q > 0.0) {
-                    w = StatsUtil.expm1(q);
+                    w = RMath.expm1(q);
                     /* ^^^^^ original code had approximation with rel.err < 2e-7 */
                     /* if t is rejected sample again at step 8 */
                     if (c * fabs(u) <= w * Math.exp(e - 0.5 * t * t)) {
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RHyper.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RHyper.java
index 32d9c6afee1cdf69754660e0d751057aefd13389..7dba16badee221972e01340629a88e35ffe0a1b1 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RHyper.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RHyper.java
@@ -100,7 +100,7 @@ public final class RHyper implements RandFunction3_Double {
         /* check parameter validity */
 
         if (!Double.isFinite(nn1in) || !Double.isFinite(nn2in) || !Double.isFinite(kkin)) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
 
         nn1in = forceint(nn1in);
@@ -108,7 +108,7 @@ public final class RHyper implements RandFunction3_Double {
         kkin = forceint(kkin);
 
         if (nn1in < 0 || nn2in < 0 || kkin < 0 || kkin > nn1in + nn2in) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
         if (nn1in >= Integer.MAX_VALUE || nn2in >= Integer.MAX_VALUE || kkin >= Integer.MAX_VALUE) {
             /*
@@ -230,7 +230,7 @@ public final class RHyper implements RandFunction3_Double {
                 nUv++;
                 if (nUv >= 10000) {
                     RError.warning(RError.SHOW_CALLER, Message.GENERIC, String.format("rhyper() branch III: giving up after %d rejections", nUv));
-                    return StatsUtil.mlError();
+                    return RMath.mlError();
                 }
 
                 if (u < p1) { /* rectangular region */
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RLogis.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RLogis.java
index 007b812f841e234635bf6880d12728f879add390..961dafe00d10a34156926a5b54379f5510e45ebc 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RLogis.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RLogis.java
@@ -18,7 +18,7 @@ public final class RLogis implements RandFunction2_Double {
     @Override
     public double evaluate(double location, double scale, RandomNumberProvider rand) {
         if (Double.isNaN(location) || !Double.isFinite(scale)) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
 
         if (scale == 0. || !Double.isFinite(location)) {
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/StatsUtil.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java
similarity index 66%
rename from com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/StatsUtil.java
rename to com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java
index 7a5348110e6dcddb2b05578952f40a334ac2f354..04e41edbf8cfb1e208b56f5ed48e6592256e29eb 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/StatsUtil.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java
@@ -18,10 +18,12 @@ import com.oracle.truffle.api.CompilerDirectives.CompilationFinal;
 import com.oracle.truffle.api.CompilerDirectives.TruffleBoundary;
 
 /**
- * Auxiliary functions and constants from GNU R. These are originally found in the source files
- * mentioned in the code.
+ * Encapsulates functions to be found in Rmath.h or in nmath directory in GnuR except for random
+ * distribution related functions, which usually have their own files.
+ *
+ * @see DPQ
  */
-public class StatsUtil {
+public class RMath {
 
     /**
      * corresponds to macro {@code ML_ERR_return_NAN} in GnuR.
@@ -30,97 +32,6 @@ public class StatsUtil {
         return Double.NaN;
     }
 
-    public static final double DBLEPSILON = 2.2204460492503131e-16;
-
-    @TruffleBoundary
-    private static void fail(String message) {
-        throw new RuntimeException(message);
-    }
-
-    // Constants and auxiliary functions originate from dpq.h and Rmath.h.
-
-    public static final double M_LN2 = 0.693147180559945309417232121458;
-
-    public static final double M_PI = 3.141592653589793238462643383280;
-
-    public static final double M_2PI = 6.283185307179586476925286766559;
-
-    public static final double M_1_SQRT_2PI = 0.398942280401432677939946059934;
-
-    public static final double M_SQRT_32 = 5.656854249492380195206754896838;
-
-    public static final double M_LOG10_2 = 0.301029995663981195213738894724;
-
-    public static final double DBL_MANT_DIG = 53;
-
-    public static final int DBL_MAX_EXP = 1024;
-
-    public static final int DBL_MIN_EXP = -1021;
-
-    public static double rdtlog(double p, boolean lowerTail, boolean logp) {
-        return lowerTail ? rdlog(p, logp) : rdlexp(p, logp);
-    }
-
-    public static double rdlog(double p, boolean logp) {
-        return logp ? p : Math.log(p);
-    }
-
-    public static double rdlexp(double x, boolean logp) {
-        return logp ? rlog1exp(x) : log1p(-x);
-    }
-
-    public static double rlog1exp(double x) {
-        return x > -M_LN2 ? Math.log(-expm1(x)) : log1p(-Math.exp(x));
-    }
-
-    public static double rdtclog(double p, boolean lowerTail, boolean logp) {
-        return lowerTail ? rdlexp(p, logp) : rdlog(p, logp);
-    }
-
-    public static double rdtqiv(double p, boolean lowerTail, boolean logp) {
-        return logp ? (lowerTail ? Math.exp(p) : -expm1(p)) : rdlval(p, lowerTail);
-    }
-
-    public static double rdtciv(double p, boolean lowerTail, boolean logp) {
-        return logp ? (lowerTail ? -expm1(p) : Math.exp(p)) : rdcval(p, lowerTail);
-    }
-
-    public static double rdlval(double p, boolean lowerTail) {
-        return lowerTail ? p : (0.5 - (p) + 0.5);
-    }
-
-    public static double rdcval(double p, boolean lowerTail) {
-        return lowerTail ? (0.5 - p + 0.5) : p;
-    }
-
-    public static double rd0(boolean logp) {
-        return logp ? Double.NEGATIVE_INFINITY : 0.;
-    }
-
-    public static double rd1(boolean logp) {
-        return logp ? 0. : 1.;
-    }
-
-    public static double rdt0(boolean lowerTail, boolean logp) {
-        return lowerTail ? rd0(logp) : rd1(logp);
-    }
-
-    public static double rdt1(boolean lowerTail, boolean logp) {
-        return lowerTail ? rd1(logp) : rd0(logp);
-    }
-
-    public static boolean rqp01check(double p, boolean logp) {
-        return (logp && p > 0) || (!logp && (p < 0 || p > 1));
-    }
-
-    public static double rdexp(double x, boolean logp) {
-        return logp ? x : Math.exp(x);
-    }
-
-    public static double rdfexp(double f, double x, boolean giveLog) {
-        return giveLog ? -0.5 * Math.log(f) + x : Math.exp(x) / Math.sqrt(f);
-    }
-
     public static double lfastchoose(double n, double k) {
         return -Math.log(n + 1.) - lbeta(n - k + 1., k + 1.);
     }
@@ -158,7 +69,7 @@ public class StatsUtil {
         double y;
         double a = Math.abs(x);
 
-        if (a < DBLEPSILON) {
+        if (a < MathConstants.DBL_EPSILON) {
             return x;
         }
         if (a > 0.697) {
@@ -173,7 +84,7 @@ public class StatsUtil {
         }
         /* Newton step for solving log(1 + y) = x for y : */
         /* WARNING: does not work for y ~ -1: bug in 1.5.0 */
-        y -= (1 + y) * (log1p(y) - x);
+        y -= (1 + y) * (RMath.log1p(y) - x);
         return y;
     }
 
@@ -216,7 +127,7 @@ public class StatsUtil {
             /*
              * Improve on speed (only); again give result accurate to IEEE double precision:
              */
-            if (Math.abs(x) < .5 * DBLEPSILON) {
+            if (Math.abs(x) < .5 * MathConstants.DBL_EPSILON) {
                 return x;
             }
 
@@ -263,4 +174,9 @@ public class StatsUtil {
         }
         return (b0 - b2) * 0.5;
     }
+
+    @TruffleBoundary
+    private static void fail(String message) {
+        throw new RuntimeException(message);
+    }
 }
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RNbinomMu.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RNbinomMu.java
index bac788a88f092635d86a66edb0f39d218c08f813..c095616fa14e46565ba778fbf3356a6f257ddc66 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RNbinomMu.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RNbinomMu.java
@@ -20,7 +20,7 @@ public final class RNbinomMu implements RandFunction2_Double {
     @Override
     public double evaluate(double size, double mu, RandomNumberProvider rand) {
         if (!Double.isFinite(mu) || Double.isNaN(size) || size <= 0 || mu < 0) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
         if (!Double.isFinite(size)) {
             size = Double.MAX_VALUE / 2.;
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RNchisq.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RNchisq.java
index da55025e1ad1d072053cd34ce4dab98b98260a8e..5e22bda3f5faf7e3fd655ed596183d3df56e6c07 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RNchisq.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RNchisq.java
@@ -22,7 +22,7 @@ public final class RNchisq implements RandFunction2_Double {
     @Override
     public double evaluate(double df, double lambda, RandomNumberProvider rand) {
         if (!Double.isFinite(df) || !Double.isFinite(lambda) || df < 0. || lambda < 0.) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
 
         if (lambda == 0.) {
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RPois.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RPois.java
index ffcc999eace1ce6645b72423f3a7ee81f7cbb191..2638fa34de84ce1797d21b9e75acdb4798ab2ba5 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RPois.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RPois.java
@@ -80,7 +80,7 @@ public final class RPois implements RandFunction1_Double {
         boolean newBigMu = false;
 
         if (!Double.isFinite(mu) || mu < 0) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
 
         if (mu <= 0.) {
@@ -219,7 +219,7 @@ public final class RPois implements RandFunction1_Double {
                  * sample t from the laplace 'hat' (if t <= -0.6744 then pk < fk for all mu >= 10.)
                  */
                 u = 2 * rand.unifRand() - 1.;
-                t = 1.8 + StatsUtil.fsign(e, u);
+                t = 1.8 + RMath.fsign(e, u);
             }
 
             if (t > -0.6744 || gotoStepF) {
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RWeibull.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RWeibull.java
index 6e1fae2faf3211f4352d2dfae59e4a0a95c3c0d7..3686f96591d283d05747f0a49deb6dbb5b4c1b93 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RWeibull.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RWeibull.java
@@ -18,7 +18,7 @@ public final class RWeibull implements RandFunction2_Double {
     @Override
     public double evaluate(double shape, double scale, RandomNumberProvider rand) {
         if (!Double.isFinite(shape) || !Double.isFinite(scale) || shape <= 0. || scale <= 0.) {
-            return scale == 0. ? 0. : StatsUtil.mlError();
+            return scale == 0. ? 0. : RMath.mlError();
         }
 
         return scale * Math.pow(-Math.log(rand.unifRand()), 1.0 / shape);
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Random2.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Random2.java
index e81ed75ed1f546923c664cb7c2d214909e9f0471..729f0d0833ed777de42e12fc552206ce7d12f58d 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Random2.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Random2.java
@@ -11,8 +11,8 @@
  */
 package com.oracle.truffle.r.library.stats;
 
-import static com.oracle.truffle.r.library.stats.StatsUtil.rdtciv;
-import static com.oracle.truffle.r.library.stats.StatsUtil.rdtqiv;
+import static com.oracle.truffle.r.library.stats.DPQ.rdtciv;
+import static com.oracle.truffle.r.library.stats.DPQ.rdtqiv;
 
 /*
  * Logic derived from GNU-R, see inline comments.
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rf.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rf.java
index 5d9426797321df6200e4dbc4a00eb9c9e551cb87..21cbb6d2bcf8609018f0fa730b73c6f59f0c8892 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rf.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rf.java
@@ -18,7 +18,7 @@ public final class Rf implements RandFunction2_Double {
     @Override
     public double evaluate(double n1, double n2, RandomNumberProvider rand) {
         if (Double.isNaN(n1) || Double.isNaN(n2) || n1 <= 0. || n2 <= 0.) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
 
         double v1;
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rnorm.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rnorm.java
index 24f2aed7d53b9d6f82635bec1bb53b441f4daa9b..1996274e87ea2c9c733dfdf4ad7ff06cb2eb0270 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rnorm.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rnorm.java
@@ -18,7 +18,7 @@ public final class Rnorm implements RandFunction2_Double {
     @Override
     public double evaluate(double mu, double sigma, RandomNumberProvider rand) {
         if (Double.isNaN(mu) || !Double.isFinite(sigma) || sigma < 0.) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
         if (sigma == 0. || !Double.isFinite(mu)) {
             return mu; /* includes mu = +/- Inf with finite sigma */
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rt.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rt.java
index 2ba6a6271371c136fe6a29a544797de96be0513f..51a266499e2e0af2b8c120ff15edcc26ae123310 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rt.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rt.java
@@ -20,7 +20,7 @@ public final class Rt implements RandFunction1_Double {
     @Override
     public double evaluate(double df, RandomNumberProvider rand) {
         if (Double.isNaN(df) || df <= 0.0) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
 
         if (!Double.isFinite(df)) {
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Runif.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Runif.java
index 57a877536fa532404eb5899b115ba11e4c812043..5701011cf53f8e9f6ee09cae6e6294696b393de8 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Runif.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Runif.java
@@ -30,7 +30,7 @@ public final class Runif implements RandFunction2_Double {
     @Override
     public double evaluate(double min, double max, RandomNumberProvider rand) {
         if (!RRuntime.isFinite(min) || !RRuntime.isFinite(max) || max < min) {
-            return StatsUtil.mlError();
+            return RMath.mlError();
         }
         if (min == max) {
             return min;
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Signrank.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Signrank.java
index 422d76cb6f7c006a3e94ad3779e76fc6d9cfb0b2..1b33cb4ddb8d4aeac8960ed8a69100feefeec1e2 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Signrank.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Signrank.java
@@ -31,12 +31,12 @@ public final class Signrank {
             if (Double.isInfinite(n)) {
                 // In GnuR these "results" seem to be generated due to the behaviour of R_forceint,
                 // and the "(int) n" cast, which ends up casting +/-infinity to integer...
-                return n < 0 ? StatsUtil.mlError() : 0;
+                return n < 0 ? RMath.mlError() : 0;
             }
 
             n = forceint(n);
             if (n < 0) {
-                return StatsUtil.mlError();
+                return RMath.mlError();
             }
 
             if (n == 0) {
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/TOMS708.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/TOMS708.java
index a61e44b42201c5c5af425fa83c95fd84718d8e8f..65056493bd7f7358f33beff9b9ef16afa8141c67 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/TOMS708.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/TOMS708.java
@@ -177,8 +177,8 @@ public class TOMS708 {
             eps = RRuntime.EPSILON; /* == DBL_EPSILON (in R, Rmath) */
 
             /* ----------------------------------------------------------------------- */
-            w = DPQ.d0(logP);
-            w1 = DPQ.d0(logP);
+            w = DPQ.rd0(logP);
+            w1 = DPQ.rd0(logP);
 
             if (a < 0.0 || b < 0.0) {
                 return result(w, w1, 1);
@@ -206,8 +206,8 @@ public class TOMS708 {
                     return result(w, w1, 6);
                 }
                 // else:
-                w = DPQ.d0(logP);
-                w1 = DPQ.d1(logP);
+                w = DPQ.rd0(logP);
+                w1 = DPQ.rd1(logP);
                 return result(w, w1, 0);
             }
             if (y == 0.0) {
@@ -215,20 +215,20 @@ public class TOMS708 {
                     return result(w, w1, 7);
                 }
                 // else:
-                w = DPQ.d1(logP);
-                w1 = DPQ.d0(logP);
+                w = DPQ.rd1(logP);
+                w1 = DPQ.rd0(logP);
                 return result(w, w1, 0);
             }
 
             if (a == 0.0) {
                 // else:
-                w = DPQ.d1(logP);
-                w1 = DPQ.d0(logP);
+                w = DPQ.rd1(logP);
+                w1 = DPQ.rd0(logP);
                 return result(w, w1, 0);
             }
             if (b == 0.0) {
-                w = DPQ.d0(logP);
-                w1 = DPQ.d1(logP);
+                w = DPQ.rd0(logP);
+                w1 = DPQ.rd1(logP);
                 return result(w, w1, 0);
             }
 
@@ -762,7 +762,7 @@ public class TOMS708 {
          */
 
         if (x == 0.) {
-            return DPQ.d0(logP);
+            return DPQ.rd0(logP);
         }
         /* ----------------------------------------------------------------------- */
         /* compute the factor x^a/(a*Beta(a,b)) */
@@ -850,7 +850,7 @@ public class TOMS708 {
             }
         }
         debugPrintf(" bpser(a=%f, b=%f, x=%f, log=%b): prelim.ans = %.14f;\n", a, b, x, logP, ans);
-        if (ans == DPQ.d0(logP) || (!logP && a <= eps * 0.1)) {
+        if (ans == DPQ.rd0(logP) || (!logP && a <= eps * 0.1)) {
             return ans;
         }
 
@@ -1059,7 +1059,7 @@ public class TOMS708 {
         /* R has M_1_SQRT_2PI , and M_LN_SQRT_2PI = ln(sqrt(2*pi)) = 0.918938.. */
 
         if (x == 0.0 || y == 0.0) {
-            return DPQ.d0(logP);
+            return DPQ.rd0(logP);
         }
         double a0 = min(a, b);
         if (a0 < 8.0) {
@@ -1081,7 +1081,7 @@ public class TOMS708 {
             double z = a * lnx + b * lny;
             if (a0 >= 1.) {
                 z -= betaln(a, b);
-                return DPQ.dExp(z, logP);
+                return DPQ.rdexp(z, logP);
             }
 
             /* ----------------------------------------------------------------------- */
@@ -1098,7 +1098,7 @@ public class TOMS708 {
 
             if (b0 <= 1.0) { /* algorithm for max(a,b) = b0 <= 1 */
 
-                double eZ = DPQ.dExp(z, logP);
+                double eZ = DPQ.rdexp(z, logP);
 
                 if (!logP && eZ == 0.0) {
                     /* exp() underflow */
diff --git a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseGammaFunctions.java b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseGammaFunctions.java
index effc12dc4a140b4275b8d7822ef084206ac59b97..5b82397a22c02c2e4d0295b89b3b570970691762 100644
--- a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseGammaFunctions.java
+++ b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseGammaFunctions.java
@@ -12,13 +12,13 @@
  */
 package com.oracle.truffle.r.nodes.builtin.base;
 
-import static com.oracle.truffle.r.library.stats.StatsUtil.DBLEPSILON;
-import static com.oracle.truffle.r.library.stats.StatsUtil.DBL_MANT_DIG;
-import static com.oracle.truffle.r.library.stats.StatsUtil.DBL_MAX_EXP;
-import static com.oracle.truffle.r.library.stats.StatsUtil.DBL_MIN_EXP;
-import static com.oracle.truffle.r.library.stats.StatsUtil.M_LOG10_2;
-import static com.oracle.truffle.r.library.stats.StatsUtil.M_PI;
-import static com.oracle.truffle.r.library.stats.StatsUtil.fmax2;
+import static com.oracle.truffle.r.library.stats.MathConstants.DBL_EPSILON;
+import static com.oracle.truffle.r.library.stats.MathConstants.DBL_MANT_DIG;
+import static com.oracle.truffle.r.library.stats.MathConstants.DBL_MAX_EXP;
+import static com.oracle.truffle.r.library.stats.MathConstants.DBL_MIN_EXP;
+import static com.oracle.truffle.r.library.stats.MathConstants.M_LOG10_2;
+import static com.oracle.truffle.r.library.stats.MathConstants.M_PI;
+import static com.oracle.truffle.r.library.stats.RMath.fmax2;
 import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.complexValue;
 import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.numericValue;
 import static com.oracle.truffle.r.runtime.RDispatch.MATH_GROUP_GENERIC;
@@ -300,7 +300,7 @@ public class BaseGammaFunctions {
             // mVal.nz = 0;
             xln = Math.log(x);
             if (kode == 1 /* && m == 1 */) { /* the R case --- for very large x: */
-                double lrg = 1 / (2. * DBLEPSILON);
+                double lrg = 1 / (2. * DBL_EPSILON);
                 if (n == 0 && x * xln > lrg) {
                     // ans[0] = -xln;
                     ans = -xln;
@@ -316,7 +316,7 @@ public class BaseGammaFunctions {
             nx = Math.min(-DBL_MIN_EXP, DBL_MAX_EXP);
             assert (nx == 1021);
             r1m5 = M_LOG10_2; // Rf_d1mach(5);
-            r1m4 = DBLEPSILON * 0.5; // Rf_d1mach(4) * 0.5;
+            r1m4 = DBL_EPSILON * 0.5; // Rf_d1mach(4) * 0.5;
             wdtol = fmax2(r1m4, 0.5e-18); /* 1.11e-16 */
 
             /* elim = approximate exponential over and underflow limit */
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 85cebcd63d1d487ea9b7df578b3df47bfaf25174..f143ca599f53f21618f57a0740fdf1daf303f033 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
@@ -74,6 +74,7 @@ import com.oracle.truffle.r.library.stats.RNbinomMu;
 import com.oracle.truffle.r.library.stats.RNchisq;
 import com.oracle.truffle.r.library.stats.RPois;
 import com.oracle.truffle.r.library.stats.RWeibull;
+import com.oracle.truffle.r.library.stats.RandGenerationFunctions;
 import com.oracle.truffle.r.library.stats.RandGenerationFunctionsFactory;
 import com.oracle.truffle.r.library.stats.Rbinom;
 import com.oracle.truffle.r.library.stats.Rf;
@@ -84,7 +85,6 @@ 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.StatsFunctionsFactory;
-import com.oracle.truffle.r.library.stats.StatsUtil;
 import com.oracle.truffle.r.library.stats.Wilcox.RWilcox;
 import com.oracle.truffle.r.library.tools.C_ParseRdNodeGen;
 import com.oracle.truffle.r.library.tools.DirChmodNodeGen;
@@ -396,7 +396,7 @@ public class CallAndExternalFunctions {
                     return getExternalModelBuiltinNode("updateform");
 
                 case "Cdqrls":
-                    return new RInternalCodeBuiltinNode(RContext.getInstance(), "stats", RInternalCode.loadSourceRelativeTo(StatsUtil.class, "lm.R"), "Cdqrls");
+                    return new RInternalCodeBuiltinNode(RContext.getInstance(), "stats", RInternalCode.loadSourceRelativeTo(RandGenerationFunctions.class, "lm.R"), "Cdqrls");
 
                 case "dnorm":
                     return StatsFunctionsFactory.Function3_1NodeGen.create(new Dnorm4());
diff --git a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/LookupAdapter.java b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/LookupAdapter.java
index 9569ccc557fd49addc15c425f2d5e320b3f2ad54..5da46c93311b5d2f68f85a60ba0b085c4b3539c8 100644
--- a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/LookupAdapter.java
+++ b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/LookupAdapter.java
@@ -29,7 +29,7 @@ import com.oracle.truffle.api.dsl.Specialization;
 import com.oracle.truffle.api.frame.VirtualFrame;
 import com.oracle.truffle.api.nodes.Node;
 import com.oracle.truffle.api.profiles.BranchProfile;
-import com.oracle.truffle.r.library.stats.StatsUtil;
+import com.oracle.truffle.r.library.stats.RandGenerationFunctions;
 import com.oracle.truffle.r.nodes.access.vector.ElementAccessMode;
 import com.oracle.truffle.r.nodes.access.vector.ExtractVectorNode;
 import com.oracle.truffle.r.nodes.builtin.RBuiltinNode;
@@ -48,9 +48,9 @@ import com.oracle.truffle.r.runtime.data.RLogical;
 import com.oracle.truffle.r.runtime.data.RMissing;
 import com.oracle.truffle.r.runtime.data.model.RAbstractStringVector;
 import com.oracle.truffle.r.runtime.ffi.DLL;
-import com.oracle.truffle.r.runtime.ffi.NativeCallInfo;
 import com.oracle.truffle.r.runtime.ffi.DLL.DLLInfo;
 import com.oracle.truffle.r.runtime.ffi.DLL.SymbolHandle;
+import com.oracle.truffle.r.runtime.ffi.NativeCallInfo;
 
 /**
  * Locator for "builtin" package function implementations. The "builtin" packages contain many
@@ -172,7 +172,7 @@ abstract class LookupAdapter extends RBuiltinNode {
     }
 
     protected static RExternalBuiltinNode getExternalModelBuiltinNode(String name) {
-        return new RInternalCodeBuiltinNode(RContext.getInstance(), "stats", RInternalCode.loadSourceRelativeTo(StatsUtil.class, "model.R"), name);
+        return new RInternalCodeBuiltinNode(RContext.getInstance(), "stats", RInternalCode.loadSourceRelativeTo(RandGenerationFunctions.class, "model.R"), name);
     }
 
     protected static final int CallNST = DLL.NativeSymbolType.Call.ordinal();
diff --git a/mx.fastr/copyrights/overrides b/mx.fastr/copyrights/overrides
index 9b4374cd86b4e87ac9df4ec8a41fd2c37c4eee29..475db3491e56345526579054797f8b1119c3077b 100644
--- a/mx.fastr/copyrights/overrides
+++ b/mx.fastr/copyrights/overrides
@@ -42,7 +42,6 @@ com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DoubleCentre
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java,gnu_r.core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rt.java,gnu_r_ihaka_core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java,gnu_r_qgamma.copyright
-com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathConstants.java,gnu_r.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QHyper.java,gnu_r_ihaka_core.copyright
 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
@@ -61,7 +60,8 @@ com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Chisq.java,g
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/SplineFunctions.java,gnu_r_splines.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/StatsFunctions.java,gnu_r_gentleman_ihaka.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RandGenerationFunctions.java,gnu_r_gentleman_ihaka.copyright
-com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/StatsUtil.java,gnu_r_ihaka.copyright
+com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java,gnu_r_ihaka.copyright
+com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathConstants.java,gnu_r_ihaka.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/TOMS708.java,gnu_r.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/SNorm.java,gnu_r_ihaka_core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/SExp.java,gnu_r_ihaka_core.copyright