Flag for enabling warnings when converting from AbstractFloat

[Joel-Dahne]

I find that one of the easiest mistakes to make when working with Arblib in Julia is to accidentally use Float64 for some internal computations. For example for x::Arb it's perfectly fine to do any of the following

2x
x / 2
x^2
π * x
x / π
x^π

But you get a wrong result for, e.g.,

2π * x
x / log(2)
x^(1 / 3) # x^(1 // 3) works though
sqrt(2) * x

One of the common things for all of the wrong examples is that they involve a conversion from AbstractFloat to Arb. The Idea I have had is to have some sort of flag or macro which, when enable, prints a warning whenever a AbstractFloat is converted to Arb. So it would behave something like

function convert(::Type{Arb}, x::AbstractFloat)
    @warn "converting a $(typeof(x)) to Arb, it's possible that this is a mistake."
    return Arb(x)
end

Does this sound interesting to you?

I don't know exactly how this would be implemented, I feel like it should be possible but I don't know exactly how. Any thoughts?

[kalmarek]

I actually like the idea! how about just using @debug macro? This is silent by definition unless an ENV flag is set

Example:

function Base.convert(::Type{Arb}, x::Union{Float16, Float32, Float64})
   @debug "Converting a float to Arb. This may result in loss of precision!"
   return Arb(x)
end

and in REPL:

julia> 2.6*Arb(1)
2.6000000000000000888178419700125232338905334472656250000000000000000000000000

julia> ENV["JULIA_DEBUG"]="Arblib"
"Arblib"

julia> 2.6*Arb(1)
┌ Debug: Converting a float to Arb. This may result in loss of precision!!
└ @ Arblib ~/.julia/dev/Arblib/src/constructors.jl:19
2.6000000000000000888178419700125232338905334472656250000000000000000000000000

Of course this will not save the user from doing Arb(2.6) directly, since the "proper" place to have this @debug statement would be in arbcall setter... Maybe we could actually move it out of arbcalls and write it by hand?

[Joel-Dahne]

A @debug macro sounds like a good solution! There's a couple of options for where this could be placed

• At convert(Arb, x::Float64)
• At Arb(x::Float64)
• At Arblib.set!(res, x::Float64)
  Only the set! method would guarantee that we catch all conversions from Float64. Most mistakes would probably be in the convert method, for example all of the examples I gave above would go under this. In some cases converting a Float64 to Arb is fine, for example when defining some arbitrary cutoff value or something like that, and hence updating the constructor or the setter could lead to a lot of false positives. But you could definitely come up with examples where someone uses the constructor explicitly but it's still a mistake, I often use Arb((0, π)) but for example Arb((0, 2π)) would be a mistake.

My current feeling is that we should update the setter. It's better to give the user false positives than to miss some true positives, in particular when working with computer assisted proofs.

Another question is which types to warn for. For sure Float16, Float32 and Float64 and I think BigFloat as well. However for Arf I'm not so sure, you can definitely make the same mistakes with Arf but on the other hand it is often used together with Arb to represent lower and upper bounds and for that reason I often convert between the two. You are also less likely to accidentally use Arf in a computation simply because it doesn't support a lot of operations. So I would go for emitting a warning for the setters from Float64 and BigFloat but not Arf.

I'll experiment a bit with this locally and see what I like!

[kalmarek]

Setter for Float64s is defined in src/arbcalls so we'd need to comment it out and add to setters.
I think the better and safest solution is to just overload Arb(x::Any) to provide fallback with @debug statement and provide another method Arb(x::Union{TypesKnownToBeExact}) without the debug. This way we're safe and the standard dispatch mechanism will allow the user to get rid of the "warnings" should he/she use a type known to be exact.

On the other hand I sometimes purposely create arbs from whole matrices of Float64 and then manipulate them to satisfy my assumptions (e.g. rows and columns satisfy given set of linear equations). That's a different use case which will probably not benefit from this approach.

[kalmarek]

I just noticed quite severe regression:

julia> ENV["JULIA_DEBUG"]=""
""

julia> @btime 1.0*$(Arb(2.6)) # @debug statement in `convert`
  108.709 ns (2 allocations: 128 bytes)
2.6000000000000000888178419700125232338905334472656250000000000000000000000000

julia> @btime 1.0*$(Arb(2.6)) # no @debug statement in `convert`
  39.899 ns (2 allocations: 128 bytes)
2.6000000000000000888178419700125232338905334472656250000000000000000000000000

which isn't very surprising: I think @debug statement needs to lookup a global constant...
so another solution is to introduce

const _debug = Ref(false)
function Base.convert(::Type{Arb}, x::Union{Float16, Float32, Float64})
    _debug[] && @warn "Converting a float to Arb. This may result in loss of precision!"
    return Arb(x)
