logrange (#821)

Author: mcabbott

Project.toml

@@ -1,6 +1,6 @@
 name = "Compat"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.13.0"
+version = "4.14.0"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"

README.md

@@ -72,6 +72,8 @@ changes in `julia`.
 
 * `allequal(f, itr)` and `allunique(f, itr)` methods. ([#47679]) (since Compat 4.13.0)
 
+* `logrange(lo, hi; length)` is like `range` but with a constant ratio, not difference. ([#39071]) (since Compat 4.14.0) Note that on Julia 1.8 and earlier, the version from Compat has slightly lower floating-point accuracy than the one in Base (Julia 1.11 and later).
+
 * `Iterators.cycle(itr, n)` is the lazy version of `repeat(vector, n)`. ([#47354]) (since Compat 4.13.0)
 
 * `@compat public foo, bar` marks `foo` and `bar` as public in Julia 1.11+ and is a no-op in Julia 1.10 and earlier. ([#50105]) (since Compat 3.47.0, 4.10.0)
@@ -163,6 +165,7 @@ Note that you should specify the correct minimum version for `Compat` in the
 [#36229]: https://github.com/JuliaLang/julia/issues/36229
 [#37978]: https://github.com/JuliaLang/julia/issues/37978
 [#39037]: https://github.com/JuliaLang/julia/issues/39037
+[#39071]: https://github.com/JuliaLang/julia/pull/39071
 [#39245]: https://github.com/JuliaLang/julia/issues/39245
 [#39285]: https://github.com/JuliaLang/julia/issues/39285
 [#39794]: https://github.com/JuliaLang/julia/issues/39794

src/Compat.jl

@@ -250,7 +250,7 @@ end
 if VERSION < v"1.8.0-DEV.487"
     export eachsplit
 
-    """
+    @doc """
         eachsplit(str::AbstractString, dlm; limit::Integer=0)
         eachsplit(str::AbstractString; limit::Integer=0)
 
@@ -275,7 +275,8 @@ if VERSION < v"1.8.0-DEV.487"
     "Ma"
     "rch"
     ```
-    """
+    """ eachsplit
+
     function eachsplit end
 
     struct SplitIterator{S<:AbstractString,F}
@@ -421,7 +422,7 @@ end
         return nothing
     end
 
-    """
+    @doc """
         pkgversion(m::Module)
 
     Return the version of the package that imported module `m`,
@@ -433,7 +434,8 @@ end
 
     To get the version of the package that imported the current module
     the form `pkgversion(@__MODULE__)` can be used.
-    """
+    """ pkgversion
+
     function pkgversion(m::Module)
         path = pkgdir(m)
         path === nothing && return nothing
@@ -457,7 +459,7 @@ if VERSION < v"1.9.0-DEV.1163"
     import Base: IteratorSize, HasLength, HasShape, OneTo
     export stack
 
-    """
+    @doc """
         stack(iter; [dims])
 
     Combine a collection of arrays (or other iterable objects) of equal size
@@ -547,10 +549,11 @@ if VERSION < v"1.9.0-DEV.1163"
     julia> hvcat(5, M...) |> size  # hvcat puts matrices next to each other
     (14, 15)
     ```
-    """
+    """ stack
+
     stack(iter; dims=:) = _stack(dims, iter)
 
-    """
+    @doc """
         stack(f, args...; [dims])
 
     Apply a function to each element of a collection, and `stack` the result.
@@ -576,7 +579,8 @@ if VERSION < v"1.9.0-DEV.1163"
      1.0  2.0  3.0   10.0   20.0   30.0  0.1   0.2   0.3
      4.0  5.0  6.0  400.0  500.0  600.0  0.04  0.05  0.06
     ```
-    """
+    """ stack(f, iter)
+
     stack(f, iter; dims=:) = _stack(dims, f(x) for x in iter)
     stack(f, xs, yzs...; dims=:) = _stack(dims, f(xy...) for xy in zip(xs, yzs...))
 
@@ -803,6 +807,226 @@ if VERSION < v"1.11.0-DEV.1579"
     Iterators.cycle(xs, n::Integer) = Iterators.flatten(Iterators.repeated(xs, n))
 end
 
+# https://github.com/JuliaLang/julia/pull/39071
+if !isdefined(Base, :logrange)  # VERSION < v"1.12.0-DEV.2" or appropriate 1.11.x after backporting
+
+    export logrange
+
+    @doc """
+        logrange(start, stop, length)
+        logrange(start, stop; length)
+
+    Construct a specialized array whose elements are spaced logarithmically
+    between the given endpoints. That is, the ratio of successive elements is
+    a constant, calculated from the length.
+
+    This is similar to `geomspace` in Python. Unlike `PowerRange` in Mathematica,
+    you specify the number of elements not the ratio.
+    Unlike `logspace` in Python and Matlab, the `start` and `stop` arguments are
+    always the first and last elements of the result, not powers applied to some base.
+
+    # Examples
+    ```
+    julia> logrange(10, 4000, length=3)
+    3-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}:
+     10.0, 200.0, 4000.0
+
+    julia> ans[2] ≈ sqrt(10 * 4000)  # middle element is the geometric mean
+    true
+
+    julia> range(10, 40, length=3)[2] ≈ (10 + 40)/2  # arithmetic mean
+    true
+
+    julia> logrange(1f0, 32f0, 11)
+    11-element Base.LogRange{Float32, Float64}:
+     1.0, 1.41421, 2.0, 2.82843, 4.0, 5.65685, 8.0, 11.3137, 16.0, 22.6274, 32.0
+
+    julia> logrange(1, 1000, length=4) ≈ 10 .^ (0:3)
+    true
+    ```
+
+    See the [`Compat.LogRange`](@ref Compat.LogRange) type for further details.
+
+    !!! compat "Julia 1.9"
+        The version of this struct in Compat.jl does not use `Base.TwicePrecision{Float64}`
+        before Julia 1.9, so it sometimes has larger floating-point errors on intermediate points.
+
+    !!! compat "Julia 1.11"
+        The printing of Compat.jl's version of the struct is also different,
+        less like `LinRange` and more like `Vector`.
+    """ logrange
+
+    logrange(start::Real, stop::Real, length::Integer) = LogRange(start, stop, Int(length))
+    logrange(start::Real, stop::Real; length::Integer) = logrange(start, stop, length)
+
+    @doc """
+        LogRange{T}(start, stop, len) <: AbstractVector{T}
+
+    A range whose elements are spaced logarithmically between `start` and `stop`,
+    with spacing controlled by `len`. Returned by [`logrange`](@ref).
+
+    Like [`LinRange`](@ref), the first and last elements will be exactly those
+    provided, but intermediate values may have small floating-point errors.
+    These are calculated using the logs of the endpoints, which are
+    stored on construction, often in higher precision than `T`.
+
+    !!! compat "Julia 1.9"
+        The version of this struct in Compat.jl does not use `Base.TwicePrecision{Float64}`
+        before Julia 1.9. Therefore it has larger floating-point errors on intermediate
+        points than shown below.
+
+    !!! compat "Julia 1.11"
+        The printing of Compat.jl's version of the struct is also different,
+        less like `LinRange` and more like `Vector`.
+
+    # Examples
+    ```
+    julia> logrange(1, 4, length=5)
+    5-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}:
+     1.0, 1.41421, 2.0, 2.82843, 4.0
+
+    julia> Base.LogRange{Float16}(1, 4, 5)
+    5-element Base.LogRange{Float16, Float64}:
+     1.0, 1.414, 2.0, 2.828, 4.0
+
+    julia> logrange(1e-310, 1e-300, 11)[1:2:end]
+    6-element Vector{Float64}:
+     1.0e-310
+     9.999999999999974e-309
+     9.999999999999981e-307
+     9.999999999999988e-305
+     9.999999999999994e-303
+     1.0e-300
+
+    julia> prevfloat(1e-308, 5) == ans[2]
+    true
+    ```
+
+    Note that integer eltype `T` is not allowed.
+    Use for instance `round.(Int, xs)`, or explicit powers of some integer base:
+
+    ```
+    julia> xs = logrange(1, 512, 4)
+    4-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}:
+     1.0, 8.0, 64.0, 512.0
+
+    julia> 2 .^ (0:3:9) |> println
+    [1, 8, 64, 512]
+    ```
+    """ LogRange
+
+    struct LogRange{T<:Real,X} <: AbstractArray{T,1}
+        start::T
+        stop::T
+        len::Int
+        extra::Tuple{X,X}
+        function LogRange{T}(start::T, stop::T, len::Int) where {T<:Real}
+            if T <: Integer
+                # LogRange{Int}(1, 512, 4) produces InexactError: Int64(7.999999999999998)
+                throw(ArgumentError("LogRange{T} does not support integer types"))
+            end
+            if iszero(start) || iszero(stop)
+                throw(DomainError((start, stop),
+                    "LogRange cannot start or stop at zero"))
+            elseif start < 0 || stop < 0
+                # log would throw, but _log_twice64_unchecked does not
+                throw(DomainError((start, stop),
+                    "LogRange does not accept negative numbers"))
+            elseif !isfinite(start) || !isfinite(stop)
+                throw(DomainError((start, stop),
+                    "LogRange is only defined for finite start & stop"))
+            elseif len < 0
+                throw(ArgumentError(string( # LazyString(
+                    "LogRange(", start, ", ", stop, ", ", len, "): can't have negative length")))
+            elseif len == 1 && start != stop
+                throw(ArgumentError(string( # LazyString(
+                    "LogRange(", start, ", ", stop, ", ", len, "): endpoints differ, while length is 1")))
+            end
+            ex = _logrange_extra(start, stop, len)
+            new{T,typeof(ex[1])}(start, stop, len, ex)
+        end
+    end
+
+    function LogRange{T}(start::Real, stop::Real, len::Integer) where {T}
+        LogRange{T}(convert(T, start), convert(T, stop), convert(Int, len))
+    end
+    function LogRange(start::Real, stop::Real, len::Integer)
+        T = float(promote_type(typeof(start), typeof(stop)))
+        LogRange{T}(convert(T, start), convert(T, stop), convert(Int, len))
+    end
+
+    Base.size(r::LogRange) = (r.len,)
+    Base.length(r::LogRange) = r.len
+
+    Base.first(r::LogRange) = r.start
+    Base.last(r::LogRange) = r.stop
+
+    function _logrange_extra(a::Real, b::Real, len::Int)
+        loga = log(1.0 * a)  # widen to at least Float64
+        logb = log(1.0 * b)
+        (loga/(len-1), logb/(len-1))
+    end
+
+    function Base.getindex(r::LogRange{T}, i::Int) where {T}
+        @inline
+        @boundscheck checkbounds(r, i)
+        i == 1 && return r.start
+        i == r.len && return r.stop
+        # Main path uses Math.exp_impl for TwicePrecision, but is not perfectly
+        # accurate, hence the special cases for endpoints above.
+        logx = (r.len-i) * r.extra[1] + (i-1) * r.extra[2]
+        x = _exp_allowing_twice64(logx)
+        return copysign(T(x), r.start)
+    end
+
+    function Base.show(io::IO, r::LogRange{T}) where {T}
+        print(io, "LogRange{", T, "}(")
+        ioc = IOContext(io, :typeinfo => T)
+        show(ioc, first(r))
+        print(io, ", ")
+        show(ioc, last(r))
+        print(io, ", ")
+        show(io, length(r))
+        print(io, ')')
+    end
+
+    # Display LogRange like LinRange -- PR widened signature of print_range to allow this
+    # function Base.show(io::IO, ::MIME"text/plain", r::LogRange)
+    #     isempty(r) && return show(io, r)
+    #     summary(io, r)
+    #     println(io, ":")
+    #     print_range(io, r, " ", ", ", "", " \u2026 ")
+    # end
+
+    _exp_allowing_twice64(x::Number) = exp(x)
+
+    if VERSION >= v"1.9.0-DEV.318"  # Julia PR #44717 allows this high-precision path:
+        
+        _exp_allowing_twice64(x::Base.TwicePrecision{Float64}) = Base.Math.exp_impl(x.hi, x.lo, Val(:ℯ))
+        
+        function _log_twice64_unchecked(x::Float64)
+            xu = reinterpret(UInt64, x)
+            if xu < (UInt64(1)<<52) # x is subnormal
+                xu = reinterpret(UInt64, x * 0x1p52) # normalize x
+                xu &= ~Base.sign_mask(Float64)
+                xu -= UInt64(52) << 52 # mess with the exponent
+            end
+            Base.TwicePrecision(Base.Math._log_ext(xu)...)
+        end
+
+        function _logrange_extra(a::Float64, b::Float64, len::Int)
+            loga = _log_twice64_unchecked(a)
+            logb = _log_twice64_unchecked(b)
+            # The reason not to do linear interpolation on log(a)..log(b) in `getindex` is
+            # that division of TwicePrecision is quite slow, so do it once on construction:
+            (loga/(len-1), logb/(len-1))
+        end
+    end
+else
+    # Ensure that Compat.LogRange is always this struct, not exported from Base
+    using Base: LogRange
+end
+
 include("deprecated.jl")
 
 end # module Compat

test/runtests.jl

@@ -780,3 +780,68 @@ end
     Base.haslength(cycle(0:3, 2)) == false  # but not sure we should test these
     Base.IteratorSize(cycle(0:3, 2)) == Base.SizeUnknown()
 end
+
+# https://github.com/JuliaLang/julia/pull/39071
+@testset "logrange" begin
+    # basic idea
+    @test logrange(2, 16, 4) ≈ [2, 4, 8, 16]
+    @test logrange(1/8, 8.0, 7) ≈ [0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]
+    @test logrange(1000, 1, 4) ≈ [1000, 100, 10, 1]
+    @test logrange(1, 10^9, 19)[1:2:end] ≈ 10 .^ (0:9)
+
+    # endpoints
+    @test logrange(0.1f0, 100, 33)[1] === 0.1f0
+    @test logrange(0.789, 123_456, 135_790)[[begin, end]] == [0.789, 123_456]
+    @test logrange(nextfloat(0f0), floatmax(Float32), typemax(Int))[end] === floatmax(Float32)
+    @test logrange(nextfloat(Float16(0)), floatmax(Float16), 66_000)[end] === floatmax(Float16)
+    @test first(logrange(pi, 2pi, 3000)) === logrange(pi, 2pi, 3000)[1] === Float64(pi)
+    if Int == Int64
+        @test logrange(0.1, 1000, 2^54)[end] === 1000.0
+    end
+
+    # empty, only, constant
+    @test first(logrange(1, 2, 0)) === 1.0
+    @test last(logrange(1, 2, 0)) === 2.0
+    @test collect(logrange(1, 2, 0)) == Float64[]
+    @test only(logrange(2pi, 2pi, 1)) === logrange(2pi, 2pi, 1)[1] === 2pi
+    @test logrange(1, 1, 3) == fill(1.0, 3)
+
+    # subnormal Float64
+    x = logrange(1e-320, 1e-300, 21) .* 1e300
+    @test x ≈ logrange(1e-20, 1, 21) rtol=1e-6
+
+    # types
+    @test eltype(logrange(1, 10, 3)) == Float64
+    @test eltype(logrange(1, 10, Int32(3))) == Float64
+    @test eltype(logrange(1, 10f0, 3)) == Float32
+    @test eltype(logrange(1f0, 10, 3)) == Float32
+    @test eltype(logrange(1, big(10), 3)) == BigFloat
+    @test logrange(big"0.3", big(pi), 50)[1] == big"0.3"
+    @test logrange(big"0.3", big(pi), 50)[end] == big(pi)
+
+    # more constructors
+    @test logrange(1,2,length=3) === Compat.LogRange(1,2,3) == Compat.LogRange{Float64}(1,2,3)
+    @test logrange(1f0, 2f0, length=3) == Compat.LogRange{Float32}(1,2,3)
+
+    # errors
+    @test_throws UndefKeywordError logrange(1, 10)  # no default length
+    @test_throws ArgumentError logrange(1, 10, -1)  # negative length
+    @test_throws ArgumentError logrange(1, 10, 1) # endpoints must not differ
+    @test_throws DomainError logrange(1, -1, 3)   # needs complex numbers
+    @test_throws DomainError logrange(-1, -2, 3)  # not supported, for now
+    @test_throws MethodError logrange(1, 2+3im, length=4)  # not supported, for now
+    @test_throws ArgumentError logrange(1, 10, 2)[true]  # bad index
+    @test_throws BoundsError logrange(1, 10, 2)[3]
+    @test_throws ArgumentError Compat.LogRange{Int}(1,4,5)  # no integer ranges
+    @test_throws MethodError Compat.LogRange(1,4, length=5)  # type does not take keyword
+    # (not sure if these should ideally be DomainError or ArgumentError)
+    @test_throws DomainError logrange(1, Inf, 3)
+    @test_throws DomainError logrange(0, 2, 3)
+    @test_throws DomainError logrange(1, NaN, 3)
+    @test_throws DomainError logrange(NaN, 2, 3)
+
+    # printing
+    @test repr(Compat.LogRange(1,2,3)) == "LogRange{Float64}(1.0, 2.0, 3)"  # like 2-arg show
+    @test_skip repr("text/plain", Compat.LogRange(1,2,3)) == "3-element Compat.LogRange{Float64, Base.TwicePrecision{Float64}}:\n 1.0, 1.41421, 2.0"
+    @test_skip repr("text/plain", Compat.LogRange(1,2,0)) == "LogRange{Float64}(1.0, 2.0, 0)"  # empty case
+end