Calculation of remainder, floating point issues

[lwhitefox]

The following remainder calculations are a bit surprising to me, as I’d expect to get a remainder of 0 (or something like ~1e-18) in each case. I guess there is some floating point issue:

julia> rem(500,0.625)
0.0

julia> rem(50,0.0625)
0.0

julia> rem(5,0.00625)
0.006249999999999723

julia> rem(0.5,0.000625)
0.0006249999999999896

julia> rem(0.05,0.0000625)
1.734723475976807e-18

In particular, for the 3rd and 4th results, the generated remainder is essentially the same as the second argument, e.g.:

julia> rem(0.5, 0.000625) ≈ 0.000625
true

If I’m testing for a 0 remainder, then the implication is that I should test for a remainder that is either 0 or the divisor, but I’m curious if this expected or if there is an issue here.

[DNF]

One would expect the remainder to either be very close to 0.0 or very close to 0.0...625, depending on whether it’s slightly bigger or smaller than a multiple. Or rather, whether the floating point roundoffs pull in one direction or the other.

[_bernhard]

This is an issue with floating point numbers in general. in contrast to integers floats are (for most values) not exact, but rounded to the nearest value that can be represented.

This is exactly what happens in your example - the divisor (and dividend) is rounded to the nearest float and depending on whether this rounding happens up or down, the last remainder is either slightly larger than the divisor (you get a result of ~2e-18) or slightly smaller (you get a result of ~divisor - 2e-18).

[natlampen]

The documentation hints that rem(x,y,r) is the same as x - y*round(x/y,r). For x = 5 and y = 0.00625 this actually returns 0. So maybe a bug is hiding somewhere?

[lwhitefox]

The comments above about floating point roundoffs are what I was alluding to as “floating point issues”. If the remainder value is e.g. 2e-18, that’s clearly roundoff.

However, @natlampen hits the nail on the head - when the “remainder” value is essentially equal to the divisor, it seems like this is an obvious result to correct. Hence I also think about a “bug hiding somewhere”.

[stevengj]

@_bernhard

This is exactly what happens in your example - the divisor (and dividend) is rounded to the nearest float

In particular, 0.00625 is not exactly representable in binary floating point, so it is rounded to a slightly larger value:

julia> big(0.00625)
0.00625000000000000034694469519536141888238489627838134765625

(In contrast, 5.0 is exactly represented.) So, you are actually computing the remainder of 5 \div 0.00625000000000000034694469519536141888238489627838134765625, which is not zero. The correct quotient is therefore 799, not 800 (=5/0.00625), and the correct remainder computed in arbitrary precision arithmetic to about 300 decimal digits (setprecision(BigFloat, 1024)) is:

julia> 5 - 799 * big(0.00625)
0.00624999999999972279118853890622631297446787357330322265625

If we round this to the nearest Float64 value, we get:

julia> Float64(5 - 799 * big(0.00625))
0.006249999999999723

julia> rem(5, 0.00625)
0.006249999999999723

So, in fact, rem is giving exactly the correct answer to the question you asked (i.e., the exact answer rounded to Float64).

It’s just that you are not asking the question you thought you were asking, because 0.00625 in floating point is not the number you think it is.

[lwhitefox]

Thanks @stevengj this is a very clear example. Again, the general issues with floating point representation are clear - in fact I tried t=displaying values following more or less the logic you did. The mistake I made was to try to use @printf and I stopped a few digits short; next time I will use big which I didn’t know of before and thoroughly displays the issues. If I’d used big I wouldn’t have posted

It would be useful to include a 1-line example like this in the docs just to help make clear what the implications of “exact” might be

[stevengj]

@lwhitefox

It would be useful to include a 1-line example like this in the docs just to help make clear what the implications of “exact” might be

Sure, but realize also that this is not specific to Julia. The same floating-point standard is followed by essentially all modern languages, and there is lots of learning material online about this (including in Julia).

(Whether big prints all of the digits, i.e. an exact decimal representation, depends on the current precision setting, since BigFloat is variable precision. I don’t remember off the top of my head if there is a function that always prints the exact decimal representation of a floating-point value, regardless of the precision and regardless of how many decimal digits are required?)