logrange.jl

# Variants for https://github.com/JuliaLang/julia/pull/39071

# In this file the iterator versions are called "logrange" 
# and take Pair, since that's where I started,
# while random-access versions are called "geomrange" so that 
# they can share the serial number, and take 3 arguments.


####################### My first take in PR

function logrange1(p::Pair{<:Real, <:Real}, n::Integer)
    lo, hi = promote(p.first, p.second)
    if lo>0 && hi>0
        (exp(x) for x in range(log(lo), log(hi), length=n))
    elseif lo<0 && hi<0
        (-exp(x) for x in range(log(-lo), log(-hi), length=n))
    else
        throw(DomainError(p, "logrange requires that first and last elements are both positive, or both negative"))
    end
end

# Equally crude AbstractVector version

function geomrange1(lo::Number, hi::Number, n::Int)
    r = range(log(lo), log(hi), length=n)
    GeomRange1{eltype(r), typeof(r)}(r, n)
end
struct GeomRange1{T,R} <: AbstractVector{T}
    range::R
    length::Int
end
Base.size(g::GeomRange1) = (g.length,)
Base.@propagate_inbounds Base.getindex(g::GeomRange1, i::Integer) = exp(g.range[i])


##################### II my attempt at an iterator

struct LogRange2{T}
    start::T
    ratio::T
    length::Int
end

Base.length(r::LogRange2) = r.length

Base.iterate(r::LogRange2) = r.start, (r.start,1)

function Base.iterate(r::LogRange2, (x,n))
    if n >= r.length
        return nothing
    else
        y = x * r.ratio
        return y, (y, n+1)
    end
end

function logrange2(p::Pair{<:Real, <:Real}, n::Integer)
    n >= 2 || throw(ArgumentError(n, "logrange must have length at least 2"))
    lo, hi = p
    (lo>0 && hi>0) || (lo<0 && hi<0) || throw(DomainError(p,
        "logrange requires that first and last elements are both positive, or both negative"))
    rat = (hi/lo)^inv(n-1)

    LogRange2(promote(lo, rat)..., n)
end

#=
lo, hi, n = 1.0, 1000.0, 4
lo, hi, n = 0.7, 3000.0, 99

rat = (hi/lo)^inv(n-1)
rat = Float64((big(hi)/big(lo))^inv(big(n-1)))

rat = (big(hi)/big(lo))^inv(big(n-1))
# lo = big(lo)

hi == let
    tmp = lo
    for i in 2:n
        tmp *= rat
    end
    tmp |> Float64
end # false
=#

################################ III, Simeon Schaub's clever function

function _log(x) # = log(x) + dx * d log(x)/dx
    l = log(x)
    x_ = exp(l)
    return l - Base.TwicePrecision(x_ - x) / x_
end

function _exp(x::Base.TwicePrecision) # = exp(x) + dx * d exp(x)/dx
    e = exp(x.hi)
    # return Float64(e - (log(e) - x) * e)
    x_ = log(e)
    return (e - (x_ - x) * e).hi
end

function logrange3(lo, hi, n)
    log_r = LinRange(_log(lo), _log(hi), n)
    return (_exp(x) for x in log_r)
end

collect(logrange3(1/8, 8., 7))
# 7-element Array{Float64,1}:
#  0.125
#  0.24999999999999994
#  0.4999999999999999
#  1.0
#  1.9999999999999998
#  4.000000000000001
#  8.0

logrange3(p::Pair, n) = logrange3(p.first, p.second, n) # my interface

# AbstractVector version

function geomrange3(lo::Number, hi::Number, n::Int)
    log_r = LinRange(_log(lo), _log(hi), n)
    GeomRange3{eltype(log_r), typeof(log_r)}(log_r)
end

struct GeomRange3{T,R} <: AbstractVector{T}
    log_r::R
end

Base.size(g::GeomRange3) = size(g.log_r)
Base.@propagate_inbounds Base.getindex(g::GeomRange3, i::Integer) = _exp(g.log_r[i])


################################ IV, an iterator set up using BigFloats

struct LogRange4{T,S}
    start::T
    ratio::S
    length::Int
end

Base.length(r::LogRange4) = r.length

Base.iterate(r::LogRange4{T,S}) where {T,S} = r.start, (convert(promote_type(T,S), r.start),1)

function Base.iterate(r::LogRange4{T}, (x,n)) where T
    if n >= r.length
        return nothing
    else
        y = x * r.ratio
        return T(y), (y, n+1)
    end
end

# using DoubleFloats
const twice = Base.TwicePrecision

