R Under development (unstable) (2023-08-05 r84874 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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. > stopifnot(require("Rmpfr")) Loading required package: 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 object is masked from 'package:gmp': outer The following objects are masked from 'package:stats': dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm The following objects are masked from 'package:base': cbind, pmax, pmin, rbind > > stopifnot(chooseMpfr(1:10, 0) == 1,# failed earlier + chooseMpfr(20, 0:20) == choose(20, 0:20), + chooseMpfr(19, 0:20) == choose(19, 0:20), + chooseMpfr (30, 4:30) * (-1)^(4:30) == + chooseMpfr.all(30, k0=4, alternating=TRUE) + ) > > cat('Time elapsed: ', proc.time(),'\n') # "stats" Time elapsed: 0.42 0.03 0.46 NA NA > > ## sumBinomMpfr() ... had embarrasing bug for a while > sBn <- Rmpfr:::sumBinomMpfr.v1 > stopifnot( + all.equal( sBn(10, sqrt), + sumBinomMpfr(10, sqrt), tol=1e-77) , + all.equal( sBn(10, log, n0=1, alternating=FALSE), + sumBinomMpfr(10, log, n0=1, alternating=FALSE), tol=1e-77) + ) > > fBin <- function(k) x^k * (1-x)^(n-k) > ## \sum_{k=0}^n (n \\ k) x^k (1-x)^{n-k} == sum(dbinom(0:n, n, prob=x)) == 1 : > for(x in runif(50)) { + n <- 1 + rpois(1, lambda=10) + cat(".") + stopifnot(all.equal(1, sumBinomMpfr(n, fBin, alternating=FALSE), + tol = 1e-15)) + };cat("\n") .................................................. > > > cat('Time elapsed: ', proc.time(),'\n') # "stats" Time elapsed: 0.75 0.04 0.81 NA NA > > if(!interactive()) warnings() > > proc.time() user system elapsed 0.75 0.04 0.81