7 ms·

Good essay. The dynamic scoping / mutable global state issue, especially with respect to rounding modes, is a total nightmare. We've been trying our hardest to

Good essay. The dynamic scoping / mutable global state issue, especially with respect to rounding modes, is a total nightmare. We've been trying our hardest to provide great support for rounding modes in Julia and it just doesn't work well. The trouble is that there are functions that need normal rounding regardless of the current rounding mode while other operations can safely respect the current rounding mode. There are still other functions for which neither is appropriate. A software-implemented sin function, for example, may completely break if used with an unexpected rounding mode, so you don't want that. But you also don't want to just completely ignore the rounding mode either. What you really want is a sin function that returns a correct result, rounded according to the rounding mode. But that means you really need five variants of the sin function and you need to dynamically dispatch between them based on the current rounding mode – i.e. global mutable state – which is, of course, quite bad for performance, not to mention the explosion of complexity. I've come to believe that the rounding mode concept as designed in IEEE 754 was only ever feasible to use when you were writing your own assembly code. The essay is not quite right about the motivation for rounding modes – they are not intended to support interval arithmetic. When we got the chance to ask him about it, Kahan was quite negative about interval arithmetic, noting that intervals tend to grow exponentially large as your computation progresses. Instead, the motivation he suggested was that you could run a computation in all the rounding modes, and if the results in all of them were reasonably close, you could be fairly certain that the result was accurate. Of course, at the higher level the essay is exactly right: the fact that this motivation isn't explained in the IEEE 754 spec is precisely what allowed this misunderstanding about the motivation to occur at all.

conistonwater 12y ago
How important is it to implement floating point arithmetic exactly as the standard specifies it? What if the same set of features is generally available, but not quite in the required way?
- StefanKarpinski 12y ago
  That depends on what you want. Some people don't care at all, but others have legitimate needs for a lot of the features described in the spec. So if you can at all, it's best to follow the spec, imo. That way languages don't do a lot of annoyingly different, similar things. And, as this essay points out, the spec is wildly successful – it is arguably the most successful spec of all time – so second-guessing it is usually a mistake. Kahan at al. have forgotten more about numerical analysis than most of us will ever know.
  - gre 12y ago
    Factor has a library for setting up floating point environments that might interest you. I think we implemented all of the functionality you might need and it's nestable with dynamic scoping. I haven't needed to use it for anything though. Code: https://github.com/slavapestov/factor/blob/master/basis/math/floats/env/env.factor https://github.com/slavapestov/factor/blob/master/basis/math... Tests: https://github.com/slavapestov/factor/blob/master/basis/math/floats/env/env-tests.factor https://github.com/slavapestov/factor/blob/master/basis/math... Docs: http://docs.factorcode.org/content/article-math.floats.env.html http://docs.factorcode.org/content/article-math.floats.env.h...
- nraynaud 12y ago
  a lot of analysis is based on ieee 754, that would mean developers would have to re-do the analysis and change the finicky algorithms. I'm personally not able to do it. I nod at the research papers I read, but I would not be able to write them.
phkahler 12y ago
>> But that means you really need five variants of the sin function and you need to dynamically dispatch between them based on the current rounding mode Just define a rounding behavior for the language and implement it that way. Don't claim full 754 support, just specify the strategy used by the language. A sin function should behave according to the design of that function and should be able to ignore any previous state of the FP hardware. I have not seen a language that directly supports setting the rounding modes, so any language libraries can do what they like - you don't need to preserve or worry about something you don't offer the option to modify.
- gre 12y ago
  Check out my other comment below. Factor has with-rounding-mode and you can nest these modes. Example: +round-up+ [ 1.0 3.0 /f double>bits .h +round-zero+ [ 1.0 3.0 /f double>bits .h ] with-rounding-mode ] with-rounding-mode ! output ! 3fd5555555555556 ! 3fd5555555555555
  - andrioni 12y ago
    Julia also supports this: with_rounding(Float64, RoundUp) do println(bits(1.0 / 3.0)) with_rounding(Float64, RoundToZero) do println(bits(1.0 / 3.0)) end end 0011111111010101010101010101010101010101010101010101010101010110 0011111111010101010101010101010101010101010101010101010101010101