function logrange4(p::Pair{<:Real, <:Real}, n::Integer)
    n >= 2 || throw(ArgumentError(n, "logrange must have length at least 2"))
    lo, hi = p
    (lo>0 && hi>0) || (lo<0 && hi<0) || throw(DomainError(p,
        "logrange requires that first and last elements are both positive, or both negative"))
    # rat = (hi/lo)^inv(n-1)
    lo = float(lo)

    bigrat = (big(hi)/big(lo))^inv(big(n-1))
    # bigrat = (Double64(hi)/Double64(lo))^inv(Double64(n-1))
    r64 = Float64(bigrat)
    rat = Base.TwicePrecision(r64, Float64(bigrat - r64))

    # # Newton's
    # m = n-1
    # x0 = (hi/lo)^inv(m) |> twice # initial, x0
    # # (m-1) * x0/m
    # # (hi/lo) / (m * x0^(m-1))
    # x0pow = prod(x0 for _ in 1:m-1)
    # # twice(hi/lo) / (m * x0pow)
    # x1 = (m-1) * x0/m + twice(hi/lo) / (m * x0pow)

    # rat = x1

    LogRange4(lo, rat, n)
end

################################ @mkitti 's first suggestion

function logrange5(p::Pair{<:Real, <:Real}, n::Integer)
           lo, hi = promote(p.first, p.second)
           invert = abs(lo) > abs(hi)
           ratio = invert ? lo/hi : hi/lo
           if ratio % 10 == 0
               exp = exp10
               log = log10
           elseif ratio % 2 == 0
               exp = exp2
               log = log2
           else
               exp = Base.exp
               log = Base.log
           end
           if ratio > 0
               if invert
                   log_ratio = -log(ratio)
                   (x > log_ratio/2 ? lo*exp(x) : hi*exp(x-log_ratio) for x in range(0, log_ratio, length=n))
               else
                   log_ratio = log(ratio)
                   (x < log_ratio/2 ? lo*exp(x) : hi*exp(x-log_ratio) for x in range(0, log_ratio, length=n))
               end
           else
               throw(DomainError(p, "logrange requires that the first and last elements are both positive, or both negative"))
           end
       end


################################ VI, my fanciest TwicePrecision iterator
# as in PR, 10-13 Jan. 
# Note that this too-fancy constructor function has broken type-stability:
# @code_warntype logrange6(1 => 10^4, 1000)

import Base: TwicePrecision, iterate, length, eltype, ndims, size

function logrange6(start::Number, stop::Number;
    length::Union{Nothing,Integer}=nothing, ratio::Union{Nothing,Number}=nothing)
    if length !== nothing
        length >= 2 || throw(ArgumentError("logrange6 must have length of at least 2"))
        _start, _stop = map(float, promote(start, stop))
        _ratio = ratio_nth_root(stop, start, Int(length)-1)
        if ratio !== nothing
            check = oftype(float(ratio), _ratio)
            ratio == check || throw(ArgumentError("ratio = $ratio does not match $check computed from length = $length"))
        end
        LogRange6(_start, _ratio, Int(length))
    elseif ratio !== nothing
        _start, _ratio = promote(start, ratio)
        num = trunc(Int, log(ratio, stop/_start))
        fin = prod(_ratio for _ in 1:num; init=start)
        _length = fin <= stop ? num+1 : num
        LogRange6(_start, _ratio, _length)
    end
end

ratio_nth_root(a::Number, b::Number, n::Int) = (a/b)^(1/n)
ratio_nth_root(a::Float32, b::Float32, n::Int) = (Float64(a)/Float64(b))^(1/n)
function ratio_nth_root(a::Real, b::Real, m::Int)
    over = TwicePrecision(a) / TwicePrecision(b)
    r1 = TwicePrecision((over.hi)^(1/m))
    r1pow = prod(r1 for _ in 1:m-1)
    r2 = (m-1)*r1/m + over / (m * r1pow)
end

struct LogRange6{T,S}
    start::T
    ratio::S
    len::Int
end

length(r::LogRange6) = r.len
size(r::LogRange6) = (r.len,)

function iterate(r::LogRange6{T,S}) where {T,S}
    y = convert(promote_type(T,S), r.start)
    r.start, (y,1)
end

function iterate(r::LogRange6{T}, (x,n)) where {T}
    if n >= r.len
        nothing
    else
        y = x * r.ratio
        convert(T,y), (y,n+1)
    end
end

eltype(::Type{<:LogRange6{T}}) where {T} = T
eltype(r::LogRange6{T}) where {T} = T

ndims(::Type{<:LogRange6}) = 1
ndims(r::LogRange6) = 1

# same interfaces:
logrange6(p::Pair, l::Int) = logrange6(p.first, p.second; length=l)
logrange6(a::Number, b::Number, l::Int) = logrange6(a, b; length=l)


