diff --git a/NEWS.md b/NEWS.md index fc3a2304ea72a..42d36b606e839 100644 --- a/NEWS.md +++ b/NEWS.md @@ -83,6 +83,7 @@ New library functions * `in!(x, s::AbstractSet)` will return whether `x` is in `s`, and insert `x` in `s` if not. * The new `Libc.mkfifo` function wraps the `mkfifo` C function on Unix platforms ([#34587]). +* `logrange(start, stop; length)` makes a range of constant ratio, instead of constant step ([#39071]) * `copyuntil(out, io, delim)` and `copyline(out, io)` copy data into an `out::IO` stream ([#48273]). * `eachrsplit(string, pattern)` iterates split substrings right to left. * `Sys.username()` can be used to return the current user's username ([#51897]). diff --git a/base/exports.jl b/base/exports.jl index fbb5b449719d6..a31ea1b0fa842 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -412,6 +412,7 @@ export isperm, issorted, last, + logrange, mapslices, max, maximum!, @@ -1095,6 +1096,7 @@ public Generator, ImmutableDict, OneTo, + LogRange, AnnotatedString, AnnotatedChar, UUID, diff --git a/base/range.jl b/base/range.jl index 9480b46d8fcf1..4f92b305564cd 100644 --- a/base/range.jl +++ b/base/range.jl @@ -70,6 +70,7 @@ Valid invocations of range are: * Call `range` with one of `stop` or `length`. `start` and `step` will be assumed to be one. See Extended Help for additional details on the returned type. +See also [`logrange`](@ref) for logarithmically spaced points. # Examples ```jldoctest @@ -252,10 +253,13 @@ end ## 1-dimensional ranges ## """ - AbstractRange{T} + AbstractRange{T} <: AbstractVector{T} -Supertype for ranges with elements of type `T`. -[`UnitRange`](@ref) and other types are subtypes of this. +Supertype for linear ranges with elements of type `T`. +[`UnitRange`](@ref), [`LinRange`](@ref) and other types are subtypes of this. + +All subtypes must define [`step`](@ref). +Thus [`LogRange`](@ref Base.LogRange) is not a subtype of `AbstractRange`. """ abstract type AbstractRange{T} <: AbstractArray{T,1} end @@ -550,6 +554,8 @@ julia> collect(LinRange(-0.1, 0.3, 5)) 0.19999999999999998 0.3 ``` + +See also [`Logrange`](@ref Base.LogRange) for logarithmically spaced points. """ struct LinRange{T,L<:Integer} <: AbstractRange{T} start::T @@ -620,7 +626,7 @@ parameters `pre` and `post` characters for each printed row, `sep` separator string between printed elements, `hdots` string for the horizontal ellipsis. """ -function print_range(io::IO, r::AbstractRange, +function print_range(io::IO, r::AbstractArray, pre::AbstractString = " ", sep::AbstractString = ", ", post::AbstractString = "", @@ -1488,3 +1494,179 @@ julia> mod(3, 0:2) # mod(3, 3) """ mod(i::Integer, r::OneTo) = mod1(i, last(r)) mod(i::Integer, r::AbstractUnitRange{<:Integer}) = mod(i-first(r), length(r)) + first(r) + + +""" + 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 +```jldoctest +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 [`LogRange`](@ref Base.LogRange) type for further details. + +See also [`range`](@ref) for linearly spaced points. + +!!! compat "Julia 1.11" + This function requires at least Julia 1.11. +""" +logrange(start::Real, stop::Real, length::Integer) = LogRange(start, stop, Int(length)) +logrange(start::Real, stop::Real; length::Integer) = logrange(start, stop, length) + + +""" + 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`. + +# Examples +```jldoctest +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: + +```jldoctest +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] +``` + +!!! compat "Julia 1.11" + This type requires at least Julia 1.11. +""" +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(LazyString( + "LogRange(", start, ", ", stop, ", ", len, "): can't have negative length"))) + elseif len == 1 && start != stop + throw(ArgumentError(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 + +size(r::LogRange) = (r.len,) +length(r::LogRange) = r.len + +first(r::LogRange) = r.start +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 _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 + +function 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 T(x) +end + +function 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 diff --git a/base/show.jl b/base/show.jl index a5bf5f73ba0ed..217b22ab39ef8 100644 --- a/base/show.jl +++ b/base/show.jl @@ -23,6 +23,13 @@ function show(io::IO, ::MIME"text/plain", r::LinRange) print_range(io, r) end +function show(io::IO, ::MIME"text/plain", r::LogRange) # display LogRange like LinRange + isempty(r) && return show(io, r) + summary(io, r) + println(io, ":") + print_range(io, r, " ", ", ", "", " \u2026 ") +end + function _isself(ft::DataType) ftname = ft.name isdefined(ftname, :mt) || return false diff --git a/base/twiceprecision.jl b/base/twiceprecision.jl index 955bfc97b16ff..0de6270cafb2d 100644 --- a/base/twiceprecision.jl +++ b/base/twiceprecision.jl @@ -785,3 +785,19 @@ _tp_prod(t::TwicePrecision) = t x.hi < y.hi || ((x.hi == y.hi) & (x.lo < y.lo)) isbetween(a, x, b) = a <= x <= b || b <= x <= a + +# These functions exist for use in LogRange: + +_exp_allowing_twice64(x::Number) = exp(x) +_exp_allowing_twice64(x::TwicePrecision{Float64}) = Math.exp_impl(x.hi, x.lo, Val(:ℯ)) + +# No error on negative x, and for NaN/Inf this returns junk: +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 &= ~sign_mask(Float64) + xu -= UInt64(52) << 52 # mess with the exponent + end + TwicePrecision(Math._log_ext(xu)...) +end diff --git a/doc/src/base/math.md b/doc/src/base/math.md index 86a7e4d8d3173..c1b18c0b66d86 100644 --- a/doc/src/base/math.md +++ b/doc/src/base/math.md @@ -39,6 +39,8 @@ Base.:(:) Base.range Base.OneTo Base.StepRangeLen +Base.logrange +Base.LogRange Base.:(==) Base.:(!=) Base.:(!==) diff --git a/test/ranges.jl b/test/ranges.jl index 549078bab45fe..4660a96dfc16a 100644 --- a/test/ranges.jl +++ b/test/ranges.jl @@ -2603,3 +2603,86 @@ end errmsg = ("deliberately unsupported for CartesianIndex", "StepRangeLen") @test_throws errmsg range(CartesianIndex(1), step=CartesianIndex(1), length=3) end + +@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) === Base.LogRange(1,2,3) == Base.LogRange{Float64}(1,2,3) + @test logrange(1f0, 2f0, length=3) == Base.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 Base.LogRange{Int}(1,4,5) # no integer ranges + @test_throws MethodError Base.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(Base.LogRange(1,2,3)) == "LogRange{Float64}(1.0, 2.0, 3)" # like 2-arg show + @test repr("text/plain", Base.LogRange(1,2,3)) == "3-element Base.LogRange{Float64, Base.TwicePrecision{Float64}}:\n 1.0, 1.41421, 2.0" + @test repr("text/plain", Base.LogRange(1,2,0)) == "LogRange{Float64}(1.0, 2.0, 0)" # empty case +end + +@testset "_log_twice64_unchecked" begin + # it roughly works + @test big(Base._log_twice64_unchecked(exp(1))) ≈ 1.0 + @test big(Base._log_twice64_unchecked(exp(123))) ≈ 123.0 + + # it gets high accuracy + @test abs(big(log(4.0)) - log(big(4.0))) < 1e-16 + @test abs(big(Base._log_twice64_unchecked(4.0)) - log(big(4.0))) < 1e-30 + + # it handles subnormals + @test abs(big(Base._log_twice64_unchecked(1e-310)) - log(big(1e-310))) < 1e-20 + + # it accepts negative, NaN, etc without complaint: + @test Base._log_twice64_unchecked(-0.0).lo isa Float64 + @test Base._log_twice64_unchecked(-1.23).lo isa Float64 + @test Base._log_twice64_unchecked(NaN).lo isa Float64 + @test Base._log_twice64_unchecked(Inf).lo isa Float64 +end