Why is 1/5 not equal to 1//5?

I am curious about the decision to have 1/5 == 1//5 return false.

I’m assuming there is a rationale (pun intended) and would like to understand why. Note that 1/2 == 1//2 returns true. This make sense in so far as some rational numbers can be represented exactly in binary, but mathematically it feels inconsistent.

An alternative approach might be something like this:

function myeq(x::T, q::Rational)::Bool where {T<:AbstractFloat}
typeof(x)(q) == x
end


In terms of performance:

julia> using BenchmarkTools

julia> @btime 1/5 == 1//5
0.875 ns (0 allocations: 0 bytes)
false

julia> @btime myeq(1/5, 1//5)
0.833 ns (0 allocations: 0 bytes)
true

1 Like

Fundamentally the reason this returns false is that they aren’t the same number. This isn’t really any different than 1/5 == Float32(1/5) returning false.

17 Likes

It’s also exactly the same reason that 2/3 != 2//3. That one seems more obvious because two thirds isn’t represented exactly in either base 10 (0.666666…) or base 2 (0.1010101010…). While one fifth is exactly representable in base 10, it’s not in base 2. It’s a repeating fraction just like two thirds… and so is similarly inexact. The rational 1//5 is exactly the number one fifth, whereas the floating point number is exactly:

julia> using Printf

julia> @printf "%.60f" 1/5
0.200000000000000011102230246251565404236316680908203125000000

18 Likes

Another reason to use strict equality is that anything (?) else creates a tangled web of paradoxes. For example, even if we made 1/5 == 1//5 we would still struggle with the famous (and non-Julia-specific)

julia> 1/10 + 2/10
0.30000000000000004

julia> 1/10 + 2/10 - 3/10
5.551115123125783e-17


If we made 1/5 == 1//5, then why shouldn’t 1/10 + 2/10 == 3/10? And there are situations such as

julia> (1/49) * 49
0.9999999999999999


Is it right to say that 1/49 == 1//49 when (1/49)*49 != (1//49)*49?

Floating point “errors” are not random (although in many contexts it’s okay to treat them as such). Floating point arithmetic is carefully defined to take the input arguments as they are (not assuming that 0.2 == 1//5, but in fact that it is exactly equal to what its representation implies: 7205759403792794 * 2^-55 – see Base.decompose(0.2)), perform the arithmetic exactly, and then round the result to the best representable value according to the floating point rounding configuration (usually “round to nearest, ties prefer even”).

Once you start accepting approximate notions of equality as “true equality,” it becomes very difficult to figure out where to draw the line any more. We have the isapprox function (and ≈ operator) for implementing approximate equality (with the option to manually specify tolerances, which you should always do in my opinion).

For generic use, chaos is likely to follow for anything except strict comparisons. This is why we do not attempt to make 0.2 == 1//5, even though 0.2 is the “written” representation of the closest floating point number to 1//5.

P.S.
In Julia and many languages, the “written” representation is the shortest string that rounds to the exact float value (so that a float can be converted losslessly back-and-forth between a string and native representation). One notable offender is MATLAB, who (probably to quiet a never-ending stream of complaints) prints 0.1 + 0.2 as 0.3 even though it still reports that 0.1 + 0.2 is not equal to 0.3. I find this behavior infuriating. EDIT: the above-linked website actually shows that many languages opt for this “cleaned-up” printing (except when using a printf construct with many digits). MATLAB is hardly alone.

6 Likes

You can also look at it the other way around — there does exist a rational number that exactly equals the Float64 that 1/5 evaluates to:

julia> r = 0b0011001100110011001100110011001100110011001100110011010//
0b10000000000000000000000000000000000000000000000000000000
0x000ccccccccccccd//0x0040000000000000

julia> r == 1/5
true

julia> @printf "%.60f" r
0.200000000000000011102230246251565404236316680908203125000000


I’ve written out the rational there using binary notation that should hopefully be evocative of writing a truncated repeating decimal — akin to writing 0.6666667 as \frac{6666667}{10000000}.

A floating point number isn’t some “fuzzy” representation of a range of numbers; it’s an exact value! It’s just that when you do some computation you might end up at a value that cannot be represented and then you round to the nearest point that does exist. I like to think of it as snapping to a grid.

15 Likes
4 Likes

Thanks all. I was aware that floating point is not exact, but I see the appeal of having == only be true if the equality is exact. And this is an option:

julia> 1/3 ≈ 1//3
true

8 Likes

Yes, you can choose non-default (even stochastic rounding with a package), but that’s for operations, for literals, i.e. the case here it’s always that same (that) rounding(?!).

We have the isapprox function (and ≈ operator) for implementing approximate equality (with the option to manually specify tolerances, which you should always do in my opinion).

Are you saying we should use the operator (or function) always (not == for float)s? Or “specify” tolerance, i.e. use the function only, since the operator doesn’t support that. I.e. is the operator considered harmful? It has the same default as the function, and I thought those defaults sensible (there also)?!

Note e.g.

julia> 2/3 ≈ Float32(2/3) and for all numbers I can think of, why not use \isapprox<TAB>?
true


Note though == is actually faster, but that doesn’t seem like an argument to use it, or the slower ≈, rather neither, e.g. < in most cases. I think equality in some form isn’t that common in (performance critical) float code. Maybe inequality in some for though, and ≉ i.e. \napprox.

I was trying to prove (all types of) approx comparison slower, with or without missing, but fail because of global scope maybe):

julia> cmp(a, b) = a ≉ b
julia> a= 1/10 + missing; b = 3/10;

julia> typeof(a)
Missing

julia> @code_native cmp(a, b)


Rounding of basic hardware arithmetic is controlled by the floating point environment, which is set via processor flags. Julia has a stub function for this, setrounding, but it’s currently (v1.9) not supported for any types except BigFloat.

I’m saying that it’s good to be explicit about what you mean by “is approximately equal” for anything but “throwaway” work, rather than relying on the default tolerances. I don’t expect the default relative tolerance of sqrt(eps(T)) will make every Float64(a/b) ≈ Float32(a/b) (since Float64 has 52 bits of mantissa and Float32 has only 23 – less than half), but it will be true for the vast majority of values. To be safe here, I’d use a relative tolerance more like 4eps(Float32).

To be honest, I’m not sure I’ve ever actually used isapprox. If you want your result to be a rational number, you probably want to do the calculation using rationals start-to-finish. Another option is to use rationalize to find a rational approximation to a float, but I’ve rarely found it useful in practice. I find most of my comparisons with floats are either inequalities or checks for equality with zero.

1 Like

Then ≈ doesn’t apply (nor ==), I mean for relative tolerance, since it’s infinite. So when do you use == or ≈ in real-world float code (I only see it useful for informal (test) code, or in tests)? I suppose you could use isapprox for the absolute tolerance, then for 0.0 too.

Now I’m curious, I think literals, i.e. for big"0.2" setrounding doesn’t apply. And neither for BigFloat(0.2). I.e. the former isn’t a calculation/operation? Or if so considered, then still only at compile time parse time so doesn’t apply, but would for the latter?

Another reason to use strict equality is that anything (?) else creates a tangled web of paradoxes.

You can use equality with interval arithmetic, but that probably also has some issues.

big"0.2" is equivalent to @big_str("0.2") (and any xyz"content" is equivalent to @xyz_str("content")). You can @edit big"0.2" to see the macro source code, but the short version is that it uses tryparse to parse the content of the string as a BigInt or BigFloat at parse time. The reason this macro doesn’t respect setrounding is because macros are evaluated at parse time rather than during the actual setrounding call (by that time it’s already a BigFloat literal).

BigFloat(0.2) just promotes 0.2::Float64 to BigFloat, so it is exactly equal to 0.2::Float64 and no rounding is required (since BigFloat covers a strict superset of Float64).

Observe the following options:

# does not work because the string is evaluated before the setrounding call
julia> setrounding(()->big"0.2",BigFloat,RoundDown)
0.2000000000000000000000000000000000000000000000000000000000000000000000000000004

# evalautes the string as part of the call, so does work
julia> setrounding(()->tryparse(BigFloat,"0.2"),BigFloat,RoundDown)
0.1999999999999999999999999999999999999999999999999999999999999999999999999999983

# basically the same as above, but simpler to write
julia> BigFloat("0.2",RoundDown)
0.1999999999999999999999999999999999999999999999999999999999999999999999999999983

# actually changes the environment before evaluating the @big_str macro
julia> old = rounding(BigFloat); setrounding(BigFloat,RoundDown); z = big"0.2"; setrounding(BigFloat,old); z
0.1999999999999999999999999999999999999999999999999999999999999999999999999999983


I’m not familiar with interval arithmetic and have certainly never used it. But it seems that (I may be wrong) “equality” with intervals is actually ranges of inequalities low <= x <= high (so is semantically more like x in interval than x == value), so I don’t think it’s really fair to call this “equality” in the spirit of this conversation. I definitely have used statements like abs(abs(x)-1) <= tol || error("x not unit magnitude") to check that x has a magnitude of approximately 1, but this is only an approximate equality check. So in that sense, I suppose I have used something like isapprox (with a manual tolerance) in practice.

1 Like

I don’t really think there’s much choice here if you want a coherent design. Suppose we made

1//5 == 1/5


Presumably we’d also want that to be true for different types, so we’d also have:

1//5 == 1f0/5f0
1//5 == big(1)/big(5)


Transitivity of equality then implies that

1f0/5f0 == 1/5 == big(1)/big(5)


But those are clearly different values as they have different decimal expansions:

julia> @printf("%0.70f\n", 1f0/5f0)
0.2000000029802322387695312500000000000000000000000000000000000000000000

julia> @printf("%0.70f\n", 1/5)
0.2000000000000000111022302462515654042363166809082031250000000000000000

julia> @printf("%0.70f\n", big(1)/big(5))
0.2000000000000000000000000000000000000000000000000000000000000000000000


Making these equal while they have different decimal values seems pretty questionable.

Another way to put this is that if equality is going to be an equivalence relation (hard to argue with that), then we need to pick equivalence classes for all numerical values. The classes Julia uses are the same ones that mathematics itself uses: two numerical values are only equal if they represent the same number. Since the rational value 1//5 represents the fraction 1/5 it can only be equal to other representations of that value and since 1/5 cannot be represented in binary floating-point, no floating-point value is equal to 1//5.

There are other equivalence classes one could pick. For example, you could make Float64 special and compare everything based on whether it maps to the same Float64 or not. In particular, that would give you 1//5 == 1/5 since Float64(1//5) is 1/5. For lower precision binary types like Float32 this wouldn’t change the notion of equality since they are exactly embedded in Float64. But for higher precision types like BigFloat, this would mean that very many values would fall into the same equivalence class. For example, we currently have:

julia> x = big(1)/big(10)
0.1000000000000000000000000000000000000000000000000000000000000000000000000000002

julia> y = x + 0.5*eps(1/10)
0.1000000000000000069388939039072283776476979255676269531250000000000000000000002

julia> x != y
true


But if we defined equality in terms of Float64, then we’d have x == y since Float64(x) == Float64(y). Maybe that’s not so bad—it’s a bit weird for two numbers with so many differing significant digits to be unequal, but maybe for binary float types that people think of as “approximate” it could be tolerable. But it becomes far less palatable if you consider “exact” types like rationals. Consider:

julia> a = 1//5
1//5

julia> b = 3602879701896397//18014398509481984
3602879701896397//18014398509481984


Obviously these are different and non-equal, right? But if we define equality in terms of Float64 equivalence classes, these values have to be considered equal since they map to the same Float64:

julia> Float64(a) == Float64(b)
true


That’s pretty unacceptable for a rational number type. What if we patch up this badness by exempting rational number types from this definition of equality? We’d still be in trouble then since we’d have a == 0.2 and b == 0.2 since those comparison are done via conversion to Float64, yet we’d have a != b. In other words, transitivity of equality would fail.

The bottom line is this: if you want a coherent notion of equality, you need to pick consistent equivalence classes of numbers and apply that equivalence relation to all different representations of numbers. If you try to mix and match different equivalence relations for different types, you end up with failures of transitivity. While there are many possible equivalence classes you could pick, the only choice that matches mathematics itself is for each numerical value to be in its own equivalence class. Which is exactly what we do.

22 Likes

Thank you. I like very much the importance that == be an equivalence relation!

4 Likes

What’s worse, I thought this would be true (and thus ≈ your go to operator solving all, I suppose though ok for most values in practice), but a bit surprisingly:

julia> 3602879701896397//18014398509481984 ≈ 1//5
false

julia> 3602879701896397/18014398509481984
0.2


Ah, this appears to be a bug. The rtol for comparing two rationals is defaulting to 0 because the methods for Base.rtoldefault are

rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))
rtoldefault(::Type{<:Real}) = 0


This second method seems relatively suspect to me. (alternatively we should have a method for rtoldefault(::Type{Rational{T}}) where {T} = one(T) // isqrt(typemax(T)) which would make this behave as expected.)

2 Likes

Seems right to me? I don’t see any way to define a nonzero default rtol for an arbitrary numeric type; 0 is the only safe fallback.

rtol = one(T) // isqrt(typemax(T)) seems like a somewhat reasonable default for Rational{T}, but I also don’t quite see in what context it is useful, so it’s hard to judge.

6 Likes