Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fuzzing test failure with TwicePrecision in test/ranges.jl #23497

Open
mbauman opened this issue Aug 29, 2017 · 4 comments
Open

Fuzzing test failure with TwicePrecision in test/ranges.jl #23497

mbauman opened this issue Aug 29, 2017 · 4 comments
Labels
ranges Everything AbstractRange

Comments

@mbauman
Copy link
Member

mbauman commented Aug 29, 2017

This fuzzing loop seems to sometimes find failures.

Found in the wild in https://julia.iblis.cnmc.tw/#/builders/1/builds/1617/steps/6/logs/stdio, reproduced locally by increasing the number of iterations.

Here are the magic numbers:

x, y = Base.TwicePrecision{Float32}(0.5767872f0, 1.9906205f-8), Base.TwicePrecision{Float32}(0.2913009f0, 1.4812992f-8)
xw, yw = asww(x), asww(y)
T = Float32
Tw = widen(T)
slopbits = 5

julia> cmp_sn2(Tw(xw/yw), astuple(x/y)..., slopbits)
false
@mbauman mbauman mentioned this issue Aug 29, 2017
Merged
@timholy
Copy link
Member

timholy commented Aug 29, 2017

In this method everything is exact until the computation of lo, which is off by a surprising margin:

julia> bits(widen(lo))
"0011111000101010110111101101001111100000000000000000000000000000"

julia> bits(((((widen(x.hi) - uh) - ul) + x.lo) - hi*y.lo)/y.hi)
"0011111000101010110111101101011101000110001010101000110110100101"

Unfortunately, it's going to be a while before I have time to figure this out (debugging this high-precision arithmetic is surprisingly time consuming). If people are desperate, just bump slopbits by 1 and it should pass.

@timholy
Copy link
Member

timholy commented Aug 29, 2017

FYI the 5 comes from 53 - 2*24, where 53 = number of precision bits of Float64 and 24 = number of precision bits of Float32. So 5 is taking the precision very, very seriously. I thought one can get away with doing that, but it's possible that the answer is no. It's also possible that some intermediate is generating a subnormal and that's the real reason.

@simonbyrne
Copy link
Contributor

I had a quick look, the problem is the line:

lo = ((((x.hi - uh) - ul) + x.lo) - hi*y.lo)/y.hi

In this case x.hi - uh == 0. The sum of the remaining part has cancellation issues which loses some precision.

@simonbyrne
Copy link
Contributor

So I found this:
https://hal.archives-ouvertes.fr/hal-01351529v2/document

I'm putting together a PR, should have something tomorrow.

@yuyichao yuyichao mentioned this issue Aug 31, 2017
Closed
This was referenced Sep 3, 2017
Merged
Merged
@brenhinkeller brenhinkeller added the ranges Everything AbstractRange label Nov 21, 2022
nsajko added a commit to nsajko/julia that referenced this issue Apr 29, 2023
@nsajko
Base: make TwicePrecision more accurate using standard algorithms
05121ef 
It seems that previously the compensated arithmetic code was all
designed ad-hoc instead of using the standard algorithms.

Some of the introduced code will also be used in `base/div.jl` to fix

I didn't check how much this improves the situation. In particular,
the example in JuliaLang#33677 still gives the same result, and I wasn't able
to evaluate JuliaLang#23497 because of how much Julia changed in the meantime.
@nsajko nsajko mentioned this issue Apr 29, 2023
Draft
nsajko added a commit to nsajko/julia that referenced this issue Apr 29, 2023
@nsajko
Base: make TwicePrecision more accurate using standard algorithms
979d8ec 
It seems that previously the compensated arithmetic code was all
designed ad-hoc instead of using the standard algorithms.

Some of the introduced code will also be used in `base/div.jl` to fix
issue JuliaLang#49450.

I didn't check how much this improves the situation. In particular,
the example in JuliaLang#33677 still gives the same result, and I wasn't able
to evaluate JuliaLang#23497 because of how much Julia changed in the meantime.
nsajko added a commit to nsajko/julia that referenced this issue May 4, 2023
@nsajko
Base: make TwicePrecision more accurate
3a64b91 
Update the arithmetic code to use algorithms published since
twiceprecision.jl got written.

Handle complex numbers.

Prevent harmful promotions.

Some of the introduced extended precision arithmetic code will also be
used in `base/div.jl` to tackle issue JuliaLang#49450 with PR JuliaLang#49561.

TODO: evaluate the accuracy improvements and write tests after JuliaLang#49616
gets merged

Results:

The example in JuliaLang#33677 still gives the same result.

I wasn't able to evaluate JuliaLang#23497 because of how much Julia changed in
the meantime.

TODO: check whether there are improvements for JuliaLang#26553

Updates JuliaLang#49589
nsajko added a commit to nsajko/julia that referenced this issue May 4, 2023
@nsajko
Base: make TwicePrecision more accurate
263de36 
Update the arithmetic code to use algorithms published since
twiceprecision.jl got written.

Handle complex numbers.

Prevent harmful promotions.

Some of the introduced extended precision arithmetic code will also be
used in `base/div.jl` to tackle issue JuliaLang#49450 with PR JuliaLang#49561.

