Skip to content

Commit

Permalink
fix isless and cmp/lexcmp for floating point
Browse files Browse the repository at this point in the history 
for now cmp() uses the total ordering, but we might change it to give a
DomainError for NaN arguments
  • Loading branch information
@JeffBezanson
JeffBezanson committed Jan 3, 2014
1 parent 756a34a commit 2343ba0
Show file tree
Hide file tree
Showing 6 changed files with 52 additions and 15 deletions.
33 changes: 31 additions & 2 deletions base/float.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,8 +158,37 @@ isequal(x::Float64, y::Float64) = fpiseq(unbox(Float64,x),unbox(Float64,y))
isless (x::Float32, y::Float32) = fpislt(unbox(Float32,x),unbox(Float32,y))
isless (x::Float64, y::Float64) = fpislt(unbox(Float64,x),unbox(Float64,y))

isless (a::Integer, b::FloatingPoint) = (a<b) | isless(float(a),b)
isless (a::FloatingPoint, b::Integer) = (a<b) | isless(a,float(b))
isless(a::FloatingPoint, b::FloatingPoint) =
(a<b) | (!isnan(a) & (isnan(b) | (signbit(a)>signbit(b))))
isless(a::Real, b::FloatingPoint) = (a<b) | isless(float(a),b)
isless(a::FloatingPoint, b::Real) = (a<b) | isless(a,float(b))

function cmp(x::FloatingPoint, y::FloatingPoint)
# TODO: for now cmp() is a total order, but we might move this logic to
# lexcmp in the future.
x<y && return -1
x>y && return 1
if isnan(x)
isnan(y) && return 0
return 1
end
isnan(y) && return -1
return signbit(y) - signbit(x)
end

function cmp(x::Real, y::FloatingPoint)
x<y && return -1
x>y && return 1
isnan(y) && return -1
return signbit(y) - signbit(x)
end

function cmp(x::FloatingPoint, y::Real)
x<y && return -1
x>y && return 1
isnan(x) && return 1
return signbit(y) - signbit(x)
end

==(x::Float64, y::Int64 ) = eqfsi64(unbox(Float64,x),unbox(Int64,y))
==(x::Float64, y::Uint64 ) = eqfui64(unbox(Float64,x),unbox(Uint64,y))
Expand Down
6 changes: 1 addition & 5 deletions base/mpfr.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ export

import
Base: (*), +, -, /, <, <=, ==, >, >=, ^, besselj, besselj0, besselj1, bessely,
bessely0, bessely1, ceil, cmp, convert, copysign, degrees2radians,
bessely0, bessely1, ceil, convert, copysign, degrees2radians,
exp, exp2, exponent, factorial, floor, hypot, isinteger, iround,
isfinite, isinf, isnan, ldexp, log, log2, log10, max, min, mod, modf,
nextfloat, prevfloat, promote_rule, radians2degrees, rem, round, show,
Expand Down Expand Up @@ -356,10 +356,6 @@ function -(x::BigFloat)
return z
end

