[Rd] sqrt(.Machine$double.xmax)^2 == Inf, but only on Windows in R

[Martin Maechler]

>>>>> Pavel Krivitsky 
>>>>>     on Wed, 28 May 2025 01:36:57 +0000 writes:

    > Dear All,
    > Thanks for looking into it, and apologies for bumping this.

    > I think the strangest thing for me is that the C and the R
    > implementations on Windows yield different results. I don't know if it
    > warrants a deeper look.

This is not so strange, really:
In many cases (and some similar to this one), we *did* want R to
be platform independent and give the same results independent of
platform. ... consequently differing on some platforms from
their platform dependent C libraries ..

Martin

    > Ultimately, I rewrote the code that relied on this behaviour. (I needed
    > something that was a finite number when squared but was still as big as
    > possible, so I divided it 1 + Machine epsilon, and it seems to work for
    > my particular situation.)

    > Best,
    > Pavel

    > On Tue, 2025-04-29 at 15:54 +0200, Tomas Kalibera wrote:
    >> 
    >> On 4/29/25 12:00, Martin Maechler wrote:
    >> > > > > > > Pavel Krivitsky via R-devel
    >> > > > > > >      on Mon, 28 Apr 2025 05:13:41 +0000 writes:
    >> >      > Hello, Under R 4.5.0 on Windows (x86-64), I get:
    >> > 
    >> >      >> sqrt(.Machine$double.xmax)^2
    >> >      > [1] Inf
    >> >      >> sqrt(.Machine$double.xmax)*sqrt(.Machine$double.xmax)
    >> >      > [1] Inf
    >> > 
    >> >      > On other hand on other platforms, including Debian Linux
    >> >      > (x86-64), I get:
    >> > 
    >> >      d> sqrt(.Machine$double.xmax)^2
    >> >      > [1] 1.797693134862315508561e+308
    >> >      d> sqrt(.Machine$double.xmax)*sqrt(.Machine$double.xmax)
    >> >      > [1] 1.797693134862315508561e+308
    >> > 
    >> >      > Windows is running inside a VirtualBox instance on the
    >> >      > same host as Linux. I don't have direct results from the
    >> >      > MacOS platforms, but based on the symptoms that had led me
    >> >      > to investigate, the behaviour is as the Linux.
    >> > 
    >> >      > Adding to the mystery, if I implement the same operation in
    >> > C, e.g.,
    >> > 
    >> >      > library(inline)
    >> >      > sqrsqrt <- cfunction(sig = c(x = "numeric"), language = "C",
    >> > "
    >> >      > double sqrtx = sqrt(Rf_asReal(x));
    >> >      > return Rf_ScalarReal(sqrtx*sqrtx);
    >> >      > ")
    >> > 
    >> >      > R on Linux yields:
    >> > 
    >> >      d> sqrsqrt(.Machine$double.xmax)
    >> >      > [1] 1.797693134862315508561e+308
    >> > 
    >> >      > i.e., the same number, whereas R on Windows yields:
    >> > 
    >> >      d> sqrsqrt(.Machine$double.xmax)
    >> >      > [1] 1.797693134862315508804e+308
    >> > 
    >> >      > which is not Inf but not the same as Linux either.
    >> > 
    >> >      > Lastly, on both platforms,
    >> > 
    >> >      d> sqrsqrt(.Machine$double.xmax) < .Machine$double.xmax
    >> >      > [1] TRUE
    >> > 
    >> >      > I am not sure if this is a bug, intended behaviour, or
    >> > something else.
    >> > 
    >> > "something else":  It is not a bug, nor intended, but should
    >> > also *not* be surprising nor a mistery:  The largest possible
    >> > double precision number is by definition "very close to
    >> > infinity" (in the space of double precision numbers)
    >> > [R 4.5.0 patched on Linux (Fedora 40; x86_64)] :
    >> > 
    >> > > (M <- .Machine$double.xmax)
    >> > [1] 1.797693e+308
    >> > > M+1 == M
    >> > [1] TRUE
    >> > > M*(1 + 2^-52)
    >> > [1] Inf
    >> > > print(1 + 2^-52, digits= 16)
    >> > [1] 1
    >> > > print(1 + 2^-52, digits= 17)
    >> > [1] 1.0000000000000002
    >> > What you see, I'd classify as quite related to R FAQ 7.31,
    >> >   := "the most frequently asked question about R
    >> >       among all the other frequently asked questions"
    >> > 
    >> > A nice reading is the README you get here
    >> >    https://github.com/ThinkR-open/seven31
    >> > which does link also to the R FAQ at
    >> >  
    >> > https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f
    >> > 
    >> > Of tangential interest only:
    >> > You mention that it is R 4.5.0 you use on Windows.
    >> > Would you (or anybody else) know if this is new behaviour or it
    >> > also happened e.g. in R 4.4.x versions on  Windows?
    >> 
    >> This depends on C math function sqrt. With sqrt implemented in mingw-
    >> w64 
    >> (mingwex), which is still the case of mingw-w64 v11, so still of 
    >> Rtools45, the result of sqrt(.Machine$double.xmax) is
    >> 
    >> "0x1p+512"
    >> 
    >> with mingw-w64 v12 (so with UCRT, and also on Linux, and also using
    >> mpfr 
    >> library) it is
    >> 
    >> "0x1.fffffffffffffp+511"
    >> 
    >> It is not required by any standard for the math functions in the C
    >> math 
    >> library to be precise (correctly rounded). But still, in this case,
    >> the 
    >> correctly rounded value would be returned once Rtools can switch to 
    >> mingw-w64 v12 (which could be no sooner than Rtools46, as mingw-w64
    >> is a 
    >> core component of the toolchain).
    >> 
    >> A bit of advertising -  I used Martin's Rmpfr package to compute this
    >> using mpfr.
    >> And there is a related blog post: 
    >> https://blog.r-project.org/2025/04/24/sensitivity-to-c-math-library-and-mingw-w64-v12
    >> 
    >> Best
    >> Tomas
    >> 
    >> 
    >> > 
    >> > 
    >> > Best regards,
    >> > Martin
    >> > 
    >> > --
    >> > Martin Maechler
    >> > ETH Zurich  and  R Core team
    >> > 
    >> > ______________________________________________
    >> > R-devel using r-project.org mailing list
    >> > https://stat.ethz.ch/mailman/listinfo/r-devel