R devel - Apr 2025 - sqrt(.Machine$double.xmax)^2 == Inf, but only on Windows in R

>>>>> 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?


Best regards,
Martin

--
Martin Maechler
ETH Zurich  and  R Core team

On Tue, 29 Apr 2025 12:00:09 +0200
Martin Maechler <maechler at stat.math.ethz.ch> wrote:
> 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?
R-4.3.1 on Windows 7 in a virtual machine gives:

dput(sqrt(.Machine$double.xmax), control = 'hex')
# 0x1p+512
sqrt(.Machine$double.xmax)^2
# [1] Inf

...which differs from R on Linux:

dput(.Machine$double.xmax, control = 'hex')
# 0x1.fffffffffffffp+1023
dput(sqrt(.Machine$double.xmax), control = 'hex')
# 0x1.fffffffffffffp+511
dput(sqrt(.Machine$double.xmax)^2, control = 'hex')
# 0x1.ffffffffffffep+1023

-- 
Best regards,
Ivan

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 at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel