Rounding and the IEEE rule

[Hans W Borchers]

Many programming languages and scientific computing systems follow the "round half to even" 

tie-breaking rule; this is also the default rounding mode in the IEEE standard for floating-point 

arithmetic (IEEE 754). So in R or Python (with NumPy), but not in, e.g., MATLAB) we have

    round(0.5) #=> 0

    round(1.5) #=> 2

    round(2.5) #=> 2

    round(3.5) #=> 4

    ...

 but in Julia 

    julia> (round(0.5), round(1.5), round(2.5), round(3.5))

    (1.0,2.0,3.0,4.0)

    julia> (round(-0.5), round(-1.5), round(-2.5), round(-3.5))

    (-1.0,-2.0,-3.0,-4.0)

Is there a special reason for Julia to follow instead the "round half away from zero" rule, for

instance compatibility with MATLAB?

[Stefan Karpinski]

We follow C, Fortran, Matlab, Python and most other programming languages here. R and NumPy's rule is pretty unusual; it has some nice statistical properties (it's apparently known as "statistician's rounding"), but is quite awkward for general programming tasks.

[John Myles White]

One question: I have the impression that the round() function is not affected by the currently chosen rounding rule in Julia. Is that right?

-- John

[Stefan Karpinski]

This isn't really related to IEEE rounding modes. Floating-point rounding modes are about choosing which of the closest representable floating-point values an operation should produce when the true value is between them. The round function is a well-defined mathematical function regardless of IEEE rounding mode. The man page for the libc round function says:

The round() functions return the integral value nearest to x rounding halfway cases away from zero, regardless of the current rounding direction.

[Kevin Squire]

Couldn't this be provided by the get_rounding/set_rounding/with_rounding framework?

[Stefan Karpinski]

It could, but it shouldn't. The rounding mode and how you round to integers are completely different matters.

[Jeff Bezanson]

In the IEEE standard, rounding results of floating-point computations
needs to be done by roundTiesToEven by default. However, for rounding
floats to integers it merely describes five different functions (one
for each commonly-used rounding method), and has no opinion on which
you should use.

I think our choice is better for general programming. Say you index an
array using the expression `A[round(x + i)]`, where i is consecutive
integers, and x is a fraction between 0 and 1. Clearly, this will
touch all array elements. But with rounding ties to even, you
*usually* get that behavior, except when x==0.5, and then you'd only
get even-index elements. That seems crazy to me.

[Andreas Lobinger]

Hello colleagues,

On Thursday, June 5, 2014 2:41:51 AM UTC+2, Stefan Karpinski wrote:

It could, but it shouldn't. The rounding mode and how you round to integers are completely different matters.

Yes.

If you deal with floating point numbers in signal processing or other system modeling cases, the stability (and/or robustness on input parameters) depend on the rounding of FP. Everyone of us learned in Numerical Analysis 1 that using FP creates you a non-linear system (a+(b+c) != (a+b)+c) but most people forget later on. That IEEE included this choose of rounding was a great step forward. A good scientific calculation environment enables this.

But still, rounding to integer (actually there is/should be ceil+floor+round) is needed for array operations.

[Hans W Borchers]

I think I now understand the difference between rounding rules and rounding to
integers as a mathematical function.

Then, how can I in Julia round a number such as 0.025 to two decimal places and
what will I get?  round() does not have a second parameter indicating the number
of digits to round to like, e.g., round() in R has.  I tried to search through
the Standard Library and did not see an obvious choice.

[Chris Foster]

Just try it, it works!

julia> round(1.2345, 2)
1.23

You can check the possible versions of round() in the repl using methods():