phkahler 12y ago
With all the state and modes and such, is there a mode that just does perhaps: round to even, clamp to infinity, never produce NaN, and never trigger an exception? Because that's basically what a lot of people want - do the math and handle the extremes in reasonable ways without throwing exceptions.
- ekimekim 12y ago
  It isn't possible to never produce a NaN. 0/0 and inf-inf, for example, have no reasonable result that can be produced. And indeed, the spec defines those operations as resulting in NaN.
- Someone 12y ago
  I know that's what a lot of people want, but I do not see how one can 'handle the extremes in reasonable ways' without producing NaNs. I have heard people advocate 0/0 = 0, but IMO that is NOT 'handling the extremes in reasonable ways'. What would you propose for 0/0?
  - Animats 12y ago
    More than that, a signaling NaN should raise an exception if used in a comparison. Otherwise, you get a bogus comparison and a bogus branch.
- stephencanon 12y ago
  In your "never produce NaN" arithmetic, what do you want sqrt(-1) to be, and why is that a better answer than NaN?
  - lmm 12y ago
    I want it to throw an exception - a language-level exception - and I'm fine with sacrificing a bit of performance (i.e. checking the flags after every operation) to do so. This is a better answer because it means I see the error where it actually happened, rather than getting a NaN out the end of my calculation and having no idea which particular computation it came from.
stephencanon 12y ago
> But that means you really need five variants of the sin function and you need to dynamically dispatch between them based on the current rounding mode. That would be one way to achieve the behavior you describe, but it is certainly not the only way, nor is it the best. You could, for example, do a computation that produces a result rounded to odd with a few extra bits (independent of the prevailing rounding mode), then rounds that result to the destination precision using the prevailing mode. This assumes that you really want correct-rounded transcendentals anyway, which at present are sufficiently slow that you might as well dispatch based on rounding mode (despite the herculean efforts of the CRLibm crew). If faithful rounding is all that is required, even simpler methods exist that are extremely performant.
- StefanKarpinski 12y ago
  Yes, this is a very fair point – it's not like our transcendental function implementations are uniformly correctly rounded anyway – being within 1 ulp is good enough. I suppose I might be happy if there were a lexical flag for each floating-point instruction (i.e. a bit of the opcode) that chose whether to use the current rounding mode or the normal one. The trouble now is that you have no choice but to use the dynamically scoped one, which is problematic often enough that one is pushed towards just opting out of changing rounding modes altogether.
plesner 12y ago
> Instead, the motivation he suggested was that you could run a computation in all the rounding modes, and if the results in all of them were reasonably close, you could be fairly certain that the result was accurate. I actually thought that was the use case I was describing, though I would expect round-positive and round-negative to be enough. Don't the other rounding modes yield results within those bounds?
- simonbyrne 12y ago
  Consider something like a-(b+c) To obtain the correct lower bound, you need to round the addition up, and the subtraction down. This is also incredibly inefficient on current processors, as the registers are flushed after each mode change.
  - plesner 12y ago
    You're right. To be absolutely sure I tried writing a little program that does the thing I describe in the post, https://gist.github.com/plesner/7a13b5c8d269315df645 https://gist.github.com/plesner/7a13b5c8d269315df645, and sure enough it breaks down completely, round default: 0.072817 round up: 0.072703 round down: 0.072931 Now I don't understand how you get a reliable result using rounding modes at all.
    - simonbyrne 12y ago
      Interval arithmetic is tricky: you have to make sure you round in the most conservative direction for each operation. As other commenters mentioned, these bounds can end up being so large as to be useless.
      - plesner 12y ago
        The way I hear this is: you have two options. - Either use interval arithmetic which is tricky and may give you useless large bounds, but those bounds are guaranteed to hold. - Or just try all the rounding modes which is less likely to give you large bounds, but now you have no guarantee that those bounds say something meaningful about your computation. Or does the second case give meaningful guarantees?
