JuliaIntervals / IntervalArithmetic.jl

Library for validated numerics using interval arithmetic
https://juliaintervals.github.io/IntervalArithmetic.jl/
Other
299 stars 71 forks source link

Document `PowerMode` #674

Open dpsanders opened 3 months ago

dpsanders commented 3 months ago

Currently to compute (1..2)^(1//7) we are converting the rational to an interval. This loses accuracy.

We should instead use the corresponding MPFR function, mpfr_rootn_si. Unfortunately Julia does not seem to do this yet; issue here.

For rationals whose denominator is a power of 2, we can alternatively use iterated sqrt to get a fast version that avoids BigFloats (but may not be correctly rounded).

As an example, tests in IntervalContractors.jl (which are originally based on ITF1788, IIRC) are failing due to this. For example,

interval(8.673020346900622e8, 8.673020346900623e8)^(1//8)

has a sup of 13.100000000000005, whereas it should be 13.100000000000001 (i.e. the result is 2 floats too large).

dpsanders commented 3 months ago

Possibly we need ExactReal(1//8) for the rational power.

dpsanders commented 3 months ago

This was implemented in a previous version. See https://github.com/JuliaIntervals/IntervalArithmetic.jl/pull/286/files

OlivierHnt commented 3 months ago

Yep it's still in the library I believe. I think you just ran into an issue with the "power mode"

The default is IntervalArithmetic.PowerMode{:fast}(), and this is controlled by overwriting IntervalArithmetic.power_mode.

OlivierHnt commented 3 months ago

Unless you think we can improve the performance of the more accurate version, then I think we can close this issue.

dpsanders commented 3 months ago

Ah, I see 🤦 Thank you! Is this documented?

OlivierHnt commented 3 months ago

No it does not seem to be in the docs.

dpsanders commented 3 months ago

I don't see any mechanism to get the current power mode?

I need this so that I can save the current power mode before fixing it to run the IntervalContractor.jl tests, then reset it to the previous value.

dpsanders commented 3 months ago

Usage:

julia> setdisplay(:full);

julia> interval(8.673020346900622e8, 8.673020346900623e8)^(1//8)
Interval{Float64}(13.099999999999998, 13.100000000000005, com, NG)

julia> IntervalArithmetic.power_mode() = IntervalArithmetic.PowerMode{:slow}()

julia> interval(8.673020346900622e8, 8.673020346900623e8)^(1//8)
Interval{Float64}(13.099999999999998, 13.100000000000001, com, NG)
dpsanders commented 3 months ago

Ah, to get the current value you just do

julia> IntervalArithmetic.power_mode()
IntervalArithmetic.PowerMode{:slow}()

!

dpsanders commented 3 months ago

When trying to change the power mode in the IntervalContractors.jl tests and then running ]test IntervalContractors, I get the warning

WARNING: Method definition power_mode() in module Main at /Users/dsanders/Dropbox/packages/IntervalContractors/test/runtests.jl:9 overwritten at /Users/dsanders/Dropbox/packages/IntervalContractors/test/runtests.jl:40.

Is this expected / normal?

OlivierHnt commented 3 months ago

Yes, that's just a warning emitted by Julia.