From 2343ba0c0c362f355f6841308d00228a13bcca06 Mon Sep 17 00:00:00 2001
From: Jeff Bezanson <jeff.bezanson@gmail.com>
Date: Fri, 3 Jan 2014 18:26:22 -0500
Subject: [PATCH] fix isless and cmp/lexcmp for floating point

for now cmp() uses the total ordering, but we might change it to give a
DomainError for NaN arguments
---
 base/float.jl     | 33 +++++++++++++++++++++++++++++++--
 base/mpfr.jl      |  6 +-----
 base/operators.jl |  8 ++++----
 base/promotion.jl |  2 --
 test/mpfr.jl      | 11 +++++++++--
 test/numbers.jl   |  7 +++++++
 6 files changed, 52 insertions(+), 15 deletions(-)

diff --git a/base/float.jl b/base/float.jl
index 12cf8518bf918..b8b0b8c13f70e 100644
--- a/base/float.jl
+++ b/base/float.jl
@@ -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))
diff --git a/base/mpfr.jl b/base/mpfr.jl
index c4302627803ea..5867fe2a5f652 100644
--- a/base/mpfr.jl
+++ b/base/mpfr.jl
@@ -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,
@@ -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])
diff --git a/base/operators.jl b/base/operators.jl
index dc38a2d083426..4a66d39e62b82 100644
--- a/base/operators.jl
+++ b/base/operators.jl
@@ -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)
 
@@ -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)
diff --git a/base/promotion.jl b/base/promotion.jl
index 538d2ebbf3e87..8992033d0b192 100644
--- a/base/promotion.jl
+++ b/base/promotion.jl
@@ -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)...)
@@ -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)
-
diff --git a/test/mpfr.jl b/test/mpfr.jl
index 1bfc656ce14eb..4b3d389b0de7b 100644
--- a/test/mpfr.jl
+++ b/test/mpfr.jl
@@ -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
@@ -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
diff --git a/test/numbers.jl b/test/numbers.jl
index eeed13e5ce1de..5a451f2986c84 100644
--- a/test/numbers.jl
+++ b/test/numbers.jl
@@ -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))