div/fld floating point issue

[simonbyrne]

Here's another, in the spirit of #3127 and #3321:

julia> t = -log1p(-0.3)
0.35667494393873234

julia> 1.4/t
3.925142552879981

julia> fld(1.4,t)
2.0

julia> div(1.4,t)
2.0

I would think that changing the trunc and floor to round would fix this?

[ViralBShah]

Is this something we need to address by 0.2?

[GunnarFarneback]

It's at least rather embarrassing. Here's another somewhat simpler and slightly worse example:

julia> 0.3/0.01
30.0

julia> fld(0.3,0.01)
28.0

What happens is that

div{T<:Real}(x::T, y::T) = convert(T,trunc((x-rem(x,y))/y))

expands to

julia> rem(0.3,0.01)
0.009999999999999983

julia> 0.3-rem(0.3,0.01)
0.29

julia> (0.3-rem(0.3,0.01))/0.01
28.999999999999996

julia> trunc((0.3-rem(0.3,0.01))/0.01)
28.0

The first two lines cast some doubt about the validity of the fix of #3127, but that's nothing to try to solve with quick fixes. From the last two lines, however, I agree with Simon that trunc should be round. In retrospect it looks obvious that the fix of #3321 should have used round in the first place.

[Yaffle]

Hello!
I am not the user of Julia language, but I searched for correct implementation of "floor division".

It seems, your implementation is not correct for:
div(18014398509481992., 7.) and div(18014398509481988., 3.).
Can you please refute this?

Thanks

[simonbyrne]

@Yaffle The accuracy of div(a,b) for floating point arguments depends on the ratio of abs(a) / abs(b). I have a proof, which I have been meaning to write up, that our div implementation is exact for Float64 when the ratio is less than 2^51. Once you get above 2^53, all bets are off, as you can't represent all possible integers anyway, but there certainly are some cases in between where we are wrong: both your examples fall in this range, so I suspect that's what is happening.

Comparing your examples with integer division (which is exact), using Julia 0.3,

julia> int(div(18014398509481992., 7.))
2573485501354571

julia> int(div(18014398509481992, 7))
2573485501354570

julia> int(div(18014398509481988., 3.))
6004799503160661

julia> int(div(18014398509481988, 3))
6004799503160662

which are both wrong. In Julia 0.4, the first one is actually correct, but this is just due to the change in round (#8750, #9344), and there are probably cases which have gone the other way.

[Yaffle]

@simonbyrne Thanks for your quick reply!

I am not an expert in floating point, but I have my own candidate (the code is in JavaScript):

var truncDiv = function (x, y) {
  var quotient = x / y;
  var remainder = x % y;
  var truncatedQuotient = Math.trunc(quotient);
  if (quotient === truncatedQuotient) {
    if (Math.abs(remainder) > Math.abs(y) - Math.abs(remainder) && Math.abs(truncatedQuotient) <= 2 / Number.EPSILON) {
      truncatedQuotient -= Math.sign(truncatedQuotient);
    }
  }
  return truncatedQuotient;
};

What do you think?

[simonbyrne]

Hmm, something like that might work, I'll have to think about it more. If you only care about integers, there is this neat result by Vincent Lefèvre (who knows much more about this stuff than I do):
http://www.vinc17.net/research/papers/rr_intdiv.pdf

Annoyingly, the actual CPU instructions (FPREM/FPREM1) also compute the last three bits of the quotient, which would be enough to determine which of the correct values we want, but there doesn't seem to be a good way to get to it: this is exposed by a libm remquo function, but the value is returned by reference which is still very slow in Julia at the moment.

[simonbyrne]

Interesting trivia: from what I can tell, computing remainders seems to be the only floating point functionality that is still handled by x87 instructions, instead of SSE.

[Yaffle]

in C++, possilbly, remainder ( http://en.cppreference.com/w/c/numeric/math/remainder ) can be used, if it is < 0, then result should be corrected.

[simonbyrne]

I think you're right.

The only potential problem I can think of with your proposed solution is x/y is in the range ( 2^52, 2^53 ) (when eps(x/y) = 1.0), and the actual numeric value of x/y has a fractional part of 1/2, in which case it can be rounded in either direction, depending on whether the quotient is odd or even.

As you suggest, using the C remainder should work, as it will match the rounding in this case.

But I'm also not sure this situation can actually happen: can you think of any such x, y for which this would occur?

[Yaffle]

But I'm also not sure this situation can actually happen: can you think of any such x, y for which this would occur?

No, I cannot, but I cannot garantee that there are no.