# AbstractVector version

function geomrange6(start::Number, stop::Number, length::Int)
    length >= 2 || throw(ArgumentError("logrange6 must have length of at least 2"))
    _start, _stop = map(float, promote(start, stop))
    _ratio = ratio_nth_root(stop, start, Int(length)-1)
    GeomRange6(_start, _ratio, Int(length))
end

struct GeomRange6{T,S} <: AbstractVector{T}
    start::T
    ratio::S
    len::Int
end

Base.length(r::GeomRange6) = r.len
Base.size(r::GeomRange6) = (r.len,)

function Base.getindex(r::GeomRange6{T}, ind::Int) where {T}
    # raw = r.start * r.ratio^(ind-1)
    raw = r.start * rec_power(r.ratio, ind-1)
    convert(T, raw)
end

# @mkitti's power function:

rec_power(x::Number, n::Int) = x^n
rec_power(tp::TwicePrecision{T}, n::Int) where {T} =
   if n == 0
       TwicePrecision{T}(1)
   elseif n == 1
       tp
   elseif n == 2
       tp*tp
   elseif n == 3
       tp*tp*tp
   else
       n2 = div(n,2)
       tpn2 = rec_power(tp, n2)
       if n2*2 < n
           tpn2 * tpn2 * tp
       else
           tpn2 * tpn2
       end
   end

################################ VII -- AbstractVector with DoubleFloats

using DoubleFloats

function geomrange7(lo::Number, hi::Number, n::Int)
    log_r = LinRange(log(Double64(lo)), log(Double64(hi)), n)
    GeomRange7{Float64, typeof(log_r)}(log_r)
end

struct GeomRange7{T,R} <: AbstractVector{T}
    range::R
end

Base.size(g::GeomRange7) = size(g.range)
Base.@propagate_inbounds Base.getindex(g::GeomRange7{T}, i::Integer) where {T} =
    convert(T, exp(g.range[i]))

################################ VIII --- Float64, with a both-ends correction idea

# This seems to hit the ends, but not every intermediate "nice" point.
# Faster than geomrange3, and I think slightly more accurate.
# With TwicePrecision this hits everything, but that has enough precision to work from one end anyway. 

struct GeomRange8{T} <: AbstractVector{T}
    start::T
    ratio::T
    error::T
    length::Int
end

Base.size(r::GeomRange8) = (r.length,)

@inline function Base.getindex(r::GeomRange8, j::Integer) 
    @boundscheck checkbounds(r, j)
    lo, rat, err, n = r.start, r.ratio, r.error, r.length
    y0 = rat^(j-1) * lo 
    dy = ((j-1) * err) / ((n-1) * rat^(n-j))
    y0 + dy
end

function geomrange8(lo::Number, hi::Number, n::Int)
    rat = (hi/lo)^(1/(n-1))
    err = hi - rat^(n-1) * lo
    GeomRange8(promote(lo, rat, err)..., n)
end

################################ Examples

collect(logrange1(1/8 => 8., 7)) # crude -- misses both ends
        geomrange1(1/8, 8.0, 7)
collect(logrange2(1/8 => 8., 7)) # ratio
collect(logrange3(1/8 => 8., 7)) # TwicePrecision-improved log/exp -- gets both ends
        geomrange3(1/8, 8.0, 7)
collect(logrange4(1/8 => 8., 7)) # big
collect(logrange5(1/8 => 8., 7)) # mkitti
collect(logrange6(1/8 => 8., 7))  # TwicePrecision + newton iterator
        geomrange6(1/8, 8.0, 7)
geomrange7(1/8, 8.0, 7) # Double64 log/exp, AbstractVector
geomrange8(1/8, 8.0, 7) # new correction scheme

collect(logrange2(3 => 1000, 17))[[begin,end]] # misses at end
collect(logrange3(3 => 1000., 17))[[begin,end]]
collect(logrange4(3 => 1000, 17))[[begin,end]]
collect(logrange5(3 => 1000, 17))[[begin,end]]
collect(logrange6(3 => 1000, 17))[[begin,end]]
geomrange7(3, 1000, 17)[[begin,end]]
geomrange8(3, 1000, 17)[[begin,end]]


collect(logrange3(7 => 123_000, 987))[[begin,end]]
collect(logrange4(7 => 123_000, 987))[[begin,end]]
collect(logrange5(7 => 123_000, 987))[[begin,end]]
collect(logrange6(7 => 123_000, 987))[[begin,end]]
geomrange7(7, 123_000, 987)[[begin,end]]
geomrange8(7, 123_000, 987)[[begin,end]]


