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

### R code from vignette source 'Rmpfr-pkg.Rnw'

###################################################
### code chunk number 1: preliminaries
###################################################
options(SweaveHooks= list(fig=function() par(mar=c(5.1, 4.1, 1.1, 2.1))),
        width = 75,
        digits = 7, # <-- here, keep R's default!
        prompt = "R> ",
        continue="   ")
Sys.setenv(LANGUAGE = "en")
if(.Platform$OS.type != "windows")
  Sys.setlocale("LC_MESSAGES","C")


###################################################
### code chunk number 2: diagnose-lib
###################################################
if(nzchar(Sys.getenv("R_MM_PKG_CHECKING"))) print( .libPaths() )
stopifnot(require("sfsmisc"))


###################################################
### code chunk number 3: exp-1
###################################################
exp(1)


###################################################
### code chunk number 4: exp-1-dig-17
###################################################
print(exp(1), digits = 17)


###################################################
### code chunk number 5: exp-1-mp
###################################################
require("Rmpfr") # after having installed the package ...
(one <- mpfr(1, 120))
exp(one)


###################################################
### code chunk number 6: factorial-1
###################################################
ns <- 1:24 ; factorial(ns)


###################################################
### code chunk number 7: factorial-full
###################################################
noquote(sprintf("%-30.0f", factorial(24)))


###################################################
### code chunk number 8: factorial-mpfr
###################################################
ns <- mpfr(1:24, 120) ; factorial(ns)


###################################################
### code chunk number 9: chooseM-ex-fake (eval = FALSE)
###################################################
## chooseMpfr.all(n = 80)


###################################################
### code chunk number 10: chooseM-run
###################################################
capture.and.write(# <- in package 'sfsmisc': ~/R/Pkgs/sfsmisc/R/misc-goodies.R
chooseMpfr.all(n = 80)
                  , 5, 2, middle = 4, i.middle = 13)


###################################################
### code chunk number 11: ex1
###################################################
(0:7) / 7  #  k/7,  for  k= 0..7  printed with R's default precision
options(digits= 16)
(0:7) / 7  #  in full  double precision accuracy
options(digits=  7) # back to default
str(.Machine[c("double.digits","double.eps", "double.neg.eps")], digits=10)
2^-(52:53)


###################################################
### code chunk number 12: n-digs
###################################################
53 * log10(2)


###################################################
### code chunk number 13: ex1
###################################################
x <- mpfr(0:7, 80)/7 # using 80 bits precision
x
7*x
7*x  - 0:7


###################################################
### code chunk number 14: Const-names
###################################################
formals(Const)$name


###################################################
### code chunk number 15: Const-ex
###################################################
Const("pi")
Const("log2")


###################################################
### code chunk number 16: pi-1000
###################################################
system.time(Pi <- Const("pi", 1000 *log2(10)))
Pi


###################################################
### code chunk number 17: pi-fn-Gauss-HB
###################################################
piMpfr <- function(prec=256, itermax = 100, verbose=TRUE) {
    m2 <- mpfr(2, prec) # '2' as mpfr number
    ## -> all derived numbers are mpfr (with precision 'prec')
    p <- m2 + sqrt(m2) # 2 + sqrt(2) = 3.414..
    y <- sqrt(sqrt(m2)) # 2^ {1/4}
    x <- (y+1/y) / m2
    it <- 0L
    repeat {
	p.old <- p
	it <- it+1L
	p <- p * (1+x) / (1+y)
	if(verbose) cat(sprintf("it=%2d, pi^ = %s, |.-.|/|.|=%e\n",
				it, formatMpfr(p, min(50, prec/log2(10))), 1-p.old/p))
	if (abs(p-p.old) <= m2^(-prec))
	    break
	if(it > itermax) {
	    warning("not converged in", it, "iterations") ; break
	}
	## else
	s <- sqrt(x)
	y <- (y*s + 1/s) / (1+y)
	x <- (s+1/s)/2
    }
    p
}

piMpfr()# indeed converges  *quadratically* fast
## with relative error
relErr <- 1 - piMpfr(256, verbose=FALSE) / Const("pi",260)
## in bits :
asNumeric(-log2(abs(relErr)))


###################################################
### code chunk number 18: Math2-def
###################################################
getGroupMembers("Math2")
showMethods("Math2", classes=c("mpfr", "mpfrArray"))


###################################################
### code chunk number 19: round-ex
###################################################
i7 <- 1/mpfr(700, 100)
c(i7, round(i7, digits = 6), signif(i7, digits = 6))


###################################################
### code chunk number 20: roundMpfr-ex
###################################################
roundMpfr(i7, precBits = 30)
roundMpfr(i7, precBits = 15)


###################################################
### code chunk number 21: asNumeric-meth
###################################################
showMethods(asNumeric)


###################################################
### code chunk number 22: format-ex
###################################################
cbind( sapply(1:7, function(d) format(i7, digits=d)) )


###################################################
### code chunk number 23: format-lrg
###################################################
x <- mpfr(2, 80) ^ ((1:4)*10000)
cbind(x) # -> show() -> print.mpfr() -> formatMpfr(.. , digits = NULL, maybe.full = FALSE)
nchar(formatMpfr(x))
nchar(formatMpfr(x, maybe.full = TRUE))