julia> methods(round)
# 15 methods for generic function "round":
round(x::Integer) at int.jl:368
round(x::Float64) at float.jl:107
round(x::Rational{T<:Integer}) at rational.jl:164
round(x::Float32) at math.jl:297
round{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}) at sparse/sparsematrix.jl:431
round(M::SymTridiagonal{T}) at linalg/tridiag.jl:32
round(M::Tridiagonal{T}) at linalg/tridiag.jl:188
iround(M::Bidiagonal{T}) at linalg/bidiag.jl:67
round{T<:Real}(::AbstractArray{T<:Real,1}) at operators.jl:354
round{T<:Real}(::AbstractArray{T<:Real,2}) at operators.jl:355
round{T<:Real}(::AbstractArray{T<:Real,N}) at operators.jl:357
round(x::Float16) at float16.jl:101
round(x::BigFloat) at mpfr.jl:639
round(x,digits::Integer) at floatfuncs.jl:85
round(x,digits::Integer,base::Integer) at floatfuncs.jl:85

It's that second last one you want. Unfortunately the list can get
quite long sometimes and it becomes hard to find what you're after.

~Chris

[Kevin Squire]

It's also documented:

 julia> help(round)

INFO: Loading help data...

Base.round(x[, digits[, base]])

   "round(x)" returns the nearest integral value of the same type as

   "x" to "x". "round(x, digits)" rounds to the specified number

   of digits after the decimal place, or before if negative, e.g.,

   "round(pi,2)" is "3.14". "round(x, digits, base)" rounds

   using a different base, defaulting to 10, e.g., "round(pi, 1, 8)"

   is "3.125".

http://julia.readthedocs.org/en/latest/stdlib/base/#Base.round

[Hans W Borchers]

Right, my fault.
Sorry, I am still wondering: round(x,n) rounds towards zero, even for non-integers
(as is documented in;

    "The round() functions return the integral value nearest to x rounding halfway

     cases away from zero, regardless of the current rounding direction.")

    julia> round(0.025, 2)
    0.03

so round(x,2) works like round(100*x)/100. (This is what I do in MATLAB.)
What I would like to have is a function round_ieee(x,n) such that

    julia> round_ieee(0.025, 2)
    0.02

Can I do this somehow with with_rounding() do ... end ?

[Simon Byrne]

On Thursday, 5 June 2014 09:47:21 UTC+1, Hans W Borchers wrote:

What I would like to have is a function round_ieee(x,n) such that

    julia> round_ieee(0.025, 2)
    0.02

Can I do this somehow with with_rounding() do ... end ?

Not generally: the IEEE notions of even and odd only pertain to binary arithmetic (as determined by the last bit of the mantissa), trying to force this notion into approximate decimal representations is going to cause all sorts of confusion: technically the Float64 represented by 0.025 is actually 1.39e-18 larger than the exact 0.025, so should it be rounded to 0.03 instead?

The IEEE754 spec itself only defines explicit rounding (as opposed to rounding after calculations) to integers (section 5.9 for those following along at home). According to the standard, IEEE754 compatible languages are required to provide methods for the following:

a) round to nearest integral: ties to even

b) round to nearest integral: ties to away (round)

c) round to integral toward zero (trunc)

d) round to integral toward positive (ceil)

e) round to integral toward negative (floor)

f) round to integral using rounding mode 

At the moment, I don't think we implement (a) or (f).

[Stefan Karpinski]

Right, calling one particular rounding style round_ieee would be very odd – the IEEE 754 spec lists various functions for rounding floats to integers and doesn't favor any of them over the others. This is unlike rounding modes where it does indicate that the default mode should be ties to even.

[Hans W Borchers]

I started this thread long time ago with a question about rounding rules and the IEEE floating point standard. I felt like being criticized for even thinking Julia could follow the "round-to-even" rule. Now I learn that Julia version 0.4 will apply this rule (as default?).

We had similar discussions about whether Julia should use the "+" or ".+" operator for adding scalars to vectors. One strong argument was that, mathematically, there is no addition of scalars to vectors. Since then Julia has switched back to using "+" again.

I would wish the replies were not that strict when questions or requests come up for a different behavior of Julia operators or functions. And because mailing lists are also an information source for users, it may be appropriate to mention such changes at the end of threads like this one.

[Andreas Noack]

