Tomas Kalibera
2025-Apr-29 13:54 UTC
[Rd] sqrt(.Machine$double.xmax)^2 == Inf, but only on Windows in R
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
Pavel Krivitsky
2025-May-28 01:36 UTC
[Rd] sqrt(.Machine$double.xmax)^2 == Inf, but only on Windows in R
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 at r-project.org?mailing list > > https://stat.ethz.ch/mailman/listinfo/r-devel