fix rational vs float comparisons (attempt 2)

base/gmp.jl

@@ -5,9 +5,9 @@ export BigInt
 import Base: *, +, -, /, <, <<, >>, >>>, <=, ==, >, >=, ^, (~), (&), (|), ($),
              binomial, cmp, convert, div, divrem, factorial, fld, gcd, gcdx, lcm, mod,
              ndigits, promote_rule, rem, show, isqrt, string, isprime, powermod,
-             widemul, sum, trailing_zeros, trailing_ones, count_ones, base, parseint,
+             sum, trailing_zeros, trailing_ones, count_ones, base, parseint,
              serialize, deserialize, bin, oct, dec, hex, isequal, invmod,
-             prevpow2, nextpow2, ndigits0z, widen
+             prevpow2, nextpow2, ndigits0z, widen, signed
 
 if Clong == Int32
     typealias ClongMax Union(Int8, Int16, Int32)
@@ -47,6 +47,8 @@ widen(::Type{Int128})  = BigInt
 widen(::Type{UInt128}) = BigInt
 widen(::Type{BigInt})  = BigInt
 
+signed(x::BigInt) = x
+
 BigInt(x::BigInt) = x
 BigInt(s::AbstractString) = parseint(BigInt,s)
 
@@ -472,9 +474,6 @@ ndigits(x::BigInt, b::Integer=10) = x.size == 0 ? 1 : ndigits0z(x,b)
 
 isprime(x::BigInt, reps=25) = ccall((:__gmpz_probab_prime_p,:libgmp), Cint, (Ptr{BigInt}, Cint), &x, reps) > 0
 
-widemul(x::Int128, y::UInt128)  = BigInt(x)*BigInt(y)
-widemul(x::UInt128, y::Int128)  = BigInt(x)*BigInt(y)
-
 prevpow2(x::BigInt) = x.size < 0 ? -prevpow2(-x) : (x <= 2 ? x : one(BigInt) << (ndigits(x, 2)-1))
 nextpow2(x::BigInt) = x.size < 0 ? -nextpow2(-x) : (x <= 2 ? x : one(BigInt) << ndigits(x-1, 2))
 

base/hashing2.jl

@@ -94,25 +94,36 @@ Special values:
 decompose(x::Integer) = x, 0, 1
 decompose(x::Rational) = num(x), 0, den(x)
 
+function decompose(x::Float16)
+    isnan(x) && return 0, 0, 0
+    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
+    n = reinterpret(UInt16, x)
+    s = (n & 0x03ff) % Int16
+    e = (n & 0x7c00 >> 10) % Int
+    s |= int16(e != 0) << 10
+    d = ifelse(signbit(x), -1, 1)
+    int(s), int(e - 25 + (e == 0)), d
+end
+
 function decompose(x::Float32)
     isnan(x) && return 0, 0, 0
     isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
-    n = reinterpret(Int32, x)
-    s = int32(n & 0x007fffff)
-    e = int32(n & 0x7f800000 >> 23)
+    n = reinterpret(UInt32, x)
+    s = (n & 0x007fffff) % Int32
+    e = (n & 0x7f800000 >> 23) % Int
     s |= int32(e != 0) << 23
-    d = ifelse(signbit(n), -1, 1)
+    d = ifelse(signbit(x), -1, 1)
     int(s), int(e - 150 + (e == 0)), d
 end
 
 function decompose(x::Float64)
     isnan(x) && return 0, 0, 0
     isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
-    n = reinterpret(Int64, x)
-    s = int64(n & 0x000fffffffffffff)
-    e = int64(n & 0x7ff0000000000000 >> 52)
+    n = reinterpret(UInt64, x)
+    s = (n & 0x000fffffffffffff) % Int64
+    e = (n & 0x7ff0000000000000 >> 52) % Int
     s |= int64(e != 0) << 52
-    d = ifelse(signbit(n), -1, 1)
+    d = ifelse(signbit(x), -1, 1)
     int(s), int(e - 1075 + (e == 0)), d
 end
 

base/int.jl

@@ -363,6 +363,17 @@ widen(::Type{UInt16}) = UInt
 widen(::Type{UInt32}) = UInt64
 widen(::Type{UInt64}) = UInt128
 
+# a few special cases,
+# Int64*UInt64 => Int128
+# |x|<=2^(k-1), |y|<=2^k-1   =>   |x*y|<=2^(2k-1)-1
+widemul(x::Signed,y::Unsigned) = widen(x)*signed(widen(y))
+widemul(x::Unsigned,y::Signed) = signed(widen(x))*widen(y)
+# multplication by Bool doesn't require widening
+widemul(x::Bool,y::Bool) = x*y
+widemul(x::Bool,y::Number) = x*y
+widemul(x::Number,y::Bool) = x*y
+
+
 ## float to integer coercion ##
 
 # requires int arithmetic defined, for the loops to work

base/rational.jl

@@ -177,30 +177,78 @@ end
 /(x::Rational, z::Complex ) = inv(z/x)
 
 ==(x::Rational, y::Rational) = (x.den == y.den) & (x.num == y.num)
