Absolutely, positively perfect ranges. Never wrong. No way. #23194

timholy · 2017-08-09T22:54:32Z

This high-precision arithmetic stuff stinks. But fortunately no one will ever, ever report a bug in this code again. And if they do, someone (probably me) is going to be swimming with the fishes. Got it?

tkelman · 2017-08-10T07:49:34Z

@nanosoldier runbenchmarks(ALL, vs=":master")

nanosoldier · 2017-08-10T11:42:25Z

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @ararslan

tkelman · 2017-08-10T12:14:24Z

test/ranges.jl

-    @test m == length(r)
-    # FIXME: these fail some small portion of the time
-    @test_skip start == first(r)
-    @test_skip stop  == last(r)


guess this was over-ambitious?

#20373 (comment)

However, I think (not sure, actually) that we should be able to insist on

Δ > 0 ? last(r) <= stop : last(r) >= stop

But we're not quite there yet. Needs more investigation.

timholy · 2017-08-10T12:30:25Z

I'm going to be traveling for a couple of days, so let me say please don't merge yet. I think another formal property we need to insist on is that

r = start:Δ:stop
range(start, Δ, length(r)) === r

and this does not satisfy that. I don't think any of the changes here are wrong, so reviews are welcome, but I think there's a little more fixing/cleanup that needs doing.

timholy · 2017-08-16T09:41:34Z

OK, this turned into a bigger deal than I wanted, but I think it's a pretty big improvement. Perhaps the most important part of this change, aside from the two bugs I fixed, was adding a whole bunch of carefully-crafted precision tests for our TwicePrecision routines. This led down a long rabbit hole, but I made a number of improvements and learned some interesting things along the way.

The most important of these is that the well-known Dekker TwicePrecision (aka, "double-double") algorithms fail miserably when one of the intermediate computations produces a subnormal. This led to a number of changes:

