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