Skip to content

Commit

Permalink
rationals: make conversion from float to rational exact.
Browse files Browse the repository at this point in the history 
This is a fairly major change, but it's in line with the recent
change to make comparisons of floats and rationals strict [#3102].
There is still some work to be done, namely:

  - handle large values correctly
  - rational-float inequalities
  • Loading branch information
@StefanKarpinski
StefanKarpinski committed Jun 11, 2013
1 parent 6d5b8d7 commit 6125749
Show file tree
Hide file tree
Showing 4 changed files with 70 additions and 41 deletions.
1 change: 1 addition & 0 deletions base/exports.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -388,6 +388,7 @@ export
prevpow2,
primes,
radians2degrees,
rationalize,
real,
realmax,
realmin,
Expand Down
8 changes: 7 additions & 1 deletion base/mpfr.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,11 +94,11 @@ end
BigFloat(x::Float32) = BigFloat(float64(x))
BigFloat(x::Rational) = BigFloat(num(x)) / BigFloat(den(x))

convert(::Type{Rational}, x::BigFloat) = convert(Rational{BigInt}, x)
convert(::Type{BigFloat}, x::Rational) = BigFloat(x) # to resolve ambiguity
convert(::Type{BigFloat}, x::Real) = BigFloat(x)
convert(::Type{FloatingPoint}, x::BigInt) = BigFloat(x)


for to in (Int8, Int16, Int32, Int64)
@eval begin
function convert(::Type{$to}, x::BigFloat)
Expand Down Expand Up @@ -138,6 +138,12 @@ promote_rule{T<:Real}(::Type{BigFloat}, ::Type{T}) = BigFloat
promote_rule{T<:FloatingPoint}(::Type{BigInt},::Type{T}) = BigFloat
promote_rule{T<:FloatingPoint}(::Type{BigFloat},::Type{T}) = BigFloat

promote_rule(::Type{Rational{BigInt}}, ::Type{BigFloat}) = Rational{BigInt}
promote_rule{T<:Integer}(::Type{Rational{T}}, ::Type{BigFloat}) = Rational{BigInt}
promote_rule{T<:FloatingPoint}(::Type{Rational{BigInt}}, ::Type{T}) = Rational{BigInt}

rationalize(x::BigFloat; tol::Real=eps(x)) = rationalize(BigInt, x, tol=tol)

# Basic arithmetic without promotion
# Unsigned addition
function +(x::BigFloat, c::Culong)
Expand Down
61 changes: 40 additions & 21 deletions base/rational.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,42 @@ end

convert{T<:Integer}(::Type{Rational{T}}, x::Rational) = Rational(convert(T,x.num),convert(T,x.den))
convert{T<:Integer}(::Type{Rational{T}}, x::Integer) = Rational(convert(T,x), convert(T,1))
function convert{T<:Integer}(::Type{Rational{T}}, x::FloatingPoint, tol::Real)

convert(::Type{Rational}, x::Rational) = x
convert(::Type{Rational}, x::Integer) = convert(Rational{typeof(x)},x)

convert(::Type{Bool}, x::Rational) = (x!=0) # to resolve ambiguity
convert{T<:Integer}(::Type{T}, x::Rational) = (isinteger(x) ? convert(T, x.num) : throw(InexactError()))
convert{T<:FloatingPoint}(::Type{T}, x::Rational) = convert(T,x.num)/convert(T,x.den)

function convert{T<:Integer}(::Type{Rational{T}}, x::FloatingPoint)
x == 0 && return zero(T)//one(T)
if !isfinite(x)
x < 0 && return -one(T)//zero(T)
x > 0 && return +one(T)//zero(T)
error(InexactError())
end
# TODO: handle values that can't be converted exactly
p = get_precision(x)-1
n = convert(T, significand(x)*2.0^p)
p -= exponent(x)
z = trailing_zeros(n)
p - z > 0 ? (n>>z)//(one(T)<<(p-z)) :
p > 0 ? (n>>p)//one(T) :
(n<<-p)//one(T)
end
convert(::Type{Rational}, x::Union(Float64,Float32)) = convert(Rational{Int}, x)

promote_rule{T<:Integer}(::Type{Rational{T}}, ::Type{T}) = Rational{T}
promote_rule{T<:Integer,S<:Integer}(::Type{Rational{T}}, ::Type{S}) = Rational{promote_type(T,S)}
promote_rule{T<:Integer,S<:Integer}(::Type{Rational{T}}, ::Type{Rational{S}}) = Rational{promote_type(T,S)}

promote_rule{T<:SmallSigned,S<:Union(Float32,Float64)}(::Type{Rational{T}}, ::Type{S}) = Rational{Int}
promote_rule{T<:SmallUnsigned,S<:Union(Float32,Float64)}(::Type{Rational{T}}, ::Type{S}) = Rational{Uint}
promote_rule{T<:Union(Int64,Int128),S<:Union(Float32,Float64)}(::Type{Rational{T}}, ::Type{S}) = Rational{T}
promote_rule{T<:Union(Uint64,Uint128),S<:Union(Float32,Float64)}(::Type{Rational{T}}, ::Type{S}) = Rational{T}

function rationalize{T<:Integer}(::Type{T}, x::FloatingPoint; tol::Real=eps(x))

This comment has been minimized.

Sign in to view

Copy link
@JeffBezanson

JeffBezanson Jun 11, 2013

Member

Maybe this could be a method of Rational.

This comment has been minimized.

Sign in to view

Copy link
@StefanKarpinski

StefanKarpinski Jun 11, 2013

Author Member

Perhaps, but it's unclear how it fits in. Now conversion from float to rational gives the exact rational value of a floating-point number. This is asking for something quite different: it finds the "simplest" rational, which when converted back to float is within a given tolerance of the original float. In particular, exact conversion from float to rational is not the same as rationalize with a tolerance of zero.

This comment has been minimized.

Sign in to view

Copy link
@JeffBezanson

JeffBezanson Jun 11, 2013

Member

I would expect setting the tolerance to zero to give something exactly equal if possible. Otherwise what does the tolerance mean? Is zero a special tolerance value that means "very small"?

This comment has been minimized.

Sign in to view

Copy link
@StefanKarpinski

StefanKarpinski Jun 11, 2013

Author Member

They're computing different things: rationalize(x, tol=0) gives a rational a//b such that a/b == x; convert(Rational, x) gives the exact rational value of the floating-point number x. These are not at all the same thing:

julia> rationalize(0.1, tol=0)
1//10

julia> convert(Rational, 0.1)
3602879701896397//36028797018963968

There are many rationals, a//b such that a/b == x, most of which are not the exact rational value represented by the floating-point value x.

This comment has been minimized.

Sign in to view

Copy link
@StefanKarpinski

StefanKarpinski Jun 11, 2013

Author Member

Perhaps rationalize needs a better name.

This comment has been minimized.

Sign in to view

Copy link
@StefanKarpinski

StefanKarpinski Jun 11, 2013

Author Member

It's also arguable that rationalize should have a default tolerance of zero, but using eps(x) seems to give more expected results.

This comment has been minimized.

Sign in to view

Copy link
@JeffBezanson

JeffBezanson Jun 11, 2013

Member

I would debate whether "compute a//b such that a/b==x" is a useful function. I think the definition should be "compute a//b such that abs(a//b - x)<=tol", in the case where a tolerance is specified.

This comment has been minimized.

Sign in to view

Copy link
@StefanKarpinski

StefanKarpinski Jun 11, 2013

Author Member

That's a fair point, but in that case you're really decomposing this into two operations:

  1. convert a float to a Rational with that exact value
  2. given a rational q, find "small" a and b such that abs(a//b-q) <= tol

Which is reasonable. I suspect that other systems such as Matlab don't express the problem this way because they don't have a first-class rational number type.

This comment has been minimized.

Sign in to view

Copy link
@StefanKarpinski

StefanKarpinski Jun 11, 2013

Author Member

I do, however, suspect that people fairly frequently want a//b such that a/b == x although maybe they shouldn't. In the 0.1 example, it's exactly what most people actually want because they probably want to recover the fact that 0.1 was given as 1/10.

This comment has been minimized.

Sign in to view

Copy link
@JeffBezanson

JeffBezanson Jun 13, 2013

Member

We have two numbers here, a//b and x. It seems reasonable to me that any tolerance should refer to the two numbers involved, not some third, non-present number a/b.

This comment has been minimized.

Sign in to view

Copy link
@StefanKarpinski

StefanKarpinski Jun 13, 2013

Author Member

Yes, that seems reasonable. I'll try to come up with something that works that way and does what people want when they want to find a rational that is close to some floating-point value. The issue is that checking abs(a//b - x) ≤ tol with the once and future promotion of (rational,float) to float is equivalent to checking that abs(a/b - x) ≤ tol. To check the actual difference (which is what a//b - x does at the moment), you have to write abs(a//b - convert(Rational,x)) ≤ tol instead.
if isnan(x); return zero(T)//zero(T); end
if x < typemin(T); return -one(T)//zero(T); end
if typemax(T) < x; return one(T)//zero(T); end
Expand All @@ -51,33 +86,17 @@ function convert{T<:Integer}(::Type{Rational{T}}, x::FloatingPoint, tol::Real)
while true
f = itrunc(y); y -= f
p, q = f*a+c, f*b+d
if p < typemin(T) || p > typemax(T) ||
q < typemin(T) || q > typemax(T)
break
end
typemin(T) <= p <= typemax(T) &&
typemin(T) <= q <= typemax(T) || break
a, b, c, d = p, q, a, b
if y == 0 || abs(a/b-x) <= tol
break
end
y = 1/y
y = inv(y)
end
return convert(T,a)//convert(T,b)
end
convert{T<:Integer}(rt::Type{Rational{T}}, x::FloatingPoint) = convert(rt,x,eps(one(x)))
convert(::Type{Rational}, x::FloatingPoint, tol::Real) = convert(Rational{Int},x,tol)
convert(::Type{Rational}, x::FloatingPoint) = convert(Rational{Int},x,eps(one(x)))
convert(::Type{Rational}, x::Integer) = convert(Rational{typeof(x)},x)
convert(::Type{Bool}, x::Rational) = (x!=0) # to resolve ambiguity

convert{T<:Rational}(::Type{T}, x::Rational) = x
convert{T<:Integer}(::Type{T}, x::Rational) =
(isinteger(x) ? convert(T, x.num) : throw(InexactError()))
convert{T<:Real}(::Type{T}, x::Rational) = convert(T,x.num)/convert(T,x.den)

promote_rule{T<:Integer}(::Type{Rational{T}}, ::Type{T}) = Rational{T}
promote_rule{T<:Integer,S<:Integer}(::Type{Rational{T}}, ::Type{S}) = Rational{promote_type(T,S)}
promote_rule{T<:Integer,S<:Integer}(::Type{Rational{T}}, ::Type{Rational{S}}) = Rational{promote_type(T,S)}
promote_rule{T<:Integer,S<:FloatingPoint}(::Type{Rational{T}}, ::Type{S}) = promote_type(T,S)
rationalize(x::Union(Float64,Float32); tol::Real=eps(x)) = rationalize(Int, x, tol=tol)

num(x::Integer) = x
den(x::Integer) = one(x)
Expand Down
41 changes: 22 additions & 19 deletions test/numbers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -333,10 +333,10 @@ end
@test !isequal(+1.0,-1.0)
@test !isequal(+Inf,-Inf)

@test isequal(-0.0f0, -0.0)
@test isequal(0.0f0, 0.0)
@test isequal(-0.0f0,-0.0)
@test isequal( 0.0f0, 0.0)
@test !isequal(-0.0f0, 0.0)
@test !isequal(0.0f0, -0.0)
@test !isequal(0.0f0 ,-0.0)

@test !isless(-Inf,-Inf)
@test isless(-Inf,-1.0)
Expand Down Expand Up @@ -414,8 +414,8 @@ end
@test isequal( 0.0, 0)
@test !isequal( 0,-0.0)
@test !isequal(-0.0, 0)
@test isless(-0.0, 0)
@test !isless( 0,-0.0)
@test isless(-0.0, 0)
@test !isless( 0,-0.0)

for x=-5:5, y=-5:5
@test (x==y)==(float64(x)==int64(y))
Expand Down Expand Up @@ -681,9 +681,12 @@ end

for a = -5:5, b = -5:5
if a == b == 0; continue; end
@test !ispow2(b) || a//b == a/b
if ispow2(b)
@test a//b == a/b
@test convert(Rational,a/b) == a//b
end
@test rationalize(a/b) == a//b
@test a//b == a//b
@test a//b == convert(Rational,a/b)
if b == 0
@test_fails integer(a//b) == integer(a/b)
else
Expand Down Expand Up @@ -1302,18 +1305,18 @@ approx_eq(a, b) = approx_eq(a, b, 1e-6)
@test int128(~0) == ~int128(0)

# issue 1552
@test isa(convert(Rational{Int8},pi),Rational{Int8})
@test convert(Rational{Int8},pi) == 22//7
@test convert(Rational{Int64},0.957762604052997) == 42499549//44373782
@test convert(Rational{Int16},0.929261477046077) == 11639//12525
@test convert(Rational{Int16},0.2264705884044309) == 77//340
@test convert(Rational{Int16},0.39999899264235683) == 2//5
@test convert(Rational{Int16},1.1264233500618559e-5) == 0//1
@test convert(Rational{Uint16},1.1264233500618559e-5) == 1//65535
@test convert(Rational{Uint16},0.6666652791223875) == 2//3
@test convert(Rational{Int8},0.9374813124660655) == 15//16
@test convert(Rational{Int8},0.003803032342443835) == 0//1
@test convert(Rational{Uint8},0.003803032342443835) == 1//255
@test isa(rationalize(Int8, pi), Rational{Int8})
@test rationalize(Int8, pi) == 22//7
@test rationalize(Int64, 0.957762604052997) == 42499549//44373782
@test rationalize(Int16, 0.929261477046077) == 11639//12525
@test rationalize(Int16, 0.2264705884044309) == 77//340
@test rationalize(Int16, 0.39999899264235683) == 2//5
@test rationalize(Int16, 1.1264233500618559e-5) == 0//1
@test rationalize(Uint16, 1.1264233500618559e-5) == 1//65535
@test rationalize(Uint16, 0.6666652791223875) == 2//3
@test rationalize(Int8, 0.9374813124660655) == 15//16
@test rationalize(Int8, 0.003803032342443835) == 0//1
@test rationalize(Uint8, 0.003803032342443835) == 1//255

# primes

Expand Down

0 comments on commit 6125749

Please sign in to comment.