R Under development (unstable) (2024-01-20 r85814 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(suppressPackageStartupMessages(require("Rmpfr")))
> ## (checking that the 32 / 64 bit  GMP message does *not* show here)
> 
> ### Try to look at the internal bit-representation of the limbs
> 
> .limbs <- function(x) {
+     stopifnot(is(x, "mpfr"))
+     lapply(x@.Data, slot, "d") # not sapply() each can have different prec. & #{limbs}
+ }
> .expo <- function(x) {
+     stopifnot(is(x, "mpfr"))
+     sapply(x@.Data, slot, "exp")
+ }
> 
> Bits <- function(x) {
+     L <- .limbs(x)# list(length n) each of "k(prec)" 32-bit ints
+     ## NB:  mpfr(2, .) and all mpfr(2^k, .) also have a 'd' ending in NA integer!
+     ##     [reason: after all, R's NA_integer_ is INT_MAX+1 = 2^31 ]
+     ## and  the  mpfr(c(NA,NaN, Inf, -Inf), .)   have *no* NA in 'd' (but all in 'exp'!
+     ## see .mpfr2list() example below
+ 
+     hasNA <- any(iNA <- sapply(lapply(L, is.na), any)) # iNA: TRUE if there's an NA
+     ## need to catch them later
+     CC <- function(ch) paste(ch, collapse="")
+     hex <- sapply(L, function(.) CC(sprintf("%x", rev(.))))
+     if(hasNA) hex[iNA] <- NA_character_
+     hex <- strsplit(hex, NULL)
+ 
+     db <- t(expand.grid(0:1,0:1,0:1,0:1, KEEP.OUT.ATTRS=FALSE)[,4:1])
+     storage.mode(db) <- "character" # "0" or "1"
+     dimnames(db) <- list(NULL, c(paste(0:9), letters[1:6]))
+     ## db is  4 x 16  matrix  with col.names "0" "1" .. "9" "a" "b" ... "f"
+ 
+     ex <- .expo(x)
+     if(is.matrix(ex)) {
+         ## 64-bit case: exponent is long == two ints
+         ## -----------  the 2nd int is in {0, -1, NA} (NA : for 0)
+         ex2 <- ex[2,]
+         ex <- ex[1,]
+     }
+     pat <- paste("(", sapply(pmax(0, ex),
+                              function(n) CC(rep.int(".", n))),
+                  ")0+$", sep="")
+     ## pat <- ifelse(iNA, NA_character_, pat)
+ 
+     getbits <- function(ch) CC(as.vector(db[,ch]))
+ 
+     ## drop trailing zeros (from r[[i]], via pat[i]), keeping correct number:
+     drop0.r <- function(i) sub(pat[i], "\\1", r[[i]])
+ 
+     if(hasNA) {
+         r <- as.list(iNA)
+         r[!iNA] <- lapply(hex[!iNA], getbits)
+         r[!iNA] <- lapply(which(!iNA), drop0.r)
+         ## FIXME this is wrong -- really have powers of 2, and want their (easy) bits :
+         r[iNA ] <- NA_character_
+         unlist(r)
+     }
+     else {
+ 	r <- lapply(hex, getbits)
+         sapply(seq_along(r), drop0.r)
+     }
+ 
+ }
> 
> x <- mpfr(r <- c(NA,NaN, Inf, -Inf), 64)
> stopifnot(identical(asNumeric(x), # mpfr has no NA, just NaN's:
+ 		    c(NaN,NaN, Inf, -Inf)),
+ 	  identical(as.character(fDec <- formatDec(x)),
+ 		    as.character(asNumeric(x))) # of different nchar() for now
+ 	  )
> formatDec(x) # should print fine (no quotes)
[1] NaN  NaN  Inf  -Inf
> 
> 
> if(FALSE) # platform dependent:
+     ## The "non-finite" mpfr value internals (in 64-bit: 'exp' has NA):
+     str(.mpfr2list(x))
> 
> 
> ## bug in Bits(): all (exact) powers of 2  will show as NA:
> 
> x <- mpfr(c(3:5,11:16, 59, 125:128, 1024:1025), 64)
> x
16 'mpfr' numbers of precision  64   bits 
 [1]    3    4    5   11   12   13   14   15   16   59  125  126  127  128 1024
[16] 1025
> data.frame(x= as.numeric(x), I(Bits(x)))
      x     Bits.x.
1     3          11
2     4        <NA>
3     5         101
4    11        1011
5    12        1100
6    13        1101
7    14        1110
8    15        1111
9    16        <NA>
10   59      111011
11  125     1111101
12  126     1111110
13  127     1111111
14  128        <NA>
15 1024        <NA>
16 1025 10000000001
> 
> x <- mpfr(c(-20:-1, 1:30),64)# w/o  0  - as its mantissa is "random" (in 64-bit)
> data.frame(x= as.numeric(x), I(Bits(x)))
     x Bits.x.
1  -20   10100
2  -19   10011
3  -18   10010
4  -17   10001
5  -16    <NA>
6  -15    1111
7  -14    1110
8  -13    1101
9  -12    1100
10 -11    1011
11 -10    1010
12  -9    1001
13  -8    <NA>
14  -7     111
15  -6     110
16  -5     101
17  -4    <NA>
18  -3      11
19  -2    <NA>
20  -1    <NA>
21   1    <NA>
22   2    <NA>
23   3      11
24   4    <NA>
25   5     101
26   6     110
27   7     111
28   8    <NA>
29   9    1001
30  10    1010
31  11    1011
32  12    1100
33  13    1101
34  14    1110
35  15    1111
36  16    <NA>
37  17   10001
38  18   10010
39  19   10011
40  20   10100
41  21   10101
42  22   10110
43  23   10111
44  24   11000
45  25   11001
46  26   11010
47  27   11011
48  28   11100
49  29   11101
50  30   11110
> b0 <- Bits(mpfr(0, 64)) # not printing it here -- they are "random" for this special case!
> 
> (half <- mpfr(0.5, 64)*(1 + mpfr(2, 64)^-16 * (-3:3)))
7 'mpfr' numbers of precision  64   bits 
[1] 0.49997711181640625  0.4999847412109375 0.49999237060546875
[4]                 0.5 0.50000762939453125  0.5000152587890625
[7] 0.50002288818359375
> Bits(half)
[1] "1111111111111101"  "111111111111111"   "1111111111111111" 
[4] NA                  "10000000000000001" "1000000000000001" 
[7] "10000000000000011"
> 
> ## pi, in varying number of bits :
> p. <- round(pi* 2^c(10,16,5*(4:8)))
> dput(p.)#-> the definition of p :
c(3217, 205887, 3294199, 105414357, 3373259426, 107944301636, 
3454217652358)
> p <- mpfr(c(3217, 205887, 3294199, 105414357,
+             3373259426, 107944301636, 3454217652358), 64)
> stopifnot(all.equal(p., p, tolerance = 1e-15))
> ## all the mantissas are those of pi, rounded differently:
> Bits(c(p, Const("pi", 64)))
[1] "110010010001"                                                    
[2] "110010010000111111"                                              
[3] "1100100100001111110111"                                          
[4] "110010010000111111011010101"                                     
[5] "11001001000011111101101010100010"                                
[6] "1100100100001111110110101010001000100"                           
[7] "110010010000111111011010101000100010000110"                      
[8] "1100100100001111110110101010001000100001011010001100001000110101"
> 
> ###--- and possibly the _internal_   sprintfMpfr() ---  see also ./tstHexBin.R
> ## TODO: use examples above for checking formatBin() <--->        ============
> spr <- Rmpfr:::sprintfMpfr
> ##=            ~~~~~~~~~~~
> (fB.04 <- formatBin(i16.04 <- mpfr(0:16,  4)))
 [1] +0b0.000p+0 +0b1.000p+0 +0b1.000p+1 +0b1.100p+1 +0b1.000p+2 +0b1.010p+2
 [7] +0b1.100p+2 +0b1.110p+2 +0b1.000p+3 +0b1.001p+3 +0b1.010p+3 +0b1.011p+3
[13] +0b1.100p+3 +0b1.101p+3 +0b1.110p+3 +0b1.111p+3 +0b1.000p+4
> (fB.60 <- formatBin(i16.60 <- mpfr(0:16, 60)))
 [1] +0b0.00000000000000000000000000000000000000000000000000000000000p+0
 [2] +0b1.00000000000000000000000000000000000000000000000000000000000p+0
 [3] +0b1.00000000000000000000000000000000000000000000000000000000000p+1
 [4] +0b1.10000000000000000000000000000000000000000000000000000000000p+1
 [5] +0b1.00000000000000000000000000000000000000000000000000000000000p+2
 [6] +0b1.01000000000000000000000000000000000000000000000000000000000p+2
 [7] +0b1.10000000000000000000000000000000000000000000000000000000000p+2
 [8] +0b1.11000000000000000000000000000000000000000000000000000000000p+2
 [9] +0b1.00000000000000000000000000000000000000000000000000000000000p+3
[10] +0b1.00100000000000000000000000000000000000000000000000000000000p+3
[11] +0b1.01000000000000000000000000000000000000000000000000000000000p+3
[12] +0b1.01100000000000000000000000000000000000000000000000000000000p+3
[13] +0b1.10000000000000000000000000000000000000000000000000000000000p+3
[14] +0b1.10100000000000000000000000000000000000000000000000000000000p+3
[15] +0b1.11000000000000000000000000000000000000000000000000000000000p+3
[16] +0b1.11100000000000000000000000000000000000000000000000000000000p+3
[17] +0b1.00000000000000000000000000000000000000000000000000000000000p+4
> stopifnot(
+     identical(sub("00p","p", spr(i16.60, bits = 10)),
+                              spr(i16.60, bits = 4)),
+     identical(spr(i16.60, bits = 4),
+               spr(i16.04, bits = 4))
+     ,
+     all.equal(i16.04, mpfr(fB.04), tolerance = 0)
+     ,
+     all.equal(i16.60, mpfr(fB.60), tolerance = 0)
+ )
> 
> ## not even this one
> two <- mpfr(2, precBits = 60)
> stopifnot(identical(two, mpfr(formatBin(two))))
> 
> 
> cat('Time elapsed: ', proc.time(),'\n') # "stats"
Time elapsed:  0.62 0.07 0.68 NA NA 
> 
> if(!interactive()) warnings()
> 
> proc.time()
   user  system elapsed 
   0.62    0.07    0.68