Commit 544f1613 authored by Martin Maechler's avatar Martin Maechler
Browse files

BernoulliQ(n) provides exact n-th Bernoulli numbers (version 0.6-1)

parent 7e7448c2
2020-06-09 Martin Maechler <maechler@stat.math.ethz.ch>
* DESCRIPTION (Version): 0.6-1
* R/Stirling-n-etc.R (BernoulliQ): Bernoulli numbers as exact fractions.
* R/zzz.R (NA_bigz_, NA_bigq_): analogous to `NA_integer_` etc
2020-06-06 Martin Maechler <maechler@stat.math.ethz.ch>
* NAMESPACE: use *same* S3method(<op>, <class>, <method_fn>) for bigz and bigq,
......
Package: gmp
Version: 0.6-0
Date: 2020-06-06
Version: 0.6-1
Date: 2020-06-10
Title: Multiple Precision Arithmetic
Author: Antoine Lucas, Immanuel Scholz, Rainer Boehme <rb-gmp@reflex-studio.de>,
Sylvain Jasson <Sylvain.Jasson@inrae.fr>,
......@@ -11,7 +11,7 @@ Description: Multiple Precision Arithmetic (big integers and rationals,
using the C library GMP (GNU Multiple Precision Arithmetic).
Depends: R (>= 3.4.0)
Imports: methods
Suggests: Rmpfr
Suggests: Rmpfr, MASS
SystemRequirements: gmp (>= 4.2.3)
License: GPL (>= 2)
BuildResaveData: no
......
......@@ -23,6 +23,7 @@ export(
"denominator", "denominator<-",
## "dim<-.bigq", "dim.bigq", "dim<-.bigz", "dim.bigz",
"div.bigq", "div.bigz", "divq.bigz",
BernoulliQ, # Rmpfr has 'Bernoulli'; want different name
"Eulerian", "Eulerian.all",
"Stirling1", "Stirling1.all", "Stirling2", "Stirling2.all",
"factorize",
......
### Really want S4 methods for all the binary operations.
### Otherwise we "never" make
### <bigz> o <bigq>
### or <bigq> o <bigz> working --
##
## But unfortunately the above seems "impossible", see
## see also setMethod() in ./matrix-prods.R
## OTOH: This *still* helps to define single-dispatch methods for asNumeric() :
## {why does it work there ??}
## This helps to define single-dispatch methods for asNumeric() :
setOldClass("bigz")#, prototype=as.bigz(integer()))
setOldClass("bigq")#, prototype=as.bigq(integer()))
......
......@@ -3,7 +3,6 @@
### The "double prec" version of this is currently in package 'nacopula'
### (MM: >>> ../../nacopula/R/special-func.R )
##' Compute Stirling numbers of the 1st kind
##'
##' s(n,k) = (-1)^{n-k} times
......@@ -32,16 +31,16 @@ Stirling1 <- function(n,k)
}
S
}
if(compute <- (nt <- length(St <- get("S1.tab", .Stirl..env))) < n) {
if(compute <- (nt <- length(St <- .Stirl..env$ S1.tab)) < n) {
## extend the "table":
length(St) <- n
for(i in (nt+1L):n) St[[i]] <- rep(as.bigz(NA), i)
for(i in (nt+1L):n) St[[i]] <- rep.bigz(NA_bigz_, i)
}
else compute <- is.na(S <- St[[n]][k])
if(compute) {
S <- S1(n,k)
## store it back:
assign("S1.tab", St, envir = .Stirl..env)
.Stirl..env$ S1.tab <- St
}
S
}
......@@ -56,11 +55,11 @@ Stirling1.all <- function(n)
{
stopifnot(length(n) == 1)
if(!n) return(as.bigz(numeric(0)))
if(get("S1.full.n", .Stirl..env) < n) {
assign("S1.full.n", n, envir = .Stirl..env)
if(.Stirl..env$ S1.full.n < n) {
.Stirl..env$ S1.full.n <- n
do.call(c, lapply(seq_len(n), Stirling1, n=n))# which fills "S1.tab"
}
else get("S1.tab", .Stirl..env)[[n]]
else .Stirl..env$ S1.tab[[n]]
}
##' Compute Stirling numbers of the 2nd kind
......@@ -99,16 +98,16 @@ Stirling2 <- function(n,k, method = c("lookup.or.store","direct"))
S2(n1, k-1) + k* S2(n1, k) else as.bigz(1) ## n = k = 1
S
}
if(compute <- (nt <- length(St <- get("S2.tab", .Stirl..env))) < n) {
if(compute <- (nt <- length(St <- .Stirl..env$ S2.tab)) < n) {
## extend the "table":
length(St) <- n
for(i in (nt+1L):n) St[[i]] <- rep(as.bigz(NA), i)
for(i in (nt+1L):n) St[[i]] <- rep.bigz(NA_bigz_, i)
}
else compute <- is.na(S <- St[[n]][k])
if(compute) {
S <- S2(n,k)
## store it back:
assign("S2.tab", St, envir = .Stirl..env)
.Stirl..env$ S2.tab <- St
}
S
})
......@@ -124,11 +123,11 @@ Stirling2.all <- function(n)
{
stopifnot(length(n) == 1)
if(!n) return(as.bigz(numeric(0)))
if(get("S2.full.n", .Stirl..env) < n) {
assign("S2.full.n", n, envir = .Stirl..env)
if(.Stirl..env$ S2.full.n < n) {
.Stirl..env$ S2.full.n <- n
do.call(c, lapply(seq_len(n), Stirling2, n=n))# which fills "S2.tab"
}
else get("S2.tab", .Stirl..env)[[n]]
else .Stirl..env$ S2.tab[[n]]
}
......@@ -174,16 +173,16 @@ Eulerian <- function(n,k, method = c("lookup.or.store","direct"))
k1*Eul(n1, k)+ (n-k)* Eul(n1, k-1) else as.bigz(1) ## n=1, k=0
r
}
if(compute <- (nt <- length(E. <- get("Eul.tab", .Stirl..env))) < n) {
if(compute <- (nt <- length(E. <- .Stirl..env$ Eul.tab)) < n) {
## extend the "table":
length(E.) <- n
for(i in (nt+1L):n) E.[[i]] <- rep(as.bigz(NA), i)
for(i in (nt+1L):n) E.[[i]] <- rep.bigz(NA_bigz_, i)
}
else compute <- is.na(E <- E.[[n]][k+1L])
if(compute) {
E <- Eul(n,k)
## store it back:
assign("Eul.tab", E., envir = .Stirl..env)
.Stirl..env$ Eul.tab <- E.
}
E
})
......@@ -199,13 +198,11 @@ Eulerian.all <- function(n)
{
stopifnot(length(n) == 1, n >= 0)
if(!n) return(as.bigz(1))
if(get("Eul.full.n", .Stirl..env) < n) {
## FIXME: do the assign() only when the lapply() below does *not* fail
## on.exit( ... ) ?
assign("Eul.full.n", n, envir = .Stirl..env)
do.call(c, lapply(0:(n-1L), Eulerian, n=n))# which fills "Eul.tab"
if(.Stirl..env$ Eul.full.n < n) {
.Stirl..env$ Eul.full.n <- n
do.call(c, lapply(0:(n-1L), Eulerian, n=n))# which fills "Eul.tab"
}
else get("Eul.tab", .Stirl..env)[[n]]
else .Stirl..env$ Eul.tab[[n]]
}
## Our environment for tables etc: no hash, as it will contain *few* objects:
......@@ -216,3 +213,49 @@ assign("S2.full.n", 0 , envir = .Stirl..env)
assign("S1.full.n", 0 , envir = .Stirl..env)
assign("Eul.tab", list(), envir = .Stirl..env) ## Eul.tab[[n]][k] == A(n, k) == < n \\ k > (Eulerian)
assign("Eul.full.n", 0 , envir = .Stirl..env)
## Bernoulli numbers (have also in 'Rmpfr' via zeta(), but they *are* rational after all!)
.Bernl..env <- new.env(parent=emptyenv(), hash=FALSE)
## Want .Bernl..env$B[n] == B_{2n}, n = 1,2,...
## ==> not storing BernoulliQ(0) = 1, nor BernoulliQ(1) = +1/2
B1 <- function(n) { # stopifnot( length(n) == 1 )
half <- as.bigq(1L, 2L)
if(n == 0L)
as.bigq(1L)
else if(n == 1L)
## B_1 = + 1/2 ("B_n+"; -1/2 according to "old" B_n- definition)
half
else if(n %% 2 == 1L)
as.bigq(0L)
else { ## n in {2, 4, 6, ..}:
n2 <- n %/% 2L
if((lB <- length(B <- .Bernl..env$B)) < n2) { ## compute Bernoulli(n)
## if(verbose) cat("n=",n,": computing from ", lB+1L, sep='', "\n")
if(!lB)
B <- as.bigq(1L, 6L) # B_2 = 1/6
## n2 >= 2; n >= 4
length(B) <- n2 # (fills missing parts with NA_bigq_
## recurse:
if(lB+1L < n2)
for(k in (lB+1L):(n2-1L)) B[k] <- B1(2*k)
k0 <- seq(length=n2-1L) # (1, 2, .., n2-1)
B. <- half - (1 + sum( chooseZ(n+1L, k0+k0) * B[k0])) / (n+1L)
B[n2] <- B.
.Bernl..env$B <- B
B.
}
else
B[n2]
}
}# end B1()
BernoulliQ <- function(n, verbose = getOption("verbose", FALSE)) {
if(!(N <- length(n))) return(as.bigq(n)) # else: length(n) >= 1 :
stopifnot(n >= 0, n == as.integer(n))
if(N == 1L)
B1(n)
else
.Call(bigrational_c, lapply(n, B1))
}
......@@ -138,7 +138,7 @@ as.bigz.bigq <- function(a, mod = NA)
{
## "FIXME": considerably faster in C / C++
if(any(ina <- is.na.bigq(a))) {
r <- as.bigz(rep.int(NA, length(a)))
r <- rep.bigz(NA_bigz_, length(a))
if(any(ii <- !ina)) {
a <- a[ii]
r[ii] <- as.bigz(numerator(a) %/% denominator(a), mod[ii])
......
* Remove full path in configure script, with example of Rmpfr configure. -*-org-*-
* Remove full path in `configure.ac` script, with example of Rmpfr configure. -*-org-*-
* Bug Fixes
** TODO Not working column assignment : A[,j] <- vec does nothing, not warn, nothing -- but A[i,] <- val works ??
** TODO m <- matrix(1:6, 2); Z <- as.bigz(m); str(Z) fails because e.g. Z[5] fails, contrary m[5]
......@@ -7,9 +7,9 @@
** DONE URGENT / BAD : Segfaults, Wrong Results, etc
** TODO completely *wrong* pmin() and pmax(<bigq>, <number>): uses Numerator in comparison ???
*** TODO: Find the *underlying* reason where pmin() goes wrong w/o warning!
** TODO Not working Arithmetic
*** DONE as.bigz(1:4) + as.bigq(7) failed, so do '*', '/', ... Fixed, using the +.bigz, etc _S3_ methods
** TODO Missing bigq %% and %/% [was 'Not working Arithmetic']
*** TODO bigq %% and %/% are not even defined (but are for 'numeric')!
*** DONE as.bigz(1:4) + as.bigq(7) failed, so do '*', '/', ... Fixed (for 0.6-0), via new add.big() etc _S3_ methods
*** DONE more checking: q <- 2.3; Fn <- gmp:::`^.bigq`; Fn(q,q) silently gave wrong result; now a good error message
*** DONE <bigz> %*% <bigq> silently gave wrong result; now works, as does crossprod(), tcrossprod()
*** DONE bigz %% and %/% should be R compatible, notably fulfilling
......@@ -21,7 +21,9 @@
see matOuter() and eqA[] in tests/arith-ex.R
* Miscellaneous
** TODO provide as.bigq.fractions(): to accept MASS::fractions() results
** TODO [MM:] is.whole() conflicts with sfsmisc::is.whole(): maybe get that as default method ?
*** Note that we use all(is.whole(.)) importantly. --> Make fast all_whole() via C ?
** gmp <--> Rmpfr:
*** Z <- as.bigz(1:7); Z[1] <- mpfr(1, 20) now gives an error {did segfault!};
but later, *with* Rmpfr we should make Z into an 'mpfr'.
......@@ -50,7 +52,7 @@ This is completely analogous to
enhance as.bigz() and as.bigq() to work for this case !
** Vectorization; other "desiderata"
*** fibnum(), etc are not vectorized -- learn from chooseZ() !
*** fibnum(), etc are *not* vectorized -- learn from chooseZ() !
*** allow +/- Inf "bigq" ?? -- +Inf :=: 1/0 -Inf :=: (-1)/0
This would allow the bigq (and bigz) arithmetic to be closer to the
......@@ -67,10 +69,10 @@ This is completely analogous to
One possibility: Using +- <LARGE>/1 for "bigq", as +- <LARGE> for "bigz"
** Matrix things
*** 3*M and M*3 are not the same: the first forgets the "nrow"
*** TODO want *different* M[1,] and M[1] for big[zq] matrices
*** TODO 3*M and M*3 are not the same: the first forgets the "nrow"
*** DONE x <- as.bigz(2:5); x[1] <- -1; x # -- a matrix !! >> no longer
*** DONE qq <- as.bigz(2:5)/7; qq[1] <- -1; qq # -- a matrix !!
*** also want *different* M[1,] and M[1] for big[zq] matrices
*** abs(.) or sign(<big..matrix>) {and probably many others!}
return vector instead of matrix, i.e.., lose the "nrow" attribute.
*** cbind(<bigz>, <bigq>) returns raw
......
\name{BernoulliQ}
\alias{BernoulliQ}
\title{Exact Bernoulli Numbers}
\description{%% ../R/Stirling-n-etc.R :
Return the \eqn{n}-th Bernoulli number \eqn{B_n}, (or \eqn{B_n^+}{Bn+},
see the reference), where \eqn{B_1 = + \frac 1 2}{B_1 = + 1/2}.
}
\usage{
BernoulliQ(n, verbose = getOption("verbose", FALSE))
}
\arguments{
\item{n}{integer \emph{vector}, \eqn{n \ge 0}{n >= 0}.}
\item{verbose}{logical indicating if computation should be traced.}
}
\value{
a big rational (class \code{\link[=bigq-class]{"bigq"}}) vector of the
Bernoulli numbers \eqn{B_n}.
}
\references{\url{https://en.wikipedia.org/wiki/Bernoulli_number}
}
\author{Martin Maechler}
\seealso{
\code{\link[Rmpfr]{Bernoulli}} in \CRANpkg{Rmpfr} in arbitrary precision
via Riemann's \eqn{\zeta}{zeta} function.
% \code{\link[DPQ]{Bern}(n)}
\code{Bern(n)} in \CRANpkg{DPQ} uses standard (double precision)
\R arithmetic for the n-th Bernoulli number.
}
\examples{
(Bn0.10 <- BernoulliQ(0:10))
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment