Gamma function overflows on Windows, but not on WSL2/Ubuntu

[DNF2]

The gamma function behaves differently on ordinary Windows 10 and under WSL2/Ubuntu.

These are both on the same computer:

WLS2:

jl> using SpecialFunctions
[ Info: Precompiling SpecialFunctions [276daf66-3868-5448-9aa4-cd146d93841b]

jl> gamma(big(3.6e13))
7.882861104028365451673058070377064242970022825665897715265241972922672893761775e+472392288679098

Windows 10:

jl> using SpecialFunctions
[ Info: Precompiling SpecialFunctions [276daf66-3868-5448-9aa4-cd146d93841b]

jl> gamma(big(3.6e13))
Inf

Version info:
WLS2:

jl> versioninfo()
Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-10875H CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)

(@temp) pkg> st
      Status `~/.julia/environments/temp/Project.toml`
  [276daf66] SpecialFunctions v1.3.0

Windows 10:

jl> versioninfo()
Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-10875H CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)

(@v1.6) pkg> st
      Status `C:\Users\dnf\.julia\environments\v1.6\Project.toml`
  [6e4b80f9] BenchmarkTools v0.7.0
  [052768ef] CUDA v3.0.3
  [e9467ef8] GLMakie v0.2.6
  [bdcacae8] LoopVectorization v0.12.12
  [5fb14364] OhMyREPL v0.5.10
  [295af30f] Revise v3.1.15
  [276daf66] SpecialFunctions v1.3.0

[stevengj]

Note that this method is simply calling mpfr_gamma, so apparently something is wrong with the MPFR library compiled for Windows?

In particular, try calling MPFR directly:

z = BigFloat()
ccall((:mpfr_gamma, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, big(3.6e13), Base.MPFR.ROUNDING_MODE[])
z

If that still has a problem, then you can file a Julia issue, since we can't fix MPFR in the SpecialFunctions package.

[simonbyrne]

You truncated the output, the actual result is:

julia> x = gamma(big(3.6e13))
7.882861104028365451673058070377064242970022825665897715265241972922672893761775e+472392288679098

In particular, if you look at the exponent, it is

julia> x.exp
1569253215571240

The problem is that the exponent is stored as a C long: on Linux (and WSL2) it is an Int64, but on Windows this is an Int32, which is unable to store a value this large, hence it overflows to an Inf. Unfortunately there isn't much we can do about this unless MPFR switches to using long long.

[DNF2]

Indeed, I can confirm that calling MPFR directly yields the same behaviour: overflow under Windows and not under Linux.