RFC: Ability to change rounding mode for all floats (cf. #2976)

[andrioni]

This adds three new functions (get_rounding, set_rounding and with_rounding) which allow users to change the rounding direction of floating point arithmetic (both Float32 and Float64). I've also changed the behavior of the old BigFloat rounding functions to agree with this new one, which I find much more intuitive and safe (see tests/rounding.jl for examples in the case of Float32/Float64).

I haven't tested this yet on Windows, and it probably shouldn't work, as it depends on the fenv.h header which is only available in C99 or POSIX systems, but I don't know how mingw works..

[vtjnash]

It appears that this should also work on windows. fenv.h exists in mingw and declares the appropriate functions and constants

[ViralBShah]

The hard coded enums are a bit scary. Is there a way to make it a bit more portable?

[ViralBShah]

BTW This is really cool stuff!

[andrioni]

Portable in what sense? The values used in fenv.h really are arbitrary (it was quite unexpected, and I had to go check the headers to find them, because they aren't documented anywhere), but are the same on OS X and Linux, at least. The way I used types is also a bit cleaner than the older approach in BigFloat, with hard-coded constant integers, and it allows a bit of finer control (as it eliminates the chance of passing arbitrary values to the function in Julia), but I'm completely open to better suggestions, as I'm also not completely comfortable right now.

[ViralBShah]

What I meant is, just like the various system specific constants that are generated in base/Makefile, can we do something similar here too, so that we can avoid using 1024, 2048, etc.?

[andrioni]

Oh, this is great, I had no idea it could be done this way.
Thanks, @vtjnash, I removed the @unix_only macro, as it hopefully should work on Windows then.

[andrioni]

Ok, I'll boot up an Ubuntu virtual machine later to see what is happening, it works fine on both OS X 10.8.3 and Arch Linux.

[bkalpert]

Very exciting! Ignorant question: Do we know that LLVM will not reorder code to nullify the effect of changing rounding mode (as per concern earlier expressed @andrioni)?

[JeffBezanson]

I believe it will work if the rounding mode is set by a called function, and not inlined.

[StefanKarpinski]

To be safe, does LLVM support some kind of instruction that acts as a barrier to reordering?

[JeffBezanson]

I don't believe so, unless you count function call. There are many LLVM list threads on this topic.

[andrioni]

I'm also a bit surprised that this did work, and it has in all my tests, so I guess the indirection provided by Julia function and/or the ccall is enough to prevent reordering. As using fenv.h directly works and Visual C++ compatibility isn't needed, this also means we'll have all floating point exception flags easily available in Julia as soon as I think of a good API.

[JeffBezanson]

Yes, it is well known that indirecting through a function call is sufficient to prevent reordering.

[bkalpert]

Minor point: Should we preserve precedent and call one mode RoundTowardZero (rather than RoundToZero)?

[StefanKarpinski]

Where's the precedent from?

[bkalpert]

My initial impression was created from fenv.h (FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO) but that is not definitive. Wikipedia on IEEE 754-1985 shows the rounding modes as:
Round to Nearest – rounds to the nearest value; if the number falls midway it is rounded to the nearest value with an even (zero) least significant bit, which occurs 50% of the time (in IEEE 754-2008 this mode is called roundTiesToEven to distinguish it from another round-to-nearest mode)
Round toward 0 – directed rounding towards zero
Round toward +∞ – directed rounding towards positive infinity
Round toward −∞ – directed rounding towards negative infinity.
Unfortunately I do not have a copy of the IEEE 754-1985 (or 2008) standard itself.

[StefanKarpinski]

I would even be in favor of shorter names: RoundUp, RoundDown, RoundNearest, RoundToZero, RoundFromZero.

[andrioni]

The names I used came from gmpy2, which in turn are inspired by those from MPFR.

[bkalpert]

The precedent is clearly varied, and I am content with either Alessandro or Stefan's choice. The capability is key!

[ViralBShah]

I prefer @StefanKarpinski 's suggested names. It doesn't seem there is a widely used set of names. Also, this pull request needs to be documented in the manual as well in the standard library before merging, or else nobody will ever know that this amazing functionality exists.

[ViralBShah]

Bump. @andrioni Can you please rebase?

I would prefer to get this merged in and bikeshed a bit later.

[ViralBShah]

@loladiro @vtjnash What about Windows?

[ViralBShah]

I believe that openlibm includes fenv.h and that should work on Windows.

[ViralBShah]

Bump again. It would be nice to get this included in the 0.2 release if possible.

[andrioni]

OK, I'm back. I'll rebase it in some hours!

[Keno]

Note that windows does not provide fsetround, you'll have to use http:https://msdn.microsoft.com/en-gb/library/e9b52ceh.aspx instead.

[andrioni]

I thought mingw would provide these functions, according to @vtjnash.

[Keno]

I thought they were macros, but I might be wrong, @vtjnash?

[vtjnash]

They are declared as functions in mingw's fenv.h

[andrioni]

Ok, sorry for the long delay, I had some bad health issues.

I'm quite out of date (a month or three), but I'll dedicate this weekend to put things on schedule again. Is this still considered for 0.2? If so, what is lacking? @ViralBShah

[andrioni]

There we go, I think everything is alright now.

[ViralBShah]

@andrioni Sorry to hear about your health, but it is great to see you back.

[ViralBShah]

@vtjnash Do you think we can merge this now, and deal with Windows issues post merging?

[vtjnash]

I'll check out whether this works on windows tonight and let you know.

[vtjnash]

this appears to be lacking documentation. per steven's excellent suggestion, it must be added before it can be merged.

[vtjnash]

this file shouldn't be part of the commit

[andrioni]

Hm, I thought I removed it, but now I see I just added the file to .gitignore. Done, thanks.

[andrioni]

I just pushed some documentation for both the manual and the standard library reference, with also some minor changes to some related BigFloat sections.

[vtjnash]

It might be good to move all of the function calls to the :openlibm library, although not critical. All this requires for working on windows is to bump openlibm to the latest version of master, which provides a complete and independent implementation of fenv.c / fenv.h.

Also, it looks like you have a build error. But otherwise, I'm fine with merging this.

[andrioni]

It's hard to believe, but apparently the build error is because of Ubuntu. On 12.04, they moved one of the subheaders, /usr/include/bits/fenv.h to /usr/include/i386-linux-gnu/bits/fenv.h, but forgot to change it in /usr/include/fenv.h, so in the end the macros don't get set, the preprocessor doesn't know about them and thus fenv_constants.jl isn't correctly created.

What should we do about it?

[StefanKarpinski]

Can we change src/fenv_constants.h so that it will pick up either location? Maybe we need to generate that file as well based on whether /usr/include/bits/fenv.h exists or not.

[andrioni]

I think we could check whether the preprocessor output parsed correctly, and if not, use the default x86/x64 values. I only found this issue on Ubuntu 12.04 (I didn't test it on other versions, I haven't got any other Ubuntu VMs handy), but not on Debian testing or on Arch Linux, so it might not be that big of an issue.

[andrioni]

The last commit uses some ugly regexes, but managed to solve it.

[andrioni]

Any comments?

[StefanKarpinski]

I'm testing locally before merging, but it looks great.

[bkalpert]

Hip, hip, hooray! One more component of reliable numerics in Julia.

[StefanKarpinski]

Yes, this is great. Kahan would/will be pleased.

[JeffBezanson]

This is great, thanks @andrioni .

I would prefer to use constants instead of types for the modes, as in:

immutable RoundingMode
  code::Int
end

const RoundToNearest = RoundingMode(0)
...

This way isa(RoundToNearest,RoundingMode) is true, which seems reasonable.

[ViralBShah]

Really glad to see this!

[simonbyrne]

Out of curiosity, what was the reason for having separate ordinary float and bigfloat functions (e.g. with_rounding vs with_bigfloat_rounding)?

[StefanKarpinski]

BigFloats support an extra rounding mode: RoundFromZero. The with_rounding_mode function only accepts rounding modes that work for everything.