end

which is more digestible

julia> @btime 1.0*$(Arb(2.6))
  48.144 ns (2 allocations: 128 bytes)
2.6000000000000000888178419700125232338905334472656250000000000000000000000000

[Joel-Dahne]

On the other hand I sometimes purposely create arbs from whole matrices of Float64 and then manipulate them to satisfy my assumptions (e.g. rows and columns satisfy given set of linear equations). That's a different use case which will probably not benefit from this approach.

I agree! For me it is also very common to construct some sort of approximation from Float64 and the compute rigorous error bounds on that. But usually my program is split into two parts, one which computes the approximation and constructs the Arb object and another one which computes the bound. In this case you could run the latter part of the program with debug logging enabled, allowing you to catch any errors.

[Joel-Dahne]

The performance regression is worrying, I didn't expect @debug to be so costly. I would prefer a solution which is zero cost if debugging is not enabled. Maybe one option could be to use a flag read at init (similar to what we do for the Flint multithreading) which enables or disables the debugging? That we it would truly be zero cost if you don't set this flag and you can pay some performance if you want to enable it.

[kalmarek]

the solution with const _debug = Ref(false) is actually less costly than reported above. I changed the constructor as suggested above like this:

const _debug = Ref(true)

#Arb
function Arb(x::Any; prec::Integer = _precision(x))
    _debug[] && @warn "Converting $(typeof(x)) to Arb. This may result in loss of precision!"
    return set!(Arb(prec = prec), x)
end
function Arb(x::T; prec::Integer = _precision(x)) where T <: Union{Rational, Integer, ArbLike, ArfLike}
    return set!(Arb(prec = prec), x)
end

julia> 2.6*Arb(2.6)
[6.7600000000000004618527782440651287048398261358993041172856528278622967320644 +/- 4.90e-77]

julia> Arblib._debug[] = true
true

julia> 2.6*Arb(2.6)
┌ Warning: Converting Float64 to Arb. This may result in loss of precision!
└ @ Arblib ~/.julia/dev/Arblib/src/constructors.jl:16
┌ Warning: Converting Float64 to Arb. This may result in loss of precision!
└ @ Arblib ~/.julia/dev/Arblib/src/constructors.jl:16
[6.7600000000000004618527782440651287048398261358993041172856528278622967320644 +/- 4.90e-77]

julia> Arblib._debug[] = false
false

julia> @btime 2.6*$(Arb(2.6))
  41.764 ns (2 allocations: 128 bytes)
[6.7600000000000004618527782440651287048398261358993041172856528278622967320644 +/- 4.90e-77]

without the @warn line:

julia> Arblib._debug[] = true
true

julia> 2.6*Arb(2.6)
[6.7600000000000004618527782440651287048398261358993041172856528278622967320644 +/- 4.90e-77]

julia> Arblib._debug[] = true
true

julia> @btime 2.6*$(Arb(2.6))
  40.423 ns (2 allocations: 128 bytes)
[6.7600000000000004618527782440651287048398261358993041172856528278622967320644 +/- 4.90e-77]

julia> 2.6*Arb(2.6)
[6.7600000000000004618527782440651287048398261358993041172856528278622967320644 +/- 4.90e-77]

julia> @btime 2.6*$(Arb(2.6))
  39.929 ns (2 allocations: 128 bytes)
[6.7600000000000004618527782440651287048398261358993041172856528278622967320644 +/- 4.90e-77]

so maybe we can pay this ~1-2ns price? Especially as branch prediction should do a pretty good job here :)

[Joel-Dahne]

Okay, lets start with this since it's the simpler solution! I'll do some local testing and try it out for a while before pushing anything probably.

I agree that this should be a trivial job for the branch predictor!

[devmotion]

An alternative with zero overhead would be to a function instead of a const Ref:

julia> using Arblib, BenchmarkTools

julia> enable_debug() = false
enable_debug (generic function with 1 method)

julia> function MyArb(x; prec::Integer = Arblib._precision(x))
           enable_debug() && @warn "Converting $(typeof(x)) to Arb. This may result in loss of precision!"
           return Arblib.set!(Arb(prec = prec), x)
       end
MyArb (generic function with 1 method)

julia> x = 2.6;

julia> y = Arb(x);

julia> @btime $x * $y;
  56.531 ns (2 allocations: 128 bytes)

julia> @btime Arb($x) * $y;
  56.727 ns (2 allocations: 128 bytes)

julia> @btime MyArb($x) * $y;
  56.654 ns (2 allocations: 128 bytes)

julia> enable_debug() = true # causes all dependent methods to be recompiled
enable_debug (generic function with 1 method)

