tests/arith-ex.R

In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

require("Rmpfr")
## includes ("gmp")# want to check "mixed arithmetic" too  __ TODO __

`%=N=%` <- function(x,y) (x == y) | (is.na(x) & is.na(y))
all.EQ <- function(x,y, tolerance = 2^-98, ...) # very small tol. for MPFR
    all.equal(x, y, tolerance=tolerance, ...)
warningI <- function(...) warning(..., immediate. = TRUE)

unlist(.Platform)

## Check that we got the "which.*" methods also from "bigq":
bcl <- c("ANY", "bigq", "bigz", "mpfr")
##if(packageVersion("gmp") >= "0.5-8") {
stopifnot(identical(bcl,
		    sort(unlist(findMethods("which.max")@signatures))),
	  identical(bcl,
		    sort(unlist(findMethods("which.min")@signatures))))
##}

options(warn = 1)# warnings *immediately*
(doExtras <- Rmpfr:::doExtras())
eps2 <-   2 * .Machine$double.eps
eps8 <-   8 * .Machine$double.eps
eps32 <- 32 * .Machine$double.eps

## must take the *larger* of the two precisions:
stopifnot(substr(format(mpfr(1, 60)/mpfr(7, 160)), 1,51) == # format() may show more digits
          "0.1428571428571428571428571428571428571428571428571")# again has extra "2" at end

(x <- mpfr(0:7, 100) / 7)
ix <- x^-1000
iX <- asNumeric(ix)

stopifnot( mpfrIs0(x - x), # badly failed on 64-bit
	  identical(-x, 0-x),# testing "- x"
          all.equal(ix, (1/x)^1000, tol= 1e-25),
          is.numeric(iX), iX[1:4] == Inf, # failed previously as we used RNDD (downward rounding)
          all.equal(log(iX[5:8]), c(559.6157879, 336.4722366, 154.1506798, 0),
                    tol = 1e-9))

## checking hexadecimal input :
stopifnot(mpfr("0xFFFFFFFFFFFFFFFFFFFF", base=16) + 1 == 2^80,
## sign(0) == 0:
          identical(sign(as(-1:1, "mpfr")), -1:1 + 0))

stopifnot(all.equal(as.numeric(x+ 1L),
                    as.numeric(x)+1L, tol = eps2),
	  as.integer(  x [x < 1]) == 0,# was *wrong* {we round()ed; previously "down"!}
	  as.integer((-x)[x < 1]) == 0,#  (ditto)
          (3 * x)/3 <= x,
          all.equal(as.numeric(x * 2L),
                    as.numeric(x + x), tol = 0))

u <- mpfr(0:17, 128)/17
two <- mpfr(2,100)
stopifnot(all.EQ(u ^ two, u ^ 2),
          identical(u ^ 2, u ^ 2L),
          all.EQ(two ^ u, 2 ^ u),
          identical(2 ^ u, 2L ^ u),
          floor  (3*u) == floor  (3/17*(0:17)),
          ceiling(u*5) == ceiling(5/17*(0:17))
          )

i7  <- mpfr(0:7,  200)/ 7
i17 <- mpfr(0:17, 300)/17
stopifnot(all.equal(as.numeric(x+1),
		    as.numeric(x)+1),
	  all.equal(round(x,2), round(asNumeric(x), 2), tol=1e-15),
	  all.equal(round(mpfr(1.152, 80), 2), 1.15), # was wrong {as.integer() bug}
	  all.equal(0:7,   7 * round ( i7, 25), tol = 2e-25),
	  all.equal(0:7,   7 * round ( i7, 50), tol = 2e-50),
	  all.equal(0:17, 17 * signif(i17,100), tol = 2e-100),
	  all.equal(0:17, 17 * signif(i17, 20), tol = 2e-20)
          )

## When we compute with 100 bits,
## we should compare relative errors with  2^-100 :
del <- abs((x+pi)-pi - x) / 2^-100
stopifnot(del <= 4) ## <= 2 already
(fd <- format(del, drop0 = TRUE))
stopifnot(all.equal(as.numeric(del),
		    as.numeric(fd), tol = 1e-15))
if(print(Sys.info()[["machine"]]) == "x86_64")
    stopifnot(fd %in% as.character(c(0:2, c(2,7)/4)))


checkPmin <- function(x, nx = as(x, "numeric")) {
    rx <- if(is(x,"mpfr")) round(x, 25) else x
    isZ <- is(x, "bigz") || is(nx, "bigz")
    M.X <- max(x, na.rm=TRUE)
    m.x <- min(x, na.rm=TRUE)
    stopifnot(all.equal(x, nx),
	      pmin(x, x, M.X) %=N=% x, x %=N=% pmax(x, m.x, x),
	      all.equal(x, pmin(x, nx, x, M.X)),
	      all.equal(x, pmax(m.x, nx, x, rx, m.x)),
	      if(isZ)TRUE else all.equal(pmin(x, 0.75), pmin(nx, 0.75)),
	      if(isZ)TRUE else all.equal(pmax(x, 0.25), pmax(nx, 0.25)))
}