-==(x::Rational, y::Integer ) = (x.den == 1) & (x.num == y)
-==(x::Integer , y::Rational) = y == x
-
-# needed to avoid ambiguity between ==(x::Real, z::Complex) and ==(x::Rational, y::Number)
-==(z::Complex , x::Rational) = isreal(z) & (real(z) == x)
-==(x::Rational, z::Complex ) = isreal(z) & (real(z) == x)
-
-==(x::FloatingPoint, q::Rational) = count_ones(q.den) <= 1 & (x == q.num/q.den) & (x*q.den == q.num)
-==(q::Rational, x::FloatingPoint) = count_ones(q.den) <= 1 & (x == q.num/q.den) & (x*q.den == q.num)
-
-# TODO: fix inequalities to be in line with equality check
 < (x::Rational, y::Rational) = x.den == y.den ? x.num < y.num :
                                widemul(x.num,y.den) < widemul(x.den,y.num)
-< (x::Rational, y::Integer ) = x.num < widemul(x.den,y)
-< (x::Rational, y::Real    ) = x.num < x.den*y
-< (x::Integer , y::Rational) = widemul(x,y.den) < y.num
-< (x::Real    , y::Rational) = x*y.den < y.num
-
 <=(x::Rational, y::Rational) = x.den == y.den ? x.num <= y.num :
                                widemul(x.num,y.den) <= widemul(x.den,y.num)
+
+
+==(x::Rational, y::Integer ) = (x.den == 1) & (x.num == y)
+==(x::Integer , y::Rational) = y == x
+< (x::Rational, y::Integer ) = x.num < widemul(x.den,y)
+< (x::Integer , y::Rational) = widemul(x,y.den) < y.num
 <=(x::Rational, y::Integer ) = x.num <= widemul(x.den,y)
-<=(x::Rational, y::Real    ) = x.num <= x.den*y
 <=(x::Integer , y::Rational) = widemul(x,y.den) <= y.num
-<=(x::Real    , y::Rational) = x*y.den <= y.num
+
+function ==(x::FloatingPoint, q::Rational)
+    if isfinite(x)
+        (count_ones(q.den) == 1) & (x*q.den == q.num)
+    else
+        x == q.num/q.den
+    end
+end
+
+==(q::Rational, x::FloatingPoint) = x == q
+
+for rel in (:<,:<=,:cmp)
+    for (Tx,Ty) in ((Rational,FloatingPoint), (FloatingPoint,Rational))
+        @eval function ($rel)(x::$Tx, y::$Ty)
+            if isnan(x) || isnan(y)
+                $(rel == :cmp ? :(throw(DomainError())) : :(return false))
+            end
+
+            xn, xp, xd = decompose(x)
+            yn, yp, yd = decompose(y)
+
+            if xd < 0
+                xn = -xn
+                xd = -xd
+            end
+            if yd < 0
+                yn = -yn
+                yd = -yd
+            end
+
+            xc, yc = widemul(xn,yd), widemul(yn,xd)
+            xs, ys = sign(xc), sign(yc)
+
+            if xs != ys
+                return ($rel)(xs,ys)
+            elseif xs == 0
+                # both are zero or ±Inf
+                return ($rel)(xn,yn)
+            end
+
+            xb, yb = ndigits0z(xc,2) + xp, ndigits0z(yc,2) + yp
+
+            if xb == yb
+                xc, yc = promote(xc,yc)
+                if xp > yp
+                    xc = (xc<<(xp-yp))
+                else
+                    yc = (yc<<(yp-xp))
+                end
+                return ($rel)(xc,yc)
+            else
+                return xc > 0 ? ($rel)(xb,yb) : ($rel)(yb,xb)
+            end
+        end
+    end
+end
+
+# needed to avoid ambiguity between ==(x::Real, z::Complex) and ==(x::Rational, y::Number)
+==(z::Complex , x::Rational) = isreal(z) & (real(z) == x)
+==(x::Rational, z::Complex ) = isreal(z) & (real(z) == x)
 
 for op in (:div, :fld, :cld)
     @eval begin

test/numbers.jl

@@ -797,6 +797,23 @@ end
 @test nextfloat(0.0) == 1//(BigInt(2)^1074)
 @test nextfloat(0.0) != 1//(BigInt(2)^1074-1)
 
+@test 1/3 < 1//3
+@test !(1//3 < 1/3)
+@test -1/3 < 1//3
+@test -1/3 > -1//3
+@test 1/3 > -1//3
+@test 1/5 > 1//5
+@test 1//3 < Inf
+@test 0//1 < Inf
+@test 1//0 == Inf
+@test -1//0 == -Inf
+@test -1//0 != Inf
+@test 1//0 != -Inf
+@test !(1//0 < Inf)
+@test !(1//3 < NaN)
+@test !(1//3 == NaN)
+@test !(1//3 > NaN)
+
 @test sqrt(2) == 1.4142135623730951
 
 @test 1+1.5 == 2.5
@@ -1977,6 +1994,10 @@ end
 @test widen(BigInt) === BigInt
 
 @test widemul(typemax(Int64),typemax(Int64)) == 85070591730234615847396907784232501249
+@test typeof(widemul(Int64(1),UInt64(1))) == Int128
+@test typeof(widemul(UInt64(1),Int64(1))) == Int128
+@test typeof(widemul(Int128(1),UInt128(1))) == BigInt
+@test typeof(widemul(UInt128(1),Int128(1))) == BigInt
 
 # .//
 @test [1,2,3] // 4 == [1//4, 2//4, 3//4]