Version:	0.7-5
Date:	2024-08-23
Title:	Multiple Precision Arithmetic
Maintainer:	Antoine Lucas <antoinelucas@gmail.com>
Description:	Multiple Precision Arithmetic (big integers and rationals, prime number tests, matrix computation), "arithmetic without limitations" using the C library GMP (GNU Multiple Precision Arithmetic).
Depends:	R (≥ 3.5.0)
Imports:	methods
Suggests:	Rmpfr, MASS, round
SystemRequirements:	gmp (>= 4.2.3)
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
BuildResaveData:	no
LazyDataNote:	not available, as we use data/*.R *and* our classes
NeedsCompilation:	yes
URL:	https://forgemia.inra.fr/sylvain.jasson/gmp
Packaged:	2024-08-23 17:22:31 UTC; antoine
Author:	Antoine Lucas [aut, cre], Immanuel Scholz [aut], Rainer Boehme [ctb], Sylvain Jasson [ctb], Martin Maechler [ctb]
Repository:	CRAN
Date/Publication:	2024-08-23 18:40:02 UTC

Apply Functions Over Matrix Margins (Rows or Columns)

Description

These are S3 methods for apply() which we re-export as S3 generic function. They “overload” the apply() function for big rationals ("bigq") and big integers ("bigz").

Usage

## S3 method for class 'bigz'
apply(X, MARGIN, FUN, ...)
## S3 method for class 'bigq'
apply(X, MARGIN, FUN, ...)

Arguments

Value

The bigz and bigq methods return a vector of class "bigz" or "bigq", respectively.

Author(s)

Antoine Lucas

See Also

apply; lapply is used by our apply() method.

Examples

 x <- as.bigz(matrix(1:12,3))
 apply(x,1,min)
 apply(x,2,max)

 x <- as.bigq(x ^ 3, d = (x + 3)^2)
 apply(x,1, min)
 apply(x,2, sum)
 ## now use the "..." to pass  na.rm=TRUE :
 x[2,3] <- NA
 apply(x,1, sum)
 apply(x,1, sum, na.rm = TRUE)

Coerce to 'numeric', not Loosing Dimensions

Description

a number-like object is coerced to type (typeof) "numeric", keeping dim (and maybe dimnames) when present.

Usage

asNumeric(x)

Arguments

Value

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

Methods

signature(x = "ANY")

the default method, which is the identity for numeric array.

signature(x = "bigq")

the method for big rationals.

signature(x = "bigq")

the method for big integers.

Note that package Rmpfr provides methods for its own number-like objects.

Author(s)

Martin Maechler

See Also

as.numeric coerces to both "numeric" and to a vector, whereas asNumeric() should keep dim (and other) attributes.

Examples

m <- matrix(1:6, 2,3)
stopifnot(identical(m, asNumeric(m)))# remains matrix

(M <- as.bigz(m) / 5) ##-> "bigq" matrix
asNumeric(M) # numeric matrix
stopifnot(all.equal(asNumeric(M), m/5))

Exact Bernoulli Numbers

Description

Return the n-th Bernoulli number B_n, (or B_n^+, see the reference), where B_1 = + \frac 1 2.

Usage

BernoulliQ(n, verbose = getOption("verbose", FALSE))

Arguments

Value

a big rational (class "bigq") vector of the Bernoulli numbers B_n.

Author(s)

Martin Maechler

References

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

See Also

Bernoulli in Rmpfr in arbitrary precision via Riemann's \zeta function. Bern(n) in DPQ uses standard (double precision) R arithmetic for the n-th Bernoulli number.

Examples

(Bn0.10 <- BernoulliQ(0:10))

Large sized rationals

Description

Class "bigq" encodes rationals encoded as ratios of arbitrary large integers (via GMP). A simple S3 class (internally a raw vector), it has been registered as formal (S4) class (via setOldClass), too.

Usage

as.bigq(n, d = 1)
## S3 method for class 'bigq'
as.character(x, b=10,...)
## S3 method for class 'bigq'
as.double(x,...)
as.bigz.bigq(a, mod=NA)
is.bigq(x)
## S3 method for class 'bigq'
is.na(x)
## S3 method for class 'bigq'
print(x, quote=FALSE, initLine = TRUE, ...)
denominator(x)
numerator(x)
NA_bigq_
c_bigq(L)

Arguments

Details

as.bigq(x) when x is numeric (aka double precision) calls the ‘GMP’ function mpq_set_d() which is documented to be exact (every finite double precision number is a rational number).

as.bigz.bigq() returns the smallest integers not less than the corresponding rationals bigq.

NA_bigq_ is computed on package load time as as.bigq(NA).

Value

An R object of (S3) class "bigq" representing the parameter value.

Author(s)

Antoine Lucas

Examples

x <- as.bigq(21,6)
x
# 7 / 2
# Wow ! result is simplified.

y <- as.bigq(5,3)

# addition works !
x + y

# You can even try multiplication, division...
x * y / 13

# and, since May 2012,
x ^ 20
stopifnot(is.bigq(x), is.bigq(x + y),
	  x ^ 20 == as.bigz(7)^20 / 2^20)

# convert to string, double
as.character(x)
as.double(x)

stopifnot( is.na(NA_bigq_) )

# Depict the "S4-class" bigq, i.e., the formal (S4) methods:
if(require("Rmpfr")) # mostly interesting there
  showMethods(class="bigq")

# an  sapply() version that works for big rationals "bigq":
sapplyQ <- function(X, FUN, ...) c_bigq(lapply(X, FUN, ...))

# dummy example showing it works (here):
qq <- as.bigq(1, 1:999)
q1 <- sapplyQ(qq, function(q) q^2)
stopifnot( identical(q1, qq^2) )

Relational Operators

Description

Binary operators which allow the comparison of values in atomic vectors.

Usage

## S3 method for class 'bigq'
sign(x)

## S3 method for class 'bigq'
e1 < e2
## S3 method for class 'bigq'
e1 <= e2
## S3 method for class 'bigq'
e1 == e2
## S3 method for class 'bigq'
e1 >= e2
## S3 method for class 'bigq'
e1 > e2
## S3 method for class 'bigq'
e1 != e2

Arguments

Examples

  x <- as.bigq(8000,21)
  x  < 2 * x
  

Basic arithmetic operators for large rationals

Description

Addition, subtraction, multiplication, division, and absolute value for large rationals, i.e. "bigq" class R objects.

Usage

add.bigq(e1, e2)
## S3 method for class 'bigq'
e1 + e2

sub.bigq(e1, e2=NULL)
## S3 method for class 'bigq'
e1 - e2

mul.bigq(e1, e2)
## S3 method for class 'bigq'
e1 * e2

div.bigq(e1, e2)
## S3 method for class 'bigq'
e1 / e2

## S3 method for class 'bigq'
e1 ^ e2

## S3 method for class 'bigq'
abs(x)

Arguments

Details

Operators can be use directly when the objects are of class "bigq": a + b, a * b, etc, and a ^ n, where n must be coercable to a biginteger ("bigz").

Value

A bigq class representing the result of the arithmetic operation.

Author(s)

Immanuel Scholz and Antoine Lucas

Examples

## 1/3 + 1 = 4/3 :
as.bigq(1,3) + 1

r <- as.bigq(12, 47)
stopifnot(r ^ 3 == r*r*r)

Large Sized Integer Values

Description

Class "bigz" encodes arbitrarily large integers (via GMP). A simple S3 class (internally a raw vector), it has been registered as formal (S4) class (via setOldClass), too.

Usage

as.bigz(a, mod = NA)
NA_bigz_
## S3 method for class 'bigz'
as.character(x, b = 10, ...)
is.bigz(x)
## S3 method for class 'bigz'
is.na(x)
## S3 method for class 'bigz'
print(x, quote=FALSE, initLine = is.null(modulus(x)), ...)
c_bigz(L)

Arguments

Details

Bigz's are integers of arbitrary, but given length (means: only restricted by the host memory). Basic arithmetic operations can be performed on bigzs as addition, subtraction, multiplication, division, modulation (remainder of division), power, multiplicative inverse, calculating of the greatest common divisor, test whether the integer is prime and other operations needed when performing standard cryptographic operations.

For a review of basic arithmetics, see add.bigz.

Comparison are supported, i.e., "==", "!=", "<", "<=", ">", and ">=".

NA_bigz_ is computed on package load time as as.bigz(NA).

Objects of class "bigz" may have a “modulus”, accessible via modulus(), currently as an attribute mod. When the object has such a modulus m, arithmetic is performed “modulo m”, mathematically “within the ring Z/mZ”. For many operations, this means

   result <- mod.bigz(result, m)  ## == result %% m

is called after performing the arithmetic operation and the result will have the attribute mod set accordingly. This however does not apply, e.g., for /, where a / b := a b^{-1} and b^{-1} is the multiplicate inverse of b with respect to ring arithmetic, or NA with a warning when the inverse does not exist. The warning can be turned off via options("gmp:warnModMismatch" = FALSE)

Powers of bigzs can only be performed, if either a modulus is going to be applied to the result bigz or if the exponent fits into an integer value. So, if you want to calculate a power in a finite group (“modulo c”), for large c do not use a ^ b %% c, but rather as.bigz(a,c) ^ b.

The following rules for the result's modulus apply when performing arithmetic operations on bigzs:

• If none of the operand has a modulus set, the result will not have a modulus.

• If both operands have a different modulus, the result will not have a modulus, except in case of mod.bigz, where the second operand's value is used.

• If only one of the operands has a modulus or both have a common (the same), it is set and used for the arithmetic operations, except in case of mod.bigz, where the second operand's value is used.

Value

An R object of (S3) class "bigz", representing the argument (x or a).

Note

    x <- as.bigz(1234567890123456789012345678901234567890)
  

will not work as R converts the number to a double, losing precision and only then convert to a "bigz" object.

Instead, use the syntax

    x <- as.bigz("1234567890123456789012345678901234567890")
  

Author(s)

Immanuel Scholz

References

The GNU MP Library, see https://gmplib.org

Examples

## 1+1=2
a <- as.bigz(1)
a + a

## Two non-small Mersenne primes:
two <- as.bigz(2)
p1 <- two^107 -1 ; isprime(p1); p1
p2 <- two^127 -1 ; isprime(p2); p2

stopifnot( is.na(NA_bigz_) )

## Calculate c = x^e mod n
## --------------------------------------------------------------------
x <- as.bigz("0x123456789abcdef") # my secret message
e <- as.bigz(3) # something smelling like a dangerous public RSA exponent
(n <- p1 * p2) #  a product of two primes
as.character(n, b=16)# as both primes were Mersenne's..

## recreate the three numbers above [for demo below]:
n. <- n; x. <- x; e. <- e # save
Rev <- function() { n <<- n.; x <<- x.; e <<- e.}

# first way to do it right
modulus(x) <- n
c <- x ^ e ; c ; Rev()

# similar second way (makes more sense if you reuse e) to do it right
modulus(e) <- n
c2 <- x ^ e
stopifnot(identical(c2, c), is.bigz(c2)) ; Rev()

# third way to do it right
c3 <- x ^ as.bigz(e, n) ; stopifnot(identical(c3, c))

# fourth way to do it right
c4 <- as.bigz(x, n) ^ e ; stopifnot(identical(c4, c))

# WRONG! (although very beautiful. Ok only for very small 'e' as here)
cc <- x ^ e %% n
cc == c

# Return result in hexa
as.character(c, b=16)

# Depict the "S4-class" bigz, i.e., the formal (S4) methods:
if(require("Rmpfr")) # mostly interesting there
  showMethods(class="bigz")

# an  sapply() version that works for big integers "bigz":
sapplyZ <- function(X, FUN, ...) c_bigz(lapply(X, FUN, ...))

# dummy example showing it works (here):
zz <- as.bigz(3)^(1000+ 1:999)
z1 <- sapplyZ(zz, function(z) z^2)
stopifnot( identical(z1, zz^2) )

Basic Arithmetic Operators for Large Integers ("bigz")

Description

Addition, substraction, multiplication, (integer) division, remainder of division, multiplicative inverse, power and logarithm functions.

Usage

add.bigz(e1, e2)
sub.bigz(e1, e2 = NULL)
mul.bigz(e1, e2)
div.bigz(e1, e2)
divq.bigz(e1,e2) ## ==  e1 %/% e2
mod.bigz(e1, e2) ## ==  e1 %%  e2
## S3 method for class 'bigz'
abs(x)

inv.bigz(a, b,...)## == (1 / a) (modulo b)
pow.bigz(e1, e2,...)## == e1 ^ e2
## S3 method for class 'bigz'
log(x, base=exp(1))
## S3 method for class 'bigz'
log2(x)
## S3 method for class 'bigz'
log10(x)

Arguments

Details

Operators can be used directly when objects are of class bigz: a + b, log(a), etc.

For details about the internal modulus state, and the rules applied for arithmetic operations on big integers with a modulus, see the bigz help page.

a / b = div(a,b) returns a rational number unless the operands have a (matching) modulus where a * b^-1 results.
a %/% b (or, equivalently, divq(a,b)) returns the quotient of simple integer division (with truncation towards zero), possibly re-adding a modulus at the end (but not using a modulus like in a / b).

r <- inv.bigz(a, m), the multiplicative inverse of a modulo m, corresponds to 1/a or a ^-1 from above when a has modulus m. Note that a not always has an inverse modulo m, in which case r will be NA with a warning that can be turned off via

options("gmp:warnNoInv" = FALSE)

.

Value

Apart from / (or div), where rational numbers (bigq) may result, these functions return an object of class "bigz", representing the result of the arithmetic operation.

Author(s)

Immanuel Scholz and Antoine Lucas

References

The GNU MP Library, see https://gmplib.org

Examples

# 1+1=2
as.bigz(1) + 1
as.bigz(2)^10
as.bigz(2)^200

# if my.large.num.string is set to a number, this returns the least byte
(my.large.num.string <- paste(sample(0:9, 200, replace=TRUE), collapse=""))
mod.bigz(as.bigz(my.large.num.string), "0xff")

# power exponents can be up to MAX_INT in size, or unlimited if a
# bigz's modulus is set.
pow.bigz(10,10000)

## Modulo 11,   7 and 8 are inverses :
as.bigz(7, mod = 11) * 8 ## ==>  1  (mod 11)
inv.bigz(7, 11)## hence, 8
a <- 1:10
(i.a <- inv.bigz(a, 11))
d <- as.bigz(7)
a %/% d  # = divq(a, d)
a %%  d  # = mod.bigz (a, d)

(ii <- inv.bigz(1:10, 16))
## with 5 warnings (one for each NA)
op <- options("gmp:warnNoInv" = FALSE)
i2 <- inv.bigz(1:10, 16) # no warnings
(i3 <- 1 / as.bigz(1:10, 16))
i4 <- as.bigz(1:10, 16) ^ -1
stopifnot(identical(ii, i2),
	  identical(as.bigz(i2, 16), i3),
	  identical(i3, i4))
options(op)# revert previous options' settings

stopifnot(inv.bigz(7, 11) == 8,
          all(as.bigz(i.a, 11) * a == 1),
          identical(a %/% d, divq.bigz(1:10, 7)),
          identical(a %%  d, mod.bigz (a, d))
 )

Exact Rational Binomial Probabilities

Description

Compute exact binomial probabilities using (big integer and) big rational arithmetic.

Usage

dbinomQ(x, size, prob, log = FALSE)

Arguments

Value

a big rational ("bigq") of the length of (recycled) x+size+prob.

Author(s)

Martin Maechler

See Also

chooseZ; R's (stats package) dbinom().

Examples

dbinomQ(0:8,8, as.bigq(1,2))
##  1/256  1/32   7/64   7/32   35/128 7/32   7/64   1/32   1/256

ph16. <- dbinomQ(0:16, size=16, prob = 1/2)  # innocous warning
ph16  <- dbinomQ(0:16, size=16, prob = as.bigq(1,2))
ph16.75 <- dbinomQ(0:16, size=16, prob = as.bigq(3,4))
ph8.75  <- dbinomQ(0:8, 8, as.bigq(3,4))
stopifnot(exprs = {
   dbinomQ(0:8,8, as.bigq(1,2)) * 2^8 == choose(8, 0:8)
   identical(ph8.75, chooseZ(8,0:8) * 3^(0:8) / 4^8)
   all.equal(ph8.75, choose (8,0:8) * 3^(0:8) / 4^8, tol=1e-15) # see exactly equal
   identical(ph16, ph16.)
   identical(ph16,
            dbinomQ(0:16, size=16, prob = as.bigz(1)/2))
   all.equal(dbinom(0:16, 16, prob=1/2), asNumeric(ph16),    tol=1e-15)
   all.equal(dbinom(0:16, 16, prob=3/4), asNumeric(ph16.75), tol=1e-15)
})

(Cumulative) Sums, Products of Large Integers and Rationals

Description

Theses are methods to ‘overload’ the sum(), cumsum() and prod() functions for big rationals and big integers.

Usage

## S3 method for class 'bigz'
cumsum(x)
## S3 method for class 'bigq'
cumsum(x)
## S3 method for class 'bigz'
sum(..., na.rm = FALSE)
## S3 method for class 'bigq'
sum(..., na.rm = FALSE)
## S3 method for class 'bigz'
prod(..., na.rm = FALSE)
## S3 method for class 'bigq'
prod(..., na.rm = FALSE)

Arguments

Value

return an element of class bigz or bigq.

Author(s)

Antoine Lucas

See Also

apply

Examples

 x <- as.bigz(1:12)
 cumsum(x)
 prod(x)
 sum(x)

 x <- as.bigq(1:12)
 cumsum(x)
 prod(x)
 sum(x)

Extract or Replace Parts of a 'bigz' or 'bigq' Object

Description

Operators acting on vectors, arrays and lists to extract or replace subsets.

Usage

## S3 method for class 'bigz'
x[i=NULL, j=NULL, drop = TRUE]
## S3 method for class 'bigq'
x[i=NULL, j=NULL, drop = TRUE]
##___ In the following, only the bigq method is mentioned (but 'bigz' is "the same"): ___
## S3 method for class 'bigq'
c(..., recursive = FALSE)
## S3 method for class 'bigq'
rep(x, times=1, length.out=NA, each=1, ...)

Arguments

Examples

  a <- as.bigz(123)
  ## indexing "outside" --> extends the vectors (filling with NA)
  a[2] <- a[1]
  a[4] <- -4

  ## create a vector of 3 a
  c(a,a,a)

  ## repeate a 5 times
  rep(a,5)

  ## with matrix: 3 x 2
  m <- matrix.bigz(1:6,3)

  m[1,] # the first row
  m[1,, drop=TRUE] # the same: drop does *not* drop
  m[1]
  m[-c(2,3),]
  m[-c(2,3)]
  m[c(TRUE,FALSE,FALSE)]

  ##_modification on matrix
  m[2,-1] <- 11

Extrema (Maxima and Minima)

Description

We provide S3 methods for min and max for big rationals (bigq) and big integers (biqz); consequently, range() works as well.

Similarly, S4 methods are provided for which.min() and which.max().

Usage

## S3 method for class 'bigz'
max(..., na.rm=FALSE)
## S3 method for class 'bigq'
max(..., na.rm=FALSE)
## S3 method for class 'bigz'
min(..., na.rm=FALSE)
## S3 method for class 'bigq'
min(..., na.rm=FALSE)

## S4 method for signature 'bigz'
which.min(x)
## S4 method for signature 'bigq'
which.max(x)

Arguments

Value

an object of class "bigz" or "bigq".

Author(s)

Antoine Lucas

See Also

max etc in base.

Examples

 x <- as.bigz(1:10)
 max(x)
 min(x)
 range(x) # works correctly via default method
 x <- x[c(7:10,6:3,1:2)]
 which.min(x) ## 9
 which.max(x) ## 4

 Q <- as.bigq(1:10, 3)
 max(Q)
 min(Q)
 (Q <- Q[c(6:3, 7:10,1:2)])
 stopifnot(which.min(Q) == which.min(asNumeric(Q)),
           which.max(Q) == which.max(asNumeric(Q)))

stopifnot(range(x) == c(1,10), 3*range(Q) == c(1,10))

Factorial and Binomial Coefficient as Big Integer

Description

Efficiently compute the factorial n! or a binomial coefficient {n\choose k} as big integer (class bigz).

Usage

factorialZ(n)
chooseZ(n, k)

Arguments

Value

a vector of big integers, i.e., of class bigz.

See Also

factorial and gamma in base R;

Examples

factorialZ(0:10)# 1 1 2 6 ... 3628800
factorialZ(0:40)# larger
factorialZ(200)

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

system.time(replicate(8, f1e4 <<- factorialZ(10000)))
nchar(as.character(f1e4))# 35660 ... (too many to even look at ..)

chooseZ(1000, 100:102)# vectorizes
chooseZ(as.bigz(2)^120, 10)
n <- c(50,80,100)
k <- c(20,30,40)
## currently with an undesirable warning: % from methods/src/eval.c  _FIXME_
stopifnot(chooseZ(n,k) == factorialZ(n) / (factorialZ(k)*factorialZ(n-k)))

Factorize a number

Description

Give all primes numbers to factor the number

Usage

factorize(n)

Arguments

Details

The factorization function uses the Pollard Rho algorithm.

Value

Vector of class bigz.

Author(s)

Antoine Lucas

References

The GNU MP Library, see https://gmplib.org

Examples

factorize(34455342)

Format Numbers Keeping Classes Distinguishable

Description

Format (generalized) numbers in a way that their classes are distinguishable. Contrary to format() which uses a common format for all elements of x, here, each entry is formatted individually.

Usage

formatN(x, ...)
## Default S3 method:
formatN(x, ...)
## S3 method for class 'integer'
formatN(x, ...)
## S3 method for class 'double'
formatN(x, ...)
## S3 method for class 'bigz'
formatN(x, ...)
## S3 method for class 'bigq'
formatN(x, ...)

Arguments

Value

a character vector of the same length as x, each entry a representation of the corresponding entry in x.

Author(s)

Martin Maechler

See Also

format, including its (sophisticated) default method; as.character.

Examples

## Note that each class is uniquely recognizable from its output:
formatN(    -2:5)# integer
formatN(0 + -2:5)# double precision
formatN(as.bigz(-2:5))
formatN(as.bigq(-2:5, 4))

Split Number into Fractional and Exponent of 2 Parts

Description

Breaks the number x into its binary significand (“fraction”) d \in [0.5, 1) and ex, the integral exponent for 2, such that x = d \cdot 2^{ex}.

If x is zero, both parts (significand and exponent) are zero.

Usage

frexpZ(x)

Arguments

Value

a list with the two components

Author(s)

Martin Maechler

See Also

log2, etc; for bigz objects built on (the C++ equivalent of) frexp(), actually GMP's ‘⁠mpz_get_d_2exp()⁠’.

Examples

frexpZ(1:10)
## and confirm :
with(frexpZ(1:10),  d * 2^exp)
x <- rpois(1000, lambda=100) * (1 + rpois(1000, lambda=16))
X <- as.bigz(x)
stopifnot(all.equal(x, with(frexpZ(x), d* 2^exp)),
          1+floor(log2(x)) == (fx <- frexpZ(x)$exp),
          fx == frexpZ(X)$exp,
          1+floor(log2(X)) == fx
)

Greatest Common Divisor (GCD) and Least Common Multiple (LCM)

Description

Compute the greatest common divisor (GCD) and least common multiple (LCM) of two (big) integers.

Usage

## S3 method for class 'bigz'
gcd(a, b)
lcm.bigz(a, b)

Arguments

Value

An element of class bigz

Author(s)

Antoine Lucas

References

The GNU MP Library, see https://gmplib.org

See Also

gcdex

Examples

gcd.bigz(210,342) # or also
lcm.bigz(210,342)
a <- 210 ; b <- 342
stopifnot(gcd.bigz(a,b) * lcm.bigz(a,b) == a * b)

## or
(a <- as.bigz("82696155787249022588"))
(b <- as.bigz("65175989479756205392"))
gcd(a,b) # 4
stopifnot(gcd(a,b) * lcm.bigz(a,b) == a * b)

Compute Bezoult Coefficient

Description

Compute g,s,t as as + bt = g = gcd(a,b). s and t are also known as Bezoult coefficients.

Usage

gcdex(a, b)

Arguments

Value

a class "bigz" vector of length 3 with (long integer) values g, s, t.

Author(s)

Antoine Lucas

References

The GNU MP Library, see https://gmplib.org

See Also

gcd.bigz

Examples

gcdex(342,654)

Base Functions in 'gmp'-ified Versions

Description

Functions from base etc which need a copy in the gmp namespace so they correctly dispatch.

Usage

outer(X, Y, FUN = "*", ...)

Arguments

See Also

outer.

Examples

twop <- as.bigz(2)^(99:103)
(mtw <- outer(twop, 0:2))
stopifnot(
   identical(dim(mtw), as.integer(c(5,3)))
 ,
   mtw[,1] == 0
 ,
   identical(as.vector(mtw[,2]), twop)
)

GMP Number Utilities

Description

gmpVersion() returns the version of the GMP library which gmp is currently linked to.

Usage

gmpVersion()

References

The GNU MP Library, see https://gmplib.org

Examples

gmpVersion()

Whole ("Integer") Numbers

Description

Check which elements of x[] are integer valued aka “whole” numbers.

Usage

is.whole(x)
## Default S3 method:
is.whole(x)
## S3 method for class 'bigz'
is.whole(x)
## S3 method for class 'bigq'
is.whole(x)

Arguments

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() method for "mpfr" numbers.

Examples

 is.integer(3) # FALSE, it's internally a double
 is.whole(3)   # TRUE
 ## integer valued complex numbers  (two FALSE) :
 is.whole(c(7, 1 + 1i, 1.2, 3.4i, 7i))
 is.whole(factorialZ(20)^(10:12)) ## "bigz" are *always* whole numbers
 q <- c(as.bigz(36)^50 / as.bigz(30)^40, 3, factorialZ(30:31), 12.25)
 is.whole(q) # F T T T F

Determine if number is (very probably) prime

Description

Determine whether the number n is prime or not, with three possible answers:

2:

n is prime,

1:

n is probably prime (without beeing certain),

0:

n is composite.

Usage

isprime(n, reps = 40)  

Arguments

Details

This function does some trial divisions, then some Miller-Rabin probabilistic primary tests. reps controls how many such tests are done, 5 to 10 is already a resonable number. More will reduce the chances of a composite being returned as “probably prime”.

Value

Author(s)

Antoine Lucas

References

The GNU MP Library, see https://gmplib.org

See Also

nextprime, factorize.

Note that for “small” n, which means something like n < 10'000'000, non-probabilistic methods (such as factorize()) are fast enough.

Examples

isprime(210)
isprime(71)

# All primes numbers from 1 to 100
t <- isprime(1:99)
(1:99)[t > 0]

table(isprime(1:10000))# 0 and 2 : surely prime or not prime

primes <- function(n) {
  ## all primes <= n
  stopifnot(length(n) == 1, n <= 1e7) # be reasonable
  p <- c(2L, as.integer(seq(3, n, by=2)))
  p[isprime(p) > 0]
}

## quite quickly, but for these small numbers
## still slower than e.g., sfsmisc::primes()
system.time(p100k <- primes(100000))

## The first couple of Mersenne primes:
p.exp <- primes(1000)
Mers <- as.bigz(2) ^ p.exp - 1
isp.M <- sapply(seq_along(Mers), function(i) isprime(Mers[i], reps=256))
cbind(p.exp, isp.M)[isp.M > 0,]
Mers[isp.M > 0]

Compute Fibonacci and Lucas numbers

Description

fibnum compute n-th Fibonacci number. fibnum2 compute (n-1)-th and n-th Fibonacci number. lucnum compute n-th lucas number. lucnum2 compute (n-1)-th and n-th lucas number.

Fibonacci numbers are define by: F_n=F_{n-1}+F_{n-2} Lucas numbers are define by: L_n=F_n+2F_{n-1}

Usage

fibnum(n)
fibnum2(n)
lucnum(n)
lucnum2(n)

Arguments

Value

Fibonacci numbers and Lucas number.

Author(s)

Antoine Lucas

References

The GNU MP Library, see https://gmplib.org

Examples

fibnum(10)
fibnum2(10)
lucnum(10)
lucnum2(10)

Matrix manipulation with gmp

Description

Overload of “all” standard tools useful for matrix manipulation adapted to large numbers.

Usage

## S3 method for class 'bigz'
matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL, mod = NA,...)

is.matrixZQ(x)

## S3 method for class 'bigz'
x %*% y
## S3 method for class 'bigq'
x %*% y
## S3 method for class 'bigq'
crossprod(x, y=NULL,...)
## S3 method for class 'bigz'
tcrossprod(x, y=NULL,...)

## S3 method for class 'bigz'
cbind(..., deparse.level=1)
## S3 method for class 'bigq'
rbind(..., deparse.level=1)
## ..... etc

Arguments

Details

The extract function ("[") is the same use for vector or matrix. Hence, x[i] returns the same values as x[i,]. This is not considered a feature and may be changed in the future (with warnings).

All matrix multiplications should work as with numeric matrices.

Special features concerning the "bigz" class: the modulus can be

Unset:

Just play with large numbers

Set with a vector of size 1:

Example: matrix.bigz(1:6,nrow=2,ncol=3,mod=7) This means you work in Z/nZ, for the whole matrix. It is the only case where the %*% and solve functions will work in Z/nZ.

Set with a vector smaller than data:

Example: matrix.bigz(1:6,nrow=2,ncol=3,mod=1:5). Then, the modulus is repeated to the end of data. This can be used to define a matrix with a different modulus at each row.

Set with same size as data:

Modulus is defined for each cell

Value

matrix(): A matrix of class "bigz" or "bigq".

is.matrixZQ(): TRUE or FALSE.

dim(), ncol(), etc: integer or NULL, as for simple matrices.

cbind(x,y,...) and rbind(x,y,...) now (2024-01, since gmp version 0.9-5), do drop deparse.level=. instead of wrongly creating an extra column or row and the "bigz" method takes all arguments into account and calls the "bigq" method in case of arguments inheriting from "bigq".

Author(s)

Antoine Lucas and Martin Maechler

See Also

Solving a linear system: solve.bigz. matrix

Examples

V <- as.bigz(v <- 3:7)
crossprod(V)# scalar product
(C <- t(V))
stopifnot(dim(C) == dim(t(v)), C == v,
          dim(t(C)) == c(length(v), 1),
          crossprod(V) == sum(V * V),
         tcrossprod(V) == outer(v,v),
          identical(C, t(t(C))),
          is.matrixZQ(C), !is.matrixZQ(V), !is.matrixZQ(5)
	)

## a matrix
x <- diag(1:4)
## invert this matrix
(xI <- solve(x))

## matrix in Z/7Z
y <- as.bigz(x,7)
## invert this matrix (result is *different* from solve(x)):
(yI <- solve(y))
stopifnot(yI %*% y == diag(4),
          y %*% yI == diag(4))

## matrix in Q
z  <- as.bigq(x)
## invert this matrix (result is the same as solve(x))
(zI <- solve(z))

stopifnot(abs(zI - xI) <= 1e-13,
          z %*% zI == diag(4),
          identical(crossprod(zI), zI %*% t(zI))
         )

A <- matrix(2^as.bigz(1:12), 3,4)
for(a in list(A, as.bigq(A, 16), factorialZ(20), as.bigq(2:9, 3:4))) {
  a.a <- crossprod(a)
  aa. <- tcrossprod(a)
  stopifnot(identical(a.a, crossprod(a,a)),
 	    identical(a.a, t(a) %*% a)
            ,
            identical(aa., tcrossprod(a,a)),
	    identical(aa., a %*% t(a))
 	   )
}# {for}

Modulus of a Big Integer

Description

The modulus of a bigz number a is “unset” when a is a regular integer, a \in Z). Or the modulus can be set to m which means a \in Z/\,m\cdot Z), i.e., all arithmetic with a is performed ‘modulo m’.

Usage

modulus(a)
modulus(a) <- value

Arguments

Examples

x <- as.bigz(24)
modulus(x) # NULL, i.e. none

# x element of Z/31Z :
modulus(x) <- 31
x+x  #  48 |-> (17 %% 31)
10*x # 240 |-> (23 %% 31)
x31 <- x

# reset modulus to "none":
modulus(x) <- NA; x; x. <- x
x <- x31
modulus(x) <- NULL; x

stopifnot(identical(x,            as.bigz(24)), identical(x, x.),
          identical(modulus(x31), as.bigz(31)))

Exported function for mpfr use

Description

Theses hidden function are provided for mpfr use. Use theses function with care.

Usage

.as.bigz(a, mod=NA)

Arguments

Value

An R object of (S3) class "bigz", representing the argument (x or a).

References

The GNU MP Library, see https://gmplib.org

Examples

.as.bigz(1)

Next Prime Number

Description

Return the next prime number, say p, with p > n.

Usage

nextprime(n)

Arguments

Details

This function uses probabilistic algorithm to identify primes. For practical purposes, it is adequate, the chance of a composite passing will be extremely small.

Value

A (probably) prime number

Author(s)

Antoine Lucas

References

The GNU MP Library, see https://gmplib.org

See Also

isprime and its references and examples.

Examples

nextprime(14)
## still very fast:
(p <- nextprime(1e7))
## to be really sure { isprime()  gives "probably prime" } :
stopifnot(identical(p, factorize(p)))

RFC 2409 Oakley Groups - Parameters for Diffie-Hellman Key Exchange

Description

RFC 2409 standardizes global unique prime numbers and generators for the purpose of secure asymmetric key exchange on the Internet.

Usage

data(Oakley1)
data(Oakley2)

Value

Oakley1 returns an object of class bigz for a 768 bit Diffie-Hellman group. The generator is stored as value with the respective prime number as modulus attribute.

Oakley2 returns an object of class bigz for a 1024 bit Diffie-Hellman group. The generator is stored as value with the respective prime number as modulus attribute.

References

The Internet Key Exchange (RFC 2409), Nov. 1998

Examples

packageDescription("gmp") # {possibly useful for debugging}

data(Oakley1)
(M1 <- modulus(Oakley1))
isprime(M1)# '1' : "probably prime"
sizeinbase(M1)#  232 digits (was 309 in older version)

Exponentiation function

Description

This function return x ^ y mod n.

This function return x ^ y mod n pow.bigz do the same when modulus is set.

Usage

powm(x, y, n)

Arguments

Value

A bigz class representing the parameter value.

Author(s)

A. L.

See Also

pow.bigz

Examples

powm(4,7,9)


x = as.bigz(4,9)
x ^ 7

Generate a random number

Description

Generate a uniformly distributed random number in the range 0 to 2^{size} -1, inclusive.

Usage

urand.bigz(nb=1,size=200, seed = 0)

Arguments

Value

A biginteger of class bigz.

Author(s)

Antoine Lucas

References

‘⁠mpz\_urandomb⁠’ from the GMP Library, see https://gmplib.org

Examples

# Integers are differents
urand.bigz()
urand.bigz()
urand.bigz()

# Integers are the same
urand.bigz(seed="234234234324323")
urand.bigz(seed="234234234324323")

# Vector
urand.bigz(nb=50,size=30)

Relational Operators

Description

Binary operators which allow the comparison of values in atomic vectors.

Usage

## S3 method for class 'bigz'
sign(x)
## S3 method for class 'bigz'
e1 == e2
## S3 method for class 'bigz'
 e1 < e2
## S3 method for class 'bigz'
e1 >= e2

Arguments

See Also

mod.bigz for arithmetic operators.

Examples

  x <- as.bigz(8000)
  x ^ 300 < 2 ^x

  sign(as.bigz(-3:3))
  sign(as.bigq(-2:2, 7))

Rounding Big Rationals ("bigq") to Decimals

Description

Rounding big rationals (of class "bigq", see as.bigq()) to decimal digits is strictly based on a (optionally choosable) definition of rounding to integer, i.e., digits = 0, the default method of which we provide as round0().

The users typically just call round(x, digits) as elsewhere, and the round() method will call round(x, digits, round0=round0).

Usage

round0(x)

roundQ(x, digits = 0, r0 = round0)

## S3 method for class 'bigq'
round(x, digits = 0)

Arguments

Value

round0() returns a vector of big integers, i.e., "bigz" classed.

roundQ(x, digits, round0) returns a vector of big rationals, "bigq", as x.

round.bigq is very simply defined as function(x, digits) roundQ(x, digits) .

Author(s)

Martin Maechler, ETH Zurich

References

The vignette “Exact Decimal Rounding via Rationals” from CRAN package round,

Wikipedia, Rounding, notably "Round half to even": https://en.wikipedia.org/wiki/Rounding#Round_half_to_even

See Also

round for (double precision) numbers in base R; roundX from CRAN package round.

Examples

qq <- as.bigq((-21:31), 10)
noquote(cbind(as.character(qq), asNumeric(qq)))
round0(qq) # Big Integer ("bigz")
## corresponds to R's own "round to even" :
stopifnot(round0(qq) == round(asNumeric(qq)))
round(qq) # == round(qq, 0): the same as round0(qq) *but* Big Rational ("bigq")

halfs <- as.bigq(1,2) + -5:12



## round0() is simply
round0 <- function (x) {
    nU <- as.bigz.bigq(xU <- x + as.bigq(1, 2)) # traditional round: .5 rounded up
    if(any(I <- is.whole.bigq(xU))) { # I <==>  x == <n>.5 : "hard case"
        I[I] <- .mod.bigz(nU[I], 2L) == 1L # rounded up is odd  ==> round *down*
        nU[I] <- nU[I] - 1L
    }
    nU
}

## 's' for simple: rounding as you learned in school:
round0s <- function(x) as.bigz.bigq(x + as.bigq(1, 2))

cbind(halfs, round0s(halfs), round0(halfs))


## roundQ() is simply
roundQ <- function(x, digits = 0, r0 = round0) {
    ## round(x * 10^d) / 10^d --  vectorizing in both (x, digits)
    p10 <- as.bigz(10) ^ digits # class: if(all(digits >= 0)) "bigz" else "bigq"
    r0(x * p10) / p10
}

Compute size of a bigz in a base

Description

Return an approximation to the number of character the integer X would have printed in base b. The approximation is never too small.

In case of powers of 2, function gives exact result.

Usage

sizeinbase(a, b=10)

Arguments

Value

integer of the same length as a: the size, i.e. number of digits, of each a[i].

Author(s)

Antoine Lucas

References

The GNU MP Library, see https://gmplib.org

Examples

sizeinbase(342434, 10)# 6 obviously

Iv <- as.bigz(2:7)^500
sizeinbase(Iv)
stopifnot(sizeinbase(Iv)       == nchar(as.character(Iv)),
          sizeinbase(Iv, b=16) == nchar(as.character(Iv, b=16)))

Solve a system of equation

Description

This generic function solves the equation a \%*\% x = b for x, where b can be either a vector or a matrix.

If a and b are rational, return is a rational matrix.

If a and b are big integers (of class bigz) solution is in Z/nZ if there is a common modulus, or a rational matrix if not.

Usage

## S3 method for class 'bigz'
solve(a, b, ...)
## S3 method for class 'bigq'
solve(a, b, ...)

Arguments

Details

It uses the Gauss and trucmuch algo ... (to be detailled).

Value

If a and b are rational, return is a rational matrix.

If a and b are big integers (of class bigz) solution is in Z/nZ if there is a common modulus, of a rational matrix if not.

Author(s)

Antoine Lucas

See Also

solve

Examples

x <- matrix(1:4,2,2)  ## standard solve :
solve(x)

q <- as.bigq(x) ## solve with rational
solve(q)

z <- as.bigz(x)
modulus(z) <- 7  ## solve in Z/7Z :
solve(z)

b <- c(1,3)
solve(q,b)
solve(z,b)

## Inversion of ("non-trivial") rational matrices :

A <- rbind(c(10, 1,  3),
           c( 4, 2, 10),
           c( 1, 8,  2))
(IA.q <- solve(as.bigq(A))) # fractions..
stopifnot(diag(3) == A %*% IA.q)# perfect

set.seed(5); B <- matrix(round(9*runif(5^2, -1,1)), 5)
B
(IB.q <- solve(as.bigq(B)))
stopifnot(diag(5) == B %*% IB.q, diag(5) == IB.q %*% B,
          identical(B, asNumeric(solve(IB.q))))

Eulerian and Stirling Numbers of First and Second Kind

Description

Compute Eulerian numbers and Stirling numbers of the first and second kind, possibly vectorized for all k “at once”.

Usage

Stirling1(n, k)
Stirling2(n, k, method = c("lookup.or.store", "direct"))
Eulerian (n, k, method = c("lookup.or.store", "direct"))

Stirling1.all(n)
Stirling2.all(n)
Eulerian.all (n)

Arguments

Details

Eulerian numbers:
A(n,k) = the number of permutations of 1,2,...,n with exactly k ascents (or exactly k descents).

Stirling numbers of the first kind:
s(n,k) = (-1)^{n-k} times the number of permutations of 1,2,...,n with exactly k cycles.

Stirling numbers of the second kind:
S^{(k)}_n is the number of ways of partitioning a set of n elements into k non-empty subsets.

Value

A(n,k), s(n,k) or S(n,k) = S^{(k)}_n, respectively.

Eulerian.all(n) is the same as sapply(0:(n-1), Eulerian, n=n) (for n > 0),
Stirling1.all(n) is the same as sapply(1:n, Stirling1, n=n), and
Stirling2.all(n) is the same as sapply(1:n, Stirling2, n=n), but more efficient.

Note

For typical double precision arithmetic,
Eulerian*(n, *) overflow (to Inf) for n \ge 172,
Stirling1*(n, *) overflow (to \pmInf) for n \ge 171, and
Stirling2*(n, *) overflow (to Inf) for n \ge 220.

Author(s)

Martin Maechler ("direct": May 1992)

References

Eulerians:

NIST Digital Library of Mathematical Functions, 26.14: https://dlmf.nist.gov/26.14

Stirling numbers:

Abramowitz and Stegun 24,1,4 (p. 824-5 ; Table 24.4, p.835); Closed Form : p.824 "C."

NIST Digital Library of Mathematical Functions, 26.8: https://dlmf.nist.gov/26.8

See Also

chooseZ for the binomial coefficients.

Examples

Stirling1(7,2)
Stirling2(7,3)

stopifnot(
 Stirling1.all(9) == c(40320, -109584, 118124, -67284, 22449, -4536, 546, -36, 1)
 ,
 Stirling2.all(9) == c(1, 255, 3025, 7770, 6951, 2646, 462, 36, 1)
 ,
 Eulerian.all(7) == c(1, 120, 1191, 2416, 1191, 120, 1)
)