Rationalize Float64 when creating intervals?

[dpsanders]

No description provided.

[dpsanders]

It may be a good idea to use rationalize when confronted with Float64. In the current version, we have:

julia> a = 1/3
0.3333333333333333

julia> @interval(a)
[3.3333333333333326e-01, 3.3333333333333331e-01] with 53 bits of precision

and the interval does not contain the true 1//3.

However,

julia> a = rationalize(1/3)
1//3

julia> @interval(a)
[3.3333333333333331e-01, 3.3333333333333337e-01] with 53 bits of precision

works correctly.

[dpsanders]

An alternative would just be to grow the interval slightly:

julia> a = 1/3
0.3333333333333333

julia> b, c = prevfloat(a), nextfloat(a)
(0.33333333333333326,0.33333333333333337)

julia> @interval(b, c)
[3.3333333333333326e-01, 3.3333333333333337e-01] with 53 bits of precision

which gives a slightly better answer, but maybe is more dangerous?

[lbenet]

I would be more in favor of the first option, though it also has some caveats:

julia> a = 1/factorial(5)
0.008333333333333333

julia> rationalize(a)  # this is fine
1//120

julia> convert(Rational{Int64},a)  # this doesn't, I don't understand why...
4803839602528529//576460752303423488

The reason I don't like the second option is that it may lead to ever growing intervals. If this is somehow controlled, in the sense that the interval obtained is thin in the extended sense (isthin yields true), the the option will indeed be a good candidate.

[dpsanders]

In convert(Rational{Int64},a), the denominator is just 2^59, i.e. it gives you the exact fraction with a power of 2 denominator that is equal to the binary representation of the float (or something like that).

rationalize seems to try to do something more clever -- maybe too clever... I don't understand the code

[dpsanders]

Somehow it should be trying to "guess" the rational by finding the repeating string in the binary expansion.

[dpsanders]

I don't see why you say that my option 2 (widen the interval slightly) will lead to ever-growing intervals. This is just when the interval is first created from a float (which, ideally, will never happen, so maybe we should just throw an error if this happens!)

[lbenet]

Regarding rationalize and convert, I'm confused, since convert calls rationalize (https://github.com/JuliaLang/julia/blob/master/base/rational.jl#L74) with zero tolerance...

Regarding option 2, the question is if it yields a thin interval or not, in the extended sense. If it doesn't, I fear that the "extra diameter" will propagate further due to interval arithmetics

[dpsanders]

Oh, apologies, you're right. But it uses tol=0 instead of the default tol=eps(x) used in rationalize.
So it gives the exact rational instead of a slightly approximate one.

I understand your worry about the "extra diameter". The problem is that we are trying to infer what the user wanted, instead of what they wrote... (Just to record the relevant discussion: JuliaLang/julia#4278)

[dpsanders]

Related (EDIT: really the same) problem:

julia> a = big(.1)
1.0000000000000001e-01 with 53 bits of precision

julia> @interval(a)
[1.0000000000000001e-01, 1.0000000000000001e-01] with 53 bits of precision

[lbenet]

Can it be ensured that @interval(prevfloat(a), nextfloat(a)) with a a Float64 is indeed a thin interval? Maybe, but still one has to be cautious:

julia> Interval(prevfloat(1/3),nextfloat(1/3))   # answer is not OK
[0.33333333333333326, 0.33333333333333337]

julia> ValidatedNumerics.isthin(ans)
false

The second case is what I fear.

[lbenet]

I edited my previous comment....

julia> @interval(prevfloat(1/3), nextfloat(1/3))
ERROR: MethodError: `prevfloat` has no method matching prevfloat(::ValidatedNumerics.Interval{Base.MPFR.BigFloat})

[dpsanders]

You were right to use Interval instead of @interval in that case!

No, it cannot be guaranteed to be a thin interval. The reason is that the float is already rounded, either down or up, so it is already one of the endpoints of the true thin interval. I do not see any way around this except to round it again.

[lbenet]

In a way, this discussion is parallel to the proposal for @bigfloat at JuliaLang/julia#4278 (comment)

[lbenet]

Maybe we should merge #11 and document the caveats, and carry on...

[dpsanders]

How about:

transf(a::Float64)   =  transf(rationalize(a))
transf(a::BigFloat)  =  transf(convert(Float64, a))  # surprising -- reduce the precision?

Note that if it gets a BigFloat, it first converts it Float64 in order to use rationalize...!

[dpsanders]

Interestingly, in the context of interval arithmetic, we are in some sense safer than the discussion referred to above, since we just need to guarantee that the resulting interval really contains the original Float64 or BigFloat (tests required...).

[lbenet]

Have you tested your proposal about transf(a::Float64)? Have you tested it?

Regarding the second one, I remember that the idea of implementing it as it is, was to avoid re-rounding a, again to avoid extra-diameter of the resulting interval. So, I am skeptic/suspicious...

[lbenet]

I was checking other parts of the code, and just realized that @round requires the type to be rounded. Is the use of @round of any help in this context?

[dpsanders]

I was checking other parts of the code, and just realized that @round
requires the type to be rounded. Is the use of @round of any help in this
context?

I'm not sure what you mean by "requires the type to be rounded". What did
you have in mind?

[lbenet]

It is used as @round(T,a_low,a_hi), where T is the type of the interval; see e.g. https://github.com/dpsanders/ValidatedNumerics.jl/blob/master/src/intervals/interval_arithmetic.jl#L39, but almost any function there also works like that.

[dpsanders]

It is used as @round(T,a_low,a_hi), where T is the type of the interval;
see e.g.
https://github.com/dpsanders/ValidatedNumerics.jl/blob/master/src/intervals/interval_arithmetic.jl#L39,
but almost any function there also works like that.

But if you feed in floating point numbers they will not be changed by
@round.

[lbenet]

You are right...

[lbenet]

I was not aware of JuliaLang/julia#8845. Is this helpful?

[dpsanders]

It's basically just nicer syntax (and perhaps slightly more performant). It doesn't change the mathematics.
Once 0.4 comes out there will be a lot of things we can rewrite!

[dpsanders]

rationalize is now used.