Power function numerical mystery

General Usage
1

I’m trying to hunt down a numerical platform deviation. Much simplification leads to this power computation, where I don’t understand the variation in the results (in the last bit of the mantissa).

julia> 15.632459f0 ^ 0.2f0
1.7330276f0

julia> 15.632459f0 .^ 0.2f0
1.7330276f0

julia> x, a = 15.632459f0, 0.2f0
(15.632459f0, 0.2f0)

julia> x ^ 0.2f0
1.7330276f0

julia> x .^ 0.2f0
1.7330275f0

julia> 15.632459f0 ^ a
1.7330276f0

julia> 15.632459f0 .^ a
1.7330275f0

julia> x ^ a
1.7330276f0

julia> x .^ a
1.7330276f0

julia> versioninfo()
Julia Version 0.6.2
Commit d386e40c17* (2017-12-13 18:08 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Xeon(R) CPU E31270 @ 3.40GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, sandybridge)

Can someone explain what’s going on here? Moreover, this variation does not happen on a Linux system with a more recent CPU:

julia> versioninfo()
Julia Version 0.6.2
Commit d386e40c17 (2017-12-13 18:08 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-7700K CPU @ 4.20GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Prescott)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, broadwell)
2

On 0.6, dot-ops fuse numeric literals into an anonymous function that gets compiled separately. I can’t reproduce, but maybe check the difference between:

f(x, a) = x ^ a
g(a) = 15.632459f0 ^ a

If f(15.632459f0, 0.2f0) and g(0.2f0) are different, then we can start digging into the LLVM/native code to see what’s different.

3

Cannot reproduce either. More precisely, I find no way to get 1.7330275f0.

4

Can you also check

h(x) = x^0.2f0
h(15.632459f0)

I suspect it is something weird with the system math library: though we use openlibm for most things, LLVM intrinsics (such as pow) still call the system library.

What OS and version of glibc are you running? You can find glibc version by ldd --version.

5

The odd results are from Windows 7 using the 64-bit binary from the download page. On Linux I get 1.7330276f0 in all cases.

6

Those all give the 1.7330275f0 result.

julia> f(x, a) = x ^ a
f (generic function with 1 method)

julia> g(a) = 15.632459f0 ^ a
g (generic function with 1 method)

julia> f(15.632459f0, 0.2f0)
1.7330275f0

julia> g(0.2f0)
1.7330275f0

julia> h(x) = x^0.2f0
h (generic function with 1 method)

julia> h(15.632459f0)
1.7330275f0
7

Windows 10 64bit, Julia 0.7 (all except second one giving 1.7330275f0):

julia> 15.632459f0 ^ 0.2f0
1.7330275f0

julia> 15.632459f0 .^ 0.2f0
1.7330276f0

julia> x, a = 15.632459f0, 0.2f0
(15.632459f0, 0.2f0)

julia> x ^ 0.2f0
1.7330275f0

julia> x .^ 0.2f0
1.7330275f0

julia> 15.632459f0 ^ a
1.7330275f0

julia> 15.632459f0 .^ a
1.7330275f0

julia> x ^ a
1.7330275f0

julia> x .^ a
1.7330275f0

Julia 0.6.2 though gives different results:

julia> 15.632459f0 ^ 0.2f0
1.7330276f0

julia> 15.632459f0 .^ 0.2f0
1.7330276f0

julia> x, a = 15.632459f0, 0.2f0
(15.632459f0, 0.2f0)

julia> x ^ 0.2f0
1.7330276f0

julia> x .^ 0.2f0
1.7330275f0

julia> 15.632459f0 ^ a
1.7330276f0

julia> 15.632459f0 .^ a
1.7330275f0

julia> x ^ a
1.7330276f0

julia> x .^ a
1.7330276f0