collect(logrange3(1 => 10^10, 21))[1:2:7] # misses intermediates
collect(logrange4(1 => 10^10, 21))[1:2:5]
collect(logrange5(1 => 10^10, 21))[1:2:5]
collect(logrange6(1 => 10^10, 21))[1:2:5]
geomrange7(1, 10^10, 21)[1:2:7]
geomrange8(1, 10^10, 21)[1:2:7]



################################ Times

# Short range:

julia> @btime sum(logrange1(1 => 10^4, 10)) # naiive generator
  360.005 ns (0 allocations: 0 bytes)
15609.350234062036

julia> @btime sum(logrange2(1 => 10^4, 10))
  1.430 ns (0 allocations: 0 bytes)
15609.350234062023

julia> @btime sum(logrange3(1 => 10^4, 10)) # TwicePrecision-improved log/exp
  1.067 μs (0 allocations: 0 bytes)
15609.350234062025

julia> @btime sum(geomrange3(1, 10^4, 10)) # ... AbstractVector version, identitcal
  1.008 μs (0 allocations: 0 bytes)
15609.350234062025

julia> @btime sum(logrange4(1 => 10^4, 10))
  8.556 μs (27 allocations: 1.16 KiB)
15609.350234062025

julia> @btime sum(logrange5(1 => 10^4, 10))
  1.699 μs (76 allocations: 1.34 KiB)
15609.350234062025

julia> @btime sum(logrange6(1 => 10^4, 10)) # TwicePrecision + newton iterator
  304.464 ns (0 allocations: 0 bytes)
15609.350234062025

julia> @btime collect(logrange6(1 => 10^4, 10));
  336.271 ns (1 allocation: 160 bytes)

julia> @btime sum(geomrange6(1, 10^4, 10)) # ... AbstractVector version, multiplies per element
  474.893 ns (0 allocations: 0 bytes)
15609.350234062025

julia> @btime sum(geomrange7(1, 10^4, 10)) # Double64 log/exp, AbstractVector
  6.397 μs (0 allocations: 0 bytes)
15609.350234062025

julia> @btime sum(geomrange8(1, 10^4, 10)) # new correction scheme
  327.383 ns (0 allocations: 0 bytes)
15609.350234062027


# Longer range, 10^4 elements:

julia> @btime sum(logrange1(1 => 10^4, 1000))
  8.766 μs (0 allocations: 0 bytes)
1.0895501856939471e6

julia> @btime sum(logrange2(1 => 10^4, 1000))
  1.093 μs (0 allocations: 0 bytes)
1.0895501856939865e6

julia> @btime sum(logrange3(1 => 10^4, 1000))
  98.705 μs (0 allocations: 0 bytes)
1.0895501856939464e6

julia> @btime sum(logrange4(1 => 10^4, 1000))
  16.040 μs (27 allocations: 1.16 KiB)
1.0895501856939467e6

julia> @btime sum(logrange5(1 => 10^4, 1000))
  128.218 μs (7006 allocations: 109.62 KiB)
1.0895501856939464e6

julia> @btime sum(logrange6(1 => 10^4, 1000)) # TwicePrecision + newton iterator
  16.685 μs (5 allocations: 128 bytes)
1.0895501856939467e6

julia> @btime sum(collect(logrange6(1 => 10^4, 1000)))
  16.952 μs (5 allocations: 8.05 KiB)
1.0895501856939462e6

julia> @btime sum(geomrange6(1, 10^4, 1000)) # ... AbstractVector version, multiplies per element
  127.034 μs (0 allocations: 0 bytes)
1.0895501856939467e6

julia> @btime sum(geomrange7(1, 10^4, 1000)) # Double64 log/exp, AbstractVector
  627.692 μs (0 allocations: 0 bytes)
1.0895501856939467e6

julia> @btime sum(geomrange8(1, 10^4, 1000)) # new correction scheme
  38.212 μs (0 allocations: 0 bytes)
1.0895501856939467e6


# Even longer -- geomrange6 is still quicker than Double64, but iterating once & collecting is faster.

julia> @btime sum(geomrange3(1, 10^4, 10^7))
  932.320 ms (0 allocations: 0 bytes)
1.0856280226250076e10

julia> @btime sum(collect(logrange6(1 => 10^4, 10^7)))
  168.196 ms (6 allocations: 76.29 MiB)
1.0856280226249672e10

julia> @btime sum(geomrange6(1, 10^4, 10^7))
  3.446 s (0 allocations: 0 bytes)
1.0856280226250078e10

julia> @btime sum(geomrange7(1, 10^4, 10^7))
  6.460 s (0 allocations: 0 bytes)
1.0856280226250078e10

julia> @btime sum(geomrange8(1, 10^4, 10^7))
  387.815 ms (0 allocations: 0 bytes)
1.0856280226250076e10