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 ζ(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

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

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

See Also

zeta is used to compute them.

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

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Rmpfr)
Loading required package: gmp

Attaching package: 'gmp'

The following objects are masked from 'package:base':

    %*%, apply, crossprod, matrix, tcrossprod

C code of R package 'Rmpfr': GMP using 64 bits per limb


Attaching package: 'Rmpfr'

The following objects are masked from 'package:stats':

    dbinom, dnorm, dpois, pnorm

The following objects are masked from 'package:base':

    cbind, pmax, pmin, rbind

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Rmpfr/Bernoulli.Rd_%03d_medium.png", width=480, height=480)
> ### Name: Bernoulli
> ### Title: Bernoulli Numbers in Arbitrary Precision
> ### Aliases: Bernoulli
> ### Keywords: arith
> 
> ### ** Examples
> 
> ## Don't show: 
> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04 LTS

locale:
[1] C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] Rmpfr_0.6-0 gmp_0.5-12 
>  .libPaths()
[1] "/home/ddbj/local/lib64/R/library"
>  packageDescription("gmp")
Package: gmp
Version: 0.5-12
Date: 2014-07-28
Title: Multiple Precision Arithmetic
Author: Antoine Lucas, Immanuel Scholz, Rainer Boehme
        <rb-gmp@reflex-studio.de>, Sylvain Jasson
        <jasson@toulouse.inra.fr>, Martin Maechler
        <maechler@stat.math.ethz.ch>
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).
Imports: methods
Suggests: Rmpfr
SystemRequirements: gmp (>= 4.2.3)
License: GPL
BuildResaveData: no
LazyDataNote: not available, as we use data/*.R *and* our classes
URL: http://mulcyber.toulouse.inra.fr/projects/gmp
Packaged: 2014-07-28 19:35:49 UTC; antoine
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2014-07-28 21:53:29
Built: R 3.3.1; x86_64-pc-linux-gnu; 2016-07-01 22:48:44 UTC; unix

-- File: /home/ddbj/local/lib64/R/library/gmp/Meta/package.rds 
> ## End(Don't show)
> Bernoulli(0:10)
11 'mpfr' numbers of precision  128   bits 
 [1]                                            1
 [2]                                          0.5
 [3]   0.1666666666666666666666666666666666666669
 [4]                                           -0
 [5] -0.03333333333333333333333333333333333333342
 [6]                                           -0
 [7]  0.02380952380952380952380952380952380952381
 [8]                                           -0
 [9] -0.03333333333333333333333333333333333333342
[10]                                           -0
[11]   0.0757575757575757575757575757575757575757
> 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
1 'mpfr' number of precision  128   bits 
[1] -9.049423963609480500529241443083535605319e27677
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>