R Under development (unstable) (2024-11-15 r87338 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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

> 
> n <- 10000 + 0:10
> stopifnot(exprs = {
+     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)
+     chooseMpfr(mpfr(1111, 2^8), 1111) == 1
+     chooseMpfr(mpfr(n,256), n  ) == 1 # was wrong in <= 2023
+     chooseMpfr(mpfr(n,256), n-1) == n #  "    "
+ })
> cat('Time elapsed: ', proc.time(),'\n') # "stats"
Time elapsed:  0.68 0.14 0.76 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:  1.12 0.15 1.21 NA NA 
> 
> if(!interactive()) warnings()
> 
> proc.time()
   user  system elapsed 
   1.12    0.15    1.21