It can be difficult to remember all threads on the mailing list that relate to changes made in Julia. Hence, it is very helpful when users make cross references such that you are doing here. A link to the exact issues on github would be even better so here they are

https://github.com/JuliaLang/julia/issues/8750

https://github.com/JuliaLang/julia/pull/9344

[cdm]

thank you for the cross reference, H.W. Borchers ...

and thank you for the link out to issues, A. Noack.

as a fellow user, i appreciate having these.

best,

cdm

[Simon Byrne]

On Friday, 26 December 2014 06:14:34 UTC-6, Hans W Borchers wrote:

I started this thread long time ago with a question about rounding rules and the IEEE floating point standard. I felt like being criticized for even thinking Julia could follow the "round-to-even" rule. Now I learn that Julia version 0.4 will apply this rule (as default?).

I apologise if you felt I was being critical, that was not at all my intention. These are complicated issues, and I don't claim to have all the solutions. In fact, we still have yet to resolve the problem with the digits argument: I have just opened an issue explaining the problem here:

https://github.com/JuliaLang/julia/issues/9464

Contributions to the discussion are certainly welcome.

-simon

[Stefan Karpinski]

I'm also sorry if you felt like you were being criticized in any way for proposing alternate behaviors. Please keep in mind that any disagreement is just that – disagreement. Even if everyone seems disagree with you, please do make your case. As you've noted, they may change their minds – or they might not. We can't come to decisions unless everyone tries to make the case for what they believe is the best behavior in a civil, reasoned fashion. To me, the fact that these behavior both changed, despite initially being unpopular, indicates just how important this kind of discussion is.

Regarding this particular change, I'm personally not really in favor of the default being the IEEE ties-to-even behavior – I don't think the numerical argument in its favor is really compelling and I think the current C-like round-ties-from-zero behavior is less surprising and should continue to be the default. But I'm not the one doing the work on this and I'm not the most informed person to make such a decision. So let's try the IEEE behavior and see how it goes. If this turns out to be problematic or annoying (I kind of suspect it will), we can still change our minds before Julia 1.0, so now is the time to give it a shot.

[Hans W Borchers]

Don't worry. I am a long-term ... user and as such got used to really harsh criticism. As others have noted, this is - besides Ruby - the most friendly mailing list I have seen (and I am sometimes not living up to that standard). What irritates me a bit are these back-and-forth decisions. I accepted ".+" for scalar plus vector operations (with a bit of teeth grinding, but listening to a good mathematical analogy), used it in some programs, only learning a few weeks later that a change agent had struck again.

I know about (some of) the problems with rounding. The perhaps mathematically most complete and correct arithmetic computing system, PARI/GP, still uses "round-to-+Inf" and 'generations' of mathematicians have lived well with that. And the "problems with the digit argument" mentioned are perhaps a reason why many systems like Octave, Matlab (up to 2014a), Mathematica, PARI/GP do not allow for a second parameter in their 'round' functions. Maybe there is really no satisfying solution here for this problem.

[Stefan Karpinski]

On Fri, Dec 26, 2014 at 6:00 PM, Hans W Borchers <hwbor...@gmail.com> wrote:

What irritates me a bit are these back-and-forth decisions. I accepted ".+" for scalar plus vector operations (with a bit of teeth grinding, but listening to a good mathematical analogy), used it in some programs, only learning a few weeks later that a change agent had struck again.

I agree it's annoying, but I think it's important to try things and recognize when an experiment was a mistake. Fortunately, we haven't done this very often.

 

I know about (some of) the problems with rounding. The perhaps mathematically most complete and correct arithmetic computing system, PARI/GP, still uses "round-to-+Inf" and 'generations' of mathematicians have lived well with that. And the "problems with the digit argument" mentioned are perhaps a reason why many systems like Octave, Matlab (up to 2014a), Mathematica, PARI/GP do not allow for a second parameter in their 'round' functions. Maybe there is really no satisfying solution here for this problem.

