statistics for accuracy and benchmarks of power function

statistics.jl

"""
```julia
using BenchmarkTools

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

julia> D = stat_ulps(nx);

julia> overview(D)
1000000
0 => 974262.0
1 => 25738.0
total: n: 1.0e6 - median: 0.013 - mean: 0.026 ± 0.158 ulps
```
"""
function ulps(x, n, p=^)
    a = p(x, n)
    b = big(x)^n
    u = abs(a - Float64(b)) / eps(a)
    return isfinite(u) ? BigInt(ceil(u)) : -1
end

function stat_ulps(nx::AbstractVector{<:Tuple{<:Integer,<:AbstractFloat}})
    D = Dict{Integer,Any}()
    try
        for (b, a) in nx
            u = ulps(a, b)
            R = get!(D, u) do
                []
            end
            t = (b, a)
            push!(R, t)
        end
    finally

    end
    for k in keys(D)
        D[k] = sort(collect(Set(D[k])); by=t -> abs(t[1]) + log2(abs(t[2]) + eps()) * eps())
    end
    return D
end

function randpower(; minex=3, maxex=typemax(Int64), alpha=1.0)
    minex = max(minex, 1)

    b = NaN
    while !(minex <= b <= maxex)
        lb = rand() ^ alpha * (log(maxex) - log(minex)) + log(minex)
        b = Int64(round(exp(lb)))
    end
    b *= rand([-1, 1])

    ab = NaN
    while !isfinite(ab)
        ab = reinterpret(Float64, rand(UInt64))
    end
    s = sign(ab)
    lab = log(abs(ab)) / b
    a = exp(lab) * s
    c = a^float(b)
    while !isfinite(c)
        a *= s
        b = b ÷ 2
        c = a^b
    end
    return b, a
end


function overview(D::Dict{<:Integer,<:Any}, maxulps=53)
    s0 = s1 = s2 = 0.0
    m = NaN
    ky = sort!(collect(keys(D)))
    n = sum(length(D[k]) for k in ky)
    println(n)
    smax = 0
    for k in ky
        l = float(length(D[k]))
        k > maxulps && (k = maxulps + 1)
        if s0 + l > n / 2 && isnan(m)
            m = k - 0.5 + (n / 2 - s0) / l
        end
        s0 += l
        s1 += l * k
        s2 += l * k * k
        if k <= maxulps
            println("$k => $l")
        else
            smax += l
        end
    end
    if smax > 0
        println(">$maxulps => $smax")
    end
    u = round(s1 / s0; digits=3)
    v = round(sqrt(s2 / s0 - u * u); digits=3)
    u, m, v = Float64.(round.((u, m, v); digits=3))
    println("total: n: $s0 - median: $m - mean: $u ± $v ulps")
end

function powerbench(;f=^, x=prevfloat(1.0), nrange=-62:62, seed=0)
    Random.seed!(seed)
    R = Dict()
    for i = nrange
        n = sign(i)*2^abs(i)
        R[n] = @benchmark $f($x,$n) samples=10000 evals=1000
    end
    for i = nrange
        n = sign(i)*(2^(abs(i)+1) - 1)
        R[n] = @benchmark $f($x,$n) samples=10000 evals=1000
    end
    for i = nrange
        n = sign(i)*(float(2^abs(i)) + 0.1)
        R[n] = @benchmark $f($x,$n) samples=10000 evals=1000
    end
    R
end

"""
produce and use the powerbench output:
```julia

R = powerbench()

A = [ones(63*4) [zeros(63);ones(63);zeros(63);ones(63)] [1:63;1:63;1:63;1:63] [zeros(63*2);1:63;1:63]];

b = [[time(R[2^i]) for i in 0:62]; [time(R[-2^i]) for i in 0:62]; [time(R[2*2^i-1]) for i in 0:62]; [time(R[-(2*2^i-1)]) for i in 0:62]];

x0, xn, x2, x1 = A \\ b

esttime(n) = x0 + xn * (n < 0) + x2 * (top_set_bit(n)-1) + x1 * (count_ones(n)-1)

```
"""

function plottimes(R::Dict; nrange=-62:62)
    p = plot()

    function between(n)
        a = sign(first(nrange))*2^abs(first(nrange))
        b = sign(last(nrange))*2^abs(last(nrange))
        a <= n <= b
    end
    nn = sort!(collect(filter(n -> isinteger(n) && between(n), keys(R))))
    xx = [sign(n) * log2(abs(n)) for n in nn]
    yy = [time(R[n]) for n in nn]
    plot!(p, xx, yy)
    nn = sort!(collect(filter(n -> !isinteger(n) && between(n), keys(R))))
    xx = [sign(n) * log2(abs(n)) for n in nn]
    yy = [time(R[n]) for n in nn]
    plot!(p, xx, yy)
    display(p)
end