x <- mpfr(0:7, 100) / 7
checkPmin(x)

nx <- (0:7)/7
(qx <- as.bigq(0:7, 7))
 x[c(2,5)] <- NA
nx[c(2,5)] <- NA
qx[c(2,5)] <- NA

Z <- as.bigz(1:7)
mZ <- mpfr(Z, 64)
stopifnot(Z == mZ, mZ == Z)

checkPmin(x, nx)
    cat("checking pmin(. bigq ): ")
    ## FIXME checkPmin(x, qx); cat("[Ok]\n")
    ##
    print( base::pmin(Z, Z, max(Z)) )# via  gmp:: rep.bigz(x, length.out = *)
    cat("checking pmin(. bigz )
 [currently with lots of pmin() and pmax(...) warnings 'incompatible methods]:\n ")
    checkPmin(Z); cat("[Ok]\n") # via gmp:: all.equal.bigz()

stopifnot(all.equal( round(x, 10),  round(nx, 10)),
          all.equal(signif(x, 10), signif(nx, 10)))

## L & x ,  x & L  failed in Rmpfr 0.2* and 0.4-2
stopifnot(identical(L <- x > 0.5, L & x),
	  identical(L, x & L),
	  identical(x > 0, x | L))

## Summary() methods {including NA
stopifnot(exprs = {
    is.na(min(x))
    is.na(max(x))
    is.na(range(x))
    is.na(sum(x))
    is.na(prod(x))
    min(x, na.rm=TRUE) == 0
    max(x, na.rm=TRUE) == 1
    range(x, na.rm=TRUE) == 0:1
    all.equal(sum (x, na.rm=TRUE)*7,  2+3+5+6+7,  tolerance = 1e-28) # 1.0975e-30
    prod(x, na.rm=TRUE) == 0
    all.equal(180, prod(x[-1], na.rm=TRUE)*7^4, tol = 1e-15) # 1.579e-16
    ##
    ## all(), any()  had memory bug [PROTECT missing, but more, somehow]
    !all(x)
    is.na( all(x[-1]) )
    any(x)
    is.na(any(x[c(2,5)]))
    ## do these *twice* {that triggered R-forge bug  #6764 }
    ! all(x, na.rm=TRUE)
      any(x, na.rm=TRUE)
    ##
    ! all(x, na.rm=TRUE)
      any(x, na.rm=TRUE)
})

##-------------- Modulo and "integer division" -------------

## R's	?Arithmetic :
##
##   ‘%%’ indicates ‘x mod y’ and ‘%/%’ indicates integer division.  It
##   is guaranteed that ‘x == (x %% y) + y * ( x %/% y )’ (up to
##   rounding error) unless ‘y == 0’ where the result of ‘%%’ is
##   ‘NA_integer_’ or ‘NaN’ (depending on the ‘typeof’ of the
##   arguments).
##
## and has 'details' about  how	 non-integer 'y' works
##
(N <- if(doExtras) 1000 else 200)
(todays.seed <- eval(parse(text=Sys.Date())))# so this is reproducible
                                        # (and constant within one day)
set.seed(todays.seed)
mm <- c(-4:4, sample(50, N-9, replace=TRUE))
for(n in seq_len(N)) {
    cat("."); if(n %% 50 == 0) cat(n,"\n")
    m <- mm[n]
    prec <- sample(52:200, 1)# "high precision" ==> can use small tol
    x <- sample(100, 50) - 20
    for(kind in c('int','real')) {
	if(kind == "real") {
	    m <- jitter(m)
	    x <- jitter(x)
	    tol.1 <- eps32 * pmax(1, 1/abs(m))
	    EQ <- function(x,y, tol = tol.1)
		isTRUE(all.equal(x, as.numeric(y), tol=tol))
	    EQ2 <- function(x,y, tol = tol.1) {
		## for the DIV - MOD identity, a small x leads to cancellation
		all((x %=N=% y) | abs(x - y) < tol*pmax(abs(x), 1)) ||
		    isTRUE(all.equal(x, as.numeric(y), tol=tol))
	    }
	} else { ## "integer"
	    EQ2 <- EQ <- function(x,y, tol) all(x %=N=% y)
	}
	i.m <- mpfr(x, prec) %% mpfr(m, prec)
	if(!EQ2(x %% m, i.m)) {
	    cat("\n -- m = ",m,"  (prec = ",prec,")\n")
	    rE <- range(rel.E <- as.numeric(1 - (x %% m)/i.m))
	    print(cbind(x, 'R.%%' = x %% m, rel.E))
	    MSG <- if(max(abs(rE)) < 1e-10) warningI else stop
	    MSG(sprintf("not all equal: range(rel.Err.) = [%g, %g]", rE[1],rE[2]))
	}
	##
	if(m != 0) {
	    ##---Check the    x	 ==  (x %% m) +	 m * ( x %/% m )  assertion ------
	    ##
	    if(EQ2(x, (x %% m) + m*( x %/% m ), tol = 1e-12)) { ## ok for R
		## --> also ok for mpfr ?
		iDm <- mpfr(x, prec) %/% mpfr(m, prec)
		rhs <- i.m + m*iDm
		if(!EQ2(x, i.m + m*iDm)) {
		    cat("\n -- m = ",m,"  (prec = ",prec,")\n")
		    print(cbind(x,' MPFR[ x%%m + m(x %/% m) ]' = as.numeric(rhs), rel.E))
		    MSG <- if(max(abs(rE)) < 1e-10) warningI else stop
		    MSG(sprintf("Identity(MOD - DIV) not all eq.: range(rel.Err.) = [%g, %g]",
				rE[1],rE[2]))
		}
	    } else {
		cat("\n hmm.. the basic	 %% <-> %/%  assertion 'fails' in *R* :\n")
		rhs <- (x %% m)	 +  m * ( x %/% m )
		rel.E <- (1 - rhs/x)
		print(cbind(x, 'x%%m + m(x %/% m)' = rhs, rel.E))
	    }
	}
    }
}

## mpfr  o  <number>  now implemented, for  '%%', too :
r <- as.double(i <- -10:20)

stopifnot(
    ## %% -------------------------------------
    mpfr(i, prec=99) %% 7  == i %% 7
    , ##
    mpfr(i, prec=99) %% 7  ==
    mpfr(i, prec=99) %% 7L
    , ##
    i %% mpfr(27, prec=99) == i %% 27
    , ##
    r %% mpfr(27, prec=99) == r %% 27
    , ## %/% -------------------------------------
    mpfr(i, prec=99) %/% 7  == i %/% 7
    , ##
    mpfr(i, prec=99) %/% 7  ==
    mpfr(i, prec=99) %/% 7L
    , ##
    mpfr(i, prec=99) %/% mpfr(27, prec=99) == i %/% 27
    , ##
    i %/% mpfr(27, prec=99) == i %/% 27
    , ##
    i %/% mpfr(27, prec=99) ==
    r %/% mpfr(27, prec=99)
    , TRUE ##
)

cat('Time elapsed: ', proc.time(),'\n') # "stats"

## Was reproducible BUG in Rmpfr-addition (on Linux, MPFR 4.x.y) -- ## but the bug was Rmpfr,
## in ../src/Ops.c, detecting if *integer*, i.e., long can be used
dn <- 1e20
dOO <- 9223372036854775808; formatC(dOO) # "9.2...e18"
(r <- dn / (dn + dOO)) # 0.915555  (double prec arithmetic)
## but *so* strange when switching to Rmpfr :  addition accidentally *subtract*!!
n <- mpfr(dn, precBits = 99)
(rM <- n / (n + dOO)) # wrongly gave " 1 'mpfr' .... 99 bits; 1.101605140483951.....
stopifnot(exprs = {
    all.equal(n + dOO, dn + dOO)
    all.equal(n / (n + dOO), r)
})

##  log(., base)  :
(ten40 <- as.bigz(10)^40)
ten40m <- mpfr(ten40)
(lt40 <- log(ten40m, 10)) # gave Error in ... : base != exp(1) is not yet implemented
##  'mpfr' ..  133 bits \\ [1] 40
stopifnot(exprs = {
    grepl("^40[.]000+$", print(format(lt40, digits = 60)))
    identical(lt40, log10(ten40m))
    identical(log(ten40m, 2), log2(ten40m))
    inherits(Pi <- Const("pi", 140), "mpfr")
    all.equal(show(log(ten40m, Pi)),
                   log(ten40m)/log(Pi), tol = 1e-40)
})



###------Standard Statistics Functions --------------------------------------------------------

x <- c(del, 1000)
stopifnot(identical(mean(x), mean(x, trim=0)))
for(tr in (0:8)/16)
    stopifnot(all.equal(mean(          x,  trim = tr),
                        mean(asNumeric(x), trim = tr), tol=1e-15))

cat('Time elapsed: ', proc.time(),'\n') # "stats"