From 991965d9e421bd5ed19b897cf92dfd8662888000 Mon Sep 17 00:00:00 2001
From: stepan <stepan.sindelar@oracle.com>
Date: Wed, 21 Dec 2016 14:35:52 +0100
Subject: [PATCH] Fixes in nmath ported code.
- R_forceint is not equivalent to Math.round,
- forceint moved to RMath
- true/false typo in TOMS708
---
.../oracle/truffle/r/library/stats/DPQ.java | 4 +--
.../oracle/truffle/r/library/stats/DPois.java | 2 +-
.../truffle/r/library/stats/Dbinom.java | 2 +-
.../oracle/truffle/r/library/stats/Geom.java | 2 +-
.../r/library/stats/MathConstants.java | 7 ----
.../truffle/r/library/stats/Pbinom.java | 4 +--
.../truffle/r/library/stats/QHyper.java | 2 +-
.../truffle/r/library/stats/RHyper.java | 2 +-
.../truffle/r/library/stats/RLogis.java | 31 ------------------
.../oracle/truffle/r/library/stats/RMath.java | 12 +++++++
.../truffle/r/library/stats/Rbinom.java | 2 +-
.../truffle/r/library/stats/Signrank.java | 2 +-
.../truffle/r/library/stats/TOMS708.java | 32 ++++++++-----------
13 files changed, 37 insertions(+), 67 deletions(-)
delete mode 100644 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RLogis.java
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 24e3dfb71a..35a19b7a93 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
@@ -43,7 +43,7 @@ public final class DPQ {
// 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 nointCheck instead
public static boolean nonint(double x) {
- return Math.abs(x - Math.round(x)) > 1e-7 * Math.max(1., Math.abs(x));
+ return Math.abs(x - RMath.forceint(x)) > 1e-7 * Math.max(1., Math.abs(x));
}
// R_D__0
@@ -156,7 +156,7 @@ public final class DPQ {
// R_P_bounds_Inf_01
public static void rpboundsinf01(double x, boolean lowerTail, boolean logP) throws EarlyReturn {
if (!Double.isFinite(x)) {
- throw new EarlyReturn(x > 0 ? rdt0(lowerTail, logP) : rdt0(lowerTail, logP));
+ throw new EarlyReturn(x > 0 ? rdt1(lowerTail, logP) : rdt0(lowerTail, logP));
}
}
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 b01d2a4d75..c60b2f0bd5 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
@@ -13,7 +13,7 @@
package com.oracle.truffle.r.library.stats;
import static com.oracle.truffle.r.library.stats.GammaFunctions.dpoisRaw;
-import static com.oracle.truffle.r.library.stats.MathConstants.forceint;
+import static com.oracle.truffle.r.library.stats.RMath.forceint;
import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
import com.oracle.truffle.r.library.stats.StatsFunctions.Function2_1;
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 509f543c81..0486d311a7 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
@@ -88,6 +88,6 @@ public final class Dbinom implements StatsFunctions.Function3_1 {
return DPQ.rd0(giveLog);
}
- return dbinomRaw(MathConstants.forceint(x), MathConstants.forceint(n), p, 1 - p, giveLog);
+ return dbinomRaw(RMath.forceint(x), RMath.forceint(n), p, 1 - p, giveLog);
}
}
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 05443e198b..29390fe211 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
@@ -27,7 +27,7 @@
*/
package com.oracle.truffle.r.library.stats;
-import static com.oracle.truffle.r.library.stats.MathConstants.forceint;
+import static com.oracle.truffle.r.library.stats.RMath.forceint;
import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandFunction1_Double;
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 0885fbeecf..61bc036c42 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
@@ -97,11 +97,4 @@ public final class MathConstants {
public static double logspaceAdd(double logx, double logy) {
return Math.max(logx, logy) + Math.log1p(Math.exp(-Math.abs(logx - logy)));
}
-
- // R_forceint
- public static double forceint(double x) {
- // Note: in GnuR this is alias for nearbyint, which may not behave exactly like Math.round,
- // especially Math.round(-0.5) == 0.0, instead of -0.0, does it matter a lot?
- return Double.isFinite(x) ? Math.round(x) : x;
- }
}
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 70fcce282d..e61ab3cb1f 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,9 +34,9 @@ public final class Pbinom implements StatsFunctions.Function3_2 {
if (DPQ.nonint(size)) {
nanProfile.enter();
DPQ.nointCheckWarning(size, "n");
- return DPQ.rd0(logP);
+ return RMath.mlError();
}
- size = Math.round(size);
+ size = RMath.forceint(size);
/* PR#8560: n=0 is a valid value */
if (size < 0 || prob < 0 || prob > 1) {
nanProfile.enter();
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 9a54383f82..e0641de73b 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
@@ -12,9 +12,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.RMath.fmax2;
import static com.oracle.truffle.r.library.stats.RMath.fmin2;
+import static com.oracle.truffle.r.library.stats.RMath.forceint;
import static com.oracle.truffle.r.library.stats.RMath.lfastchoose;
import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
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 4cd77bcd0f..197316d6df 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
@@ -13,7 +13,7 @@
package com.oracle.truffle.r.library.stats;
import static com.oracle.truffle.r.library.stats.MathConstants.M_LN_SQRT_2PI;
-import static com.oracle.truffle.r.library.stats.MathConstants.forceint;
+import static com.oracle.truffle.r.library.stats.RMath.forceint;
import com.oracle.truffle.api.CompilerDirectives.TruffleBoundary;
import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandFunction3_Double;
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
deleted file mode 100644
index f21be8e24c..0000000000
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RLogis.java
+++ /dev/null
@@ -1,31 +0,0 @@
-/*
- * 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) 1998--2008, The R Core Team
- * Copyright (c) 2016, 2016, Oracle and/or its affiliates
- *
- * All rights reserved.
- */
-package com.oracle.truffle.r.library.stats;
-
-import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandFunction2_Double;
-import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandomNumberProvider;
-
-public final class RLogis extends RandFunction2_Double {
- @Override
- public double execute(double location, double scale, RandomNumberProvider rand) {
- if (Double.isNaN(location) || !Double.isFinite(scale)) {
- return RMath.mlError();
- }
-
- if (scale == 0. || !Double.isFinite(location)) {
- return location;
- } else {
- double u = rand.unifRand();
- return location + scale * Math.log(u / (1. - u));
- }
- }
-}
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java
index 02c6164331..49fc18d2d6 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RMath.java
@@ -37,6 +37,18 @@ public class RMath {
return -Math.log(n + 1.) - lbeta(n - k + 1., k + 1.);
}
+ /**
+ * Implementation of {@code R_forceint}, which is not equal to {@code Math.round}, because it
+ * returns {@code double} and so it can handle values that do not fit into long.
+ */
+ public static double forceint(double x) {
+ // Note: in GnuR this is alias for nearbyint
+ if (Double.isNaN(x)) {
+ return 0;
+ }
+ return Math.floor(x + 0.5);
+ }
+
public static double fsign(double x, double y) {
if (Double.isNaN(x) || Double.isNaN(y)) {
return x + y;
diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rbinom.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rbinom.java
index edd2864e79..37964392e6 100644
--- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rbinom.java
+++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rbinom.java
@@ -30,7 +30,7 @@ public final class Rbinom extends RandFunction2_Double {
if (!Double.isFinite(nin)) {
return RRuntime.INT_NA;
}
- double r = MathConstants.forceint(nin);
+ double r = RMath.forceint(nin);
if (r != nin) {
return RRuntime.INT_NA;
}
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 9a976a351a..f13e125ea9 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
@@ -11,7 +11,7 @@
package com.oracle.truffle.r.library.stats;
-import static com.oracle.truffle.r.library.stats.MathConstants.forceint;
+import static com.oracle.truffle.r.library.stats.RMath.forceint;
import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandFunction1_Double;
import com.oracle.truffle.r.library.stats.RandGenerationFunctions.RandomNumberProvider;
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 65056493bd..b9f92d9cb1 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
@@ -40,7 +40,7 @@ public class TOMS708 {
// R_Log1_Exp
public static double log1Exp(double x) {
- return ((x) > -M_LN2 ? log(-rexpm1(x)) : log1p(-exp(x)));
+ return ((x) > -M_LN2 ? log(-rexpm1(x)) : RMath.log1p(-exp(x)));
}
private static double sin(double v) {
@@ -59,10 +59,6 @@ public class TOMS708 {
return Math.sqrt(v);
}
- private static double log1p(double v) {
- return Math.log1p(v);
- }
-
private static double exp(double v) {
return Math.exp(v);
}
@@ -357,7 +353,7 @@ public class TOMS708 {
w1 = Double.NEGATIVE_INFINITY; // = 0 on log-scale
}
bgrat.w = w1;
- bgrat.bgrat(b0, a0, y0, x0, 15 * eps, false);
+ bgrat.bgrat(b0, a0, y0, x0, 15 * eps, true);
w1 = bgrat.w;
if (bgrat.ierr != 0) {
ierr = 8;
@@ -594,7 +590,7 @@ public class TOMS708 {
// exp(a*lnx) underflows for large (a * lnx); e.g. large a ==> using log_r := log(r):
// r = r1 * exp(a * lnx) * exp(bm1 * 0.5 * lnx);
// log(r)=log(b) + log1p(gam1(b)) + b * log(z) + (a * lnx) + (bm1 * 0.5 * lnx),
- double logR = Math.log(b) + log1p(gam1(b)) + b * Math.log(z) + nu * lnx;
+ double logR = Math.log(b) + RMath.log1p(gam1(b)) + b * Math.log(z) + nu * lnx;
// FIXME work with log_u = log(u) also when log_p=FALSE (??)
// u is 'factored out' from the expansion {and multiplied back, at the end}:
double logU = logR - (algdiv(b, a) + b * Math.log(nu)); // algdiv(b,a) =
@@ -613,7 +609,7 @@ public class TOMS708 {
double l = // := *w/u .. but with care: such that it also works when u underflows to 0:
logW ? ((w == Double.NEGATIVE_INFINITY) ? 0. : Math.exp(w - logU)) : ((w == 0.) ? 0. : Math.exp(Math.log(w) - logU));
- debugPrintf(" bgrat(a=%f, b=%f, x=%f, *)\n -> u=%f, l='w/u'=%f, ", a, b, x, u, l);
+ debugPrintf(" bgrat(a=%f, b=%f, x=%f, *)\n -> u=%f, l='w/u'=%f, \n", a, b, x, u, l);
double qR = gratR(b, z, logR, eps); // = q/r of former grat1(b,z, r, &p, &q)
double v = 0.25 / (nu * nu);
@@ -713,7 +709,7 @@ public class TOMS708 {
} while (fabs(c) > tol);
if (logP) {
- ans += log1p(a * s);
+ ans += RMath.log1p(a * s);
} else {
ans *= a * s + 1.0;
}
@@ -831,7 +827,7 @@ public class TOMS708 {
if (logP) {
/* FIXME? potential for improving log(t) */
- ans = z + log(a0 / a) + log1p(gam1(b0)) - log(t);
+ ans = z + log(a0 / a) + RMath.log1p(gam1(b0)) - log(t);
} else {
ans = exp(z) * (a0 / a) * (gam1(b0) + 1.0) / t;
}
@@ -871,14 +867,14 @@ public class TOMS708 {
} while (n < 1e7 && fabs(w) > tol);
if (fabs(w) > tol) { // the series did not converge (in time)
// warn only when the result seems to matter:
- if ((logP && !(a * sum > -1. && fabs(log1p(a * sum)) < eps * fabs(ans))) || (!logP && fabs(a * sum + 1) != 1.)) {
+ if ((logP && !(a * sum > -1. && fabs(RMath.log1p(a * sum)) < eps * fabs(ans))) || (!logP && fabs(a * sum + 1) != 1.)) {
emitWarning(" bpser(a=%f, b=%f, x=%f,...) did not converge (n=1e7, |w|/tol=%f > 1; A=%f)", a, b, x, fabs(w) / tol, ans);
}
}
debugPrintf(" -> n=%d iterations, |w|=%f %s %f=tol:=eps/a ==> a*sum=%f\n", (int) n, fabs(w), (fabs(w) > tol) ? ">!!>" : "<=", tol, a * sum);
if (logP) {
if (a * sum > -1.0) {
- ans += log1p(a * sum);
+ ans += RMath.log1p(a * sum);
} else {
ans = ML_NEGINF;
}
@@ -1115,7 +1111,7 @@ public class TOMS708 {
double c = (gam1(a) + 1.0) * (gam1(b) + 1.0) / z;
/* FIXME? log(a0*c)= log(a0)+ log(c) and that is improvable */
- return (logP ? eZ + log(a0 * c) - log1p(a0 / b0) : eZ * (a0 * c) / (a0 / b0 + 1.0));
+ return (logP ? eZ + log(a0 * c) - RMath.log1p(a0 / b0) : eZ * (a0 * c) / (a0 / b0 + 1.0));
}
/* else : ALGORITHM FOR 1 < b0 < 8 */
@@ -1141,7 +1137,7 @@ public class TOMS708 {
t = gam1(apb) + 1.0;
}
- return (logP ? log(a0) + z + log1p(gam1(b0)) - log(t) : a0 * exp(z) * (gam1(b0) + 1.0) / t);
+ return (logP ? log(a0) + z + RMath.log1p(gam1(b0)) - log(t) : a0 * exp(z) * (gam1(b0) + 1.0) / t);
} else {
/* ----------------------------------------------------------------------- */
@@ -1248,9 +1244,9 @@ public class TOMS708 {
z = gam1(apb) + 1.0;
}
// L50:
- double c = giveLog ? log1p(gam1(a)) + log1p(gam1(b)) - log(z) : (gam1(a) + 1.0) * (gam1(b) + 1.0) / z;
+ double c = giveLog ? RMath.log1p(gam1(a)) + RMath.log1p(gam1(b)) - log(z) : (gam1(a) + 1.0) * (gam1(b) + 1.0) / z;
debugPrintf(" brcmp1(mu,a,b,*): a0 < 1, b0 <= 1; c=%.15f\n", c);
- return giveLog ? ans + log(a0) + c - log1p(a0 / b0) : ans * (a0 * c) / (a0 / b0 + 1.0);
+ return giveLog ? ans + log(a0) + c - RMath.log1p(a0 / b0) : ans * (a0 * c) / (a0 / b0 + 1.0);
}
// else: algorithm for a0 < 1 < b0 < 8
// L60:
@@ -1278,7 +1274,7 @@ public class TOMS708 {
}
debugPrintf(" brcmp1(mu,a,b,*): a0 < 1 < b0 < 8; t=%.15f\n", t);
// L72:
- return giveLog ? log(a0) + esum(mu, z, true) + log1p(gam1(b0)) - log(t) : a0 * esum(mu, z, false) * (gam1(b0) + 1.0) / t;
+ return giveLog ? log(a0) + esum(mu, z, true) + RMath.log1p(gam1(b0)) - log(t) : a0 * esum(mu, z, false) * (gam1(b0) + 1.0) / t;
} else {
@@ -1302,7 +1298,7 @@ public class TOMS708 {
y0 = 1.0 / (h + 1.0);
lambda = a - (a + b) * x;
}
- double lx0 = -log1p(b / a); // in both cases
+ double lx0 = -RMath.log1p(b / a); // in both cases
debugPrintf(" brcmp1(mu,a,b,*): a,b >= 8; x0=%.15f, lx0=log(x0)=%.15f\n", x0, lx0);
// L110:
--
GitLab