###################################################
### code chunk number 24: Math-group
###################################################
getGroupMembers("Math")


###################################################
### code chunk number 25: Matrix-ex
###################################################
head(x <- mpfr(0:7, 64)/7) ; mx <-  x
dim(mx) <- c(4,2)


###################################################
### code chunk number 26: mpfrArr-ex
###################################################
dim(aa <- mpfrArray(1:24, precBits = 80, dim = 2:4))


###################################################
### code chunk number 27: pr-mpfrArr-fake (eval = FALSE)
###################################################
## aa


###################################################
### code chunk number 28: pr-mpfrArr-do
###################################################
capture.and.write(aa, 11, 4)


###################################################
### code chunk number 29: crossprod
###################################################
mx[ 1:3, ] + c(1,10,100)
crossprod(mx)


###################################################
### code chunk number 30: apply-mat
###################################################
apply(7 * mx, 2, sum)


###################################################
### code chunk number 31: Ei-curve
###################################################
getOption("SweaveHooks")[["fig"]]()
curve(Ei,  0, 5, n=2001);  abline(h=0,v=0, lty=3)


###################################################
### code chunk number 32: Li2-1
###################################################
if(mpfrVersion() >= "2.4.0")  ## Li2() is not available in older MPFR versions
  all.equal(Li2(1), Const("pi", 128)^2/6, tol = 1e-30)


###################################################
### code chunk number 33: Li2-curve
###################################################
getOption("SweaveHooks")[["fig"]]()
if(mpfrVersion() >= "2.4.0")
   curve(Li2, -2, 13,   n=2000); abline(h=0,v=0, lty=3)


###################################################
### code chunk number 34: erf-curves
###################################################
getOption("SweaveHooks")[["fig"]]()
curve(erf, -3,3, col = "red", ylim = c(-1,2))
curve(erfc, add = TRUE, col = "blue")
abline(h=0, v=0, lty=3); abline(v=c(-1,1), lty=3, lwd=.8, col="gray")
legend(-3,1, c("erf(x)", "erfc(x)"), col = c("red","blue"), lty=1)


###################################################
### code chunk number 35: integrateR-dnorm
###################################################
integrateR(dnorm,0,2000)
integrateR(dnorm,0,2000, rel.tol=1e-15)
integrateR(dnorm,0,2000, rel.tol=1e-15, verbose=TRUE)


###################################################
### code chunk number 36: integ-exp-double
###################################################
(Ie.d <- integrateR(exp,            0     , 1, rel.tol=1e-15, verbose=TRUE))


###################################################
### code chunk number 37: integ-exp-mpfr
###################################################
(Ie.m <- integrateR(exp, mpfr(0,200), 1, rel.tol=1e-25, verbose=TRUE))
(I.true <- exp(mpfr(1, 200)) - 1)
## with absolute errors
as.numeric(c(I.true - Ie.d$value,
             I.true - Ie.m$value))


###################################################
### code chunk number 38: integ-poly-double
###################################################
if(require("polynom")) {
    x <- polynomial(0:1)
    p <- (x-2)^4 - 3*(x-3)^2
    Fp <- as.function(p)
    print(pI <- integral(p)) # formally
    print(Itrue <- predict(pI, 5) - predict(pI, 0)) ## == 20
} else {
    Fp <- function(x) (x-2)^4 - 3*(x-3)^2
    Itrue <- 20
}
(Id <- integrateR(Fp, 0,      5))
(Im <- integrateR(Fp, 0, mpfr(5, 256),
                  rel.tol = 1e-70, verbose=TRUE))
## and the numerical errors, are indeed of the expected size:
256 * log10(2) # - expect ~ 77 digit accuracy for mpfr(*., 256)
as.numeric(Itrue - c(Im$value, Id$value))


###################################################
### code chunk number 39: pnorm-extr
###################################################
pnorm(-1234, log.p=TRUE)


###################################################
### code chunk number 40: pnorm-extr
###################################################
(p123 <- Rmpfr::pnorm(mpfr(-123, 66), log.p=TRUE)) # is based on
(ec123 <- erfc(123 * sqrt(mpfr(0.5, 66+4))) / 2) # 1.95....e-3288

(p333 <- Rmpfr::pnorm(mpfr(-333, 66), log.p=TRUE))
exp(p333)
stopifnot(p123 == log(roundMpfr(ec123, 66)), ## '==' as we implemented our pnorm() 
          all.equal(p333, -55451.22709, tol=1e-8))


###################################################
### code chunk number 41: mpfr-erange
###################################################
(old_erng <- .mpfr_erange() )


###################################################
### code chunk number 42: double-erange
###################################################
unlist( .Machine[c("double.min.exp", "double.max.exp")] )


###################################################
### code chunk number 43: really-min
###################################################
2^(-1022 - 52)


###################################################
### code chunk number 44: mpfr-all-eranges
###################################################
.mpfr_erange(.mpfr_erange_kinds) ## and then set
# use very slightly smaller than extreme values:
(myERng <- (1-2^-52) * .mpfr_erange(c("min.emin","max.emax")))
.mpfr_erange_set(value = myERng) # and to see what happened:
.mpfr_erange()


###################################################
### code chunk number 45: GMP-numbs
###################################################
.mpfr_gmp_numbbits()