mpfr-class: Class "mpfr" of Multiple Precision Floating Point Numbers

Description

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

For 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: This method is identical to the (exported) formatMpfr() function, see its help page for details.

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"): The show method for "mpfr" classed objects calls the S3 method print.mpfr(object) with several optional arguments, itself based on the format() method which calls formatMpfr().

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.

Author

Martin Maechler

Examples

Run this code

## 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  :
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)
}

Run the code above in your browser using DataLab

Description

Usage

Arguments

Objects from the Class

Slots

Methods

Author

See Also

Examples