Strange interaction of rounding mode and `rationalize` on 0.4

[dpsanders]

I have just found the following strange interaction of the rounding mode and rationalize on 0.4 (commit 9ad3aa0, 9 days old master) -- the result depends rather strongly on rounding mode:

julia> set_rounding(Float64, RoundNearest)
0

julia> rationalize(0.1)
1//10

julia> set_rounding(Float64, RoundUp)
0

julia> rationalize(0.1)
1//11

julia> set_rounding(Float64, RoundDown)
0

julia> rationalize(0.1)
1//9

cc @simonbyrne

[simonbyrne]

Yeah, this is a problem with rounding modes. Rounding modes only apply to basic arithmetic operations (+, -, *, / and sqrt). Most other functions are written assuming round-to-nearest mode, and so will give funny results. In order to do this correctly, we would have to check and modify the rounding mode on each function call, but that would impose significant overhead.

This might improve in future AVX-512, which is supposed to allow for operation-level rounding modes, but until that point, I think we should just document that "there be dragons".

The irony in this case is that the problem here is the div function with floating point arguments: ideally, we would implement this as

function div(x::Float64,y::Float64)
    q = with_rounding(Float64,RoundToZero) do
        x / y
    end
    trunc(q)
end

But that is significantly slower than our current (more complicated) method.

[dpsanders]

The current behaviour seems to be a regression from 0.3, where the result is 1//10 for all of the above rounding examples. However, the code from 0.3, when run on 0.4, also seems to give the above incorrect results. What is it that has changed between 0.3 and 0.4? As far as I can see, the version of rationalize from 0.3.7 uses only / and convert.

Can I define my own version of rationalize that will work also on 0.4 with non-nearest rounding? (Even if it is a bit slower.) I need this for the ValidatedNumerics package, which is where I noticed the problem. There is a branch always_round_up which tries to get around the problem of constantly changing the rounding mode by just fixing [arbitrarily] RoundUp and using prevfloat to get the corresponding RoundDown behaviour. Does this sound reasonable / possible?

[dpsanders]

My mistake, the code from 0.3 does seem to work correctly on 0.4 with all rounding modes. I guess I will just go with that...

[simonbyrne]

You almost certainly don't want to change the global rounding mode without setting it back: so many functions expect round-to-nearest that you are going to get very odd results. Debugging using non-standard rounding modes is not fun (see, for example, #9847)

I had a quick skim though your package: it looks quite neat. My suggestion, if you're using 0.4, would be to avoid using with_rounding at the higher level, and instead exploit the ability to dispatch on RoundingMode by defining different methods, e.g.

foo(x, ::RoundingMode{:Up}) = ...
foo(x, ::RoundingMode{:Down}) = ...
foo(x, ::RoundingMode{:Nearest}) = ...

For the basic arithmetic operations, you can then define these using with_rounding, but for other functions you can define specific methods to call.

[JeffBezanson]

Is there a better way to write floating point div? We're currently using a pretty naive generic algorithm.

[simonbyrne]

Assuming round-to-nearest is set, I think we're provably accurate as long as abs(x/y) < 2^51, but we could improve that up to 2^53 (we can't go any higher as we can't represent all integers beyond that point). @Yaffle had some suggestions in https://github.com/JuliaLang/julia/issues/4156#issuecomment-72633282.

[dpsanders]

@simonbyrne Thanks for your comments, and for taking a look at the package.

The idea is that any functions used while doing interval arithmetic will be those specialized to be used for intervals (a work in progress... trig functions are particularly gnarly), which should be designed to respect rounding. I like the idea of dispatching on rounding mode, but currently we're trying to maintain 0.3 compatibility. No doubt once 0.4 comes out, we'll quietly drop that...

The point of just fixing RoundUp once and for all is to simply avoid the constant changing of rounding mode. On one particular micro-benchmark, this seems to be around a factor of 20 faster (!). Another possibility would be to stick with RoundNearest and just widen the interval in both directions (both prevfloat and nextfloat), though this is unnecessarily wide.