lmm 12y ago
Honestly this sounds like a perfect use case for monads - you write your computation with a dependency on the rounding mode, and you can either ripple that all the way up (to your whole program if need be), or specify an explicit rounding mode at a level where it makes sense. It would also allow sensible (if poorly-performing) rounding for subtraction cases. Has anyone implemented it like that?
- ufo 12y ago
  I think monads might possibly somewhat alleviate the problem of writing 5 different versions of sin by hand but I think you would still need to do some dynamic dispatching to handle the rounding modes, which according to Stephan is going to be a big performance hit.

conistonwater 12y ago

How important is it to implement floating point arithmetic exactly as the standard specifies it? What if the same set of features is generally available, but not quite in the required way?

StefanKarpinski 12y ago
That depends on what you want. Some people don't care at all, but others have legitimate needs for a lot of the features described in the spec. So if you can at all, it's best to follow the spec, imo. That way languages don't do a lot of annoyingly different, similar things. And, as this essay points out, the spec is wildly successful – it is arguably the most successful spec of all time – so second-guessing it is usually a mistake. Kahan at al. have forgotten more about numerical analysis than most of us will ever know.
- gre 12y ago
  Factor has a library for setting up floating point environments that might interest you. I think we implemented all of the functionality you might need and it's nestable with dynamic scoping. I haven't needed to use it for anything though. Code: https://github.com/slavapestov/factor/blob/master/basis/math/floats/env/env.factor https://github.com/slavapestov/factor/blob/master/basis/math... Tests: https://github.com/slavapestov/factor/blob/master/basis/math/floats/env/env-tests.factor https://github.com/slavapestov/factor/blob/master/basis/math... Docs: http://docs.factorcode.org/content/article-math.floats.env.html http://docs.factorcode.org/content/article-math.floats.env.h...
nraynaud 12y ago
a lot of analysis is based on ieee 754, that would mean developers would have to re-do the analysis and change the finicky algorithms. I'm personally not able to do it. I nod at the research papers I read, but I would not be able to write them.

phkahler 12y ago

>> But that means you really need five variants of the sin function and you need to dynamically dispatch between them based on the current rounding mode Just define a rounding behavior for the language and implement it that way. Don't claim full 754 support, just specify the strategy used by the language. A sin function should behave according to the design of that function and should be able to ignore any previous state of the FP hardware. I have not seen a language that directly supports setting the rounding modes, so any language libraries can do what they like - you don't need to preserve or worry about something you don't offer the option to modify.

gre 12y ago
Check out my other comment below. Factor has with-rounding-mode and you can nest these modes. Example: +round-up+ [ 1.0 3.0 /f double>bits .h +round-zero+ [ 1.0 3.0 /f double>bits .h ] with-rounding-mode ] with-rounding-mode ! output ! 3fd5555555555556 ! 3fd5555555555555
- andrioni 12y ago
  Julia also supports this: with_rounding(Float64, RoundUp) do println(bits(1.0 / 3.0)) with_rounding(Float64, RoundToZero) do println(bits(1.0 / 3.0)) end end 0011111111010101010101010101010101010101010101010101010101010110 0011111111010101010101010101010101010101010101010101010101010101

phkahler 12y ago

With all the state and modes and such, is there a mode that just does perhaps: round to even, clamp to infinity, never produce NaN, and never trigger an exception? Because that's basically what a lot of people want - do the math and handle the extremes in reasonable ways without throwing exceptions.

