dynamic rounding control and constant folding

[Simon Kornblith]

https://github.com/JuliaLang/julia/pull/8364 brought this to mind. Rather than continue the diversion there, I figured I'd start a thread here. Consider:

```julia
function f()
    set_rounding(Float64, RoundUp)
    println(1.0+5.0e-324)
    set_rounding(Float64, RoundNearest)
    nothing
end
```

This prints 1.0 because LLVM apparently performs constant folding with round to nearest during compilation, but if you paste the method body at the REPL, it gets interpreted, so it prints 1.0000000000000002. Does this bother anyone else?

Simon

[Chris Foster]

Oh, that's nasty, definitely disturbing.

Does anyone here know why dynamic scoping of IEEE floating point is
even desirable? It seems like it would be sufficient to support only
lexical scoping for interesting purposes like interval arithmetic
(actually, I don't know about any other interesting uses for the
rounding mode, which probably means I'm a bit naive.).

With lexical scoping the compiler would at least have a chance of
doing the right thing without throwing away a bunch of important
optimizations.

From a bit of reading, it seems like C99 just provides you with the
ability to turn of constant folding at lexical scope (along with some
other things) using FENV_ACCESS, or at a whole program level with
optimization flags. Ugh.

~Chris

[Simon Byrne]

To follow up from my comment on the thread, it turns out that dynamic-mode rounding isn't actually required: in the standard it is just a "should", which, 

• indicates that among several possibilities, one is recommended as particularly suitable, without mentioning or excluding others; or that a certain course of action is preferred but not necessarily required; or that (in the negative form) a certain course of action is deprecated but not prohibited (“should” means “is recommended to”).

The main argument for dynamic-mode rounding seems to be for debugging: Kahan gives some examples:

* http://www.cs.berkeley.edu/~wkahan/Stnfrd50.pdf

* http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf  (slides 57-60)

* http://www.eecs.berkeley.edu/~wkahan/Mindless.pdf  (section 11)

Now, I'm not sure that anyone else actually does this, but it is an interesting idea. However the problem is that the needs for this are fundamentally different than the needs for something like interval arithmetic. As a simple example, what should something like 

y = cos(x)

do when called with a RoundDown attribute? There are a

(1) run exactly the same code you would usually, but with the RoundDown attribute set on the FPU/SSE.

(2) call a completely different function, which would ideally return the nearest floating point number that is less than cos(x).

(1) is conceptually simple to do, but doesn't really make any sense outside of the debugging context. (2) is arguably more useful in other cases, but fiendishly difficult and computationally expensive to do correctly (see, for example, CRlibm). 

-Simon

[Toivo Henningsson]

It seems to me that dynamically scoped rounding is pretty much the same as forcibly adding an implicit rounding mode argument to any function that does floating point computations, and forcing the function to relate to this argument. Almost no functions are written to behave in a sensible manner wrt an externally supplied rounding mode.

Instead, I think that those few functions that are written with rounding mode in mind should take the rounding mode as an explicit argument, and, most of all, it should not be possible to manipulate the results of a function that was written without consideration to rounding modes by changing the rounding mode before calling it, since that would expose implementation details in an unintended way, making the function a leaky abstraction.

I guess this points in the direction of lexically scoped rounding control, or to have versions of the basic floating point operations that take an explicit rounding mode argument.

[Chris Foster]

On Thu, Oct 2, 2014 at 11:55 PM, Toivo Henningsson <toiv...@gmail.com> wrote:
> It seems to me that dynamically scoped rounding is pretty much the same as
> forcibly adding an implicit rounding mode argument to any function that does
> floating point computations, and forcing the function to relate to this
> argument. Almost no functions are written to behave in a sensible manner wrt
> an externally supplied rounding mode.
>
> Instead, I think that those few functions that are written with rounding
> mode in mind should take the rounding mode as an explicit argument, and,
> most of all, it should not be possible to manipulate the results of a
> function that was written without consideration to rounding modes by
> changing the rounding mode before calling it, since that would expose
> implementation details in an unintended way, making the function a leaky
> abstraction.

That makes a lot of sense. As soon as a function f combines the
builtin floating point operations in a nontrivial way, setting the
rounding mode outside of f doesn't get you any guarantee that the
result f(x) is correctly rounded.

I don't think it's practical, but now I can't stop thinking about
multiple dispatch based on the current rounding mode.

[Patrick O'Leary]

On Wednesday, October 1, 2014 7:01:44 PM UTC-5, Simon Kornblith wrote:

This prints 1.0 because LLVM apparently performs constant folding with round to nearest during compilation, but if you paste the method body at the REPL, it gets interpreted, so it prints 1.0000000000000002. Does this bother anyone else?

Julia code is never interpreted--there is no interpreter. This is a function, so global scope issues shouldn't be a problem. Do different optimization passes run for code entered at the REPL?

[Simon Kornblith]

There is an interpreter that is sometimes used for top-level code, but can't handle some language constructs like loops.

begin

    set_rounding(Float64, RoundUp)
    println(1.0+5.0e-324)
    set_rounding(Float64, RoundNearest)
    nothing
end

prints 1.0000000000000002 but

begin
    for i = 1:1; end

    set_rounding(Float64, RoundUp)
    println(1.0+5.0e-324)
    set_rounding(Float64, RoundNearest)
    nothing
end

prints 1.0.

Simon

[Patrick O'Leary]

On Saturday, October 4, 2014 7:39:36 AM UTC-5, Simon Kornblith wrote:

On Saturday, October 4, 2014 7:49:13 AM UTC-4, Patrick O'Leary wrote:

On Wednesday, October 1, 2014 7:01:44 PM UTC-5, Simon Kornblith wrote:

This prints 1.0 because LLVM apparently performs constant folding with round to nearest during compilation, but if you paste the method body at the REPL, it gets interpreted, so it prints 1.0000000000000002. Does this bother anyone else?

Julia code is never interpreted--there is no interpreter. This is a function, so global scope issues shouldn't be a problem. Do different optimization passes run for code entered at the REPL?

There is an interpreter that is sometimes used for top-level code, but can't handle some language constructs like loops.

I admit I've never seen that before. I could have sworn I was just parroting Jeff here, but it seems I'm remembering something that he never wrote. My apologies.

Obvious follow-up: given that there are potential correctness, or at least consistency, issues, does a partial-language interpreter make sense? Is there a performance advantage for simple statements? Is it needed for bootstrapping?

[Chris Foster]

On Sun, Oct 5, 2014 at 7:36 AM, Patrick O'Leary

> Obvious follow-up: given that there are potential correctness, or at least
> consistency, issues, does a partial-language interpreter make sense? Is
> there a performance advantage for simple statements? Is it needed for
> bootstrapping?

I think the issue here is orthogonal to interpreter vs jit compiler.
The problem is that the compiler uses the RoundNearest rounding mode
during the constant folding optimization pass, but this might be
inconsistent with the dynamic floating point rounding mode set at
runtime using set_rounding.

~Chris

[Patrick O'Leary]

Ahh, I think I'm finally getting there (I've never played with rounding modes, so I'm learning here.)  So LLVM does its constant fold, but since the surrounding code isn't actually executing the FPU is still stuck in the default rounding mode when the fold occurs? That just sounds fundamentally broken. Which may have been the point y'all were trying to make. Thanks for patiently explaining.

[Simon Byrne]

The problem seems to be that LLVM doesn't actually do anything to support rounding modes (e.g. via intrinsics): our current get/set_rounding functions just call openlibm functions, which LLVM can't reason about. I guess this is something that needs to be done on their end.

As far as an interface goes, we could just add an extra argument, e.g.

+(::RoundingMode, x, y) = ...

and/or abuse Unicode combining characters, e.g.

+⃖(x,y) = +(RoundDown, x, y)

+⃗(x,y) = +(RoundUp, x, y)

Unfortunately there's no inward pointing arrow for RoundToZero...

-Simon

[Stefan Karpinski]

I'm pretty convinced at this point that the whole rounding mode business was never well enough designed to actually be used by anything but hand-written assembly code.

[Toivo Henningsson]

Including to be used by Julia behind the scenes to expose it under a saner model? I agree that the whole global flags business is kind of a disaster, but beyond that, it seems that the real trouble is that we're at the mercy of llvm here, which doesn't seem to have better support for rounding modes than C does.