[Rd] Floating point issue

[Martin Maechler]

>>>>> Brodie Gaslam via R-devel 
>>>>>     on Sun, 10 Jul 2022 21:10:29 -0400 writes:

    > The results are not exactly the same.  Notice that on Bill's system the 
    > bit pattern of 10^25 and 10000000000000000905969664 are the same, but 
    > not so on yours.  So there is a mismatch happening on parsing between 
    > your M1 mac and other's systems.

    > This is the main thing I wanted to point out.  But since I'm here I'm 
    > going to add some additional lazily researched speculation.

    > As Dirk points out, M1 does not have long double, and if you look at 
    > what I think is responsible for parsing of numbers like the ones we're 
    > discussing here, we see[1]:

    > double R_strtod5(const char *str, char **endptr, char dec,
    >                  Rboolean NA, int exact)
    > {
    >      LDOUBLE ans = 0.0;

    > IIRC long double on systems that implement it as 80 bits (most that have 
    > x87 coprocessors) has 63-4 bits of precision, vs 53 for 64 bit long 
    > double.  Roughly speaking, that's 19-20 digits of base 10 precision for 
    > long double, vs 15-16 for 64 bit double.  Then:

    >> substr(rep("10000000000000000905969664", 2),  c(1, 1), c(16, 20))
    > [1] "1000000000000000"     "10000000000000000905"

    > Again, I have not carefully researched this, but it seems likely that 
    > parsing is producing different a different outcome in this case because 
    > the intermediate values generated can be kept at higher precision on 
    > systems with 80 bit long doubles prior to coercing to double for the 
    > final result.

    > IIRC, if you need invertible deparsing/parsing I think you can use:

    > deparse(1e25, control=c('hexNumeric'))

    > Although I don't have an 80-bit-less system to test on (and I am too 
    > lazy to recompile R without long double to test).

    > Best,
    > B.

    > [1]: https://github.com/r-devel/r-svn/blob/master/src/main/util.c#L1993

An excellent analysis, thank you, Brodie and Bill Dunlap.

10 days ago I also had a colleague at ETH, involved in teaching stats,
wanting to submit a bug report about a model fit in nlme which
gave different results (well, a confidence interval of a
correlation which was not at all significant (and estimated
close to zero) changed from +/-0.45 to +/-0.85 ...
... and when I asked for sessionInfo() etc,
the result was that on the  M1  the result differed from all
other platforms.
I did not go into further analysis, but I know that nlme does
*not* use any LDOUBLE | long double  but still, the M1 was less
accurate than any other platform... and I also could not see this
(the wider conf.ints) when I used a version of R configured with

   configure --disable-long-double  

[If anyone is interested I can provide the repr.ex. R code.]

>From my current experiences, I dare to say that the M1 with all
its speed is just a tad less reliable numerically than the
Intel/AMD floating point implementations..

Martin




    > On 7/10/22 5:38 PM, Antoine Fabri wrote:
    >> Thanks, I get the exact same results as yours
    >> 
    >> ``` r
    >> bitC <- function(x) noquote(vapply(as.double(x), function(x) { # split one
    >> double
    >> b <- substr(as.character(rev(numToBits(x))), 2L, 2L)
    >> paste0(c(b[1L], " ", b[2:12], " | ", b[13:64]), collapse = "")
    >> }, ""))
    >> bitC(10^25)
    >> #> [1] 0 10001010010 | 0000100010110010101000101100001010000000001010010001
    >> bitC(10000000000000000905969664)
    >> #> [1] 0 10001010010 | 0000100010110010101000101100001010000000001010010010
    >> bitC(10000000000000000905969664 - 10^25)
    >> #> [1] 0 10000011110 | 0000000000000000000000000000000000000000000000000000
    >> bitC(1e25)
    >> #> [1] 0 10001010010 | 0000100010110010101000101100001010000000001010010001
    >> ```
    >> 
    >> <sup>Created on 2022-07-10 by the [reprex package](
    >> https://reprex.tidyverse.org) (v2.0.1)</sup>
    >> 
    >> Le dim. 10 juil. 2022 à 22:23, Bill Dunlap <williamwdunlap using gmail.com> a
    >> écrit :
    >> 
    >>> The following function, 'bitC' from ?numToBits, displays the bits in a
    >>> double precision number, separated into the sign bit, the 11 exponent bits,
    >>> and the 52 bits in the mantissa.  I've shown the results with your numbers
    >>> from R-2.4.0 on my Windows 11 Lenovo laptop: what do you get?
    >>> 
    >>>> bitC <- function(x) noquote(vapply(as.double(x), function(x) { # split
    >>> one double
    >>> +     b <- substr(as.character(rev(numToBits(x))), 2L, 2L)
    >>> +     paste0(c(b[1L], " ", b[2:12], " | ", b[13:64]), collapse = "")
    >>> +   }, ""))
    >>>> bitC(10^25)
    >>> # [1] 0 10001010010 | 0000100010110010101000101100001010000000001010010001
    >>>> bitC(10000000000000000905969664)
    >>> # [1] 0 10001010010 | 0000100010110010101000101100001010000000001010010001
    >>>> bitC(10000000000000000905969664 - 10^25)
    >>> # [1] 0 00000000000 | 0000000000000000000000000000000000000000000000000000
    >>>> bitC(1e25)
    >>> # [1] 0 10001010010 | 0000100010110010101000101100001010000000001010010001
    >>> 
    >>> -Bill
    >>> 
    >>> On Sun, Jul 10, 2022 at 7:00 AM Antoine Fabri <antoine.fabri using gmail.com>
    >>> wrote:
    >>> 
    >>>> Dear r-devel,
    >>>> 
    >>>> For some numbers, the printed value is not equivalent to the input :
    >>>> 
    >>>> options(scipen = 999)
    >>>> ## GOOD
    >>>> 1e24
    >>>> #> [1]  999999999999999983222784
    >>>> 1e24 == 999999999999999983222784
    >>>> #> [1] TRUE
    >>>> 
    >>>> ## BAD
    >>>> 1e25
    >>>> #> [1] 10000000000000000905969664
    >>>> 1e25 == 10000000000000000905969664
    >>>> #> [1] FALSE
    >>>> 
    >>>> ## STILL BAD
    >>>> 10000000000000000905969664
    >>>> #> [1] 10000000000000003053453312
    >>>> 
    >>>> ## GOOD AGAIN
    >>>> 10000000000000003053453312
    >>>> #> [1] 10000000000000003053453312
    >>>> 
    >>>> # Additionally
    >>>> 10000000000000000000000000 == 1e25
    >>>> #> [1] FALSE
    >>>> 
    >>>> Are these bugs ?
    >>>> 
    >>>> [[alternative HTML version deleted]]
    >>>> 
    >>>> ______________________________________________
    >>>> R-devel using r-project.org mailing list
    >>>> https://stat.ethz.ch/mailman/listinfo/r-devel
    >>>> 
    >>> 
    >> 
    >> [[alternative HTML version deleted]]
    >> 
    >> ______________________________________________
    >> R-devel using r-project.org mailing list
    >> https://stat.ethz.ch/mailman/listinfo/r-devel

    > ______________________________________________
    > R-devel using r-project.org mailing list
    > https://stat.ethz.ch/mailman/listinfo/r-devel