function cmp(x::BigFloat, y::BigFloat)
ccall((:mpfr_cmp, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y)
end

function sqrt(x::BigFloat)
z = BigFloat()
ccall((:mpfr_sqrt, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
Expand Down
8 changes: 4 additions & 4 deletions base/operators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,10 +28,13 @@ isless(x::Real, y::Real) = x<y

ifelse(c::Bool, x, y) = Intrinsics.select_value(c, x, y)

cmp(x,y) = isless(x,y) ? -1 : isless(y,x) ? 1 : 0
cmp(x,y) = isless(x,y) ? -1 : ifelse(isless(y,x), 1, 0)
lexcmp(x,y) = cmp(x,y)
lexless(x,y) = lexcmp(x,y)<0

# cmp returns -1, 0, +1 indicating ordering
cmp(x::Real, y::Real) = int(sign(x-y))

max(x,y) = ifelse(y < x, x, y)
min(x,y) = ifelse(x < y, x, y)

Expand Down Expand Up @@ -112,9 +115,6 @@ mod1{T<:Real}(x::T, y::T) = y-mod(y-x,y)
rem1{T<:Real}(x::T, y::T) = rem(x-1,y)+1
fld1{T<:Real}(x::T, y::T) = fld(x-1,y)+1

# cmp returns -1, 0, +1 indicating ordering
cmp{T<:Real}(x::T, y::T) = int(sign(x-y))

# transposed multiply
Ac_mul_B (a,b) = ctranspose(a)*b
A_mul_Bc (a,b) = a*ctranspose(b)
Expand Down
2 changes: 0 additions & 2 deletions base/promotion.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,6 @@ mod(x::Real, y::Real) = mod(promote(x,y)...)
mod1(x::Real, y::Real) = mod1(promote(x,y)...)
rem1(x::Real, y::Real) = rem1(promote(x,y)...)
fld1(x::Real, y::Real) = fld1(promote(x,y)...)
cmp(x::Real, y::Real) = cmp(promote(x,y)...)

max(x::Real, y::Real) = max(promote(x,y)...)
min(x::Real, y::Real) = min(promote(x,y)...)
Expand All @@ -202,4 +201,3 @@ no_op_err(name, T) = error(name," not defined for ",T)
max{T<:Real}(x::T, y::T) = ifelse(y < x, x, y)
min{T<:Real}(x::T, y::T) = ifelse(x < y, x, y)
minmax{T<:Real}(x::T, y::T) = x < y ? (x, y) : (y, x)

11 changes: 9 additions & 2 deletions test/mpfr.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,6 +222,15 @@ imi = BigFloat(-Inf)
@test !isequal(BigFloat(0.0),BigFloat(-0.0))
@test isequal(z, BigFloat(NaN))

# total ordering
@test isless(big(-0.0), big(0.0))
@test cmp(big(-0.0), big(0.0)) == -1
@test cmp(big(0.0), big(-0.0)) == 1
@test isless(big(1.0), big(NaN))
@test cmp(big(1.0), big(NaN)) == -1
@test cmp(big(NaN), big(NaN)) == 0
@test cmp(big(NaN), big(1.0)) == 1

# signbit
@test signbit(BigFloat(-1.0)) == 1
@test signbit(BigFloat(1.0)) == 0
Expand Down Expand Up @@ -741,5 +750,3 @@ tol = 1e-3
@test_approx_eq_eps f+(1//1) g tol

@test_approx_eq_eps f+one(Rational{BigInt}) g tol

# new tests
7 changes: 7 additions & 0 deletions test/numbers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -466,6 +466,13 @@ end
@test isless(-0.0, 0)
@test !isless( 0,-0.0)

@test isless(-0.0, 0.0f0)
@test lexcmp(-0.0, 0.0f0) == -1
@test lexcmp(0.0, -0.0f0) == 1
@test lexcmp(NaN, 1) == 1
@test lexcmp(1, NaN) == -1
@test lexcmp(NaN, NaN) == 0

for x=-5:5, y=-5:5
@test (x==y)==(float64(x)==int64(y))
@test (x!=y)==(float64(x)!=int64(y))
Expand Down

8 comments on commit 2343ba0

@JeffBezanson
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ref #5290

@StefanKarpinski
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I really think you need to write up what the logic behind these comparison operators is. I lost track when lexcmp was introduced, which I still am not entirely happy with.

@JeffBezanson
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, as I said in the comments, this very well may not be the final version. lexcmp was introduced much later than cmp, and we might have handled floats differently if they'd been introduced at the same time.

We had cmp fall back to isless, meaning cmp uses the standard total order. I kept this behavior in this commit, for now.

There are accurate descriptions of all of these in the help text. We have

  • Canonical partial order, <
  • Canonical total order, cmp, isless
  • Lexicographic total order, lexcmp, lexless

Our total order on floats (isless) is not IEEE-compliant, because of -NaN. This doesn't bother me too much and I don't think it's crucial to change it.

But, we have a few options:

  1. Make cmp mean lexicographic total order, and get rid of lexcmp
  2. Define cmp in terms of <, and give a DomainError on NaNs, but allow lexcmp to work on NaNs
  3. Keep isless as it is now, but have lexcmp and lexless implement the IEEE order

@jiahao
Copy link
Member

@jiahao jiahao commented on 2343ba0 Jan 4, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Once a decision is made, the description should go in the manual together with an explanation of the differences between ==, isequal and is/===.

@gitfoxi
Copy link
Contributor

@gitfoxi gitfoxi commented on 2343ba0 Jan 5, 2014 via email

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jiahao
Copy link
Member

@jiahao jiahao commented on 2343ba0 Jan 5, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@gitfoxi I appreciate the thoughts on performance, but correctness has to come first. Obviously being fast and correct is ideal, but being fast and wrong is just pointless.

Some definition of comparison is a single instruction and for common CPUs it's the IEEE one

I believe it's the arithmetic comparison that can compile down to a single instruction, not the lexical one which is being corrected here.

@JeffBezanson
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not as bad as all that. Our < et. al. are what you typically use in numerical code, and they are of the single-instruction variety. In sorting, you have to deal with NaNs specially anyway, so we have careful custom sorting for float types. In fact even the typical float < is not optimal, and we use integer comparisons instead. But the way this sorting algorithm behaves, at the end of the day, corresponds to isless.

@StefanKarpinski
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some definition of comparison is a single instruction and for common CPUs it's the IEEE one. The performance of a comparison-centered algorithm, such as sorting, will depend on it.

You cannot actually use the standard IEEE 754 float comparison operations to correctly sort floating-point numbers. If you use IEEE ordering to sort floating-point numbers, most comparison-based algorithms will be incorrect because they depend on the sorted values being totally ordered. The presence of unordered elements like NaN completely screws things up because !(x < y) does not imply y <= x. Moreover, your negative and positive zeros will be arbitrarily mixed together, which is not desirable.

To correctly sort floating-point numbers you have to use a total ordering. One such ordering is specified by IEEE 754: negative NaNs, then negative values from -Inf to -0.0, then positive values from 0.0 to Inf, then positive NaNs. Matlab uses a different ordering: negative values from -Inf to -0.0, then positive values from 0.0 to Inf, and then all NaNs, regardless of their sign or payload – in fact, the sign and payload of NaNs will be destroyed by sorting. Neither of these total orderings is implemented in hardware that I'm aware of. Moreover, expressing either of these total orderings is pretty complicated. We implement isless(x::Float64, y::Float64) using an LLVM intrinsic like so:

    HANDLE(fpislt,2) {
        Value *xi = JL_INT(x);
        Value *yi = JL_INT(y);
        x = FP(x);
        fy = FP(y);
        return builder.CreateOr(
            builder.CreateAnd(
                builder.CreateFCmpORD(x, x),
                builder.CreateFCmpUNO(fy, fy)
            ),
            builder.CreateAnd(
                builder.CreateFCmpORD(x, fy),
                builder.CreateOr(
                    builder.CreateAnd(
                        builder.CreateICmpSGE(xi, ConstantInt::get(xi->getType(), 0)),
                        builder.CreateICmpSLT(xi, yi)
                    ),
                    builder.CreateAnd(
                        builder.CreateICmpSLT(xi, ConstantInt::get(xi->getType(), 0)),
                        builder.CreateICmpUGT(xi, yi)
                    )
                )
            )
        );
    }

This translates approximately to

(!isnan(x) & isnan(y)) | !isnan(x) & !isnan(y) & ((xi >= 0) & (xi < yi) | (xi < 0) & (xi > yi))

where xi and yi are x and y respectively reinterpreted as integers. Although this is actually faster than one might expect, fortunately to do actual sorting of floating-point numbers, we don't use this operation at all. Instead, we do the following:

  1. Move all the NaNs to the end of an array in a way that's stable for both the non-NaN and NaN parts.
  2. Pivot the non-NaN numbers by their sign, separating negative and positive floats.
  3. Use integer comparisons to sort the negative and positive regions separately.

This is all implemented in base/sort.jl and composes transparently with arbitrary sorting algorithms. The end-result is that floating-point arrays are sorted as if they used the isless operation with a stable sorting algorithm. The stability matters because NaNs are all equivalent according to isless but they are actually programmatically distinguishable. Moreover, we want sorting an untyped array that happens to only hold floats to end up in the same order as a typed array of floats. Since sorting an untyped array would use a stable sort and the isless operation, the NaNs end up at the end of the untyped array in the same order they appeared in the original array. That's why we do step 1 above that way. We verify that our sorting is stable for NaNs by generating NaNs with random payloads and sorting them.

The end result of all of this is that our floating-point sorting is fast, correct, and consistent, regardless of the type of array that's sorted.

Please sign in to comment.