ekimekim 12y ago
It isn't possible to never produce a NaN. 0/0 and inf-inf, for example, have no reasonable result that can be produced. And indeed, the spec defines those operations as resulting in NaN.
Someone 12y ago
I know that's what a lot of people want, but I do not see how one can 'handle the extremes in reasonable ways' without producing NaNs. I have heard people advocate 0/0 = 0, but IMO that is NOT 'handling the extremes in reasonable ways'. What would you propose for 0/0?
- Animats 12y ago
  More than that, a signaling NaN should raise an exception if used in a comparison. Otherwise, you get a bogus comparison and a bogus branch.
stephencanon 12y ago
In your "never produce NaN" arithmetic, what do you want sqrt(-1) to be, and why is that a better answer than NaN?
- lmm 12y ago
  I want it to throw an exception - a language-level exception - and I'm fine with sacrificing a bit of performance (i.e. checking the flags after every operation) to do so. This is a better answer because it means I see the error where it actually happened, rather than getting a NaN out the end of my calculation and having no idea which particular computation it came from.

stephencanon 12y ago

> But that means you really need five variants of the sin function and you need to dynamically dispatch between them based on the current rounding mode. That would be one way to achieve the behavior you describe, but it is certainly not the only way, nor is it the best. You could, for example, do a computation that produces a result rounded to odd with a few extra bits (independent of the prevailing rounding mode), then rounds that result to the destination precision using the prevailing mode. This assumes that you really want correct-rounded transcendentals anyway, which at present are sufficiently slow that you might as well dispatch based on rounding mode (despite the herculean efforts of the CRLibm crew). If faithful rounding is all that is required, even simpler methods exist that are extremely performant.

StefanKarpinski 12y ago
Yes, this is a very fair point – it's not like our transcendental function implementations are uniformly correctly rounded anyway – being within 1 ulp is good enough. I suppose I might be happy if there were a lexical flag for each floating-point instruction (i.e. a bit of the opcode) that chose whether to use the current rounding mode or the normal one. The trouble now is that you have no choice but to use the dynamically scoped one, which is problematic often enough that one is pushed towards just opting out of changing rounding modes altogether.

plesner 12y ago

> Instead, the motivation he suggested was that you could run a computation in all the rounding modes, and if the results in all of them were reasonably close, you could be fairly certain that the result was accurate. I actually thought that was the use case I was describing, though I would expect round-positive and round-negative to be enough. Don't the other rounding modes yield results within those bounds?

simonbyrne 12y ago
Consider something like a-(b+c) To obtain the correct lower bound, you need to round the addition up, and the subtraction down. This is also incredibly inefficient on current processors, as the registers are flushed after each mode change.
- plesner 12y ago
  You're right. To be absolutely sure I tried writing a little program that does the thing I describe in the post, https://gist.github.com/plesner/7a13b5c8d269315df645 https://gist.github.com/plesner/7a13b5c8d269315df645, and sure enough it breaks down completely, round default: 0.072817 round up: 0.072703 round down: 0.072931 Now I don't understand how you get a reliable result using rounding modes at all.
  - simonbyrne 12y ago
    Interval arithmetic is tricky: you have to make sure you round in the most conservative direction for each operation. As other commenters mentioned, these bounds can end up being so large as to be useless.
    - plesner 12y ago
      The way I hear this is: you have two options. - Either use interval arithmetic which is tricky and may give you useless large bounds, but those bounds are guaranteed to hold. - Or just try all the rounding modes which is less likely to give you large bounds, but now you have no guarantee that those bounds say something meaningful about your computation. Or does the second case give meaningful guarantees?

lmm 12y ago

Honestly this sounds like a perfect use case for monads - you write your computation with a dependency on the rounding mode, and you can either ripple that all the way up (to your whole program if need be), or specify an explicit rounding mode at a level where it makes sense. It would also allow sensible (if poorly-performing) rounding for subtraction cases. Has anyone implemented it like that?

ufo 12y ago
I think monads might possibly somewhat alleviate the problem of writing 5 different versions of sin by hand but I think you would still need to do some dynamic dispatching to handle the rounding modes, which according to Stephan is going to be a big performance hit.