RFC: basic floating point exception functionality

[simonbyrne]

This adds basic functionality for checking and clearing floating point exception flags (#2976). Usage is currently:

using Base.Rounding
get_floatexcept(Float64)
is_floatexcept(Float64,FEUnderflow)

Still to do:

• Settle on a consistent interface
• Incorporate BigFloat flags: I've now added these, but the exceptions don't exactly match the IEEE ones
• A way to store and restore flags
• At the moment, almost every operation seems to flag FEInexact (input, printing, etc.), even for exact numbers (e.g. x=1.0): I'm not sure what the correct behaviour is here.

cc: @andrioni

[andrioni]

I'll give it a try. About the FEInexact, it might be because of random internal operations or even grisu, it's pretty hard to avoid it.

[andrioni]

function test2()
       clear_floatexcept()
       before = get_floatexcept()
       ret = 1.0 / 7
       after = get_floatexcept()
       return ret, before, after
end

This returns no exceptions on both cases, which is very weird. I'm not sure if LLVM is reordering the operations, or anything else is happening. If anyone is good with LLVM IR, here goes it:

define { double, %FloatExceptionSet, %FloatExceptionSet } @julia_test215504() {
top:
  call void @julia_clear_floatexcept15323(%FloatExceptionSet { i32 61 }), !dbg !1356
  %0 = call i32 inttoptr (i64 140646004314944 to i32 (i32)*)(i32 61), !dbg !1357
  %1 = insertvalue %FloatExceptionSet undef, i32 %0, 0, !dbg !1357, !julia_type !1358
  %2 = call i32 inttoptr (i64 140646004314944 to i32 (i32)*)(i32 61), !dbg !1359
  %3 = insertvalue %FloatExceptionSet undef, i32 %2, 0, !dbg !1359, !julia_type !1358
  %4 = insertvalue { double, %FloatExceptionSet, %FloatExceptionSet } { double 0x3FC2492492492492, %FloatExceptionSet undef, %FloatExceptionSet undef }, %FloatExceptionSet %1, 1, !dbg !1360, !julia_type !1361
  %5 = insertvalue { double, %FloatExceptionSet, %FloatExceptionSet } %4, %FloatExceptionSet %3, 2, !dbg !1360, !julia_type !1361
  ret { double, %FloatExceptionSet, %FloatExceptionSet } %5, !dbg !1360
}

[andrioni]

Oh, I guess LLVM is constant-folding, so it doesn't raise the exception at runtime.

[andrioni]

Yeah, indeed that was the case:

function test2(x)
       clear_floatexcept()
       before = get_floatexcept()
       ret = 1.0 / x
       after = get_floatexcept()
       return ret, before, after
end

julia> test2(2)
(0.5,{},{})

julia> test2(7)
(0.14285714285714285,{},{FEInexact})

[simonbyrne]

Okay, I've added functions to raise floating point flags.

The C functions to save and restore flags are fegetexceptflag and fesetexceptflag, respectively, however they use an fexcept_t object, which according to the standard is "implementation defined", so we need some way to get its typedef from the header, similar to how the FE_* constants are obtained. Any suggestions for how to do this?

[andrioni]

In the worst case, you can just grab a sizeof(fexcept_t) and allocate a buffer big enough for it.

[simonbyrne]

@andrioni That does seem like a reasonable approach, but I'm guessing that can't be done using the preprocessor as in fenv_constants.h. Unfortunately, I know very little C, so if you know of a good way to implement it, I'd be grateful.

[vtjnash]

You can stick a function in src/sys.c which returns the required information. You should find several examples of this if you search for _sizeof_

Or you could use the implementations of these functions from openlibm, which vary have various definitions of fexcept_t from uint16 to uint16, and go with the largest type (sparc64) to cover all of them.

[simonbyrne]

@vtjnash Fantastic, thanks. I now have the save/restore working (for IEEE types).

[ViralBShah]

About grisu and FEInexact, we could save the state before entry and restore it upon exit? Is it worthwhile to even attempt this?

I think the interface is fine to start with.

[simonbyrne]

I've been thinking a little more about this, and it would be much more useful if this could be incorporated into the general julia exception handling machinery: the IEEE standard covers this in section 8.3.

Delayed exception handling would be easy to implement: we could define a method like:

function watch_floatexcept(f,T,fe)
    state = save_floatexcept(T,fe)
    clear_floatexcept(T,fe)
    f()
    is_floatexcept(T,fe) && throw(FloatError(fe))
    restore_floatexcept(T,state,fe)
end

which could then be used:

try
    watch_floatexcept(Float64,FEUnderflow) do
        # do something
    end
catch e
    # handle exception
end

(this could also be done with a macro).

Immediate handling would be much more difficult: we would need something similar to how DomainErrors are handled for functions like sin. I don't understand the internals of julia sufficiently well to determine how this might be applied more widely.

[simonbyrne]

A possible option is to trap SIGFPE, though this is well beyond my C knowledge.

[ihnorton]

Have a look at the existing fpe handler (throws DivisionError)

[simonbyrne]

@ihnorton Thanks: so does that mean we're already trapping it? That should make it easier to implement, though again, beyond my abilities.

[ViralBShah]

Is this something we should include in 0.3? It is independent new functionality.

[quinnj]

I'm not a wizard on floating point exceptions, but if I can help on grisu issues, just let me know.

[dpsanders]

What is the status on this -- would be great to have the clear_floatexcept() and get_floatexcept() functionality!

[dpsanders]

I have extracted the relevant parts of the first commit into a module:

https://gist.github.com/dpsanders/321546d749370403e0bc

@simonbyrne Please let me know if there is any problem with this, of course. [And could you please check that I got the value of JL_FE_INEXACT right (0x20 = 32) ? It seems to be correct.]

Maybe some version of this could be registered as a package if it is not going to be included in Base.

EDIT: I tried to check out your branch (to a new directory) and compile, but got some errors.

[simonbyrne]

I don't know. My ideal preference would be to make use of floating point exception masking/unmasking, as I outlined here: #5234 (comment), but I realise that this is unlikely to happen soon.

This interface is a bit unwieldy, but it does have the advantage of at least being straightforward to implement.

I'm not sure about JL_FE_INEXACT: technically this is system dependent. If we're not going to put this PR in Base, we should at least export this variable.

[ViralBShah]

Should we close this as stale?