diff --git a/NEWS.md b/NEWS.md index 7655f66a..86404d43 100644 --- a/NEWS.md +++ b/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). diff --git a/R/ops64.R b/R/ops64.R index 058f10e4..ee43e73e 100644 --- a/R/ops64.R +++ b/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 diff --git a/man/ops64.Rd b/man/ops64.Rd index 3c1de089..ee3c4a97 100644 --- a/man/ops64.Rd +++ b/man/ops64.Rd @@ -70,6 +70,10 @@ or other objects for which methods have been written for - especially 'integer64 \description{ Binary operators for integer64 vectors. } +\details{ +\code{\link{^}} with 'integer' or 'integer64' exponent now calculates the power precisely and returns +an overflow warning, if an overflow appears. +} \examples{ as.integer64(1:12) - 1 options(integer64_semantics="new") diff --git a/src/integer64.c b/src/integer64.c index 08dc44cc..8fc1d727 100644 --- a/src/integer64.c +++ b/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_; diff --git a/src/integer64.h b/src/integer64.h index 6289fb5c..ecf6e2f6 100644 --- a/src/integer64.h +++ b/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; \ diff --git a/tests/testthat/test-ops64.R b/tests/testthat/test-ops64.R index c9ac0b99..2e486574 100644 --- a/tests/testthat/test-ops64.R +++ b/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 + ) +})