julia> MyArb(x) * y
┌ Warning: Converting Float64 to Arb. This may result in loss of precision!
└ @ Main REPL[27]:2
[6.7600000000000004618527782440651287048398261358993041172856528278622967320644 +/- 4.90e-77]

This is approach is e.g. also used by TimerOutputs (https://github.com/KristofferC/TimerOutputs.jl#overhead). In comparison, the const approach has some overhead:

julia> const _debug = Ref(false);

julia> function MyArb(x; prec::Integer = Arblib._precision(x))
           _debug[] && @warn "Converting $(typeof(x)) to Arb. This may result in loss of precision!"
           return Arblib.set!(Arb(prec = prec), x)
       end
MyArb (generic function with 1 method)

julia> @btime MyArb($x) * $y;
  59.130 ns (2 allocations: 128 bytes)

Not sure how important the additional performance is for you 🙂 I can make a PR with an enable_debug() switch if you prefer this approach (in the other case as well but maybe you are already working on it?).

[Joel-Dahne]

That looks like a great idea! If I understand it correctly it would mean that you pay the cost of recompilation when switching the flag but otherwise it should be zero overheard since the check would be optimized out?

If you could make a PR that would be amazing! The reason that it has not been moving forward for a while is because some other parts of my current project took longer than expected so I have not yet gotten to checking all the code, but eventually I'll have to do this. The main question that remains to answer is exactly what methods to add this check to. This is likely problem dependent but I think adding it to all Float64 methods in Arb would make sense, so all C functions ending with _d. But as long as we have a god way to toggle this debugging it would be easy to add or remove the check to functions at a later point.

[devmotion]

If I understand it correctly it would mean that you pay the cost of recompilation when switching the flag but otherwise it should be zero overheard since the check would be optimized out?

Yes, exactly. A simple example that illustrates how it is optimized and recompiled by the compiler:

julia> enable_debug() = true
enable_debug (generic function with 1 method)

julia> f() = enable_debug() ? 1 : 1.0
f (generic function with 1 method)

julia> @code_llvm f()
;  @ REPL[2]:1 within `f'
define i64 @julia_f_256() {
top:
  ret i64 1
}

julia> enable_debug() = false
enable_debug (generic function with 1 method)

julia> @code_llvm f()
;  @ REPL[2]:1 within `f'
define double @julia_f_277() {
top:
  ret double 1.000000e+00
}

[devmotion]

The same with const debug:

julia> const debug = Ref(true)
Base.RefValue{Bool}(true)

julia> g() = debug[] ? 1 : 1.0
g (generic function with 1 method)

julia> @code_llvm g()
;  @ REPL[7]:1 within `g'
define { {}*, i8 } @julia_g_290([8 x i8]* noalias nocapture align 8 dereferenceable(8) %0) {
top:
; ┌ @ refvalue.jl:56 within `getindex'
; │┌ @ Base.jl:33 within `getproperty'
    %1 = load i8, i8* inttoptr (i64 139923040448560 to i8*), align 16
    %2 = and i8 %1, 1
    %.not = icmp eq i8 %2, 0
; └└
  %spec.select = select i1 %.not, { {}*, i8 } { {}* inttoptr (i64 139923021809568 to {}*), i8 -127 }, { {}*, i8 } { {}* inttoptr (i64 139923014266976 to {}*), i8 -126 }
  ret { {}*, i8 } %spec.select
}

julia> @code_warntype g()
Variables
  #self#::Core.Const(g)

Body::Union{Float64, Int64}
1 ─ %1 = Base.getindex(Main.debug)::Bool
└──      goto #3 if not %1
2 ─      return 1
3 ─      return 1.0

julia> debug[] = false
false

julia> @code_llvm g()
;  @ REPL[7]:1 within `g'
define { {}*, i8 } @julia_g_427([8 x i8]* noalias nocapture align 8 dereferenceable(8) %0) {
top:
; ┌ @ refvalue.jl:56 within `getindex'
; │┌ @ Base.jl:33 within `getproperty'
    %1 = load i8, i8* inttoptr (i64 139923040448560 to i8*), align 16
    %2 = and i8 %1, 1
    %.not = icmp eq i8 %2, 0
; └└
  %spec.select = select i1 %.not, { {}*, i8 } { {}* inttoptr (i64 139923021809568 to {}*), i8 -127 }, { {}*, i8 } { {}* inttoptr (i64 139923014266976 to {}*), i8 -126 }
  ret { {}*, i8 } %spec.select
}

julia> @code_warntype g()
Variables
  #self#::Core.Const(g)

Body::Union{Float64, Int64}
1 ─ %1 = Base.getindex(Main.debug)::Bool
└──      goto #3 if not %1
2 ─      return 1
3 ─      return 1.0

[Joel-Dahne]

Great! Then we get truly zero overhead when the debugging is disabled!