I would also be ok with round-to-Inf behavior since to me that's what the mathematical round function does. Of course, it's also nice that round(-x) == -round(x) for all x. Also, we've traditionally tried to keep things close to C when there aren't other languages (Matlab, Python) that serve as better models.

[Tamas Papp]

My rule of thumb is that if I need to worry about the IEEE rounding mode
making a difference in my results, then my problem is really
ill-conditioned and rounding the least of my worries. So I usually just
ignore it.

Same applies to rounding for data displays: if I am worried about 3.5
being rounded to 3 or 4, maybe I should just display that extra
significant digit.

Best,

Tamas

[Tobias Knopp]

Hans,

yes this can be frustrating. But Julia is in flux and I have to say that I am very happy that changes are reverted if they turn out to be not practically (e.g. ".+" which I also was not happy with).

One important thing to note: If you don't like this back and forth it should be better to just use a stable release. I am using 0.3 since summer for my research because it is important for me that thing just work. Master certainly contains several goodies I am eager to try but for serious coding it is a no go when one has to always adapt to recent changes on a development branch. Once 0.4 is released I will port thing in a single pass.

Cheers

Tobi

[Hans W Borchers]

Just a remark: "round-towards-+Inf" is not the same as "round-away-from-zero" and 

does not fulfill "round(-x) = -round(x)".

    gp> round([1.5, 2.5])

    %1 = [2, 3]

    gp> round([-1.5, -2.5])

    %2 = [-1, -2]

I am astonished that PARI/GP is rarely mentioned when comparing Julia to other 

scientific computing systems -- not even in the long list of issue #8750 or on 

the 'Gang of Forty' slide of the Vision talk at MIT.

But PARI/GP is, besides Mathematica, the most competent mathematical computing 

system and has set the standard in number theory computations for over 30 years.

BTW:

@Tamas Some mathematicians and statisticians **do** care about rounding rules.

[Tamas Papp]

I did not mean to imply that no one cares about rounding rules, or that
they are unimportant -- sorry if I was not clear. I only said that in my
(necessarily limited) experience, when IEEE rounding rules start to
matter then I am facing a different, usually more fundamental, problem.

However, these discussions are necessarily very abstract. Having a
concrete issue or use case where you find one rounding mode preferable
to another would help focus the discussion.

IMO functions like round are usually for displaying numbers, not for
computation, so am fine with rounding away from 0. Especially since for
decimal displays converted from an underlying binary (float), it may not
matter for practical purposes. Eg

julia> round(1.0000000005, 9)
1.000000001

may look like rounding away from 0, but because of inexact
representation, rounding is clearcut anyway:

julia> rationalize(1.0000000005)-1
1//1999999834

Best,

Tamas

[Hans W Borchers]

> However, these discussions are necessarily very abstract. Having a 

> concrete issue or use case where you find one rounding mode preferable 

> to another would help focus the discussion.

@Tamas For an example where it might have the potential to cause problems, see

    Potential future problem: rounding mode changed in core Julia #13

    https://github.com/tlycken/Interpolations.jl/issues/13

[Tamas Papp]

Hi Hans,

If I understand corretly, #13 is a potential problem arising from a
change in behavior, not an argument in favor of one rounding mode vs
another per se.

The only reasonably cogent argument I know is mentioned in #8750
(unbiasedness). I am glad that the issue is now closed, since that means
no more electrons are wasted on it, but I am still curious about
practical examples where rounding mode makes a difference in numerical
code. All I have seen rely on ill-conditioning, which is a separate
issue.

Best,

Tamas

On Sat, Dec 27 2014, Hans W Borchers <hwbor...@gmail.com> wrote:

[Hans W Borchers]

Tamas,

I agree with you that 'unbiased rounding' is a strong argument in favor of the

"round-to-even" rule. It appears natural that e.g. R falls for this rule.

The interpolation case seems to favor a rule that moves all ties in the same

direction, being it upwards or downwards, and independent of the sign. Because

otherwise the length of interpolation intervals is not so easily predictable.

This thread is getting too long now. If you want your arguments to be taken

