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.

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

power uses Float64 exponents for integers #53967

Merged
merged 25 commits into from
Apr 22, 2024

Conversation

KlausC
Copy link
Contributor

@KlausC KlausC commented Apr 5, 2024

Improve performance of ^(::Float64, n::Integer) in the case of abs(n) > 2^13.

While pow_body is unreliable for abs(n) > 2^25 this implementation provides errors of a few ULPs, while runtime is capped to that of the Float64 implementation.

Fixes #53881
See also #53886.

@dkarrasch dkarrasch added the maths Mathematical functions label Apr 5, 2024
@KlausC
Copy link
Contributor Author

KlausC commented Apr 6, 2024

@oscardssmith I would be interested in your opinion!

@oscardssmith
Copy link
Member

this looks really good to me. I need to spend a little time benchmarking to make sure all the details look good, but I'm theory, this looks great!

@KlausC KlausC marked this pull request as draft April 9, 2024 09:05
base/math.jl Outdated Show resolved Hide resolved
base/math.jl Outdated Show resolved Hide resolved
base/math.jl Outdated Show resolved Hide resolved
base/math.jl Outdated Show resolved Hide resolved
base/math.jl Outdated Show resolved Hide resolved
@KlausC
Copy link
Contributor Author

KlausC commented Apr 13, 2024

Results of benchmark tests:
The abscissa is sign(n) * log2(abs(n) for integer exponents n.
The ordinate is in time per x^n for x close to but > 1.0.
The first graph show execution times in ns for ´x^n´ power_by_squaring, the second one for n -> exp(log(x)*n).
Hardware: Intel® Core™ i7-3610QM × 8

grafik

The result of the merged version with parameters nmin = -2^12, nmax = 3*2^13:

grafik

@KlausC

This comment was marked as off-topic.

@KlausC KlausC marked this pull request as ready for review April 13, 2024 11:44
@oscardssmith
Copy link
Member

note that the subnormal performance issue is only present on Intel cpus, AMD deals with subnormal much better.

@KlausC
Copy link
Contributor Author

KlausC commented Apr 13, 2024

This screenshot shows the runtimes in v1.10.2.

image

@KlausC
Copy link
Contributor Author

KlausC commented Apr 14, 2024

Statistics of accuracy evaluations. Note the very rare cases of 2 ULPs for low exponents (4 or 8), which can be explained by subnormal effects!
Evaluation functions defined in script.

julia> nx = [randpower(maxex = 2^63-1) for _ in 1:10^6];

julia> D = stat_ulps(nx);

julia> overview(D)
1000000
0 => 947081.0
1 => 52919.0
total: n: 1.0e6 - median: 0.028 - mean: 0.053 ± 0.224 ulps

julia> nx = [randpower(maxex = 2^12) for _ in 1:10^6];

julia> D = stat_ulps(nx);

julia> overview(D)
1000000
0 => 767308.0
1 => 232689.0
2 => 3.0
total: n: 1.0e6 - median: 0.152 - mean: 0.233 ± 0.422 ulps

julia> nx = [randpower(maxex = 20) for _ in 1:10^6];

julia> D = stat_ulps(nx);

julia> overview(D)
1000000
0 => 807276.0
1 => 192718.0
2 => 6.0
total: n: 1.0e6 - median: 0.119 - mean: 0.193 ± 0.394 ulps

julia> D[2]
6-element Vector{Any}:
 (4, -1.3443718769635203e-77)
 (-4, -8.590649363573623e76)
 (-4, 9.096283725387172e76)
 (8, -3.621102224996767e-39)
 (-8, -3.2458640793370833e38)
 (-8, 2.727732223112257e38)

# for example the first case in D[2]:

julia> n, x = (4, -1.3443718769635203e-77)
(4, -1.3443718769635203e-77)

julia> ulps(x, n)
2

julia> x^n
3.266462489987237e-308

@KlausC
Copy link
Contributor Author

KlausC commented Apr 16, 2024

The failed tests seem unrelated.
I am all set. Can someone review and merge if appropriate,

@KlausC
Copy link
Contributor Author

KlausC commented Apr 18, 2024

I have no idea, why test compiler/codegen uses ^ as a guinea-pig and also don't see how I should change the code in order to make this test succeed.

base/math.jl Outdated Show resolved Hide resolved
base/math.jl Outdated Show resolved Hide resolved
test/math.jl Outdated Show resolved Hide resolved
@oscardssmith
Copy link
Member

oscardssmith commented Apr 21, 2024

Other than the final 2 comments, I think is ready to merge. Also, thank you so much for fixing this! It's great to have a second person who understands this code.

@oscardssmith oscardssmith merged commit fe49d56 into JuliaLang:master Apr 22, 2024
4 of 7 checks passed
@KlausC KlausC deleted the krc/pow-float-path branch April 22, 2024 07:49
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
maths Mathematical functions
Projects
None yet
Development

Successfully merging this pull request may close these issues.

BUG: ^(::Float64, ::Union{Int, Float64}) incorrect for very large negative exponents
3 participants