Search code examples
rfloating-accuracygmp

conversion bigq to mpfr with Rmpfr package

The help documentation of the Rmpfr R package claims that the .bigq2mpfr() function uses the minimal precision necessary for correct representation when the precB argument is NULL :

Description:

     Coerce from and to big integers (‘bigz’) and ‘mpfr’ numbers.

     Further, coerce from big rationals (‘bigq’) to ‘mpfr’ numbers.

Usage:

     .bigz2mpfr(x, precB = NULL)
     .bigq2mpfr(x, precB = NULL)
     .mpfr2bigz(x, mod = NA)

Arguments:

       x: an R object of class ‘bigz’, ‘bigq’ or ‘mpfr’ respectively.

   precB: precision in bits for the result.  The default, ‘NULL’, means
          to use the _minimal_ precision necessary for correct
          representation.

However when converting 31/3 one gets a bad approximation:

> x <- as.bigq(31,3)
> .bigq2mpfr(x)
1 'mpfr' number of precision  8   bits 
[1] 10.31

By looking inside the .bigq2mpfr() function we see the detailed procedure:

N <- numerator(x)
D <- denominator(x)
if (is.null(precB)) {
    eN <- frexpZ(N)$exp
    eD <- frexpZ(D)$exp
    precB <- eN + eD + 1L
}
.bigz2mpfr(N, precB)/.bigz2mpfr(D, precB)

Firstly I do not understand why precB is taken as follows. The exp output of the frexpZ() is the exponent in binary decomposition:

> frexpZ(N)
$d
[1] 0.96875

$exp
[1] 5

> 0.96875*2^5
[1] 31

Here we get precB=8 and the result is then identical to:

> mpfr(31, precBits=8)/mpfr(3, precBits=8)
1 'mpfr' number of precision  8   bits 
[1] 10.31

I am under the impression one should rather replace precB with 2^precB but I'd like to get some advices about that:

> mpfr(31, precBits=8)/mpfr(3, precBits=2^8)
1 'mpfr' number of precision  256   bits 
[1] 10.33333333333333333333333333333333333333333333333333333333333333333333333333329
> mpfr(31, precBits=8)/mpfr(3, precBits=2^9)
1 'mpfr' number of precision  512   bits 
[1] 10.3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333329
> mpfr(31, precBits=8)/mpfr(3, precBits=2^7)
1 'mpfr' number of precision  128   bits 
[1] 10.33333333333333333333333333333333333332
Solution

  • This has been corrected in a newer version of the package:

    > x <- as.bigq(31,3)
> .bigq2mpfr(x)
1 'mpfr' number of precision  128   bits 
[1] 10.33333333333333333333333333333333333332