compactify_rational inaccurate
Closed this issue · 10 comments
Due to the implementation using Float64 arithmetic, the function does not what it promises in corner cases.
julia> compactify_rational(FastRational(2^52//(2^52-1)), FastRational(1//(2^52)))
1//1
julia>
The return value is outside the given tolerance interval and the midpoint is a proper FastRational{Int}. ... same for Rational.
It would be beneficial to provide an implementation based on integer or rational arithmetic. Would you accept a PR which replaces compactify_rational?
it is not intended for use with that fine an interval -- we need to do something though, so we are both comfortable.
The algorithm I have (and cite) is the only one of which I am aware that delivers the least denominator (and corresponding or least of a few numerators) given the constraints. Are you aware of a way to do this using integers or rationals (although I imagine the floating point algorithm would be faster than a rational version, while an integer version may well be the fastest)?
I think, I have an integer version in my memory. Still missing to handle integer overflows, like you exercised in fast rational arithmetic.
We could use Base.rationalize when the process would overflow, as it is designed to be robust that way.
If you are manipulating integers, knowing before the fact if an overflow is possible (not if it is assured, that is more complex) is a very fast test [particularly if we force the integers to be strictly positive, by special casing zero, and twice reflecting negative values].
if you have a way forward -- PR it
the userfacing function is now named compactify
(it runs very quickly relative to rationalize, and is not intended for use with the narrowest of radii)
For my own edification, what is the role/characterization of au, bu and v, and
what meaning, mathematical semantics do each of the triples lo3, hi3 evince?
# compute Q(midpoint - rad) and Q(midpoint + rad) approximately
# preserving integer type
function compact_reduce(midpoint::Q, radius::Q) where {Q<:Rationals}
midnum, midden = midpoint.num, midpoint.den
radnum, radden = radius.num, radius.den
u = midpoint ÷ max(midden, midnum)
au = midnum * u
bu = midden * u
v = (bu ÷ radden) * radnum
lo3 = compact_integer(au - v, bu, au - v + 1, bu)
hi3 = compact_integer(au + v - 1, bu, au + v, bu)
tuple(lo3..., hi3...)
end
The whole function is to calculate approximations to midpoint - rad and midpoint + rad without using widened arithmetic operations.
The line with u = is wrong and should read u = typemax(basetype(Q)) ÷ max(midden, midnum).
That gives a multiplier u, which brings au or bu close to the limits of the base type. Given v as defined yields (au - v) // bu <= midpoint - rad < (au - v + 1) // bu, where the outer bounds have the same type as midpoint, while the exact midpoint - rad typically requires widened integer size.
Unfortunately we lose exactness, if we don't use widened integers for the mid - rad arguments. Therefore I preferred the low/high approach with normal integers.
Unfortunately we lose exactness, if we don't use widened integers for the mid - rad arguments. Therefore I preferred the low/high approach with normal integers.
I understand that preference. There is an important use case where loosing exactness is precisely what is needed -- keeping rationals within their original IntNN determined domain through a sequence of arithmetic operations that otherwise would overflow the type. I have been looking for a good way to intersperse that sort of inexact yet determinably close reduction without relying on it after each op, and with logic that reduces effectively given a problem's character. For example, exploring constrained tilings of bounded, regularly subdivided spaces [e.g. a checkerboard].