Julia: floating point `div` differs in optimized code

[JeffBezanson]

I ran into this trying to get tests to pass on jb/functions. Also happens on master:

julia> div(0.3,0.01)
29.0

julia> f()=div(0.3,0.01)
f (generic function with 1 method)

julia> f()
30.0

[vchuravy]

That looks like constant folding to me:

julia> @code_llvm f()

define double @julia_f_21524() {
top:
  %0 = call double @llvm.rint.f64(double 3.000000e+01)
  ret double %0
}

compared to

julia> f(x, y) = x ÷ y
f (generic function with 2 methods)

julia> f(0.3, 0.01)
29.0

julia> @code_llvm f(0.3, 0.01)

define double @julia_f_21540(double, double) {
top:
  %2 = frem double %0, %1
  %3 = fsub double %0, %2
  %4 = fdiv double %3, %1
  %5 = call double @llvm.rint.f64(double %4)
  ret double %5
}

[yuyichao]

julia> rem(0.3, 0.01)
0.009999999999999983

julia> f() = rem(0.3, 0.01)
f (generic function with 1 method)

julia> f()
0.0

[yuyichao]

And seems like a LLVM bug (unless we are optimizing llvmcall ourselves....)

julia> function k()
           Base.llvmcall("""
                         %1 = frem double 3.000000e-01, 1.000000e-02
                         ret double %1
                         """, Float64, Tuple{})
       end
k (generic function with 1 method)

julia> @code_llvm k()

define double @julia_k_22866() #0 {
top:
  ret double 0.000000e+00
}

julia> k()
0.0

julia> function k2(a, b)
           Base.llvmcall("""
                         %3 = frem double %0, %1
                         ret double %3
                         """, Float64, Tuple{Float64,Float64}, a, b)
       end
k2 (generic function with 1 method)

julia> @code_llvm k2(0.3, 0.01)

define double @julia_k2_22867(double, double) #0 {
top:
  %2 = frem double %0, %1
  ret double %2
}

julia> k2(0.3, 0.01)
0.009999999999999983

[Keno]

Also seems to happen on llvm-svn.

[vchuravy]

Possibly related: https://llvm.org/bugs/show_bug.cgi?id=3316

[KristofferC]

Still happens.

[simonbyrne]

I've reopened the upstream issue with an example.

[simonbyrne]

Okay, I've put together a fix here: https://github.com/simonbyrne/llvm/commit/2718254b0627fbf07034b9f0e5ad245536bd8dfa

@Keno how do I go about actually getting this into LLVM?

[yuyichao]

http://llvm.org/docs/Phabricator.html#requesting-a-review-via-the-web-interface

[Keno]

Paste your diff here: https://reviews.llvm.org/differential/diff/create/ (Ideally with git dif -U99999 to get all the context), then fill out the form. It'll ask you to pick reviewers, but as long you put it up, I can find some people to review it for you.

[simonbyrne]

Submitted here: https://reviews.llvm.org/D29346

[Keno]

Found you some reviewers. Lets hope they know what they're doing.

[simonbyrne]

Now fixed upstream: https://github.com/llvm-mirror/llvm/commit/872b505b043539bb5d8dd6211266b651740f7e69

[tkelman]

Do we want to carry that patch, rebased for LLVM 3.9 if necessary?

Have our tests just been working around this?

[simonbyrne]

I think we've just been ignoring it, since it is a pretty weird edge case (calling rem with constant floating point arguments)

[jayvn]

When will the fix to llvm be out then ?
I ran into

julia> div(1,0.2)
4.0

Should be 5.
Which is why I'm here. As a workaround I tried
Int(a/b), but Int(2.9999998) throws InexactError().

julia> Int(2.99999999)
ERROR: InexactError()
Stacktrace:
 [1] convert(::Type{Int64}, ::Float64) at ./float.jl:679
 [2] Int64(::Float64) at ./sysimg.jl:77

julia> versioninfo()
Julia Version 0.6.1
Commit 0d7248e2ff (2017-10-24 22:15 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-7500 CPU @ 3.40GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Prescott)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, broadwell)

Will the Int(2.999998) bug be fixed as well with the patch ?

[GunnarFarneback]

0.2 can't be represented exactly as a binary floating point value, so it's reasonable for div(1, 0.2) to be either 4.0 or 5.0 depending on whether 0.2 turns out to represented by a value that is slightly smaller or slightly larger than 0.2. That's not what this issue is about though, but whether optimizations cause inconsistent results for the same computation.

2.999998 is not an integer so Int(2.999998) throwing InexactError is working as intended. What you probably want is round(a / b) or round(Int, a / b).

Further discussion about the nature of floating point arithmetics should go on the Discourse forum since it's off-topic for this particular issue.

[StefanKarpinski]

0.2 can't be represented exactly as a binary floating point value, so it's reasonable for div(1, 0.2) to be either 4.0 or 5.0 depending on whether 0.2 turns out to represented by a value that is slightly smaller or slightly larger than 0.2.

This is not random, however: there is only one Float64 value that 0.2 can correctly represent and it is slightly larger than 1/5. However, five times that number is slightly larger than 1.

[simonbyrne]

Also, fix is now in master, so can be closed.

[tkelman]

Master of Julia, or only master of llvm? Is it tested, if fixed here?

[simonbyrne]

Ah, I forgot I had compiled with LLVM-svn.

[KristofferC]

Fixed in release now.