a compensated division routine that ensures that div12 won't produce intermediate subnormals
the realization that TwicePrecision basically makes no sense for Float16: TwicePrecision{Float16} arithmetic routines produce subnormals at the drop of a hat (there just aren't enough exponent bits in Float16), and so it's far, far better to just use Float32 or Float64 if you want something that produces higher-precision intermediate results.

As a consequence, this switches our StepRangeLen types to (by default) only use TwicePrecision for Float64 ranges; for Float32 and Float16, the ref and step fields will be Float64.

I've punted on making any of the changes discussed starting at #20373 (comment). True perfection will have to wait for a future pull request.

I removed the "backport" label, since this is by now far too intrusive of a change to backport; I can create a separate patch for fixing #23178 on 0.6 after this merges.

timholy · 2017-08-16T12:44:55Z

Travis and appveyor failures seem unrelated.

StefanKarpinski · 2017-08-16T19:13:57Z

Funny that, I just cooked up a compensated multiplication procedure the other night when I was thinking about this during a bout of insomnia:

minabs(x, y) = ifelse(abs(x) ≤ abs(y), x, y)
maxabs(x, y) = ifelse(abs(x) > abs(y), x, y)

function compensated_addition(x::T, y::T) where T<:AbstractFloat
    z = x + y
    z, minabs(x, y) - (z - maxabs(x, y))
end

function compensated_multiplication(x::Float64, y::Float64)
    m, p = Base.decompose(x)
    n, q = Base.decompose(y)
    x*y, (UInt64(m)*UInt64(n) & 0x000fffffffffffff)*2.0^(p + q)
end

Chaining compensated values is a real pain, it's almost as if one wants a type for that 😸

timholy · 2017-08-16T21:31:00Z

Your compensated_addition is exactly add12 in this PR (and what's mistakenly called add2 in current master).

I'm not sure I quite understand what's happening in the multiplication routine, since I think the second isn't a correction of the first. If that's your intention, I think it's well-accepted that at this point the best way to achieve that with floating-point is to use the fact that fma has sufficient precision and is incredibly fast (on supported hardware), see this PR. (The ifelse is necessary only to get annoying things like Inf, NaN, and signed zeros correct.)

StefanKarpinski · 2017-08-16T22:27:52Z

I'm not sure I quite understand what's happening in the multiplication routine, since I think the second isn't a correction of the first.

It does work though:

julia> x, y = 0.1, 0.3
(0.1, 0.3)

julia> compensated_multiplication(x, y)
(0.03, 1.6653345369377347e-18)

julia> big(x)*big(y) - x*y
1.665334536937734749005689281744683170958705837282325806780747257107577752321959e-18

How it works: x == m*2.0^p and y == n*2.0^q so x*y == (m*n)*2.0^(p + q). Floating-point multiplication gives you the first 52 bits of m*n, you just need the last 52 bits, which you can get by doing integer multiplication, which gives the right answer mod 2^64, and you can then use & 0x000fffffffffffff to do a fast mod 2^52 operation and get those trailing bits. Then just multiply by the appropriate power of two and you've got floating-point correction to the product.

I guess using fma would work too, but I wasn't familiar with that so this is what I came up with :)

simonbyrne

I'm still a bit confused why we need Veltkamp splitting. What was wrong with the old truncbits method?

simonbyrne · 2017-08-16T22:51:09Z

base/twiceprecision.jl

+Arithmetic", §4.4.1.
+"""
+function splitprec(x::T) where {T<:AbstractFloat}
+    xs, xe = frexp(x)


frexp/ldexp is probably not the most efficient way to do this. It might be simpler to just check if its magnitude is above some threshold, and scale it down if necessary

simonbyrne · 2017-08-16T22:53:58Z

base/float.jl

@@ -825,28 +829,7 @@ fpinttype(::Type{Float64}) = UInt64
 fpinttype(::Type{Float32}) = UInt32


I'm still not a huge fan of this name: I keep reading this as "f pint type". It seems like there should be a more general name for this (i.e. give me the unsigned integer type with the same size as this type).

We could do typeof(reinterpret(Unsigned, zero(F))) but that seems a bit silly.

I changed it to uinttype.

simonbyrne · 2017-08-16T23:18:37Z

base/float.jl

@@ -1,5 +1,7 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license

+const IEEEFloat = Union{Float16, Float32, Float64}


I realise that it is beyond the scope of this PR, but it would be nice if this were an abstract type so I could reuse all this for Float128 once #22649 is fixed.

I just moved this from Base.Math. I'll leave that kind of change to you 😄. Should be pretty simple, though.

simonbyrne · 2017-08-16T23:40:16Z

@StefanKarpinski I think that only works if your calculation rounds down (or if you've changed the rounding mode). e.g. try

julia> h,l = compensated_multiplication(0.11,0.3)
(0.033, 8.604228440844963e-19)

julia> big(0.11)*big(0.3) -h
-2.609024107869117876158510242052582464154129416271767419321925274289242224767804e-18

StefanKarpinski · 2017-08-17T02:42:43Z

Yeah, might need some sign fiddling in some cases.

timholy · 2017-08-17T07:59:52Z

Clever, thanks for explaining.

timholy · 2017-08-17T10:31:52Z

I'm still a bit confused why we need Veltkamp splitting. What was wrong with the old truncbits method?

We don't really; I think it was left over from an older version that still used some splitting. (With the old truncbits method, the products hi1*hi2, hi1*lo2, and lo1*hi2 were all exact, but lo1*lo2 was not, so Veltkamp splitting is better.) I just deleted the method, though. The only actual use of splitprec was the convert-an-integer-exactly variant.

Thanks for the review!

simonbyrne · 2017-08-17T16:33:53Z

base/float.jl

@@ -711,12 +713,12 @@ end

 Test whether a floating point number is subnormal.
 """
-function issubnormal end
+function issubnormal(x::T) where {T<:IEEEFloat}
+    y = reinterpret(uinttype(T), x)


can use reinterpret(Unsigned, x)

timholy · 2017-08-18T09:00:14Z

The last commit revealed a strange bug on 32-bit in -5.0:1.0:5.0 === -5.0:1.0:5.0. Not sure I understand why deserialization/hashing is being invoked. @vtjnash? (ref. d584142 which added a Int32 type annotation in deserialize_type)

timholy · 2017-08-20T14:55:03Z

This seems good to go. I'll merge soon barring concerns.

mbauman · 2017-08-29T15:38:56Z

test/ranges.jl

+    @test cmp_sn2(Tw(xw+yw), astuple(x+y)..., slopbits)
+    @test cmp_sn2(Tw(xw-yw), astuple(x-y)..., slopbits)
+    @test cmp_sn2(Tw(xw*yw), astuple(x*y)..., slopbits)
+    @test cmp_sn2(Tw(xw/yw), astuple(x/y)..., slopbits)


Never, never any bugs, right? ;) Sorry to be the bearer of bad news, but I imagine you'll be interested in #23497.

By definition, if you think you've seen a bug in this code, you're just imagining things. Really. Honest. No bugs.

Ahhh! (sploosh)

timholy added the backport pending 0.6 label Aug 9, 2017

timholy mentioned this pull request Aug 10, 2017

float ranges: if a+n*s == b, length(a:s:b) should be n+1 #20373

Closed

tkelman reviewed Aug 10, 2017

View reviewed changes

timholy force-pushed the teh/fix_23178 branch from f371fa9 to 890ef13 Compare August 16, 2017 09:17

timholy removed the backport pending 0.6 label Aug 16, 2017

timholy requested review from simonbyrne and StefanKarpinski August 16, 2017 09:41

simonbyrne reviewed Aug 16, 2017

View reviewed changes

timholy force-pushed the teh/fix_23178 branch from 890ef13 to 6274ac3 Compare August 17, 2017 10:31

simonbyrne reviewed Aug 17, 2017

View reviewed changes

timholy force-pushed the teh/fix_23178 branch from 6274ac3 to 5d07a22 Compare August 17, 2017 17:24

timholy mentioned this pull request Aug 17, 2017

map(Float64, ::StepRangeLen{BigFloat,BigFloat,BigFloat}) has changed since v0.5? #23300

Closed

This was referenced Aug 18, 2017

Test strange deserialzation bug on 32bit #23330

Closed

Fix #23300 #23348

Merged

timholy force-pushed the teh/fix_23178 branch from c94ed9f to 2efb0e6 Compare August 19, 2017 14:11

timholy added 2 commits August 20, 2017 03:02

Better TwicePrecision utilities

67a8605

Fix overflow in linspace construction. Fixes #23178

b07b94d

timholy added 2 commits August 20, 2017 03:02

Hit linspace endpoints exactly and fix range tests. Fixes #20521.

a903c87

Ranges: use Float64 rather than TwicePrecision for Float16, Float32

a698230

timholy force-pushed the teh/fix_23178 branch from 2efb0e6 to a698230 Compare August 20, 2017 10:56

timholy merged commit e6c8b6b into master Aug 22, 2017

timholy deleted the teh/fix_23178 branch August 22, 2017 11:59

mbauman reviewed Aug 29, 2017

View reviewed changes

timholy mentioned this pull request Sep 15, 2019

Define interval arithmetics using broadcasting JuliaMath/IntervalSets.jl#55

Open

timholy mentioned this pull request Oct 25, 2019

Revisit the algorithm used to compute ranges with floating points #33677

Open

timholy mentioned this pull request Jan 31, 2023

More stringent range testing #48465

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Absolutely, positively perfect ranges. Never wrong. No way. #23194

Absolutely, positively perfect ranges. Never wrong. No way. #23194

timholy commented Aug 9, 2017

tkelman commented Aug 10, 2017

nanosoldier commented Aug 10, 2017

tkelman Aug 10, 2017

timholy Aug 10, 2017

timholy commented Aug 10, 2017

timholy commented Aug 16, 2017

timholy commented Aug 16, 2017

StefanKarpinski commented Aug 16, 2017 •

edited

Loading

timholy commented Aug 16, 2017

StefanKarpinski commented Aug 16, 2017

simonbyrne left a comment

simonbyrne Aug 16, 2017

simonbyrne Aug 16, 2017

timholy Aug 17, 2017

simonbyrne Aug 16, 2017

timholy Aug 17, 2017

simonbyrne commented Aug 16, 2017

StefanKarpinski commented Aug 17, 2017

timholy commented Aug 17, 2017

timholy commented Aug 17, 2017

simonbyrne Aug 17, 2017

timholy commented Aug 18, 2017

timholy commented Aug 20, 2017

mbauman Aug 29, 2017

timholy Aug 29, 2017 •

edited

Loading

		@@ -825,28 +829,7 @@ fpinttype(::Type{Float64}) = UInt64
		fpinttype(::Type{Float32}) = UInt32

		@@ -1,5 +1,7 @@
		# This file is a part of Julia. License is MIT: https://julialang.org/license

		const IEEEFloat = Union{Float16, Float32, Float64}

Absolutely, positively perfect ranges. Never wrong. No way. #23194

Absolutely, positively perfect ranges. Never wrong. No way. #23194

Conversation

timholy commented Aug 9, 2017

tkelman commented Aug 10, 2017

nanosoldier commented Aug 10, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

timholy commented Aug 10, 2017

timholy commented Aug 16, 2017

timholy commented Aug 16, 2017

StefanKarpinski commented Aug 16, 2017 • edited Loading

timholy commented Aug 16, 2017

StefanKarpinski commented Aug 16, 2017

simonbyrne left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

simonbyrne commented Aug 16, 2017

StefanKarpinski commented Aug 17, 2017

timholy commented Aug 17, 2017

timholy commented Aug 17, 2017

Choose a reason for hiding this comment

timholy commented Aug 18, 2017

timholy commented Aug 20, 2017

Choose a reason for hiding this comment

timholy Aug 29, 2017 • edited Loading

Choose a reason for hiding this comment

StefanKarpinski commented Aug 16, 2017 •

edited

Loading

timholy Aug 29, 2017 •

edited

Loading