power uses Float64 exponents for integers

[KlausC]

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.

[KlausC]

@oscardssmith I would be interested in your opinion!

[oscardssmith]

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]

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

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

[oscardssmith]

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

[KlausC]

This screenshot shows the runtimes in v1.10.2.

[KlausC]

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]

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

[KlausC]

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.

[oscardssmith]

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.