Handle log2, log10 more robustly (#355)

* Gemini: log2, log10 handling improvement

* Gemini: Inline LOG264/LOG1064 in logbase to preserve original function ordering

* Gemini: Add explanatory comments for logarithm routines and macros

* Improve comments

* delint

* Fine-tune NEWS

* Gemini: Add test coverage for general base vector and scalar cases in log.integer64

Author: MichaelChirico

NEWS.md

@@ -18,6 +18,8 @@
 
 1. `as.list()` doesn't overflow the protection stack for integer64 inputs (#345). Thanks @NicChr for the report and fix.
 
+1. `log()` operations, especially for `log(x, base=10)` and `log(x, base=2)`, have better precision and consistency (#180).
+
 ## NOTES
 
 1. The R version dependency has been bumped from 3.5.0 (2018) to 3.6.0 (2019).

src/integer64.c

@@ -467,6 +467,7 @@ SEXP sqrt_integer64(SEXP e1_, SEXP ret_){
   return ret_;
 }
 
+// Natural logarithm: called when base is NULL (i.e. log(x)).
 SEXP log_integer64(SEXP e1_, SEXP ret_){
   long long i, n = LENGTH(ret_);
   long long * e1 = (long long *) REAL(e1_);
@@ -479,6 +480,8 @@ SEXP log_integer64(SEXP e1_, SEXP ret_){
   return ret_;
 }
 
+// Vector base: called when length(base) > 1 (e.g. log(x, base=c(2, 10))).
+// Recycles x and base; dispatches to LOG264/LOG1064 if base element is 2 or 10.
 SEXP logvect_integer64(SEXP e1_, SEXP e2_, SEXP ret_){
   long long i, n = LENGTH(ret_);
   long long i1, n1 = LENGTH(e1_);
@@ -488,20 +491,39 @@ SEXP logvect_integer64(SEXP e1_, SEXP e2_, SEXP ret_){
   double * ret = REAL(ret_);
   Rboolean naflag = FALSE;
   mod_iterate(n1, n2, i1, i2) {
-    LOGVECT64(e1[i], e2[i], ret[i], naflag)
+    if (e2[i2] == 2.0) {
+      LOG264(e1[i1], ret[i], naflag)
+    } else if (e2[i2] == 10.0) {
+      LOG1064(e1[i1], ret[i], naflag)
+    } else {
+      LOGVECT64(e1[i1], e2[i2], ret[i], naflag)
+    }
   }
   if (naflag) warning(INTEGER64_NAN_CREATED_WARNING);
   return ret_;
 }
 
+// Scalar base: called when length(base) == 1 (e.g. log(x, base=2), log(x, base=10)).
+// Dispatches to LOG264/LOG1064 (matching base R's logbase), or precomputes logl(base) once.
 SEXP logbase_integer64(SEXP e1_, SEXP base_, SEXP ret_){
   long long i, n = LENGTH(ret_);
   long long * e1 = (long long *) REAL(e1_);
-  long double logbase = (long double) log(asReal(base_));
+  double base = asReal(base_);
   double * ret = REAL(ret_);
-  Rboolean naflag = (asReal(base_)>0) ? FALSE : TRUE;
-  for(i=0; i<n; i++) {
-    LOGBASE64(e1[i], logbase, ret[i], naflag)
+  Rboolean naflag = (base > 0) ? FALSE : TRUE;
+  if (base == 2.0) {
+    for(i=0; i<n; i++) {
+      LOG264(e1[i], ret[i], naflag)
+    }
+  } else if (base == 10.0) {
+    for(i=0; i<n; i++) {
+      LOG1064(e1[i], ret[i], naflag)
+    }
+  } else {
+    long double logbase = (long double) logl((long double)base);
+    for(i=0; i<n; i++) {
+      LOGBASE64(e1[i], logbase, ret[i], naflag)
+    }
   }
   if (naflag) warning(INTEGER64_NAN_CREATED_WARNING);
   return ret_;

src/integer64.h

@@ -225,7 +225,7 @@ else {                                                   \
     else { \
         if (e1 < 0) \
             naflag = TRUE; \
-        ret = (double) sqrt((long double)e1); \
+        ret = (double) sqrtl((long double)e1); \
     }
 
 #define LOG64(e1, ret, naflag) \
@@ -237,11 +237,12 @@ else {                                                   \
             naflag = TRUE; \
     }
 
+// NB: cast to double _after_ dividing in 'long double' for max precision.
 #define LOGVECT64(e1, e2, ret, naflag) \
     if (e1 == NA_INTEGER64) \
         ret = NA_REAL; \
     else { \
-        ret = (double) logl((long double)e1) / log(e2); \
+        ret = (double) (logl((long double)e1) / logl((long double)e2)); \
         if (isnan(ret)) \
             naflag = TRUE; \
     }
@@ -250,7 +251,7 @@ else {                                                   \
     if (e1 == NA_INTEGER64) \
         ret = NA_REAL; \
     else { \
-        ret = (double) logl((long double)e1) / e2; \
+        ret = (double) (logl((long double)e1) / e2); \
         if (isnan(ret)) \
             naflag = TRUE; \
     }

tests/testthat/test-integer64.R

@@ -336,8 +336,13 @@ test_that("arithmetic & basic math works", {
   expect_identical(sqrt(as.integer64(c(0L, 1L, 4L, 9L))), as.numeric(0:3))
   expect_identical(log(x), log(as.numeric(x)))
   expect_identical(log(as.integer64(c(1L, 2L, 4L, 8L)), base=2L), as.numeric(0:3))
+  expect_identical(log(as.integer64(c(1L, 2L, 4L, 8L)), base=c(2L, 2L, 2L, 2L)), as.numeric(0:3))
+  expect_identical(log(as.integer64(c(8L, 1000L)), base=c(2L, 10L)), c(3.0, 3.0))
+  expect_identical(log(as.integer64(c(8L, 27L, 1000L, 25L)), base=c(2L, 3L, 10L, 5L)), c(3.0, 3.0, 3.0, 2.0))
+  expect_identical(log(as.integer64(c(1L, 3L, 9L, 27L)), base=3L), as.numeric(0:3))
   expect_identical(log2(as.integer64(c(1L, 2L, 4L, 8L))), as.numeric(0:3))
-  expect_identical(log10(as.integer64(c(1L, 10L, 100L, 1000L))), as.numeric(0:3), tolerance=1e-7)
+  expect_identical(log10(as.integer64(c(1L, 10L, 100L, 1000L))), as.numeric(0:3))
+  expect_identical(log(as.integer64(c(1L, 10L, 100L, 1000L)), base=10L), as.numeric(0:3))
 
   expect_identical(trunc(x), x)
   expect_identical(floor(x), x)
@@ -346,6 +351,19 @@ test_that("arithmetic & basic math works", {
   expect_identical(round(x), x)
 
   expect_identical(round(x, -1L), as.integer64(rep(c(0L, 10L), each=5L)))
+  # Halfway round-to-even tests
+  expect_identical(
+    round(as.integer64(c(-25L, -15L, -5L, 5L, 15L, 25L)), -1L),
+    as.integer64(c(-20L, -20L, 0L, 0L, 20L, 20L))
+  )
+  expect_identical(
+    round(as.integer64(c(-250L, -150L, -50L, 50L, 150L, 250L)), -2L),
+    as.integer64(c(-200L, -200L, 0L, 0L, 200L, 200L))
+  )
+  expect_identical(
+    round(as.integer64(c(-2500L, -1500L, -500L, 500L, 1500L, 2500L)), -3L),
+    as.integer64(c(-2000L, -2000L, 0L, 0L, 2000L, 2000L))
+  )
 })
 
 test_that("basic statistics work", {
@@ -757,6 +775,9 @@ local({
   with_parameters_test_that(
     "Old \\dontshow{} tests in ?format.integer64 continue working",
     {
+      # base::round() for negative digits uses floating-point approximations in fround.c
+      # which can fail round-half-to-even on platforms with x87 FPUs (e.g. i686) (#180).
+      skip_if(s < 0L && round(5.0 * 10.0^(-s - 1L), s) != 0.0, "base::round() does not round to even on this platform")
       r <- as.integer64(round(as.integer(i), s))
       r64 <- round(as.integer64(i), s)
       expect_identical(r, r64)