TODO: evaluate the accuracy improvements and write tests after JuliaLang#49616
gets merged

Results:

The example in JuliaLang#33677 still gives the same result.

I wasn't able to evaluate JuliaLang#23497 because of how much Julia changed in
the meantime.

TODO: check whether there are improvements for JuliaLang#26553

Updates JuliaLang#49589
nsajko added a commit to nsajko/julia that referenced this issue May 4, 2023
@nsajko
Base: make TwicePrecision more accurate
5814884 
Update the arithmetic code to use algorithms published since
twiceprecision.jl got written.

Handle complex numbers.

Prevent harmful promotions.

Some of the introduced extended precision arithmetic code will also be
used in `base/div.jl` to tackle issue JuliaLang#49450 with PR JuliaLang#49561.

TODO: evaluate the accuracy improvements and write tests after JuliaLang#49616
gets merged

Results:

The example in JuliaLang#33677 still gives the same result.

I wasn't able to evaluate JuliaLang#23497 because of how much Julia changed in
the meantime.

TODO: check whether there are improvements for JuliaLang#26553

Updates JuliaLang#49589
nsajko added a commit to nsajko/julia that referenced this issue May 4, 2023
@nsajko
Base: make TwicePrecision more accurate
40294ae 
Update the arithmetic code to use algorithms published since
twiceprecision.jl got written.

Handle complex numbers.

Prevent some harmful promotions.

Some of the introduced extended precision arithmetic code will also be
used in `base/div.jl` to tackle issue JuliaLang#49450 with PR JuliaLang#49561.

TODO: write tests

TODO: some tests fail sometimes:
```
Test Failed at /home/nsajko/tmp/julia/test/ranges.jl:223
  Expression: cmp_sn2(Tw(xw / yw), astuple(x / y)..., slopbits)
```

Results:

The example in JuliaLang#33677 still gives the same result.

I wasn't able to evaluate JuliaLang#23497 because of how much Julia changed in
the meantime.

Proof that JuliaLang#26553 is fixed:
```julia-repl
julia> const TwicePrecision = Base.TwicePrecision
Base.TwicePrecision

julia> a = TwicePrecision{Float64}(0.42857142857142855, 2.3790493384824782e-17)
Base.TwicePrecision{Float64}(0.42857142857142855, 2.3790493384824782e-17)

julia> b = TwicePrecision{Float64}(3.4, 8.881784197001253e-17)
Base.TwicePrecision{Float64}(3.4, 8.881784197001253e-17)

julia> a_minus_b_inexact = Tuple(a - b)
(-2.9714285714285715, 1.0150610510858574e-16)

julia> BigFloat(a) == sum(map(BigFloat, Tuple(a)))
true

julia> BigFloat(b) == sum(map(BigFloat, Tuple(b)))
true

julia> a_minus_b = BigFloat(a) - BigFloat(b)
-2.971428571428571428571428571428577679589762353999797347402693647078382455095635

julia> Float64(a_minus_b) == first(a_minus_b_inexact)
true

julia> Float64(a_minus_b - Float64(a_minus_b)) == a_minus_b_inexact[begin + 1]
true
```

Fixes JuliaLang#26553

Updates JuliaLang#49589
nsajko added a commit to nsajko/julia that referenced this issue May 4, 2023
@nsajko
Base: make TwicePrecision more accurate
97914c8 
Update the arithmetic code to use algorithms published since
twiceprecision.jl got written.

Handle complex numbers.

Prevent some harmful promotions.

Some of the introduced extended precision arithmetic code will also be
used in `base/div.jl` to tackle issue JuliaLang#49450 with PR JuliaLang#49561.

Results:

The example in JuliaLang#33677 still gives the same result.

I wasn't able to evaluate JuliaLang#23497 because of how much Julia changed in
the meantime.

Proof that JuliaLang#26553 is fixed:
```julia-repl
julia> const TwicePrecision = Base.TwicePrecision
Base.TwicePrecision

julia> a = TwicePrecision{Float64}(0.42857142857142855, 2.3790493384824782e-17)
Base.TwicePrecision{Float64}(0.42857142857142855, 2.3790493384824782e-17)

julia> b = TwicePrecision{Float64}(3.4, 8.881784197001253e-17)
Base.TwicePrecision{Float64}(3.4, 8.881784197001253e-17)

julia> a_minus_b_inexact = Tuple(a - b)
(-2.9714285714285715, 1.0150610510858574e-16)

julia> BigFloat(a) == sum(map(BigFloat, Tuple(a)))
true

julia> BigFloat(b) == sum(map(BigFloat, Tuple(b)))
true

julia> a_minus_b = BigFloat(a) - BigFloat(b)
-2.971428571428571428571428571428577679589762353999797347402693647078382455095635

julia> Float64(a_minus_b) == first(a_minus_b_inexact)
true

julia> Float64(a_minus_b - Float64(a_minus_b)) == a_minus_b_inexact[begin + 1]
true
```

Fixes JuliaLang#26553

Updates JuliaLang#49589
@giordano giordano mentioned this issue Oct 7, 2024
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
ranges Everything AbstractRange
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@mbauman @timholy @simonbyrne @brenhinkeller