into account, I suggest you join the discussion at Simon Byrne's issue #9464.

Thanks to all for your patience.

For those interested, have a look at  PARI/GP home .

[Stuart Brorson]

On Sat, 27 Dec 2014, Tamas Papp wrote:

> The only reasonably cogent argument I know is mentioned in #8750
> (unbiasedness). I am glad that the issue is now closed, since that means
> no more electrons are wasted on it, but I am still curious about
> practical examples where rounding mode makes a difference in numerical
> code. All I have seen rely on ill-conditioning, which is a separate
> issue.

My understanding is that the ability to round up or down is useful for
examining your computation's sensitivity to round off errors. This
use case assumes you can manually set the rounding mode to one of the
flavors supported by IEEE-754. For example, you'd set rounding to
"round towards +inf", do your computation and get your result. Then
you'd set rounding to "round towards -inf", and do the same thing.
Then compare the two results. If they're "close enough", then your
computation is robust. If the results differ wildly, then you have a
problem -- probably your computation is ill conditioned, as you say.

FWIW, this is similar to (but not the same as) "interval arithmetic",
another way to quantify your computation's sensitivity to round-off
errors.

http://en.wikipedia.org/wiki/Interval_arithmetic

I don't know if Julia supports the ability to change rounding mode on
the fly. However, I believe that the Intel C compiler does support
changing rounding modes at compile time, albeit using compiler flags
-- not as convenient as issuing a library call or typing a command at
a command prompt.

Stuart

[Stefan Karpinski]

On Sat, Dec 27, 2014 at 5:17 AM, Hans W Borchers <hwbor...@gmail.com> wrote:

Just a remark: "round-towards-+Inf" is not the same as "round-away-from-zero" and 

does not fulfill "round(-x) = -round(x)".

Right, I realize it wasn't clear from my email, but that was an argument against the current C-like round-away-from-zero behavior and for mathy round-to-inf behavior.

 

I am astonished that PARI/GP is rarely mentioned when comparing Julia to other 

scientific computing systems -- not even in the long list of issue #8750 or on 

the 'Gang of Forty' slide of the Vision talk at MIT.

I have to confess that I've never, before it being mentioned here, heard of PARI/GP.

[Jason Riedy]

And Tamas Papp writes:
> My rule of thumb is that if I need to worry about the IEEE rounding mode
> making a difference in my results, then my problem is really
> ill-conditioned and rounding the least of my worries. So I usually just
> ignore it.

You know of ill-conditioning, so you are not the concern.
Round-to-even is meant to allay bias problems for the vast
majority of programmers, like those responsible for the Vancouver
stock exchange index in 1982. They used truncation which, over
time, roughly halved the index. Rounding upward in the same
scenario would have doubled it...

The choices in 754 were made with a nearly revolting amount of
deliberation. Always makes me sad to see them shrugged away.

[Stefan Karpinski]

On Sun, Dec 28, 2014 at 9:09 AM, Jason Riedy <ja...@lovesgoodfood.com> wrote:

The choices in 754 were made with a nearly revolting amount of
deliberation.  Always makes me sad to see them shrugged away.

I can assure you, we have never shrugged away anything in IEEE. We've actually asked Kahan in person and by email about a number of things. And the fact that we've been discussing this so extensively can't really be considered a shrug.

Round-to-even is meant to allay bias problems for the vast
majority of programmers, like those responsible for the Vancouver
stock exchange index in 1982.  They used truncation which, over
time, roughly halved the index.  Rounding upward in the same
scenario would have doubled it...

I just don't see this particular argument as very compelling – that's one example where ties to even happens to fix a problem, accidentally as it happens. There seem to be just as many cases where other rounding modes are better. My main argument is two-fold: a) none of the choices of rounding modes for ties are obviously better across the board and b) following the C behavior seems the least surprising to the largest number of people.

But that's just my opinion. We are going with IEEE rounding as a default. Not only are we not shrugging off IEEE rules, we're going with them.