Class "mpfr" of Multiple Precision Floating Point Numbers
"mpfr" is the class of Multiple Precision
Floatingpoint numbers with Reliable arithmetic.
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.
hypot(x,y, rnd.mode = c("N","D","U","Z","A"))x,y |
an object of class |
rnd.mode |
a 1-letter string specifying how rounding
should happen at C-level conversion to MPFR, see |
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.
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.
signature(x = "mpfr"): ...
signature(y = "mpfr", x = "ANY"), and
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).
signature(a = "ANY", b = "mpfrArray"), is
log(abs(B(a,b))) where B(a,b) is the
Beta function, beta(a,b).
signature(a = "mpfr", b = "ANY"),
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 >= 0 is required.
Here, we provide (non-NaN) for all numeric a, b.
When either a, b, or a+b is a negative
integer, Γ(.) 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).
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.
signature(e1 = "mpfr", e2 = "ANY"): ...
signature(e1 = "ANY", e2 = "mpfr"): ...
signature(e1 = "mpfr", e2 = "missing"): ...
signature(e1 = "mpfr", e2 = "mpfr"): ...
signature(e1 = "mpfr", e2 = "integer"): ...
signature(e1 = "mpfr", e2 = "numeric"): ...
signature(e1 = "integer", e2 = "mpfr"): ...
signature(e1 = "numeric", e2 = "mpfr"): ...
signature(e1 = "mpfr", e2 = "mpfr"): ...
signature(e1 = "mpfr", e2 = "integer"): ...
signature(e1 = "mpfr", e2 = "numeric"): ...
signature(e1 = "integer", e2 = "mpfr"): ...
signature(e1 = "numeric", e2 = "mpfr"): ...
signature(e1 = "mpfr", e2 = "mpfr"): ...
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.
signature(x = "mpfr"): works via
signature(x = "mpfr"): a simple wrapper of
the quantile.default method from stats.
signature(object = "mpfr"): modeled after
summary.default, ensuring to provide the full "mpfr"
range of numbers.
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, 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.
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.
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.
signature(x = "mpfr"): ...
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.
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.
signature(x = "mpfr"): ...
signature(x = "mpfr"): ...
signature(x = "mpfr"): ...
signature(x = "mpfr"): ...
signature(x = "mpfr"): ...
signature(object = "mpfr"): ...
signature(x = "mpfr"): ...
signature(z = "mpfr"): simply return z
or 0 (as "mpfr" numbers of correct precision), as mpfr
numbers are ‘real’ numbers.
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.
signature(target = "mpfr", current = "mpfr"),
signature(target = "mpfr", current = "ANY"), and
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.
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.
signature(from = "mpfr", to = "bigz"): coerces
to biginteger, see bigz in package gmp.
signature(from = "mpfr", to = "numeric"): ...
signature(from = "mpfr", to = "character"): ...
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.
signature(x = "mpfr"): makes x into an n x 1 mpfrMatrix.
signature(x = "mpfr"): gives the index of
the first minimum, see which.min.
signature(x = "mpfr"): gives the index of
the first maximum, see which.max.
Many more methods (“functions”) automagically work for
"mpfr" number vectors (and matrices, see the
mpfrMatrix class doc),
notably
sort, order, quantile,
rank.
Martin Maechler
The "mpfrMatrix" class, which extends the
"mpfr" one.
## 30 digit precision
str(x <- mpfr(c(2:3, pi), prec = 30 * log2(10)))
x^2
x[1] / x[2] # 0.66666... ~ 30 digits
## indexing - as with numeric vectors
stopifnot(identical(x[2], x[[2]]),
## indexing "outside" gives NA (well: "mpfr-NaN" for now):
is.na(x[5]),
## whereas "[[" cannot index outside:
is(try(x[[5]]), "try-error"),
## and only select *one* element:
is(try(x[[2:3]]), "try-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)
}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.