On this note: am I correct that with round-to-nearest, there is no way to tell in which direction the rounding occurred (up or down)? (Although I understand that whether rounding occurred can [sometimes...] be checked with the JL_FE_INEXACT exception from your #6170 ?)

[simonbyrne]

am I correct that with round-to-nearest, there is no way to tell in which direction the rounding occurred (up or down)

Not directly, no, though you could use double-double operations which would tell you how much rounding occurred.

As for the trig functions, you could try wrapping crlibm?

[dpsanders]

OK thanks, I'll have a go with double-double.

Actually I already have a crlibm branch where I started that task...! (And apparently that requires that the rounding mode be left in RoundNearest.)

[simonbyrne]

Okay, I think we can close this as "can't fix".

[aaljuffali]

The most accurate way is to input the number in string format and rationalize it directly without going to floating point. I am new to Julia but I have once wrote a function in Python to that. It is mostly string manipulation and integer arithmetic.

[pao]

@aaljuffali There's one intermediate step, but that's pretty close to what happens when you use it with BigFloat:

julia> rationalize(3.1415926535897932384626)
165707065//52746197

julia> rationalize(big"3.1415926535897932384626")
8111581125742254341//2581996464905601587

[aaljuffali]

The reason for rationalizing is to have exact representation. Big or FloatsXX cannot represent some floating point numbers unless you go with a true Decimal FP representation. It is available in other languages. In Python it is called "Decimal" part of the standard library. It is available for C for free. I wish Julia will have it one day.

[tkelman]

@aaljuffali you're in luck: https://github.com/stevengj/DecFP.jl

[simonbyrne]

@aaljuffali There is also Decimals.jl, which is perhaps closer in spirit to the Python library, being arbitrary precision.

[aaljuffali]

Thank you all. I have installed Decimals.jl and wrote a str2rat function for who ever needs to use it for exact representation . You might need to adapt it to your application and do some error checking. This is to give you a head start. I am not responsible for any errors!!!

using Decimals

-------------------------------------------

Function to convert a string to rational.

-------------------------------------------

function str2rat(s::ASCIIString)

decimalNumber = decimal(s) signInt = decimalNumber.s mantissaInt = decimalNumber.c exponentInt = decimalNumber.q

if signInt == 0 ansRational = mantissaInt // 1 else ansRational = -mantissaInt // 1 end # signInt

if exponentInt < 0 ansRational = ansRational * (1 // 10^-exponentInt) else ansRational = ansRational * (10^exponentInt // 1) end # exponentInt

return ansRational

end

[dpsanders]

@aaljuffali Your solution is a good one if starting from a string. However, in my ValidatedNumerics package, the idea is that you can write

@interval(2.3*3.2)

and it will return an interval containing the true result of the operation 2.3*3.2

For this, I currently use rationalize applied to a Float64. If you have an alternative suggestion, I would be happy to hear it. I do not wish to convert the Float64s to strings, since I have the feeling that this would be slow (although I have not done a proper benchmark).

[aaljuffali]

Here is a way that might be satisfactory:

1- assume you have two numbers to multiply A and B 2- Choose a factor T so that it is 10^x where x is the number of decimal places you need. 3- Do C=round(A_T) and D=round(B_T). Round each to the nearest integer. 4- Multiply the two integers. 5- Divide the result by T^2.

for example:

julia> 2.3*3.2 7.359999999999999

julia> A=2.3; B=3.2; T=10.0; C=round(A_T); D=round(B_T); C*D/T^2 7.36

You can also create a Rational number like: int(C)//int(T) and int(D)//int(T) to preserve them from further misrepresentation.

[simonbyrne]

The problem with this is that in Julia the conversion is done at parse time, in other words, by the time Julia sees it, 2.3 is already 2.29999999999999982236431605997495353221893310546875 (which is why in v0.4, you will get 2.3 < 23//10 being true). So in essence, you're trying to undo the parser.

Perhaps we do need some way to represent decimals in native Julia.

[kmsquire]

Perhaps we do need some way to represent decimals in native Julia.

d"2.3"? This is what DecFP.jl does.

Any other alternative would require modifying the parser. Might be a reasonable thing to do--e.g., following the Float32 convention, 2.3d0 could be parsed as a decimal number. (Of course, you wouldn't be able to parse this as 2.3 * d0 then.)

[ScottPJones]

Perhaps we do need some way to represent decimals in native Julia.

@simonbyrne I've tried to suggest that a couple of times already, and have been shot down... @kmsquire I don't know enough about the internals yet (pardon me for not having had the time to thoroughly research this before commenting!) but would it be possible for the parser to leave floating point literals alone, as strings, until it really needs to decide that it needs a binary floating point value or something else? Having the Scheme parser convert it to a Float64 immediately seems a bit premature, IMO. I realize this might be a somewhat breaking change (at least in the internals), but I think it might be worthwhile [even after having lived through the Tuplecalypse :wink:]

[simonbyrne]

@ScottPJones I'm sure that there are plenty more -calypse's to come (we're still only on v0.4 after all).

If you want to learn more about Julia's internals, a good place to start is Jeff's talk from JuliaCon 2014. The trouble with allowing operations on literals is that you could end up with things like x = 0.1; foo(x) and foo(0.1) giving different answers. One possible alternative might be to delay conversion until after the macro expansion stage, similar to how big integers are parsed, which would mean that macros are able to access the literals, but not functions.

[aaljuffali]

@simonbyrne You are 100% correct. Try to convince Dr. Sanders that he will not be able to have exact number representation for some numbers using binary floating point. This is why I was trying to move him to use Rationals to preserve the number throughout the calculations.

[aaljuffali]

@kmsquire DecFP.jl looks promising but it it requires Julia V0.4 and once V0.4 is released I will look at it. It needs more documentation the way it stands now.

[dpsanders]

@simonbyrne Exactly, the idea is to work around the decision (or limitation) of the parser to immediately convert 2.3 into its floating-point version. Currently I convert 2.3 into its rational representation, mainly for mathematical elegance (rather than speed).

Basically, currently the following is done in ValidatedNumerics.jl (approximately): @interval(2.3*3.4) is replaced by make_interval(23//10)*make_interval(34//10), which is converted into interval(23)/interval(10)*interval(34)/interval(10), in which each object is an interval. Directed rounding is used in the division of intervals, so that we get

julia> @floatinterval(2.3*3.2)
[7.359999999999997, 7.360000000000002]

which indeed contains the correct result (7.36).

It would be great to pass through the literals at least to the macro level.

@kmsquire I much prefer to write 2.3*3.2 than d"2.3" * d"3.2". Having said that, you can also do

julia> @floatinterval("2.3"*"3.2")
[7.3599999999999985, 7.360000000000002]

which gives a slightly different answer, but still contains the correct result.

@aaljuffali I do not need "convincing" of anything -- I understand perfectly well how binary floating point works, thank you. I hope my explanation of what I'm actually doing helps. Hopefully there will be more decent documentation of the ValidatedNumerics package soon...

[ScottPJones]

@dpsanders Maybe you could put some (gentle) pressure on the core Julia team :grinning: about passing literals to the macro level, and also having decimal arithmetic as part of the base language (of course, after it's been well tested in @stevengj 's DecFP.jl package). I've already been pushing them in this direction, it would be good for them to see it's not just "that PITA Scott" who wants this!

[StefanKarpinski]

This isn't about "pressure" or lack of willingness. Having literals be of a different type is simply not feasible for reasons that have been explained over and over again. If someone can propose a way to do this that doesn't introduce type instabilities all over the place, that would be wonderful, but no one has ever come up with a workable proposal. @simonbyrne's idea of delaying it until after macro expansion is on the right track but needs significant fleshing out.

[ScottPJones]

@StefanKarpinski I thought that's what I was saying, delaying it until after macro expansion... we'd discussed that on another thread earlier. Why would that cause type instabilities? By having decimal arithmetic as part of the base language, I meant just having Decimal32, Decimal64, Decimal128, and BigDecimal types available in base, and the syntax that @stevengj as already added for DecFP.jl for decimal literals.

[StefanKarpinski]

The most straightforward way to do this is just to parse a floating-point literal like 1.23 as @float64_str "1.23" and then providing a default implementation of @float64_str that parses the floating-point value and produces a Float64. I'm not sure if there are bootstrapping issues with this approach. It also doesn't in any way require having any decimal type in base Julia. One potential concern is that this might blow up the size of ASTs, but I don't think there are actually enough floating-literals to make that an issue.

[ScottPJones]

Having the decimal type in base Julia is a separate issue... it's absolutely a requirement in many areas, and doesn't add much code at all. I'll give exact numbers next week.

[StefanKarpinski]

How can it possibly be "absolutely a requirement" to have a decimal type in base? Why can't it be a package?

[ScottPJones]

It is absolutely a requirement for many types of applications... you misread/misquoted me. It could be a package... but deserves to be in base (much more than a lot of stuff that is in there, IMO)

[aaljuffali]

@StefanKarpinski The minute you convert to a binary Float64 by using the internal parser or your own you will loose the true representation of sum numbers like 2.3. See example below:

julia> a = 2.3 2.3 julia> bits(a) "0100000000000010011001100110011001100110011001100110011001100110"

look at the end of binary float64. You will see repeated 1's and 0's and 2.3 is actually 2.299... We have not even began any calculation. The more calculation you do the more significant digits you loose. This issue is not tied to Julia but to all programing languages out their except the ones that use decimal floating point. I hope this explains it.

[StefanKarpinski]

Is that how that works? I've always wondered.

[pao]

I think we need to move on from this issue; we're no longer talking about either rounding modes or rationalize.

@aaljuffali These issues are well understood, sarcasm aside. (If you do a little deep code spelunking, you'll see that Stefan wrote the bits function that you used, and is probably just getting frustrated at the direction of conversation.) Those issues are also not terribly important for a large class of problems, as long as you know that it is happening. And conversely, there are rational numbers exactly representable in fixed-precision binary floating point that cannot be represented exactly in fixed-precision decimal floating point.

[tkelman]

By having decimal arithmetic as part of the base language, I meant just having Decimal32, Decimal64, Decimal128, and BigDecimal types available in base

Please open a new issue for this suggestion instead of continuing to discuss it here, so we can elaborate in more detail on the many reasons to oppose bringing these features into base as a required dependency for the entire language.

Allowing the parser to be extensible so this functionality can be more cleanly implemented in a package is yet another separate issue and should be discussed in one. This idea is much less likely to encounter resistance, however

I'm not sure if there are bootstrapping issues with this approach.

will need to be discussed (not here).

[dpsanders]

@simonbyrne I finally got round to implementing your suggestion of dispatching on the different rounding-mode types:

https://github.com/dpsanders/ValidatedNumerics.jl/blob/crlibm/crlibm_examples/CRlibm.jl

There seems to be a strange issue whereby the method being defined with the next rounding mode replaces the previous one, instead of being added to the method list. Do you have any suggestion as to what's going on?

EDIT: Defining them explicitly (see bottom of file), instead of using metaprogramming, does work.

EDIT2: Seems to be fixed now, sorry for the noise. Metaprogramming is always slippery.

[0joshuaolson1]

It would be helpful to know where these separate issues are being discussed.

Pretty please?

[pao]

@MrMormon, what have you tried searching for?