Package 'Rmpfr'

R MPFR - Multiple Precision Floating-Point Reliable

Description

Rmpfr provides S4 classes and methods for arithmetic including transcendental ("special") functions for arbitrary precision floating point numbers, here often called “mpfr - numbers”. To this end, it interfaces to the LGPL'ed MPFR (Multiple Precision Floating-Point Reliable) Library which itself is based on the GMP (GNU Multiple Precision) Library.

Details

Package:	Rmpfr
Title:	Interface R to MPFR - Multiple Precision Floating-Point Reliable
Version:	0.9-6
Date:	2024-08-20
DateNote:	Previous CRAN version 0.9-5 on 2024-01-20
Type:	Package
Authors@R:	c(person("Martin","Maechler", role = c("aut","cre"), email = "[email protected]", comment = c(ORCID="0000-0002-8685-9910")) , person(c("Richard", "M."), "Heiberger", role = "ctb", email="[email protected]", comment = "formatHex(), *Bin, *Dec") , person(c("John", "C."), "Nash", role = "ctb", email="[email protected]", comment = "hjkMpfr(), origin of unirootR()") , person(c("Hans", "W."), "Borchers", role = "ctb", email="[email protected]", comment = "optimizeR(*, \"GoldenRatio\"); origin of hjkMpfr()") , person("Mikael", "Jagan", role = "ctb", comment = c("safer convert.c", ORCID = "0000-0002-3542-2938")) )
Description:	Arithmetic (via S4 classes and methods) for arbitrary precision floating point numbers, including transcendental ("special") functions. To this end, the package interfaces to the 'LGPL' licensed 'MPFR' (Multiple Precision Floating-Point Reliable) Library which itself is based on the 'GMP' (GNU Multiple Precision) Library.
SystemRequirements:	gmp (>= 4.2.3), mpfr (>= 3.1.0), pdfcrop (part of TexLive) is required to rebuild the vignettes.
SystemRequirementsNote:	'MPFR' (MP Floating-Point Reliable Library, https://www.mpfr.org/) and 'GMP' (GNU Multiple Precision library, https://gmplib.org/), see >> README.md
Depends:	gmp (>= 0.6-1), R (>= 3.6.0)
Imports:	stats, utils, methods
Suggests:	DPQmpfr, MASS, Bessel, polynom, sfsmisc (>= 1.1-14)
SuggestsNote:	MASS, polynom, sfsmisc: only for vignette;
Enhances:	dfoptim, pracma, DPQ
EnhancesNote:	mentioned in Rd xrefs | used in example
URL:	https://rmpfr.r-forge.r-project.org/
BugReports:	https://r-forge.r-project.org/tracker/?group_id=386
License:	GPL (>= 2)
Encoding:	UTF-8
Repository:	https://r-forge.r-universe.dev
RemoteUrl:	https://github.com/r-forge/rmpfr
RemoteRef:	HEAD
RemoteSha:	4f94e1132cf56bb50f5756683e935bb0a532d3e1
Author:	Martin Maechler [aut, cre] (<https://orcid.org/0000-0002-8685-9910>), Richard M. Heiberger [ctb] (formatHex(), *Bin, *Dec), John C. Nash [ctb] (hjkMpfr(), origin of unirootR()), Hans W. Borchers [ctb] (optimizeR(*, "GoldenRatio"); origin of hjkMpfr()), Mikael Jagan [ctb] (safer convert.c, <https://orcid.org/0000-0002-3542-2938>)
Maintainer:	Martin Maechler <[email protected]>

Index of help topics:

.bigq2mpfr              Conversion Utilities gmp <-> Rmpfr
Bernoulli               Bernoulli Numbers in Arbitrary Precision
Bessel_mpfr             Bessel functions of Integer Order in multiple
                        precisions
Mnumber-class           Class "Mnumber" and "mNumber" of "mpfr" and
                        regular numbers and arrays from them
Rmpfr-package           R MPFR - Multiple Precision Floating-Point
                        Reliable
array_or_vector-class   Auxiliary Class "array_or_vector"
asNumeric-methods       Methods for 'asNumeric(<mpfr>)'
atomicVector-class      Virtual Class "atomicVector" of Atomic Vectors
c.mpfr                  MPFR Number Utilities
cbind                   "mpfr" '...' - Methods for Functions cbind(),
                        rbind()
chooseMpfr              Binomial Coefficients and Pochhammer Symbol aka
                        Rising Factorial
determinant.mpfrMatrix
                        Functions for mpfrMatrix Objects
factorialMpfr           Factorial 'n!'  in Arbitrary Precision
formatHex               Flexibly Format Numbers in Binary, Hex and
                        Decimal Format
formatMpfr              Formatting MPFR (multiprecision) Numbers
frexpMpfr               Base-2 Representation and Multiplication of
                        Mpfr Numbers
getPrec                 Rmpfr - Utilities for Precision Setting,
                        Printing, etc
hjkMpfr                 Hooke-Jeeves Derivative-Free Minimization R
                        (working for MPFR)
igamma                  Incomplete Gamma Function
integrateR              One-Dimensional Numerical Integration - in pure
                        R
is.whole.mpfr           Whole ("Integer") Numbers
log1mexp                Compute f(a) = log(1 +/- exp(-a)) Numerically
                        Optimally
matmult                 (MPFR) Matrix (Vector) Multiplication
mpfr                    Create "mpfr" Numbers (Objects)
mpfr-class              Class "mpfr" of Multiple Precision Floating
                        Point Numbers
mpfrArray               Construct "mpfrArray" almost as by 'array()'
mpfrMatrix-class        Classes "mpfrMatrix" and "mpfrArray"
num2bigq                BigQ / BigRational Approximation of Numbers
optimizeR               High Precision One-Dimensional Optimization
outer                   Base Functions etc, as an Rmpfr version
pbetaI                  Accurate Incomplete Beta / Beta Probabilities
                        For Integer Shapes
pmax                    Parallel Maxima and Minima
pnorm                   Distribution Functions with MPFR Arithmetic
qnormI                  Gaussian / Normal Quantiles 'qnorm()' via
                        Inversion
roundMpfr               Rounding to Binary bits, "mpfr-internally"
sapplyMpfr              Apply a Function over a "mpfr" Vector
seqMpfr                 "mpfr" Sequence Generation
str.mpfr                Compactly Show STRucture of Rmpfr Number Object
sumBinomMpfr            (Alternating) Binomial Sums via Rmpfr
unirootR                One Dimensional Root (Zero) Finding - in pure R
zeta                    Special Mathematical Functions (MPFR)

The following (help pages) index does not really mention that we provide many methods for mathematical functions, including gamma, digamma, etc, namely, all of R's (S4) Math group (with the only exception of trigamma), see the list in the examples. Additionally also pnorm, the “error function”, and more, see the list in zeta, and further note the first vignette (below).

Partial index:

mpfr	Create "mpfr" Numbers (Objects)
mpfrArray	Construct "mpfrArray" almost as by array()
mpfr-class	Class "mpfr" of Multiple Precision Floating Point Numbers
mpfrMatrix-class	Classes "mpfrMatrix" and "mpfrArray"
Bernoulli	Bernoulli Numbers in Arbitrary Precision
Bessel_mpfr	Bessel functions of Integer Order in multiple precisions
c.mpfr	MPFR Number Utilities
cbind	"mpfr" ... - Methods for Functions cbind(), rbind()
chooseMpfr	Binomial Coefficients and Pochhammer Symbol aka
	Rising Factorial
factorialMpfr	Factorial 'n!' in Arbitrary Precision
formatMpfr	Formatting MPFR (multiprecision) Numbers
getPrec	Rmpfr - Utilities for Precision Setting, Printing, etc
roundMpfr	Rounding to Binary bits, "mpfr-internally"
seqMpfr	"mpfr" Sequence Generation
sumBinomMpfr	(Alternating) Binomial Sums via Rmpfr
zeta	Special Mathematical Functions (MPFR)
integrateR	One-Dimensional Numerical Integration - in pure R
unirootR	One Dimensional Root (Zero) Finding - in pure R
optimizeR	High Precisione One-Dimensional Optimization
hjkMpfr	Hooke-Jeeves Derivative-Free Minimization R (working for MPFR)

Further information is available in the following vignettes:

Rmpfr-pkg	Arbitrarily Accurate Computation with R: The 'Rmpfr' package (source, pdf)
log1mexp-note	Acccurately Computing log(1 - exp(.)) -- Assessed by Rmpfr (source, pdf)

Author(s)

Martin Maechler

References

MPFR (MP Floating-Point Reliable Library), https://www.mpfr.org/

GMP (GNU Multiple Precision library), https://gmplib.org/

and see the vignettes mentioned above.

See Also

The R package gmp for big integer gmp and rational numbers (bigrational) on which Rmpfr depends.

Examples

## Using  "mpfr" numbers instead of regular numbers...
n1.25 <- mpfr(5, precBits = 256)/4
n1.25

## and then "everything" just works with the desired chosen precision:hig
n1.25 ^ c(1:7, 20, 30) ## fully precise; compare with
print(1.25 ^ 30, digits=19)

exp(n1.25)

## Show all math functions which work with "MPFR" numbers (1 exception: trigamma)
getGroupMembers("Math")

## We provide *many* arithmetic, special function, and other methods:
showMethods(classes = "mpfr")
showMethods(classes = "mpfrArray")

---

Auxiliary Class "array_or_vector"

Description

"array_or_vector" is the class union of c("array", "matrix", "vector") and exists for its use in signatures of method definitions.

Details

Using "array_or_vector" instead of just "vector" in a signature makes an important difference: E.g., if we had setMethod(crossprod, c(x="mpfr", y="vector"), function(x,y) CPR(x,y)), a call crossprod(x, matrix(1:6, 2,3)) would extend into a call of CPR(x, as(y, "vector")) such that CPR()'s second argument would simply be a vector instead of the desired $2\times3$ matrix.

Objects from the Class

A virtual Class: No objects may be created from it.

Examples

showClass("array_or_vector")

---

Methods for asNumeric(<mpfr>)

Description

Methods for function asNumeric (in package gmp).

Usage

## S4 method for signature 'mpfrArray'
asNumeric(x)

Arguments

x	a “number-like” object, here, a mpfr or typically mpfrArrayone.

Value

an R object of type (typeof) "numeric", a matrix or array if x had non-NULL dimension dim().

Methods

signature(x = "mpfrArray")

this method also dispatches for mpfrMatrix and returns a numeric array.

signature(x = "mpfr")

for non-array/matrix, asNumeric(x) is basically the same as as.numeric(x).

Author(s)

Martin Maechler

See Also

our lower level (non-generic) toNum(). Further, asNumeric (package gmp), standard R's as.numeric().

Examples

x <- (0:7)/8 # (exact)
X <- mpfr(x, 99)
stopifnot(identical(asNumeric(x), x),
	  identical(asNumeric(X), x))

m <- matrix(1:6, 3,2)
(M <- mpfr(m, 99) / 5) ##-> "mpfrMatrix"
asNumeric(M) # numeric matrix
stopifnot(all.equal(asNumeric(M), m/5),
          identical(asNumeric(m), m))# remains matrix

---

Virtual Class "atomicVector" of Atomic Vectors

Description

The class "atomicVector" is a virtual class containing all atomic vector classes of base R, as also implicitly defined via is.atomic.

Objects from the Class

A virtual Class: No objects may be created from it.

Methods

In the Matrix package, the "atomicVector" is used in signatures where typically “old-style” "matrix" objects can be used and can be substituted by simple vectors.

Extends

The atomic classes "logical", "integer", "double", "numeric", "complex", "raw" and "character" are extended directly. Note that "numeric" already contains "integer" and "double", but we want all of them to be direct subclasses of "atomicVector".

Author(s)

Martin Maechler

See Also

is.atomic, integer, numeric, complex, etc.

Examples

showClass("atomicVector")

---

Bernoulli Numbers in Arbitrary Precision

Description

Computes the Bernoulli numbers in the desired (binary) precision. The computation happens via the zeta function and the formula

$B_k = -k \zeta(1 - k),$

and hence the only non-zero odd Bernoulli number is $B_1 = +1/2$. (Another tradition defines it, equally sensibly, as $-1/2$.)

Usage

Bernoulli(k, precBits = 128)

Arguments

k	non-negative integer vector
precBits	the precision in bits desired.

Value

an mpfr class vector of the same length as k, with i-th component the k[i]-th Bernoulli number.

Author(s)

Martin Maechler

References

https://en.wikipedia.org/wiki/Bernoulli_number

See Also

zeta is used to compute them.

The next version of package gmp is to contain BernoulliQ(), providing exact Bernoulli numbers as big rationals (class "bigq").

Examples

Bernoulli(0:10)
plot(as.numeric(Bernoulli(0:15)), type = "h")

curve(-x*zeta(1-x), -.2, 15.03, n=300,
      main = expression(-x %.% zeta(1-x)))
legend("top", paste(c("even","odd  "), "Bernoulli numbers"),
       pch=c(1,3), col=2, pt.cex=2, inset=1/64)
abline(h=0,v=0, lty=3, col="gray")
k <- 0:15; k[1] <- 1e-4
points(k, -k*zeta(1-k), col=2, cex=2, pch=1+2*(k%%2))

## They pretty much explode for larger k :
k2 <- 2*(1:120)
plot(k2, abs(as.numeric(Bernoulli(k2))), log = "y")
title("Bernoulli numbers exponential growth")

Bernoulli(10000)# - 9.0494239636 * 10^27677

---

Bessel functions of Integer Order in multiple precisions

Description

Bessel functions of integer orders, provided via arbitrary precision algorithms from the MPFR library.

Note that the computation can be very slow when n and x are large (and of similar magnitude).

Usage

Ai(x)
j0(x)
j1(x)
jn(n, x, rnd.mode = c("N","D","U","Z","A"))
y0(x)
y1(x)
yn(n, x, rnd.mode = c("N","D","U","Z","A"))

Arguments

x	a numeric or mpfr vector.
n	non-negative integer (vector).
rnd.mode	a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see mpfr.

Value

Computes multiple precision versions of the Bessel functions of integer order, $J_n(x)$ and $Y_n(x)$, and—when using MPFR library 3.0.0 or newer—also of the Airy function $Ai(x)$. Note that currently Ai(x) is very slow to compute for large x.

See Also

besselJ, and besselY compute the same bessel functions but for arbitrary real order and only precision of a bit more than ten digits.

Examples

x <- (0:100)/8 # (have exact binary representation)
stopifnot(exprs = {
   all.equal(besselY(x, 0), bY0 <- y0(x))
   all.equal(besselJ(x, 1), bJ1 <- j1(x))
   all.equal(yn(0,x), bY0)
   all.equal(jn(1,x), bJ1)
})

mpfrVersion() # now typically 4.1.0
if(mpfrVersion() >= "3.0.0") { ## Ai() not available previously
  print( aix <- Ai(x) )
  plot(x, aix, log="y", type="l", col=2)
  stopifnot(
    all.equal(Ai (0) , 1/(3^(2/3) * gamma(2/3)))
    , # see https://dlmf.nist.gov/9.2.ii
    all.equal(Ai(100), mpfr("2.6344821520881844895505525695264981561e-291"), tol=1e-37)
  )
  two3rd <- 2/mpfr(3, 144)
  print( all.equal(Ai(0), 1/(3^two3rd * gamma(two3rd)), tol=0) ) # 1.7....e-40
  if(Rmpfr:::doExtras()) withAutoprint({ # slowish:
     system.time(ai1k <- Ai(1000)) # 1.4 sec (on 2017 lynne)
     stopifnot(all.equal(print(log10(ai1k)),
                         -9157.031193409585185582, tol=2e-16)) # seen 8.8..e-17 | 1.1..e-16
  })
} # ver >= 3.0

---

"mpfr" '...' - Methods for Functions cbind(), rbind()

Description

cbind and rbind methods for signature ... (see dotsMethods are provided for class Mnumber, i.e., for binding numeric vectors and class "mpfr" vectors and matrices ("mpfrMatrix") together.

Usage

cbind(..., deparse.level = 1)
rbind(..., deparse.level = 1)

Arguments

...	matrix-/vector-like R objects to be bound together, see the base documentation, cbind.
deparse.level	integer determining under which circumstances column and row names are built from the actual arguments' ‘expression’, see cbind.

Value

typically a ‘matrix-like’ object, here typically of class "mpfrMatrix".

Methods

... = "Mnumber"

is used to (c|r)bind multiprecision “numbers” (inheriting from class "mpfr") together, maybe combined with simple numeric vectors.

... = "ANY"

reverts to cbind and rbind from package base.

Author(s)

Martin Maechler

See Also

cbind2, cbind, Documentation in base R's methods package

Examples

cbind(1, mpfr(6:3, 70)/7, 3:0)

---

Binomial Coefficients and Pochhammer Symbol aka Rising Factorial

Description

Compute binomial coefficients, chooseMpfr(a,n) being mathematically the same as choose(a,n), but using high precision (MPFR) arithmetic.

chooseMpfr.all(n) means the vector choose(n, 1:n), using enough bits for exact computation via MPFR. However, chooseMpfr.all() is now deprecated in favor of chooseZ from package gmp, as that is now vectorized.

pochMpfr() computes the Pochhammer symbol or “rising factorial”, also called the “Pochhammer function”, “Pochhammer polynomial”, “ascending factorial”, “rising sequential product” or “upper factorial”,

$x^{(n)}=x(x+1)(x+2)\cdots(x+n-1)= \frac{(x+n-1)!}{(x-1)!} = \frac{\Gamma(x+n)}{\Gamma(x)}.$

Usage

chooseMpfr (a, n, rnd.mode = c("N","D","U","Z","A"))
chooseMpfr.all(n, precBits=NULL, k0=1, alternating=FALSE)
pochMpfr(a, n, rnd.mode = c("N","D","U","Z","A"))

Arguments

a	a numeric or mpfr vector.
n	an integer vector; if not of length one, n and a are recycled to the same length.
rnd.mode	a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see mpfr.
precBits	integer or NULL for increasing the default precision of the result.
k0	integer scalar
alternating	logical, for chooseMpfr.all(), indicating if alternating sign coefficients should be returned, see below.

Value

For

chooseMpfr(), pochMpfr():

an mpfr vector of length max(length(a), length(n));

chooseMpfr.all(n, k0):

a mpfr vector of length n-k0+1, of binomial coefficients $C_{n,m}$ or, if alternating is true, $(-1)^m\cdot C_{n,m}$ for $m \in$ k0:n.

Note

Currently this works via a (C level) for(i in 1:n)-loop which really slow for large n, say $10^6$, with computational cost $O(n^2)$. In such cases, if you need high precision choose(a,n) (or Pochhammer(a,n)) for large n, preferably work with the corresponding factorial(mpfr(..)), or gamma(mpfr(..)) terms.

See Also

choose(n,m) (base R) computes the binomial coefficient $C_{n,m}$ which can also be expressed via Pochhammer symbol as $C_{n,m} = (n-m+1)^{(m)}/m!$.

chooseZ from package gmp; for now, factorialMpfr.

For (alternating) binomial sums, directly use sumBinomMpfr, as that is potentially more efficient.

Examples

pochMpfr(100, 4) == 100*101*102*103 # TRUE
a <- 100:110
pochMpfr(a, 10) # exact (but too high precision)
x <- mpfr(a, 70)# should be enough
(px <- pochMpfr(x, 10)) # the same as above (needing only 70 bits)
stopifnot(pochMpfr(a, 10) == px,
          px[1] ==prod(mpfr(100:109, 100)))# used to fail

(c1 <- chooseMpfr(1000:997, 60)) # -> automatic "correct" precision
stopifnot(all.equal(c1, choose(1000:997, 60), tolerance=1e-12))

## --- Experimenting & Checking
n.set <- c(1:10, 20, 50:55, 100:105, 200:203, 300:303, 500:503,
           699:702, 999:1001)
if(!Rmpfr:::doExtras()) { ## speed up: smaller set
  n. <- n.set[-(1:10)]
  n.set <- c(1:10, n.[ c(TRUE, diff(n.) > 1)])
}
C1 <- C2 <- numeric(length(n.set))
for(i.n in seq_along(n.set)) {
  cat(n <- n.set[i.n],":")
  C1[i.n] <- system.time(c.c <- chooseMpfr.all(n) )[1]
  C2[i.n] <- system.time(c.2 <- chooseMpfr(n, 1:n))[1]
  stopifnot(is.whole(c.c), c.c == c.2,
            if(n > 60) TRUE else all.equal(c.c, choose(n, 1:n), tolerance = 1e-15))
  cat(" [Ok]\n")
}
matplot(n.set, cbind(C1,C2), type="b", log="xy",
        xlab = "n", ylab = "system.time(.)  [s]")
legend("topleft", c("chooseMpfr.all(n)", "chooseMpfr(n, 1:n)"),
       pch=as.character(1:2), col=1:2, lty=1:2, bty="n")

## Currently, chooseMpfr.all() is faster only for large n (>= 300)
## That would change if we used C-code for the *.all() version

## If you want to measure more:
measureMore <- TRUE
measureMore <- FALSE
if(measureMore) { ## takes ~ 2 minutes (on "lynne", Intel i7-7700T, ~2019)
  n.s <- 2^(5:20)
  r <- lapply(n.s, function(n) {
      N <- ceiling(10000/n)
      cat(sprintf("n=%9g => N=%d: ",n,N))
      ct <- system.time(C <- replicate(N, chooseMpfr(n, n/2)))
      cat("[Ok]\n")
      list(C=C, ct=ct/N)
  })
  print(ct.n <- t(sapply(r, `[[`, "ct")))
  hasSfS <- requireNamespace("sfsmisc")
  plot(ct.n[,"user.self"] ~ n.s, xlab=quote(n), ylab="system.time(.) [s]",
       main = "CPU Time for  chooseMpfr(n, n/2)",
       log ="xy", type = "b", axes = !hasSfS)
  if(hasSfS) for(side in 1:2) sfsmisc::eaxis(side)
  summary(fm <- lm(log(ct.n[,"user.self"]) ~ log(n.s), subset = n.s >= 10^4))
  ## --> slope ~= 2  ==>  It's  O(n^2)
  nn <- 2^seq(11,21, by=1/16) ; Lcol <- adjustcolor(2, 1/2)
  bet <- coef(fm)
  lines(nn, exp(predict(fm, list(n.s = nn))), col=Lcol, lwd=3)
  text(500000,1, substitute(AA %*% n^EE,
                            list(AA = signif(exp(bet[1]),3),
                                 EE = signif(    bet[2], 3))), col=2)
} # measure more

---

Factorial 'n!' in Arbitrary Precision

Description

Efficiently compute $n!$ in arbitrary precision, using the MPFR-internal implementation. This is mathematically (but not numerically) the same as $\Gamma(n+1)$.

factorialZ (package gmp) should typically be used instead of factorialMpfr() nowadays. Hence, factorialMpfr now is somewhat deprecated.

Usage

factorialMpfr(n, precBits = max(2, ceiling(lgamma(n+1)/log(2))),
              rnd.mode = c("N","D","U","Z","A"))

Arguments

n	non-negative integer (vector).
precBits	desired precision in bits (“binary digits”); the default sets the precision high enough for the result to be exact.
rnd.mode	a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see mpfr.

Value

a number of (S4) class mpfr.

See Also

factorial and gamma in base R.

factorialZ (package gmp), to replace factorialMpfr, see above.

chooseMpfr() and pochMpfr() (on the same page).

Examples

factorialMpfr(200)

n <- 1000:1010
f1000 <- factorialMpfr(n)
stopifnot(1e-15 > abs(as.numeric(1 - lfactorial(n)/log(f1000))))

## Note that---astonishingly--- measurements show only
## *small* efficiency gain of ~ 10% : over using the previous "technique"
system.time(replicate(8, f1e4 <- factorialMpfr(10000)))
system.time(replicate(8, f.1e4 <- factorial(mpfr(10000,
                            prec=1+lfactorial(10000)/log(2)))))

---

Flexibly Format Numbers in Binary, Hex and Decimal Format

Description

Show numbers in binary, hex and decimal format. The resulting character-like objects can be back-transformed to "mpfr" numbers via mpfr().

Usage

formatHex(x, precBits = min(getPrec(x)), style = "+", expAlign = TRUE)

formatBin(x, precBits = min(getPrec(x)), scientific = TRUE,
          left.pad = "_", right.pad = left.pad, style = "+", expAlign = TRUE)
formatDec(x, precBits = min(getPrec(x)), digits = decdigits,
          nsmall = NULL, scientific = FALSE, style = "+",
          decimalPointAlign = TRUE, ...)

Arguments

x	a numeric or mpfr R object.
precBits	integer, the number of bits of precision, typically derived from x, see getPrec. Numeric, i.e., double precision numbers have 53 bits. For more detail, see mpfr.
style	a single character, to be used in sprintf's format (fmt), immediately after the " sets a sign in the output, i.e., "+" or "-", where as style = " " may seem more standard.
expAlign	logical indicating if for scientific (“exponential”) representations the exponents should be aligned to the same width, i.e., zero-padded to the same number of digits.
scientific	logical indicating that formatBin should display the binary representation in scientific notation (mpfr(3, 5) is displayed as +0b1.1000p+1). When FALSE, formatBin will display the binary representation in regular format shifted to align binary points (mpfr(3, 5) is displayed +0b11.000).
...	additional optional arguments. formatHex, formatBin: precBits is the only ... argument acted on. Other ... arguments are ignored. formatDec: precBits is acted on. Any argument accepted by format (except nsmall) is acted on. Other ... arguments are ignored.
left.pad, right.pad	characters (one-character strings) that will be used for left- and right-padding of the formatted string when scientific=FALSE. Do not change these unless for display-only purpose !!
nsmall	only used when scientific is false, then passed to format(). If NULL, the default is computed from the range of the non-zero values of x.
digits	integer; the number of decimal digits displayed is the larger of this argument and the internally generated value that is a function of precBits. This is related to but different than digits in format.
decimalPointAlign	logical indicating if padding should be used to ensure that the resulting strings align on the decimal point (".").

Details

For the hexadecimal representation, when the precision is not larger than double precision, sprintf() is used directly, otherwise formatMpfr() is used and post processed. For the binary representation, the hexadecimal value is calculated and then edited by substitution of the binary representation of the hex characters coded in the HextoBin vector. For binary with scientific=FALSE, the result of the scientific=TRUE version is edited to align binary points. For the decimal representation, the hexadecimal value is calculated with the specified precision and then sent to the format function for scientific=FALSE or to the sprintf function for scientific=TRUE.

Value

a character vector (or matrix) like x, say r, containing the formatted represention of x, with a class (unless left.pad or right.pad were not "_"). In that case, formatHex() and formatBin() return class "Ncharacter"; for that, mpfr(.) has a method and will basically return x, i.e., work as inverse function.

Since Rmpfr version 0.6-2, the S3 class "Ncharacter" extends "character". (class(.) is now of length two and class(.)[2] is "character".). There are simple [ and print methods; modifying or setting dim works as well.

Author(s)

Richard M. Heiberger [email protected], with minor tweaking by Martin M.

References

R FAQ 7.31: Why doesn't R think these numbers are equal? system.file("../../doc/FAQ")

See Also

mpfr, sprintf

Examples

FourBits <- mpfr(matrix(0:31, 8, 4, dimnames = list(0:7, c(0,8,16,24))),
                 precBits=4) ## 4 significant bits
FourBits

formatHex(FourBits)
formatBin(FourBits, style = " ")
formatBin(FourBits, scientific=FALSE)
formatDec(FourBits)

## as "Ncharacter"  'inherits from' "character", this now works too :
data.frame(Dec = c( formatDec(FourBits) ), formatHex(FourBits),
           Bin = formatBin(FourBits, style = " "))

FBB <- formatBin(FourBits) ; clB <- class(FBB)
(nFBB <- mpfr(FBB))
stopifnot(class(FBB)[1] == "Ncharacter",
          all.equal(nFBB, FourBits, tol=0))

FBH <- formatHex(FourBits) ; clH <- class(FBH)
(nFBH <- mpfr(FBH))
stopifnot(class(FBH)[1] == "Ncharacter",
          all.equal(nFBH, FourBits, tol=0))

## Compare the different "formattings"  (details will change, i.e. improve!)%% FIXME
M <- mpfr(c(-Inf, -1.25, 1/(-Inf), NA, 0, .5, 1:2, Inf), 3)
data.frame(fH = formatHex(M), f16 = format(M, base=16),
           fB = formatBin(M), f2  = format(M, base= 2),
           fD = formatDec(M), f10 = format(M), # base = 10 is default
           fSci= format(M, scientific=TRUE),
           fFix= format(M, scientific=FALSE))

## Other methods ("[", t()) also work :
stopifnot(dim(F1 <- FBB[, 1, drop=FALSE]) == c(8,1), identical(class(  F1), clB),
          dim(t(F1)) == c(1,8),                      identical(class(t(F1)),clB),
          is.null(dim(F.2 <- FBB[,2])),              identical(class( F.2), clB),
          dim(F22 <- FBH[1:2, 3:4]) == c(2,2), identical(class(F22),  clH),
          identical(class(FBH[2,3]), clH), is.null(dim(FBH[2,3])),
          identical(FBH[2,3:4], F22[2,]),
          identical(FBH[2,3], unname(FBH[,3][2])),
          TRUE)

TenFrac <- matrix((1:10)/10, dimnames=list(1:10, expression(1/x)))
TenFrac9 <- mpfr(TenFrac, precBits=9) ## 9 significant bits
TenFrac9
formatHex(TenFrac9)
formatBin(TenFrac9)
formatBin(TenFrac9, scientific=FALSE)
formatDec(TenFrac9)
formatDec(TenFrac9, precBits=10)

---

Formatting MPFR (multiprecision) Numbers

Description

Flexible formatting of “multiprecision numbers”, i.e., objects of class mpfr. formatMpfr() is also the mpfr method of the generic format function.

The formatN() methods for mpfr numbers renders them differently than their double precision equivalents, by appending "_M".

Function .mpfr2str() is the low level work horse for formatMpfr() and hence all print()ing of "mpfr" objects.

Usage

formatMpfr(x, digits = NULL, trim = FALSE, scientific = NA,
           maybe.full = (!is.null(digits) && is.na(scientific)) || isFALSE(scientific),
           base = 10, showNeg0 = TRUE, max.digits = Inf,
           big.mark = "", big.interval = 3L,
           small.mark = "", small.interval = 5L,
           decimal.mark = ".",
           exponent.char = if(base <= 14) "e" else if(base <= 36) "E" else "|e",
           exponent.plus = TRUE,
           zero.print = NULL, drop0trailing = FALSE, ...)

## S3 method for class 'mpfr'
formatN(x, drop0trailing = TRUE, ...)

.mpfr2str(x, digits = NULL, maybe.full = !is.null(digits), base = 10L)

Arguments

x	an MPFR number (vector or array).
digits	how many significant digits (in the base chosen!) are to be used in the result. The default, NULL, uses enough digits to represent the full precision, often one or two digits more than “you” would expect. For bases 2,4,8,16, or 32, MPFR requires digits at least 2. For such bases, digits = 1 is changed into 2, with a message.
trim	logical; if FALSE, numbers are right-justified to a common width: if TRUE the leading blanks for justification are suppressed.
scientific	either a logical specifying whether MPFR numbers should be encoded in scientific format (“exponential representation”), or an integer penalty (see options("scipen")). Missing values correspond to the current default penalty.
maybe.full	logical, passed to .mpfr2str().
base	an integer in $2,3,..,62$; the base (“basis”) in which the numbers should be represented. Apart from the default base 10, binary (base = 2) or hexadecimal (base = 16) are particularly interesting.
showNeg0	logical indicating if “negative” zeros should be shown with a "-". The default, TRUE is intentially different from format(<numeric>).
exponent.char	the “exponent” character to be used in scientific notation. The default takes into account that for base $B \ge 15$, "e" is part of the (mantissa) digits and the same is true for "E" when $B \ge 37$.
exponent.plus	logical indicating if "+" should be for positive exponents in exponential (aka “scientific”) representation. This used to be hardcoded to FALSE; the new default is compatible to R's format()ing of numbers and helps to note visually when exponents are in use.
max.digits	a (large) positive number to limit the number of (mantissa) digits, notably when digits is NULL (as by default). Otherwise, a numeric digits is preferred to setting max.digits (which should not be smaller than digits).
big.mark, big.interval, small.mark, small.interval, decimal.mark, zero.print, drop0trailing	used for prettying decimal sequences, these are passed to prettyNum and that help page explains the details.
...	further arguments passed to or from other methods.

Value

a character vector or array, say cx, of the same length as x. Since Rmpfr version 0.5-3 (2013-09), if x is an mpfrArray, then cx is a character array with the same dim and dimnames as x.

Note that in scientific notation, the integer exponent is always in decimal, i.e., base 10 (even when base is not 10), but of course meaning base powers, e.g., in base 32, "u.giE3"is the same as "ugi0" which is $32^3$ times "u.gi". This is in contrast, e.g., with sprintf("%a", x) where the powers after "p" are powers of $2$.

Note

Currently, formatMpfr(x, scientific = FALSE) does not work correctly, e.g., for x <- Const("pi", 128) * 2^c(-200,200), i.e., it uses the scientific / exponential-style format. This is considered bogous and hopefully will change.

Author(s)

Martin Maechler

References

The MPFR manual's description of ‘⁠mpfr_get_str()⁠’ which is the C-internal workhorse for .mpfr2str() (on which formatMpfr() builds).

See Also

mpfr for creation and the mpfr class description with its many methods. The format generic, and the prettyNum utility on which formatMpfr is based as well. The S3 generic function formatN from package gmp.

.mpfr_formatinfo(x) provides the (cheap) non-string parts of .mpfr2str(x); the (base 2) exp exponents are also available via .mpfr2exp(x).

Examples

## Printing of MPFR numbers  uses formatMpfr() internally.
 ## Note how each components uses the "necessary" number of digits:
 ( x3 <- c(Const("pi", 168), mpfr(pi, 140), 3.14) )
 format(x3[3], 15)
 format(x3[3], 15, drop0 = TRUE)# "3.14" .. dropping the trailing zeros
 x3[4] <- 2^30
 x3[4] # automatically drops trailing zeros
 format(x3[1], dig = 41, small.mark = "'") # (41 - 1 = ) 40 digits after "."

 rbind(formatN(           x3,  digits = 15),
       formatN(as.numeric(x3), digits = 15))

 (Zero <- mpfr(c(0,1/-Inf), 20)) # 0 and "-0"
 xx <- c(Zero, 1:2, Const("pi", 120), -100*pi, -.00987)
 format(xx, digits = 2)
 format(xx, digits = 1, showNeg0 = FALSE)# "-0" no longer shown

## Output in other bases :
formatMpfr(mpfr(10^6, 40), base=32, drop0trailing=TRUE)
## "ugi0"
mpfr("ugi0", base=32) #-> 1'000'000


## This now works: The large number shows "as" large integer:
x <- Const("pi", 128) * 2^c(-200,200)
formatMpfr(x, scientific = FALSE) # was 1.955...e-60  5.048...e+60

i32 <- mpfr(1:32, precBits = 64)
format(i32,   base=  2, drop0trailing=TRUE)
format(i32,   base= 16, drop0trailing=TRUE)
format(1/i32, base=  2, drop0trailing=TRUE)# using scientific notation for [17..32]
format(1/i32, base= 32)
format(1/i32, base= 62, drop0trailing=TRUE)
format(mpfr(2, 64)^-(1:16), base=16, drop0trailing=TRUE)

---

Base-2 Representation and Multiplication of Mpfr Numbers

Description

MPFR - versions of the C99 (and POSIX) standard C (and C++) mathlib functions frexp() and ldexp().

frexpMpfr(x) computes base-2 exponent e and “mantissa”, or fraction r, such that $x = r * 2^e$, where $r \in [0.5, 1)$ (unless when x is in c(0, -Inf, Inf, NaN) where r == x and e is 0), and $e$ is integer valued.

ldexpMpfr(f, E) is the inverse of frexpMpfr(): Given fraction or mantissa f and integer exponent E, it returns $x = f * 2^E$. Viewed differently, it's the fastest way to multiply or divide MPFR numbers with $2^E$.

Usage

frexpMpfr(x,    rnd.mode = c("N", "D", "U", "Z", "A"))
ldexpMpfr(f, E, rnd.mode = c("N", "D", "U", "Z", "A"))

Arguments

x	numeric (coerced to double) vector.
f	numeric fraction (vector), in $[0.5, 1)$.
E	integer valued, exponent of 2, i.e., typically in (-1024-50):1024, otherwise the result will underflow to 0 or overflow to +/- Inf.
rnd.mode	a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see mpfr.

Value

frexpMpfr returns a list with named components r (of class mpfr) and e (integer valued, of type integer is small enough, "double" otherwise).

Author(s)

Martin Maechler

References

On unix-alikes, typically man frexp and man ldexp

See Also

Somewhat related, .mpfr2exp(). frexp() and ldexp() in package DPQ.

Examples

set.seed(47)
x <- c(0, 2^(-3:3), (-1:1)/0,
       sort(rlnorm(2^12, 10, 20) * sample(c(-1,1), 512, replace=TRUE)))
head(xM <- mpfr(x, 128), 11)
str(rFM <- frexpMpfr(xM))
d.fr <- with(rFM, data.frame(x=x, r=asNumeric(r), e=e))
head(d.fr , 16)
tail(d.fr)
ar <- abs(rFM$r)
stopifnot(0.5 <= ar[is.finite(x) & x != 0], ar[is.finite(x)] < 1,
          is.integer(rFM$e))
ldx <- with(rFM, ldexpMpfr(r, e))
(iN <- which(is.na(x))) # 10
stopifnot(exprs = {
  all.equal(xM, ldx, tol = 2^-124) # allow 4 bits loss, but apart from the NA, even:
  identical(xM[-iN], ldx[-iN])
  is.na(xM [iN])
  is.na(ldx[iN])
})

---

Conversion Utilities gmp <-> Rmpfr

Description

Coerce from and to big integers (bigz) and mpfr numbers.

Further, coerce from big rationals (bigq) to mpfr numbers.

Usage

.bigz2mpfr(x, precB = NULL, rnd.mode = c('N','D','U','Z','A'))
.bigq2mpfr(x, precB = NULL, rnd.mode = c('N','D','U','Z','A'))
.mpfr2bigz(x, mod = NA)
.mpfr2bigq(x)

Arguments

x	an R object of class bigz, bigq or mpfr respectively.
precB	precision in bits for the result. The default, NULL, means to use the minimal precision necessary for correct representation.
rnd.mode	a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see details of mpfr.
mod	a possible modulus, see as.bigz in package gmp.

Details

Note that we also provide the natural (S4) coercions, as(x, "mpfr") for x inheriting from class "bigz" or "bigq".

Value

a numeric vector of the same length as x, of the desired class.

See Also

mpfr(), as.bigz and as.bigq in package gmp.

Examples

S <- gmp::Stirling2(50,10)
 show(S)
 SS <- S * as.bigz(1:3)^128
 stopifnot(all.equal(log2(SS[2]) - log2(S), 128, tolerance=1e-15),
           identical(SS, .mpfr2bigz(.bigz2mpfr(SS))))

 .bigz2mpfr(S)            # 148 bit precision
 .bigz2mpfr(S, precB=256) # 256 bit

 ## rational --> mpfr:
 sq <- SS / as.bigz(2)^100
 MP <- as(sq, "mpfr")
 stopifnot(identical(MP, .bigq2mpfr(sq)),
           SS == MP * as(2, "mpfr")^100)

 ## New since 2024-01-20:   mpfr --> big rational "bigq"
 Pi <- Const("pi", 128)
 m <- Pi* 2^(-5:5)
 (m <- c(m, mpfr(2, 128)^(-5:5)))

 ## 1 x large num/denom, then 2^(-5:5) as frac
 tail( Q <- .mpfr2bigq(m) , 12)
 getDenom <- Rmpfr:::getDenom
 stopifnot(is.whole(m * (d.m <- getDenom(m))))
 stopifnot(all.equal(m, mpfr(Q, 130), tolerance = 2^-130)) # I see even
           all.equal(m, mpfr(Q, 130), tolerance = 0) # TRUE

 m <- m * mpfr(2, 128)^200 # quite a bit larger
 tail( Q. <- .mpfr2bigq(m) , 12) # large integers ..
 stopifnot(is.whole(m * (d.m <- getDenom(m))))
 stopifnot(all.equal(m, mpfr(Q., 130), tolerance = 2^-130)) # I see even
           all.equal(m, mpfr(Q., 130), tolerance = 0) # TRUE

 m2 <- m * mpfr(2, 128)^20000 ## really huge
 stopifnot(is.whole(m2 * (d.m2 <- getDenom(m2))))
 denominator(Q2 <- .mpfr2bigq(m2)) ## all 1 ! (all m2 ~~ 2^20200 )
 stopifnot(all.equal(m2, mpfr(Q2, 130), tolerance = 2^-130)) # I see even
           all.equal(m2, mpfr(Q2, 130), tolerance = 0) # TRUE

---

Hooke-Jeeves Derivative-Free Minimization R (working for MPFR)

Description

An implementation of the Hooke-Jeeves algorithm for derivative-free optimization.

This is a slight adaption hjk() from package dfoptim.

Usage

hjkMpfr(par, fn, control = list(), ...)

Arguments

par	Starting vector of parameter values. The initial vector may lie on the boundary. If lower[i]=upper[i] for some i, the i-th component of the solution vector will simply be kept fixed.
fn	Nonlinear objective function that is to be optimized. A scalar function that takes a real vector as argument and returns a scalar that is the value of the function at that point.
control	list of control parameters. See Details for more information.
...	Additional arguments passed to fn.

Details

Argument control is a list specifing changes to default values of algorithm control parameters. Note that parameter names may be abbreviated as long as they are unique.

The list items are as follows:

tol

Convergence tolerance. Iteration is terminated when the step length of the main loop becomes smaller than tol. This does not imply that the optimum is found with the same accuracy. Default is 1.e-06.

maxfeval

Maximum number of objective function evaluations allowed. Default is Inf, that is no restriction at all.

maximize

A logical indicating whether the objective function is to be maximized (TRUE) or minimized (FALSE). Default is FALSE.

target

A real number restricting the absolute function value. The procedure stops if this value is exceeded. Default is Inf, that is no restriction.

info

A logical variable indicating whether the step number, number of function calls, best function value, and the first component of the solution vector will be printed to the console. Default is FALSE.

If the minimization process threatens to go into an infinite loop, set either maxfeval or target.

Value

A list with the following components:

par	Best estimate of the parameter vector found by the algorithm.
value	value of the objective function at termination.
convergence	indicates convergence (TRUE) or not (FALSE).
feval	number of times the objective fn was evaluated.
niter	number of iterations (“steps”) in the main loop.

Note

This algorithm is based on the Matlab code of Prof. C. T. Kelley, given in his book “Iterative methods for optimization”. It has been implemented for package dfoptim with the permission of Prof. Kelley.

This version does not (yet) implement a cache for storing function values that have already been computed as searching the cache makes it slower.

Author(s)

Hans W Borchers [email protected]; for Rmpfr: John Nash, June 2012. Modifications by Martin Maechler.

References

C.T. Kelley (1999), Iterative Methods for Optimization, SIAM.

Quarteroni, Sacco, and Saleri (2007), Numerical Mathematics, Springer.

See Also

Standard R's optim; optimizeR provides one-dimensional minimization methods that work with mpfr-class numbers.

Examples

## simple smooth example:
ff <- function(x) sum((x - c(2:4))^2)
str(rr <- hjkMpfr(rep(mpfr(0,128), 3), ff, control=list(info=TRUE)))

doX <- Rmpfr:::doExtras(); cat("doExtras: ", doX, "\n") # slow parts only if(doX)

## Hooke-Jeeves solves high-dim. Rosenbrock function  {but slowly!}
rosenbrock <- function(x) {
    n <- length(x)
    sum (100*((x1 <- x[1:(n-1)])^2 - x[2:n])^2 + (x1 - 1)^2)
}
par0 <- rep(0, 10)
str(rb.db <- hjkMpfr(rep(0, 10), rosenbrock, control=list(info=TRUE)))
if(doX) {
## rosenbrook() is quite slow with mpfr-numbers:
str(rb.M. <- hjkMpfr(mpfr(numeric(10), prec=128), rosenbrock,
                     control = list(tol = 1e-8, info=TRUE)))
}


##  Hooke-Jeeves does not work well on non-smooth functions
nsf <- function(x) {
  f1 <- x[1]^2 + x[2]^2
  f2 <- x[1]^2 + x[2]^2 + 10 * (-4*x[1] - x[2] + 4)
  f3 <- x[1]^2 + x[2]^2 + 10 * (-x[1] - 2*x[2] + 6)
  max(f1, f2, f3)
}
par0 <- c(1, 1) # true min 7.2 at (1.2, 2.4)
h.d <- hjkMpfr(par0,            nsf) # fmin=8 at xmin=(2,2)
if(doX) {
## and this is not at all better (but slower!)
h.M <- hjkMpfr(mpfr(c(1,1), 128), nsf, control = list(tol = 1e-15))
}
## --> demo(hjkMpfr) # -> Fletcher's chebyquad function m = n -- residuals

---

Incomplete Gamma Function

Description

For MPFR version >= 3.2.0, the following MPFR library function is provided: mpfr_gamma_inc(a,x), the R interface of which is igamma(a,x), where igamma(a,x) is the “upper” incomplete gamma function

$\Gamma(a,x) :=: \Gamma(a) - \gamma(a,x),$

where

$\gamma(a,x) := \int_0^x t^{a-1} e^{-t} dt,$

and hence

$\Gamma(a,x) := \int_x^\infty t^{a-1} e^{-t} dt,$

and

$\Gamma(a) := \gamma(a, \infty).$

As R's pgamma(x,a) is

$\code{pgamma(x, a)} := \gamma(a,x) / \Gamma(a),$

we get

        igamma(a,x) ==  gamma(a) * pgamma(x, a, lower.tail=FALSE)

Usage

igamma(a, x, rnd.mode = c("N", "D", "U", "Z", "A"))

Arguments

a, x	an object of class mpfr or numeric.
rnd.mode	a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see mpfr.

Value

a numeric vector of “common length”, recyling along a and x.

Author(s)

R interface: Martin Maechler

References

NIST Digital Library of Mathematical Functions, section 8.2. https://dlmf.nist.gov/8.2.i

Wikipedia (2019). Incomplete gamma function; https://en.wikipedia.org/wiki/Incomplete_gamma_function

See Also

R's gamma (function) and pgamma (probability distribution).

Examples

## show how close pgamma() is :
x <- c(seq(0,20, by=1/4), 21:50, seq(55, 100, by=5))
if(mpfrVersion() >= "3.2.0") { print(
all.equal(igamma(Const("pi", 80), x),
          pgamma(x, pi, lower.tail=FALSE) * gamma(pi),
          tol=0, formatFUN = function(., ...) format(., digits = 7)) #-> 2.75e-16 (was 3.13e-16)
)
## and ensure *some* closeness:
stopifnot(exprs = {
   all.equal(igamma(Const("pi", 80), x),
             pgamma(x, pi, lower.tail=FALSE) * gamma(pi),
             tol = 1e-15)
})
} # only if MPFR version >= 3.2.0

---

One-Dimensional Numerical Integration - in pure R

Description

Numerical integration of one-dimensional functions in pure R, with care so it also works for "mpfr"-numbers.

Currently, only classical Romberg integration of order ord is available.

Usage

integrateR(f, lower, upper, ..., ord = NULL,
           rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol,
           max.ord = 19, verbose = FALSE)

Arguments

f	an R function taking a numeric or "mpfr" first argument and returning a numeric (or "mpfr") vector of the same length. Returning a non-finite element will generate an error.
lower, upper	the limits of integration. Currently must be finite. Do use "mpfr"-numbers to get higher than double precision, see the examples.
...	additional arguments to be passed to f.
ord	integer, the order of Romberg integration to be used. If this is NULL, as per default, and either rel.tol or abs.tol are specified, the order is increased until convergence.
rel.tol	relative accuracy requested. The default is 1.2e-4, about 4 digits only, see the Note.
abs.tol	absolute accuracy requested.
max.ord	only used, when neither ord or one of rel.tol, abs.tol are specified: Stop Romberg iterations after the order reaches max.ord; may prevent infinite loops or memory explosion.
verbose	logical or integer, indicating if and how much information should be printed during computation.

Details

Note that arguments after ... must be matched exactly.

For convergence, both relative and absolute changes must be smaller than rel.tol and abs.tol, respectively.

rel.tol cannot be less than max(50*.Machine$double.eps, 0.5e-28) if abs.tol <= 0.

Value

A list of class "integrateR" (as from standard R's integrate()) with a print method and components

value	the final estimate of the integral.
abs.error	estimate of the modulus of the absolute error.
subdivisions	for Romberg, the number of function evaluations.
message	"OK" or a character string giving the error message.
call	the matched call.

Note

f must accept a vector of inputs and produce a vector of function evaluations at those points. The Vectorize function may be helpful to convert f to this form.

If you want to use higher accuracy, you must set lower or upper to "mpfr" numbers (and typically lower the relative tolerance, rel.tol), see also the examples.

Note that the default tolerances (rel.tol, abs.tol) are not very accurate, but the same as for integrate, which however often returns considerably more accurate results than requested. This is typically not the case for integrateR().

Note

We use practically the same print S3 method as print.integrate, provided by R, with a difference when the message component is not "Ok".

Author(s)

Martin Maechler

References

Bauer, F.L. (1961) Algorithm 60 – Romberg Integration, Communications of the ACM 4(6), p.255.

See Also

R's standard, integrate, is much more adaptive, also allowing infinite integration boundaries, and typically considerably faster for a given accuracy.

Examples

## See more  from  ?integrate
## this is in the region where  integrate() can get problems:
integrateR(dnorm,0,2000)
integrateR(dnorm,0,2000, rel.tol=1e-15)
(Id <- integrateR(dnorm,0,2000, rel.tol=1e-15, verbose=TRUE))
Id$value == 0.5 # exactly

## Demonstrating that 'subdivisions' is correct:
Exp <- function(x) { .N <<- .N+ length(x); exp(x) }
.N <- 0; str(integrateR(Exp, 0,1, rel.tol=1e-10), digits=15); .N

### Using high-precision functions -----

## Polynomials are very nice:
integrateR(function(x) (x-2)^4 - 3*(x-3)^2, 0, 5, verbose=TRUE)
# n= 1, 2^n=        2 | I =            46.04, abs.err =      98.9583
# n= 2, 2^n=        4 | I =               20, abs.err =      26.0417
# n= 3, 2^n=        8 | I =               20, abs.err =  7.10543e-15
## 20 with absolute error < 7.1e-15
## Now, using higher accuracy:
I <- integrateR(function(x) (x-2)^4 - 3*(x-3)^2, 0, mpfr(5,128),
                rel.tol = 1e-20, verbose=TRUE)
I ; I$value  ## all fine

## with floats:
integrateR(exp,      0     , 1, rel.tol=1e-15, verbose=TRUE)
## with "mpfr":
(I <- integrateR(exp, mpfr(0,200), 1, rel.tol=1e-25, verbose=TRUE))
(I.true <- exp(mpfr(1, 200)) - 1)
## true absolute error:
stopifnot(print(as.numeric(I.true - I$value)) < 4e-25)

## Want absolute tolerance check only (=> set 'rel.tol' very high, e.g. 1):
(Ia <- integrateR(exp, mpfr(0,200), 1, abs.tol = 1e-6, rel.tol=1, verbose=TRUE))

## Set 'ord' (but no  '*.tol') --> Using 'ord'; no convergence checking
(I <- integrateR(exp, mpfr(0,200), 1,  ord = 13, verbose=TRUE))

---

Whole ("Integer") Numbers

Description

Check which elements of x[] are integer valued aka “whole” numbers,including MPFR numbers (class mpfr).

Usage

## S3 method for class 'mpfr'
is.whole(x)

Arguments

x	any R vector, here of class mpfr.

Value

logical vector of the same length as x, indicating where x[.] is integer valued.

Author(s)

Martin Maechler

See Also

is.integer(x) (base package) checks for the internal mode or class, not if x[i] are integer valued.

The is.whole() methods in package gmp.

Examples

is.integer(3) # FALSE, it's internally a double
 is.whole(3)   # TRUE
 x <- c(as(2,"mpfr") ^ 100, 3, 3.2, 1000000, 2^40)
 is.whole(x) # one FALSE, only

---

Compute f(a) = $\mathrm{log}$(1 +/- $\mathrm{exp}$(-a)) Numerically Optimally

Description

Compute f(a) = log(1 - exp(-a)), respectively g(x) = log(1 + exp(x)) quickly numerically accurately.

Usage

log1mexp(a, cutoff = log(2))
log1pexp(x, c0 = -37, c1 = 18, c2 = 33.3)

Arguments

a	numeric (or mpfr) vector of positive values.
x	numeric vector, may also be an "mpfr" object.
cutoff	positive number; log(2) is “optimal”, but the exact value is unimportant, and anything in $[0.5, 1]$ is fine.
c0, c1, c2	cutoffs for log1pexp; see below.

Value

$log1mexp(a) := f(a) = \log(1 - \exp(-a)) = \mathrm{log1p}(-\exp(-a)) = \log(-\mathrm{expm1}(-a))$

or, respectively,

$log1pexp(x) := g(x) = \log(1 + \exp(x)) = \mathrm{log1p}(\exp(x))$

computed accurately and quickly.

Author(s)

Martin Maechler, May 2002; log1pexp() in 2012

References

Martin Mächler (2012). Accurately Computing $\log(1-\exp(-|a|))$; https://CRAN.R-project.org/package=Rmpfr/vignettes/log1mexp-note.pdf.

Examples

fExpr <- expression(
          DEF   = log(1 - exp(-a)),
          expm1 = log(-expm1(-a)),
          log1p = log1p(-exp(-a)),
          F     = log1mexp(a))
a. <- 2^seq(-58, 10, length = 256)
a <- a. ; str(fa <- do.call(cbind, as.list(fExpr)))
head(fa)# expm1() works here
tail(fa)# log1p() works here

## graphically:
lwd <- 1.5*(5:2); col <- adjustcolor(1:4, 0.4)
op <- par(mfcol=c(1,2), mgp = c(1.25, .6, 0), mar = .1+c(3,2,1,1))
  matplot(a, fa, type = "l", log = "x", col=col, lwd=lwd)
  legend("topleft", fExpr, col=col, lwd=lwd, lty=1:4, bty="n")
  # expm1() & log1mexp() work here

  matplot(a, -fa, type = "l", log = "xy", col=col, lwd=lwd)
  legend("left", paste("-",fExpr), col=col, lwd=lwd, lty=1:4, bty="n")
  # log1p() & log1mexp() work here
par(op)

aM <- 2^seqMpfr(-58, 10, length=length(a.)) # => default prec = 128
a <- aM; dim(faM <- do.call(cbind, as.list(fExpr))) # 256 x 4, "same" as 'fa'
## Here, for small 'a' log1p() and even 'DEF' is still good enough
l_f <- asNumeric(log(-faM))
all.equal(l_f[,"F"], l_f[,"log1p"], tol=0) # see TRUE (Lnx 64-bit)
io <- a. < 80 # for these, the differences are small
all.equal(l_f[io,"F"], l_f[io,"expm1"], tol=0) # see 6.662e-9
all.equal(l_f[io,"F"], l_f[io, "DEF" ], tol=0)
stopifnot(exprs = {
  all.equal(l_f[,"F"], l_f[,"log1p"],     tol= 1e-15)
  all.equal(l_f[io,"F"], l_f[io,"expm1"], tol= 1e-7)
  all.equal(l_f[io,"F"], l_f[io, "DEF" ], tol= 1e-7)
})
## For 128-bit prec, if we go down to 2^-130, "log1p" is no longer ok:
aM2 <- 2^seqMpfr(-130, 10, by = 1/2)
a <- aM2; fa2 <- do.call(cbind, as.list(fExpr))
head(asNumeric(fa2), 12)
tail(asNumeric(fa2), 12)

matplot(a, log(-fa2[,1:3]) -log(-fa2[,"F"]),  type="l", log="x",
        lty=1:3, lwd=2*(3:1)-1, col=adjustcolor(2:4, 1/3))
legend("top", colnames(fa2)[1:3], lty=1:3, lwd=2*(3:1)-1, col=adjustcolor(2:4, 1/3))

cols <- adjustcolor(2:4, 1/3); lwd <- 2*(3:1)-1
matplot(a, 1e-40+abs(log(-fa2[,1:3]) -log(-fa2[,"F"])),  type="o", log="xy",
        main = "log1mexp(a) -- approximation rel.errors, mpfr(*, prec=128)",
        pch=21:23, cex=.6, bg=5:7, lty=1:2, lwd=lwd, col=cols)
legend("top", colnames(fa2)[1:3], bty="n", lty=1:2, lwd=lwd, col=cols,
        pch=21:23, pt.cex=.6, pt.bg=5:7)


## -------------------------- log1pexp() [simpler] --------------------

curve(log1pexp, -10, 10, asp=1)
abline(0,1, h=0,v=0, lty=3, col="gray")

## Cutoff c1 for log1pexp() -- not often "needed":
curve(log1p(exp(x)) - log1pexp(x), 16, 20, n=2049)
## need for *some* cutoff:
x <- seq(700, 720, by=2)
cbind(x, log1p(exp(x)), log1pexp(x))

## Cutoff c2 for log1pexp():
curve((x+exp(-x)) - x, 20, 40, n=1025)
curve((x+exp(-x)) - x, 33.1, 33.5, n=1025)

---

(MPFR) Matrix (Vector) Multiplication

Description

Matrix / vector multiplication of mpfr (and “simple” numeric) matrices and vectors.

matmult (x,y, fPrec = 2) or crossprod(x,y, fPrec = 2) use higher precision in underlying computations.

Usage

matmult(x, y, ...)

Arguments

x, y

numeric or mpfrMatrix-classed R objects, i.e. semantically numeric matrices or vectors.

...

arguments passed to the hidden underlying .matmult.R() work horse which is also underlying the %*%, crossprod(), and tcrossprod() methods, see the mpfrMatrix class documentation: fPrec a multiplication factor, a positive number determining the number of bits precBits used for the underlying multiplication and summation arithmetic. The default is fPrec = 1. Setting fPrec = 2 doubles the precision which has been recommended, e.g., by John Nash. precBits the number of bits used for the underlying multiplication and summation arithmetic; by default precBits = fPrec * max(getPrec(x), getPrec(y)) which typically uses the same accuracy as regular mpfr-arithmetic would use.

Value

a (base R) matrix or mpfrMatrix, depending on the classes of x and y.

Note

Using matmult(x,y) instead of x %*% y, makes sense mainly if you use non-default fPrec or precBits arguments.

The crossprod(), and tcrossprod() function have the identical optional arguments fPrec or precBits.

Author(s)

Martin Maechler

See Also

%*%, crossprod, tcrossprod.

Examples

## FIXME: add example

## 1)  matmult()  <-->  %*%

## 2)  crossprod() , tcrossprod()  %% <--> ./mpfrMatrix-class.Rd  examples  (!)

---

Class "Mnumber" and "mNumber" of "mpfr" and regular numbers and arrays from them

Description

Classes "Mnumber" "mNumber" are class unions of "mpfr" and regular numbers and arrays from them.
Its purpose is for method dispatch, notably defining a cbind(...) method where ... contains objects of one of the member classes of "Mnumber".

Classes "mNumber" is considerably smaller is it does not contain "matrix" and "array" since these also extend "character" which is not really desirable for generalized numbers. It extends the simple "numericVector" class by mpfr* classes.

Methods

%*%

signature(x = "mpfrMatrix", y = "Mnumber"): ...

crossprod

signature(x = "mpfr", y = "Mnumber"): ...

tcrossprod

signature(x = "Mnumber", y = "mpfr"): ...

etc. These are documented with the classes mpfr and or mpfrMatrix.

See Also

the array_or_vector sub class; cbind-methods.

Examples

## "Mnumber" encompasses (i.e., "extends") quite a few
##  "vector / array - like" classes:
showClass("Mnumber")
stopifnot(extends("mpfrMatrix", "Mnumber"),
          extends("array",      "Mnumber"))

Mnsub <- names(getClass("Mnumber")@subclasses)
(mNsub <- names(getClass("mNumber")@subclasses))
## mNumber has *one* subclass which is not in Mnumber:
setdiff(mNsub, Mnsub)# namely "numericVector"
## The following are only subclasses of "Mnumber", but not of "mNumber":
setdiff(Mnsub, mNsub)

---

Create "mpfr" Numbers (Objects)

Description

Create multiple (i.e. typically high) precision numbers, to be used in arithmetic and mathematical computations with R.

Usage

mpfr(x, precBits, ...)
## Default S3 method:
mpfr(x, precBits, base = 10,
     rnd.mode = c("N","D","U","Z","A"), scientific = NA, ...)

Const(name = c("pi", "gamma", "catalan", "log2"), prec = 120L,
      rnd.mode = c("N","D","U","Z","A"))

is.mpfr(x)

Arguments

x	a numeric, mpfr, bigz, bigq, or character vector or array.
precBits, prec	a number, the maximal precision to be used, in bits; i.e. 53 corresponds to double precision. Must be at least 2. If missing, getPrec(x) determines a default precision.
base	(only when x is character) the base with respect to which x[i] represent numbers; base $b$ must fulfill $2 \le b \le 62$.
rnd.mode	a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see details.
scientific	(used only when x is the result of formatBin(), i.e., of class "Bcharacter":) logical indicating that the binary representation of x is in scientific notation. When TRUE, mpfr() will substitute 0 for _; when NA, mpfr() will guess, and use TRUE when finding a "p" in x; see also formatBin.
name	a string specifying the mpfrlib - internal constant computation. "gamma" is Euler's gamma ($\gamma$), and "catalan" Catalan's constant.
...	potentially further arguments passed to and from methods.

Details

The "mpfr" method of mpfr() is a simple wrapper around roundMpfr().

MPFR supports the following rounding modes,

GMP_RNDN:

round to nearest (roundTiesToEven in IEEE 754-2008).

GMP_RNDZ:

round toward zero (roundTowardZero in IEEE 754-2008).

GMP_RNDU:

round toward plus infinity (“Up”, roundTowardPositive in IEEE 754-2008).

GMP_RNDD:

round toward minus infinity (“Down”, roundTowardNegative in IEEE 754-2008).

GMP_RNDA:

round away from zero (new since MPFR 3.0.0).

The ‘round to nearest’ ("N") mode, the default here, works as in the IEEE 754 standard: in case the number to be rounded lies exactly in the middle of two representable numbers, it is rounded to the one with the least significant bit set to zero. For example, the number 5/2, which is represented by (10.1) in binary, is rounded to (10.0)=2 with a precision of two bits, and not to (11.0)=3. This rule avoids the "drift" phenomenon mentioned by Knuth in volume 2 of The Art of Computer Programming (Section 4.2.2).

When x is character, mpfr() will detect the precision of the input object.

Value

an object of (S4) class mpfr, or for mpfr(x) when x is an array, mpfrMatrix, or mpfrArray which the user should just as a normal numeric vector or array.

is.mpfr() returns TRUE or FALSE.

Author(s)

Martin Maechler

References

The MPFR team. (202x). GNU MPFR – The Multiple Precision Floating-Point Reliable Library; see https://www.mpfr.org/mpfr-current/#doc or directly https://www.mpfr.org/mpfr-current/mpfr.pdf.

See Also

The class documentation mpfr contains more details. Use asNumeric() from gmp to transform back to double precision ("numeric").

Examples

mpfr(pi, 120) ## the double-precision pi "translated" to 120-bit precision

pi. <- Const("pi", prec = 260) # pi "computed" to correct 260-bit precision
pi. # nicely prints 80 digits [260 * log10(2) ~= 78.3 ~ 80]

Const("gamma",   128L) # 0.5772...
Const("catalan", 128L) # 0.9159...

x <- mpfr(0:7, 100)/7 # a more precise version of  k/7, k=0,..,7
x
1 / x

## character input :
mpfr("2.718281828459045235360287471352662497757") - exp(mpfr(1, 150))
## ~= -4 * 10^-40
## Also works for  NA, NaN, ... :
cx <- c("1234567890123456", 345, "NA", "NaN", "Inf", "-Inf")
mpfr(cx)

## with some 'base' choices :
print(mpfr("111.1111", base=2)) * 2^4

mpfr("af21.01020300a0b0c", base=16)
##  68 bit prec.  44833.00393694653820642

mpfr("ugi0", base = 32) == 10^6   ## TRUE

## --- Large integers from package 'gmp':
Z <- as.bigz(7)^(1:200)
head(Z, 40)
## mfpr(Z) by default chooses the correct *maximal* default precision:
mZ. <- mpfr(Z)
## more efficiently chooses precision individually
m.Z <- mpfr(Z, precBits = frexpZ(Z)$exp)
## the precBits chosen are large enough to keep full precision:
stopifnot(identical(cZ <- as.character(Z),
                    as(mZ.,"character")),
          identical(cZ, as(m.Z,"character")))

## compare mpfr-arithmetic with exact rational one:
stopifnot(all.equal(mpfr(as.bigq(355,113), 99),
                    mpfr(355, 99) / 113,	tol = 2^-98))

## look at different "rounding modes":
sapply(c("N", "D","U","Z","A"), function(RND)
       mpfr(c(-1,1)/5, 20, rnd.mode = RND), simplify=FALSE)

symnum(sapply(c("N", "D","U","Z","A"),
              function(RND) mpfr(0.2, prec = 5:15, rnd.mode = RND) < 0.2 ))

---

Class "mpfr" of Multiple Precision Floating Point Numbers

Description

"mpfr" is the class of Multiple Precision Floatingpoint numbers with Reliable arithmetic.

sFor the high-level user, "mpfr" objects should behave as standard R's numeric vectors. They would just print differently and use the prespecified (typically high) precision instead of the double precision of ‘traditional’ R numbers (with class(.) == "numeric" and typeof(.) == "double").

hypot(x,y) computes the hypothenuse length $z$ in a rectangular triangle with “leg” side lengths $x$ and $y$, i.e.,

$z = hypot(x,y) = \sqrt{x^2 + y^2},$

in a numerically stable way.

Usage

hypot(x,y, rnd.mode = c("N","D","U","Z","A"))

Arguments

x, y	an object of class mpfr.
rnd.mode	a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see mpfr.

Objects from the Class

Objects are typically created by mpfr(<number>, precBits).

summary(<mpfr>) returns an object of class "summaryMpfr" which contains "mpfr" but has its own print method.

Slots

Internally, "mpfr" objects just contain standard R lists where each list element is of class "mpfr1", representing one MPFR number, in a structure with four slots, very much parallelizing the C struc in the mpfr C library to which the Rmpfr package interfaces.

An object of class "mpfr1" has slots

prec:

"integer" specifying the maxmimal precision in bits.

exp:

"integer" specifying the base-2 exponent of the number.

sign:

"integer", typically -1 or 1, specifying the sign (i.e. sign(.)) of the number.

d:

an "integer" vector (of 32-bit “limbs”) which corresponds to the full mantissa of the number.

Methods

abs

signature(x = "mpfr"): ...

atan2

signature(y = "mpfr", x = "ANY"), and

atan2

signature(x = "ANY", y = "mpfr"): compute the arc-tangent of two arguments: atan2(y, x) returns the angle between the x-axis and the vector from the origin to $(x, y)$, i.e., for positive arguments atan2(y, x) == atan(y/x).

lbeta

signature(a = "ANY", b = "mpfrArray"), is $\log(|B(a,b)|)$ where $B(a,b)$ is the Beta function, beta(a,b).

beta

signature(a = "mpfr", b = "ANY"),

beta

signature(a = "mpfr", b = "mpfr"), ..., etc: Compute the beta function $B(a,b)$, using high precision, building on internal gamma or lgamma. See the help for R's base function beta for more. Currently, there, $a,b \ge 0$ is required. Here, we provide (non-NaN) for all numeric a, b.

When either $a$, $b$, or $a+b$ is a negative integer, $\Gamma(.)$ has a pole there and is undefined (NaN). However the Beta function can be defined there as “limit”, in some cases. Following other software such as SAGE, Maple or Mathematica, we provide finite values in these cases. However, note that these are not proper limits (two-dimensional in $(a,b)$), but useful for some applications. E.g., $B(a,b)$ is defined as zero when $a+b$ is a negative integer, but neither $a$ nor $b$ is. Further, if $a > b > 0$ are integers, $B(-a,b)= B(b,-a)$ can be seen as $(-1)^b * B(a-b+1,b)$.

dim<-

signature(x = "mpfr"): Setting a dimension dim on an "mpfr" object makes it into an object of class "mpfrArray" or (more specifically) "mpfrMatrix" for a length-2 dimension, see their help page; note that t(x) (below) is a special case of this.

Ops

signature(e1 = "mpfr", e2 = "ANY"): ...

Ops

signature(e1 = "ANY", e2 = "mpfr"): ...

Arith

signature(e1 = "mpfr", e2 = "missing"): ...

Arith

signature(e1 = "mpfr", e2 = "mpfr"): ...

Arith

signature(e1 = "mpfr", e2 = "integer"): ...

Arith

signature(e1 = "mpfr", e2 = "numeric"): ...

Arith

signature(e1 = "integer", e2 = "mpfr"): ...

Arith

signature(e1 = "numeric", e2 = "mpfr"): ...

Compare

signature(e1 = "mpfr", e2 = "mpfr"): ...

Compare

signature(e1 = "mpfr", e2 = "integer"): ...

Compare

signature(e1 = "mpfr", e2 = "numeric"): ...

Compare

signature(e1 = "integer", e2 = "mpfr"): ...

Compare

signature(e1 = "numeric", e2 = "mpfr"): ...

Logic

signature(e1 = "mpfr", e2 = "mpfr"): ...

Summary

signature(x = "mpfr"): The S4 Summary group functions, max, min, range, prod, sum, any, and all are all defined for MPFR numbers. mean(x, trim) for non-0 trim works analogously to mean.default.

median

signature(x = "mpfr"): works via

quantile

signature(x = "mpfr"): a simple wrapper of the quantile.default method from stats.

summary

signature(object = "mpfr"): modeled after summary.default, ensuring to provide the full "mpfr" range of numbers.

Math

signature(x = "mpfr"): All the S4 Math group functions are defined, using multiple precision (MPFR) arithmetic, from getGroupMembers("Math"), these are (in alphabetical order):

abs, sign, sqrt, ceiling, floor, trunc, cummax, cummin, cumprod, cumsum, exp, expm1, log, log10, log2, log1p, cos, cosh, sin, sinh, tan, tanh, acos, acosh, asin, asinh, atan, atanh, cospi, sinpi, tanpi, gamma, lgamma, digamma, and trigamma.

Currently, trigamma is not provided by the MPFR library and hence not yet implemented.
Further, the cum*() methods are not yet implemented.

factorial

signature(x = "mpfr"): this will round the result when x is integer valued. Note however that factorialMpfr(n) for integer n is slightly more efficient, using the MPFR function ‘⁠mpfr_fac_ui⁠’.

Math2

signature(x = "mpfr"): round(x, digits) and signif(x, digits) methods. Note that these do not change the formal precision ('prec' slot), and you may often want to apply roundMpfr() in addition or preference.

as.numeric

signature(x = "mpfr"): ...

as.vector

signature(x = "mpfrArray"): as for standard arrays, this “drops” the dim (and dimnames), i.e., transforms x into an ‘MPFR’ number vector, i.e., class mpfr.

[[

signature(x = "mpfr", i = "ANY"), and

[

signature(x = "mpfr", i = "ANY", j = "missing", drop = "missing"): subsetting aka “indexing” happens as for numeric vectors.

format

signature(x = "mpfr"), further arguments digits = NULL, scientific = NA, etc: returns character vector of same length as x; when digits is NULL, with enough digits to recreate x accurately. For details, see formatMpfr.

is.finite

signature(x = "mpfr"): ...

is.infinite

signature(x = "mpfr"): ...

is.na

signature(x = "mpfr"): ...

is.nan

signature(x = "mpfr"): ...

log

signature(x = "mpfr"): ...

show

signature(object = "mpfr"): ...

sign

signature(x = "mpfr"): ...

Re, Im

signature(z = "mpfr"): simply return z or 0 (as "mpfr" numbers of correct precision), as mpfr numbers are ‘real’ numbers.

Arg, Mod, Conj

signature(z = "mpfr"): these are trivial for our ‘real’ mpfr numbers, but defined to work correctly when used in R code that also allows complex number input.

all.equal

signature(target = "mpfr", current = "mpfr"),

all.equal

signature(target = "mpfr", current = "ANY"), and

all.equal

signature(target = "ANY", current = "mpfr"): methods for numerical (approximate) equality, all.equal of multiple precision numbers. Note that the default tolerance (argument) is taken to correspond to the (smaller of the two) precisions when both main arguments are of class "mpfr", and hence can be considerably less than double precision machine epsilon .Machine$double.eps.

coerce

signature(from = "numeric", to = "mpfr"): as(., "mpfr") coercion methods are available for character strings, numeric, integer, logical, and even raw. Note however, that mpfr(., precBits, base) is more flexible.

coerce

signature(from = "mpfr", to = "bigz"): coerces to biginteger, see bigz in package gmp.

coerce

signature(from = "mpfr", to = "numeric"): ...

coerce

signature(from = "mpfr", to = "character"): ...

unique

signature(x = "mpfr"), and corresponding S3 method (such that unique(<mpfr>) works inside base functions), see unique.

Note that duplicated() works for "mpfr" objects without the need for a specific method.

t

signature(x = "mpfr"): makes x into an $n \times 1$ mpfrMatrix.

which.min

signature(x = "mpfr"): gives the index of the first minimum, see which.min.

which.max

signature(x = "mpfr"): gives the index of the first maximum, see which.max.

Note

Many more methods (“functions”) automagically work for "mpfr" number vectors (and matrices, see the mpfrMatrix class doc), notably sort, order, quantile, rank.

Author(s)

Martin Maechler

See Also

The "mpfrMatrix" class, which extends the "mpfr" one.

roundMpfr to change precision of an "mpfr" object which is typically desirable instead of or in addition to signif() or round(); is.whole() from gmp, etc.

Special mathematical functions such as some Bessel ones, e.g., jn; further, zeta(.) $(= \zeta(.))$, Ei() etc. Bernoulli numbers and the Pochhammer function pochMpfr.

Examples

## 30 digit precision
(x <- mpfr(c(2:3, pi), prec = 30 * log2(10)))
str(x) # str() displays *compact*ly => not full precision
x^2
x[1] / x[2] # 0.66666... ~ 30 digits

## indexing - as with numeric vectors
stopifnot(exprs = {
   identical(x[2], x[[2]])
   ## indexing "outside" gives NA (well: "mpfr-NaN" for now):
   is.na(x[5])
   ## whereas "[[" cannot index outside:
   inherits(tryCatch(x[[5]], error=identity), "error")
   ## and only select *one* element:
   inherits(tryCatch(x[[2:3]], error=identity), "error")
})

## factorial() & lfactorial would work automagically via [l]gamma(),
## but factorial() additionally has an "mpfr" method which rounds
f200 <- factorial(mpfr(200, prec = 1500)) # need high prec.!
f200
as.numeric(log2(f200))# 1245.38 -- need precBits >~ 1246 for full precision

##--> see  factorialMpfr() for more such computations.

##--- "Underflow" **much** later -- exponents have 30(+1) bits themselves:

mpfr.min.exp2 <- - (2^30 + 1)
two <- mpfr(2, 55)
stopifnot(two ^ mpfr.min.exp2 == 0)
## whereas
two ^ (mpfr.min.exp2 * (1 - 1e-15))
## 2.38256490488795107e-323228497   ["typically"]

##--- "Assert" that {sort}, {order}, {quantile}, {rank}, all work :

p <- mpfr(rpois(32, lambda=500), precBits=128)^10
np <- as.numeric(log(p))
(sp <- summary(p))# using the print.summaryMpfr() method
stopifnot(all(diff(sort(p)) >= 0),
   identical(order(p), order(np)),
   identical(rank (p), rank (np)),
   all.equal(sapply(1:9, function(Typ) quantile(np, type=Typ, names=FALSE)),
      sapply(lapply(1:9, function(Typ) quantile( p, type=Typ, names=FALSE)),
	     function(x) as.numeric(log(x))),
      tol = 1e-3),# quantiles: interpolated in orig. <--> log scale
 TRUE)

m0 <- mpfr(numeric(), 99)
xy <- expand.grid(x = -2:2, y = -2:2) ; x <- xy[,"x"] ; y <- xy[,"y"]
a2. <- atan2(y,x)

stopifnot(identical(which.min(m0), integer(0)),
	  identical(which.max(m0), integer(0)),
          all.equal(a2., atan2(as(y,"mpfr"), x)),
	  max(m0) == mpfr(-Inf, 53), # (53 is not a feature, but ok)
	  min(m0) == mpfr(+Inf, 53),
	  sum(m0) == 0, prod(m0) == 1)

## unique(), now even base::factor()  "works" on <mpfr> :
set.seed(17)
p <- rlnorm(20) * mpfr(10, 100)^-999
pp <- sample(p, 50, replace=TRUE)
str(unique(pp)) # length 18 .. (from originally 20)
## Class 'mpfr' [package "Rmpfr"] of length 18 and precision 100
## 5.56520587824e-999 4.41636588227e-1000 ..
facp <- factor(pp)
str(facp) # the factor *levels* are a bit verbose :
# Factor w/ 18 levels "new(\"mpfr1\", ...........)" ...
# At least *some* factor methods work :
stopifnot(exprs = {
  is.factor(facp)
  identical(unname(table(facp)),
            unname(table(asNumeric(pp * mpfr(10,100)^1000))))
})

## ((unfortunately, the expressions are wrong; should integer "L"))
#
## More useful: levels with which to *invert* factor() :
## -- this is not quite ok:
## simplified from 'utils' :
deparse1 <- function(x, ...) paste(deparse(x, 500L, ...), collapse = " ")
if(FALSE) {
 str(pp.levs <- vapply(unclass(sort(unique(pp))), deparse1, ""))
 facp2 <- factor(pp, levels = pp.levs)
}

---