Calculate power for int64 exponent precisely

NEWS.md

@@ -20,6 +20,8 @@
 
 1. `log()` operations, especially for `log(x, base=10)` and `log(x, base=2)`, have better precision and consistency (#180).
 
+1. `^.integer64` with integer or integer64 exponent now calculates the power precisely and returns an overflow warning, if an overflow appears (#288).
+
 ## NOTES
 
 1. The R version dependency has been bumped from 3.5.0 (2018) to 3.6.0 (2019).

R/ops64.R

@@ -15,6 +15,10 @@
 #' @param e1,e2,x numeric or complex vectors or objects which can be coerced to such,
 #'   or other objects for which methods have been written for - especially 'integer64' vectors.
 #'
+#' @details
+#' [`^`] with 'integer' or 'integer64' exponent now calculates the power precisely and returns
+#'   an overflow warning, if an overflow appears.
+#'
 #' @returns
 #'   - [`&`], [`|`], [`!`], [`!=`], [`==`], [`<`], [`<=`], [`>`], [`>=`] return a logical vector
 #'   - [`/`] returns a double vector

src/integer64.c

@@ -366,10 +366,11 @@ SEXP power_integer64_integer64(SEXP e1_, SEXP e2_, SEXP ret_){
   long long * e1 = (long long *) REAL(e1_);
   long long * e2 = (long long *) REAL(e2_);
   long long * ret = (long long *) REAL(ret_);
-  long double longret;
   Rboolean naflag = FALSE;
     mod_iterate(n1, n2, i1, i2) {
-        POW64(e1[i1],e2[i2],ret[i],naflag, longret)
+      if (pow64_overflow(e1[i1], e2[i2], &ret[i])) {
+        naflag = TRUE;
+      }
     }
     if (naflag)warning(INTEGER64_OVERFLOW_WARNING);
   return ret_;

src/integer64.h

@@ -118,6 +118,64 @@ static inline bool mul64_overflow(long long a, long long b, long long *res) {
 #endif
 }
 
+static inline bool pow64_overflow(long long base, long long exp, long long *res) {
+  // special cases: n^0, 1^m, NA^m, n^NA, n^1, n^2, 0^m, -1^m, n^-m
+  if (exp == 0) {
+    *res = 1;
+    return false;
+  }
+  if (base == 1) {
+    *res = 1;
+    return false;
+  }
+  if (base == NA_INTEGER64 || exp == NA_INTEGER64) {
+    *res = NA_INTEGER64;
+    return false;
+  }
+  if (exp == 1) {
+    *res = base;
+    return false;
+  }
+  if (exp == 2) {
+    if (mul64_overflow(base, base, res)) {
+      *res = NA_INTEGER64;
+      return true;
+    }
+    return false;
+  }
+  if (base == 0) {
+    if (exp < 0) {
+      *res = NA_INTEGER64;
+      return true;
+    }
+    *res = 0;
+    return false;
+  }
+  if (base == -1) {
+    *res = exp & 1 ? -1 : 1;
+    return false;
+  }
+  if (exp < 0) {
+    *res = 0;
+    return false;
+  }
+  long long r = 1;
+  while (1) {
+    if (exp & 1 && mul64_overflow(r, base, &r)) {
+      *res = NA_INTEGER64;
+      return true;
+    }
+    exp >>= 1;
+    if (!exp) break;
+    if (mul64_overflow(base, base, &base)) {
+      *res = NA_INTEGER64;
+      return true;
+    }
+  }
+  *res = r;
+  return false;
+}
+
 #define PLUS64(e1,e2,ret,naflag) \
     if (e1 == NA_INTEGER64 || e2 == NA_INTEGER64) \
         ret = NA_INTEGER64; \
@@ -154,18 +212,6 @@ static inline bool mul64_overflow(long long a, long long b, long long *res) {
           ret = llroundl(longret); \
     }
 
-#define POW64(e1,e2,ret,naflag, longret) \
-    if (e1 == NA_INTEGER64 || e2 == NA_INTEGER64) \
-        ret = NA_INTEGER64; \
-    else { \
-        longret = pow(e1, (long double) e2); \
-        if (isnan(longret)){ \
-          naflag = TRUE; \
-          ret = NA_INTEGER64; \
-        }else \
-          ret = llroundl(longret); \
-    }
-
 #define POW64REAL(e1,e2,ret,naflag,longret) \
     if (e1 == NA_INTEGER64 || ISNAN(e2)) \
         ret = NA_INTEGER64; \

tests/testthat/test-ops64.R

@@ -383,3 +383,98 @@ test_that("Edge cases for character/factor comparisons work", {
   # nolint next: expect_comparison_linter. Checking '==' method
   expect_true(as.integer64("999999999999999999") == as.factor("999999999999999999"))
 })
+
+test_that("power with integer64", {
+  overflow_warning = "NAs produced by integer64 overflow"
+
+  # within integer range
+  x = as.integer(sqrt(.Machine$integer.max))
+  x = seq(-x, x)
+  expect_identical(as.integer64(x)^2L, as.integer64(x^2L))
+  expect_identical(as.integer64(x)^2.0, as.integer64(x^2.0))
+
+  # within integer64 range, which fails with double exponent
+  expect_identical(as.integer64(2147483650)^2L, as.integer64("4611686027017322500"))
+  expect_identical(as.integer64(-2147483650)^2L, as.integer64("4611686027017322500"))
+  expect_identical(as.integer64(94906267L)^2L, as.integer64(94906267L) * as.integer64(94906267L))
+
+  # integer64 base with integer64 exponent
+  expect_identical(as.integer64(10L) ^ as.integer64(3L), as.integer64(1000L))
+  expect_identical(as.integer64(2147483650) ^ as.integer64(2L), as.integer64("4611686027017322500"))
+  expect_identical(as.integer64(2L) ^ 62L, as.integer64("4611686018427387904"))
+  expect_identical(as.integer64(2L) ^ as.integer64(62L), as.integer64("4611686018427387904"))
+  expect_identical(as.integer64(3L) ^ 39L, as.integer64("4052555153018976267"))
+
+  # Special bases: 0, 1, -1
+  expect_identical(as.integer64(0L) ^ 0L, as.integer64(1L))
+  expect_identical(as.integer64(0L) ^ 1L, as.integer64(0L))
+  expect_identical(as.integer64(0L) ^ 5L, as.integer64(0L))
+  expect_warning(expect_identical(as.integer64(0L) ^ (-1L), NA_integer64_), overflow_warning)
+
+  expect_identical(as.integer64(1L) ^ 0L, as.integer64(1L))
+  expect_identical(as.integer64(1L) ^ 100L, as.integer64(1L))
+  expect_identical(as.integer64(1L) ^ (-5L), as.integer64(1L))
+  expect_identical(as.integer64(1L) ^ lim.integer64()[2L], as.integer64(1L))
+  expect_identical(as.integer64(1L) ^ NA_integer64_, as.integer64(1L))
+  expect_identical(as.integer64(1L) ^ NA_integer_, as.integer64(1L))
+
+  expect_identical(as.integer64(-1L) ^ 0L, as.integer64(1L))
+  expect_identical(as.integer64(-1L) ^ 1L, as.integer64(-1L))
+  expect_identical(as.integer64(-1L) ^ 2L, as.integer64(1L))
+  expect_identical(as.integer64(-1L) ^ 3L, as.integer64(-1L))
+  expect_identical(as.integer64(-1L) ^ (-1L), as.integer64(-1L))
+  expect_identical(as.integer64(-1L) ^ (-2L), as.integer64(1L))
+  expect_identical(as.integer64(-1L) ^ (-3L), as.integer64(-1L))
+  expect_identical(as.integer64(-1L) ^ lim.integer64()[1L], as.integer64(-1L))
+  expect_identical(as.integer64(-1L) ^ lim.integer64()[2L], as.integer64(-1L))
+
+  # Negative exponents with |base| >= 2
+  expect_identical(as.integer64(2L) ^ (-1L), as.integer64(0L))
+  expect_identical(as.integer64(2L) ^ (-5L), as.integer64(0L))
+  expect_identical(as.integer64(-2L) ^ (-1L), as.integer64(0L))
+  expect_identical(as.integer64(-2L) ^ (-2L), as.integer64(0L))
+
+  # Boundary base values with exponent 1
+  expect_identical(lim.integer64()[1L] ^ 1L, lim.integer64()[1L])
+  expect_identical(lim.integer64()[2L] ^ 1L, lim.integer64()[2L])
+
+  # Overflow detection with both odd and even exponents
+  expect_warning(expect_identical(as.integer64(2147483650) ^ 3L, NA_integer64_), overflow_warning)
+  expect_warning(expect_identical(as.integer64(100000L) ^ 8L, NA_integer64_), overflow_warning)
+  expect_warning(expect_identical(as.integer64(2L) ^ 63L, NA_integer64_), overflow_warning)
+  expect_warning(expect_identical(as.integer64(2L) ^ 64L, NA_integer64_), overflow_warning)
+  expect_warning(expect_identical(as.integer64(3L) ^ 40L, NA_integer64_), overflow_warning)
+  expect_warning(expect_identical(as.integer64(3L) ^ 64L, NA_integer64_), overflow_warning)
+  expect_warning(expect_identical(lim.integer64()[1L] ^ 2L, NA_integer64_), overflow_warning)
+  expect_warning(expect_identical(lim.integer64()[2L] ^ 2L, NA_integer64_), overflow_warning)
+
+  # Missing values and empty inputs
+  expect_no_warning(expect_identical(NA_integer64_ ^ 0L, as.integer64(1L)))
+  expect_no_warning(expect_identical(NA_integer64_ ^ 2L, NA_integer64_))
+  expect_identical(NA_integer64_ ^ 2L, NA_integer64_)
+  expect_no_warning(expect_identical(NA_integer64_ ^ (-1L), NA_integer64_))
+  expect_identical(NA_integer64_ ^ (-1L), NA_integer64_)
+  expect_no_warning(expect_identical(NA_integer64_ ^ NA_integer64_, NA_integer64_))
+  expect_identical(NA_integer64_ ^ NA_integer64_, NA_integer64_)
+  expect_no_warning(expect_identical(as.integer64(0L) ^ NA_integer64_, NA_integer64_))
+  expect_identical(as.integer64(0L) ^ NA_integer64_, NA_integer64_)
+  expect_no_warning(expect_identical(as.integer64(-1L) ^ NA_integer64_, NA_integer64_))
+  expect_identical(as.integer64(-1L) ^ NA_integer64_, NA_integer64_)
+  expect_no_warning(expect_identical(as.integer64(2L) ^ NA_integer_, NA_integer64_))
+  expect_identical(as.integer64(2L) ^ NA_integer_, NA_integer64_)
+  expect_no_warning(expect_identical(as.integer64(2L) ^ NA_integer64_, NA_integer64_))
+  expect_identical(as.integer64(2L) ^ NA_integer64_, NA_integer64_)
+  expect_no_warning(expect_identical(as.integer64(2L) ^ NA_real_, NA_integer64_))
+  expect_identical(as.integer64(2L) ^ NA_real_, NA_integer64_)
+  expect_identical(integer64() ^ 2L, integer64())
+  expect_identical(as.integer64(2L) ^ integer(), integer64())
+  expect_identical(integer64() ^ integer(), integer64())
+
+  # Vectorization and recycling
+  expect_identical(c(as.integer64(2L), as.integer64(3L)) ^ 2L, as.integer64(c(4L, 9L)))
+  expect_identical(as.integer64(2L) ^ c(1L, 2L, 3L), as.integer64(c(2L, 4L, 8L)))
+  expect_warning(
+    expect_identical(c(as.integer64(2L), as.integer64(2L)) ^ c(2L, 64L), c(as.integer64(4L), NA_integer64_)),
+    overflow_warning
+  )
+})