From 7382e9190d9dcd673d04d2e505d4a39c4c04e14f Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Mon, 2 Jan 2017 14:22:48 -0600 Subject: [PATCH 1/2] Make LinSpace generic and endpoint-preserving. FloatRange->StepRangeLen. Fixes #14420. --- NEWS.md | 28 +++ base/abstractarray.jl | 6 +- base/coreimg.jl | 1 + base/dates/types.jl | 3 + base/deprecated.jl | 6 +- base/exports.jl | 2 +- base/float.jl | 44 ++-- base/mpfr.jl | 5 + base/operators.jl | 44 ++-- base/range.jl | 497 +++++++++++++++++----------------------- base/rational.jl | 4 + base/sysimg.jl | 2 + base/traits.jl | 22 ++ base/twiceprecision.jl | 508 +++++++++++++++++++++++++++++++++++++++++ doc/src/stdlib/math.md | 1 + test/core.jl | 2 +- test/ranges.jl | 146 +++++++----- 17 files changed, 928 insertions(+), 393 deletions(-) create mode 100644 base/traits.jl create mode 100644 base/twiceprecision.jl diff --git a/NEWS.md b/NEWS.md index a899202425e93..f6dd9408fe99b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -130,6 +130,33 @@ This section lists changes that do not have deprecation warnings. (since it is shorthand for `NTuple{N,T} where T`). To get the old behavior of matching any tuple, use `NTuple{N,Any}` ([#18457]). + * `FloatRange` has been replaced by `StepRangeLen`, and the internal + representation of `LinSpace` has changed. Aside from changes in + the internal field names, this leads to several differences in + behavior ([#18777]): + + + Both `StepRangeLen` and `LinSpace` can represent ranges of + arbitrary object types---they are no longer limited to + floating-point numbers. + + + For ranges that produce `Float64`, `Float32`, or `Float16` + numbers, `StepRangeLen` can be used to produce values with + little or no roundoff error due to internal arithmetic that is + typically twice the precision of the output result. + + + To take advantage of this precision, `linspace(start, stop, + len)` now returns a range of type `StepRangeLen` rather than + `LinSpace` when `start` and `stop` are + `FloatNN`. `LinSpace(start, stop, len)` always returns a + `LinSpace`. + + + `StepRangeLen(a, step, len)` constructs a single-precision range + using the values and types of `a` and `step` as given, whereas + `range(a, step, len)` will attempt to match inputs `a::FloatNN` + and `step::FloatNN` to rationals and construct a `StepRangeLen` + that internally uses twice-precision arithmetic. These two + outcomes exhibit differences in both precision and speed. + Library improvements -------------------- @@ -891,6 +918,7 @@ Language tooling improvements [#18628]: https://github.com/JuliaLang/julia/issues/18628 [#18644]: https://github.com/JuliaLang/julia/issues/18644 [#18690]: https://github.com/JuliaLang/julia/issues/18690 +[#18777]: https://github.com/JuliaLang/julia/issues/18777 [#18839]: https://github.com/JuliaLang/julia/issues/18839 [#18931]: https://github.com/JuliaLang/julia/issues/18931 [#18965]: https://github.com/JuliaLang/julia/issues/18965 diff --git a/base/abstractarray.jl b/base/abstractarray.jl index fbdf2f51133a0..1f408445398b3 100644 --- a/base/abstractarray.jl +++ b/base/abstractarray.jl @@ -781,11 +781,9 @@ full(x::AbstractArray) = x map{T<:Real}(::Type{T}, r::StepRange) = T(r.start):T(r.step):T(last(r)) map{T<:Real}(::Type{T}, r::UnitRange) = T(r.start):T(last(r)) -map{T<:AbstractFloat}(::Type{T}, r::FloatRange) = FloatRange(T(r.start), T(r.step), r.len, T(r.divisor)) +map{T<:AbstractFloat}(::Type{T}, r::StepRangeLen) = convert(StepRangeLen{T}, r) function map{T<:AbstractFloat}(::Type{T}, r::LinSpace) - new_len = T(r.len) - new_len == r.len || error("$r: too long for $T") - LinSpace(T(r.start), T(r.stop), new_len, T(r.divisor)) + LinSpace(T(r.start), T(r.stop), length(r)) end ## unsafe/pointer conversions ## diff --git a/base/coreimg.jl b/base/coreimg.jl index 2e2aa989048a6..3f1ea2b1ead62 100644 --- a/base/coreimg.jl +++ b/base/coreimg.jl @@ -23,6 +23,7 @@ include("options.jl") # core operations & types include("promotion.jl") include("tuple.jl") +include("traits.jl") include("range.jl") include("expr.jl") include("error.jl") diff --git a/base/dates/types.jl b/base/dates/types.jl index db85b10cdd7ad..347d9059388ec 100644 --- a/base/dates/types.jl +++ b/base/dates/types.jl @@ -310,3 +310,6 @@ import Base: sleep, Timer, timedwait sleep(time::Period) = sleep(toms(time) / 1000) Timer(time::Period, repeat::Period=Second(0)) = Timer(toms(time) / 1000, toms(repeat) / 1000) timedwait(testcb::Function, time::Period) = timedwait(testcb, toms(time) / 1000) + +(::Type{Base.TypeOrder}){T<:AbstractTime}(::Type{T}) = Base.HasOrder() +(::Type{Base.TypeArithmetic}){T<:AbstractTime}(::Type{T}) = Base.ArithmeticOverflows() diff --git a/base/deprecated.jl b/base/deprecated.jl index 78eed7327b74c..84ff810f81483 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -1471,7 +1471,7 @@ end # Deprecate manually vectorized `big` methods in favor of compact broadcast syntax @deprecate big(r::UnitRange) big.(r) @deprecate big(r::StepRange) big.(r) -@deprecate big(r::FloatRange) big.(r) +@deprecate big(r::StepRangeLen) big.(r) @deprecate big(r::LinSpace) big.(r) @deprecate big{T<:Integer,N}(x::AbstractArray{T,N}) big.(x) @deprecate big{T<:AbstractFloat,N}(x::AbstractArray{T,N}) big.(x) @@ -1842,4 +1842,8 @@ eval(Dates, quote end end) +# FloatRange replaced by StepRangeLen + +@deprecate FloatRange{T}(start::T, step, len, den) Base.floatrange(T, start, step, len, den) + # End deprecations scheduled for 0.6 diff --git a/base/exports.jl b/base/exports.jl index 52e5786a37aef..736038da9e5cc 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -62,7 +62,7 @@ export ExponentialBackOff, Factorization, FileMonitor, - FloatRange, + StepRangeLen, Future, Hermitian, UniformScaling, diff --git a/base/float.jl b/base/float.jl index bbfe7ee333807..5b80d1033a3c0 100644 --- a/base/float.jl +++ b/base/float.jl @@ -752,6 +752,29 @@ fpinttype(::Type{Float16}) = UInt16 # maximum float exponent without bias @pure exponent_raw_max{T<:AbstractFloat}(::Type{T}) = Int(exponent_mask(T) >> significand_bits(T)) +## TwicePrecision utilities +# The numeric constants are half the number of bits in the mantissa +for (F, T, n) in ((Float16, UInt16, 5), (Float32, UInt32, 12), (Float64, UInt64, 26)) + @eval begin + function truncbits(x::$F, nb) + @_inline_meta + truncmask(x, typemax($T) << nb) + end + function truncmask(x::$F, mask) + @_inline_meta + reinterpret($F, mask & reinterpret($T, x)) + end + function splitprec(x::$F) + @_inline_meta + hi = truncmask(x, typemax($T) << $n) + hi, x-hi + end + end +end + +truncbits(x, nb) = x +truncmask(x, mask) = x + ## Array operations on floating point numbers ## float{T<:AbstractFloat}(A::AbstractArray{T}) = A @@ -763,17 +786,11 @@ function float{T}(A::AbstractArray{T}) convert(AbstractArray{typeof(float(zero(T)))}, A) end -for fn in (:float,) - @eval begin - $fn(r::StepRange) = $fn(r.start):$fn(r.step):$fn(last(r)) - $fn(r::UnitRange) = $fn(r.start):$fn(last(r)) - $fn(r::FloatRange) = FloatRange($fn(r.start), $fn(r.step), r.len, $fn(r.divisor)) - function $fn(r::LinSpace) - new_len = $fn(r.len) - new_len == r.len || error(string(r, ": too long for ", $fn)) - LinSpace($fn(r.start), $fn(r.stop), new_len, $fn(r.divisor)) - end - end +float(r::StepRange) = float(r.start):float(r.step):float(last(r)) +float(r::UnitRange) = float(r.start):float(last(r)) +float(r::StepRangeLen) = StepRangeLen(float(r.ref), float(r.step), length(r), r.offset) +function float(r::LinSpace) + LinSpace(float(r.start), float(r.stop), length(r)) end # big, broadcast over arrays @@ -781,8 +798,7 @@ end function big end # no prior definitions of big in sysimg.jl, necessitating this broadcast(::typeof(big), r::UnitRange) = big(r.start):big(last(r)) broadcast(::typeof(big), r::StepRange) = big(r.start):big(r.step):big(last(r)) -broadcast(::typeof(big), r::FloatRange) = FloatRange(big(r.start), big(r.step), r.len, big(r.divisor)) +broadcast(::typeof(big), r::StepRangeLen) = StepRangeLen(big(r.ref), big(r.step), length(r), r.offset) function broadcast(::typeof(big), r::LinSpace) - big(r.len) == r.len || throw(ArgumentError(string(r, ": too long for ", big))) - LinSpace(big(r.start), big(r.stop), big(r.len), big(r.divisor)) + LinSpace(big(r.start), big(r.stop), length(r)) end diff --git a/base/mpfr.jl b/base/mpfr.jl index 76ba6fcee3d41..807fd36966b88 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -978,4 +978,9 @@ function Base.deepcopy_internal(x::BigFloat, stackdict::ObjectIdDict) return y end +function lerpi(j::Integer, d::Integer, a::BigFloat, b::BigFloat) + t = BigFloat(j)/d + fma(t, b, fma(-t, a, a)) +end + end #module diff --git a/base/operators.jl b/base/operators.jl index acc39a95bb508..07f829d34a09f 100644 --- a/base/operators.jl +++ b/base/operators.jl @@ -976,45 +976,29 @@ for f in (:+, :-) range($f(first(r1),first(r2)), $f(step(r1),step(r2)), r1l) end - function $f{T<:AbstractFloat}(r1::FloatRange{T}, r2::FloatRange{T}) + function $f{T}(r1::LinSpace{T}, r2::LinSpace{T}) len = r1.len (len == r2.len || throw(DimensionMismatch("argument dimensions must match"))) - divisor1, divisor2 = r1.divisor, r2.divisor - if divisor1 == divisor2 - FloatRange{T}($f(r1.start,r2.start), $f(r1.step,r2.step), - len, divisor1) - else - d1 = Int(divisor1) - d2 = Int(divisor2) - d = lcm(d1,d2) - s1 = div(d,d1) - s2 = div(d,d2) - FloatRange{T}($f(r1.start*s1, r2.start*s2), - $f(r1.step*s1, r2.step*s2), len, d) - end - end - - function $f{T<:AbstractFloat}(r1::LinSpace{T}, r2::LinSpace{T}) - len = r1.len - (len == r2.len || - throw(DimensionMismatch("argument dimensions must match"))) - divisor1, divisor2 = r1.divisor, r2.divisor - if divisor1 == divisor2 - LinSpace{T}($f(r1.start, r2.start), $f(r1.stop, r2.stop), - len, divisor1) - else - linspace(convert(T, $f(first(r1), first(r2))), - convert(T, $f(last(r1), last(r2))), len) - end + linspace(convert(T, $f(first(r1), first(r2))), + convert(T, $f(last(r1), last(r2))), len) end - $f(r1::Union{FloatRange, OrdinalRange, LinSpace}, - r2::Union{FloatRange, OrdinalRange, LinSpace}) = + $f(r1::Union{StepRangeLen, OrdinalRange, LinSpace}, + r2::Union{StepRangeLen, OrdinalRange, LinSpace}) = $f(promote_noncircular(r1, r2)...) end end +function +{T,S}(r1::StepRangeLen{T,S}, r2::StepRangeLen{T,S}) + len = length(r1) + (len == length(r2) || + throw(DimensionMismatch("argument dimensions must match"))) + StepRangeLen(first(r1)+first(r2), step(r1)+step(r2), len) +end + +-(r1::StepRangeLen, r2::StepRangeLen) = +(r1, -r2) + # Pair immutable Pair{A,B} diff --git a/base/range.jl b/base/range.jl index d58238d93ee0f..ee89fa65ee208 100644 --- a/base/range.jl +++ b/base/range.jl @@ -1,5 +1,60 @@ # This file is a part of Julia. License is MIT: http://julialang.org/license +colon(a::Real, b::Real) = colon(promote(a,b)...) + +colon{T<:Real}(start::T, stop::T) = UnitRange{T}(start, stop) + +range(a::Real, len::Integer) = UnitRange{typeof(a)}(a, oftype(a, a+len-1)) + +colon{T}(start::T, stop::T) = colon(start, oftype(stop-start, 1), stop) + +range(a, len::Integer) = range(a, oftype(a-a, 1), len) + +# first promote start and stop, leaving step alone +colon{A<:Real,C<:Real}(start::A, step, stop::C) = colon(convert(promote_type(A,C),start), step, convert(promote_type(A,C),stop)) +colon{T<:Real}(start::T, step::Real, stop::T) = colon(promote(start, step, stop)...) + +""" + colon(start, [step], stop) + +Called by `:` syntax for constructing ranges. +""" +colon{T<:Real}(start::T, step::T, stop::T) = _colon(TypeOrder(T), TypeArithmetic(T), start, step, stop) +_colon{T}(::HasOrder, ::Any, start::T, step, stop::T) = StepRange(start, step, stop) +# for T<:Union{Float16,Float32,Float64} see twiceprecision.jl +_colon{T}(::HasOrder, ::ArithmeticRounds, start::T, step, stop::T) = StepRangeLen(start, step, floor(Int, (stop-start)/step)+1) +_colon{T}(::Any, ::Any, start::T, step, stop::T) = StepRangeLen(start, step, floor(Int, (stop-start)/step)+1) + +""" + :(start, [step], stop) + +Range operator. `a:b` constructs a range from `a` to `b` with a step size of 1, and `a:s:b` +is similar but uses a step size of `s`. These syntaxes call the function `colon`. The colon +is also used in indexing to select whole dimensions. +""" +colon{T}(start::T, step, stop::T) = StepRange(start, step, stop) + +""" + range(start, [step], length) + +Construct a range by length, given a starting value and optional step (defaults to 1). +""" +range{T}(a::T, step, len::Integer) = _range(TypeOrder(T), TypeArithmetic(T), a, step, len) +_range{T,S}(::HasOrder, ::ArithmeticOverflows, a::T, step::S, len::Integer) = StepRange{T,S}(a, step, convert(T, a+step*(len-1))) +_range{T,S}(::Any, ::Any, a::T, step::S, len::Integer) = StepRangeLen{typeof(a+0*step),T,S}(a, step, len) + +# AbstractFloat specializations +colon{T<:AbstractFloat}(a::T, b::T) = colon(a, T(1), b) +range(a::AbstractFloat, len::Integer) = range(a, oftype(a, 1), len) + +colon{T<:Real}(a::T, b::AbstractFloat, c::T) = colon(promote(a,b,c)...) +colon{T<:AbstractFloat}(a::T, b::AbstractFloat, c::T) = colon(promote(a,b,c)...) +colon{T<:AbstractFloat}(a::T, b::Real, c::T) = colon(promote(a,b,c)...) + +range(a::AbstractFloat, st::AbstractFloat, len::Integer) = range(promote(a, st)..., len) +range(a::Real, st::AbstractFloat, len::Integer) = range(float(a), st, len) +range(a::AbstractFloat, st::Real, len::Integer) = range(a, float(st), len) + ## 1-dimensional ranges ## abstract Range{T} <: AbstractArray{T,1} @@ -96,192 +151,74 @@ immutable OneTo{T<:Integer} <: AbstractUnitRange{T} end OneTo{T<:Integer}(stop::T) = OneTo{T}(stop) -colon(a::Real, b::Real) = colon(promote(a,b)...) - -colon{T<:Real}(start::T, stop::T) = UnitRange{T}(start, stop) - -range(a::Real, len::Integer) = UnitRange{typeof(a)}(a, oftype(a, a+len-1)) - -colon{T}(start::T, stop::T) = StepRange(start, one(stop-start), stop) - -range{T}(a::T, len::Integer) = - StepRange{T, typeof(a-a)}(a, one(a-a), a+oftype(a-a,(len-1))) - -# first promote start and stop, leaving step alone -# this is for non-numeric ranges where step can be quite different -colon{A<:Real,C<:Real}(a::A, b, c::C) = colon(convert(promote_type(A,C),a), b, convert(promote_type(A,C),c)) - -""" - colon(start, [step], stop) - -Called by `:` syntax for constructing ranges. -""" -colon{T<:Real}(start::T, step, stop::T) = StepRange(start, step, stop) - -""" - :(start, [step], stop) - -Range operator. `a:b` constructs a range from `a` to `b` with a step size of 1, and `a:s:b` -is similar but uses a step size of `s`. These syntaxes call the function `colon`. The colon -is also used in indexing to select whole dimensions. -""" -colon{T<:Real}(start::T, step::T, stop::T) = StepRange(start, step, stop) -colon{T<:Real}(start::T, step::Real, stop::T) = StepRange(promote(start, step, stop)...) - -colon{T}(start::T, step, stop::T) = StepRange(start, step, stop) +## Step ranges parametrized by length """ - range(start, [step], length) - -Construct a range by length, given a starting value and optional step (defaults to 1). + StepRangeLen{T,R,S}(ref::R, step::S, len, [offset=1]) + +A range `r` where `r[i]` produces values of type `T`, parametrized by +a `ref`erence value, a `step`, and the `len`gth. By default `ref` is +the starting value `r[1]`, but alternatively you can supply it as the +value of `r[offset]` for some other index `1 <= offset <= len`. In +conjunction with `TwicePrecision` this can be used to implement ranges +that are free of roundoff error. """ -range{T,S}(a::T, step::S, len::Integer) = StepRange{T,S}(a, step, convert(T, a+step*(len-1))) - -## floating point ranges - -immutable FloatRange{T<:AbstractFloat} <: Range{T} - start::T - step::T - len::T - divisor::T -end -FloatRange(a::AbstractFloat, s::AbstractFloat, l::Real, d::AbstractFloat) = - FloatRange{promote_type(typeof(a),typeof(s),typeof(d))}(a,s,l,d) - -# float rationalization helper -function rat(x) - y = x - a = d = 1 - b = c = 0 - m = maxintfloat(Float32) - while abs(y) <= m - f = trunc(Int,y) - y -= f - a, c = f*a + c, a - b, d = f*b + d, b - max(abs(a),abs(b)) <= convert(Int,m) || return c, d - oftype(x,a)/oftype(x,b) == x && break - y = inv(y) +immutable StepRangeLen{T,R,S} <: Range{T} + ref::R # reference value (might be smallest-magnitude value in the range) + step::S # step value + len::Int # length of the range + offset::Int # the index of ref + + function StepRangeLen(ref::R, step::S, len::Integer, offset::Integer = 1) + len >= 0 || throw(ArgumentError("length cannot be negative, got $len")) + 1 <= offset <= max(1,len) || throw(ArgumentError("StepRangeLen: offset must be in [1,$len], got $offset")) + new(ref, step, len, offset) end - return a, b -end - -function colon{T<:AbstractFloat}(start::T, step::T, stop::T) - step == 0 && throw(ArgumentError("range step cannot be zero")) - start == stop && return FloatRange{T}(start,step,1,1) - (0 < step) != (start < stop) && return FloatRange{T}(start,step,0,1) - - # float range "lifting" - r = (stop-start)/step - n = round(r) - lo = prevfloat((prevfloat(stop)-nextfloat(start))/n) - hi = nextfloat((nextfloat(stop)-prevfloat(start))/n) - if lo <= step <= hi - a0, b = rat(start) - a = convert(T,a0) - if a/convert(T,b) == start - c0, d = rat(step) - c = convert(T,c0) - if c/convert(T,d) == step - e = lcm(b,d) - a *= div(e,b) - c *= div(e,d) - eT = convert(T,e) - if (a+n*c)/eT == stop - return FloatRange{T}(a, c, n+1, eT) - end - end - end - end - FloatRange{T}(start, step, floor(r)+1, one(step)) end -colon{T<:AbstractFloat}(a::T, b::T) = colon(a, one(a), b) - -colon{T<:Real}(a::T, b::AbstractFloat, c::T) = colon(promote(a,b,c)...) -colon{T<:AbstractFloat}(a::T, b::AbstractFloat, c::T) = colon(promote(a,b,c)...) -colon{T<:AbstractFloat}(a::T, b::Real, c::T) = colon(promote(a,b,c)...) - -range(a::AbstractFloat, len::Integer) = FloatRange(a,one(a),len,one(a)) -range(a::AbstractFloat, st::AbstractFloat, len::Integer) = FloatRange(a,st,len,one(a)) -range(a::Real, st::AbstractFloat, len::Integer) = FloatRange(float(a), st, len, one(st)) -range(a::AbstractFloat, st::Real, len::Integer) = FloatRange(a, float(st), len, one(a)) +StepRangeLen{R,S}(ref::R, step::S, len::Integer, offset::Integer = 1) = + StepRangeLen{typeof(ref+0*step),R,S}(ref, step, len, offset) ## linspace and logspace -immutable LinSpace{T<:AbstractFloat} <: Range{T} +immutable LinSpace{T} <: Range{T} start::T stop::T - len::T - divisor::T -end - -function linspace{T<:AbstractFloat}(start::T, stop::T, len::T) - len == round(len) || throw(InexactError()) - 0 <= len || error("linspace($start, $stop, $len): negative length") - if len == 0 - n = convert(T, 2) - if isinf(n*start) || isinf(n*stop) - start /= n; stop /= n; n = one(T) + len::Int + lendiv::Int + + function LinSpace(start,stop,len) + len >= 0 || throw(ArgumentError("linspace($start, $stop, $len): negative length")) + if len == 1 + start == stop || throw(ArgumentError("linspace($start, $stop, $len): endpoints differ")) + return new(start, stop, 1, 1) end - return LinSpace(-start, -stop, -one(T), n) + new(start,stop,len,max(len-1,1)) end - if len == 1 - start == stop || error("linspace($start, $stop, $len): endpoints differ") - return LinSpace(-start, -start, zero(T), one(T)) - end - n = convert(T, len - 1) - len - n == 1 || error("linspace($start, $stop, $len): too long for $T") - a0, b = rat(start) - a = convert(T,a0) - if a/convert(T,b) == start - c0, d = rat(stop) - c = convert(T,c0) - if c/convert(T,d) == stop - e = lcm(b,d) - a *= div(e,b) - c *= div(e,d) - s = convert(T,n*e) - if isinf(a*n) || isinf(c*n) - s, p = frexp(s) - p2 = oftype(s,2)^p - a /= p2; c /= p2 - end - if a*n/s == start && c*n/s == stop - return LinSpace(a, c, len, s) - end - end - end - a, c, s = start, stop, n - if isinf(a*n) || isinf(c*n) - s, p = frexp(s) - p2 = oftype(s,2)^p - a /= p2; c /= p2 - end - if a*n/s == start && c*n/s == stop - return LinSpace(a, c, len, s) - end - return LinSpace(start, stop, len, n) end -function linspace{T<:AbstractFloat}(start::T, stop::T, len::Real) - T_len = convert(T, len) - T_len == len || throw(InexactError()) - linspace(start, stop, T_len) + +function LinSpace(start, stop, len::Integer) + T = typeof((stop-start)/len) + LinSpace{T}(start, stop, len) end """ - linspace(start::Real, stop::Real, n::Real=50) + linspace(start, stop, n=50) Construct a range of `n` linearly spaced elements from `start` to `stop`. ```jldoctest julia> linspace(1.3,2.9,9) -9-element LinSpace{Float64}: - 1.3,1.5,1.7,1.9,2.1,2.3,2.5,2.7,2.9 +1.3:0.2:2.9 ``` """ -linspace(start::Real, stop::Real, len::Real=50) = - linspace(promote(AbstractFloat(start), AbstractFloat(stop))..., len) +linspace(start, stop, len::Real=50) = linspace(start, stop, Int(len)) + +linspace(start::Real, stop::Real, len::Integer) = linspace(promote(start, stop)..., len) +linspace{T<:Integer}(start::T, stop::T, len::Integer) = linspace(Float64, start, stop, len, 1) +# for Float16, Float32, and Float64 see twiceprecision.jl +linspace{T<:Real}(start::T, stop::T, len::Integer) = LinSpace{T}(start, stop, len) +linspace{T}(start::T, stop::T, len::Integer) = LinSpace{T}(start, stop, len) function show(io::IO, r::LinSpace) print(io, "linspace(") @@ -329,7 +266,7 @@ function print_range(io::IO, r::Range, # screen, with the middle columns summarized by horz, vert, or diag ellipsis maxpossiblecols = div(screenwidth, 1+sepsize) # assume each element is at least 1 char + 1 separator colsr = n <= maxpossiblecols ? (1:n) : [1:div(maxpossiblecols,2)+1; (n-div(maxpossiblecols,2)):n] - rowmatrix = r[colsr]' # treat the range as a one-row matrix for print_matrix_row + rowmatrix = reshape(r[colsr], 1, length(colsr)) # treat the range as a one-row matrix for print_matrix_row A = alignment(io, rowmatrix, 1:m, 1:length(rowmatrix), screenwidth, screenwidth, sepsize) # how much space range takes if n <= length(A) # cols fit screen, so print out all elements print(io, pre) # put in pre chars @@ -374,7 +311,7 @@ size(r::Range) = (length(r),) isempty(r::StepRange) = (r.start != r.stop) & ((r.step > zero(r.step)) != (r.stop > r.start)) isempty(r::AbstractUnitRange) = first(r) > last(r) -isempty(r::FloatRange) = length(r) == 0 +isempty(r::StepRangeLen) = length(r) == 0 isempty(r::LinSpace) = length(r) == 0 """ @@ -397,8 +334,8 @@ julia> step(linspace(2.5,10.9,85)) """ step(r::StepRange) = r.step step(r::AbstractUnitRange) = 1 -step(r::FloatRange) = r.step/r.divisor -step{T}(r::LinSpace{T}) = ifelse(r.len <= 0, convert(T,NaN), (r.stop-r.start)/r.divisor) +step(r::StepRangeLen) = r.step +step(r::LinSpace) = (last(r)-first(r))/r.lendiv unsafe_length(r::Range) = length(r) # generic fallback @@ -411,8 +348,8 @@ unsafe_length(r::AbstractUnitRange) = Integer(last(r) - first(r) + 1) unsafe_length(r::OneTo) = r.stop length(r::AbstractUnitRange) = unsafe_length(r) length(r::OneTo) = unsafe_length(r) -length(r::FloatRange) = Integer(r.len) -length(r::LinSpace) = Integer(r.len + signbit(r.len - 1)) +length(r::StepRangeLen) = r.len +length(r::LinSpace) = r.len function length{T<:Union{Int,UInt,Int64,UInt64}}(r::StepRange{T}) isempty(r) && return zero(T) @@ -451,12 +388,12 @@ end first{T}(r::OrdinalRange{T}) = convert(T, r.start) first{T}(r::OneTo{T}) = one(T) -first{T}(r::FloatRange{T}) = convert(T, r.start/r.divisor) -first{T}(r::LinSpace{T}) = convert(T, (r.len-1)*r.start/r.divisor) +first(r::StepRangeLen) = unsafe_getindex(r, 1) +first(r::LinSpace) = r.start last{T}(r::OrdinalRange{T}) = convert(T, r.stop) -last{T}(r::FloatRange{T}) = convert(T, (r.start + (r.len-1)*r.step)/r.divisor) -last{T}(r::LinSpace{T}) = convert(T, (r.len-1)*r.stop/r.divisor) +last(r::StepRangeLen) = unsafe_getindex(r, length(r)) +last(r::LinSpace) = r.stop minimum(r::AbstractUnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : first(r) maximum(r::AbstractUnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : last(r) @@ -469,15 +406,12 @@ copy(r::Range) = r ## iteration -start(r::FloatRange) = 0 -done(r::FloatRange, i::Int) = length(r) <= i -next{T}(r::FloatRange{T}, i::Int) = - (convert(T, (r.start + i*r.step)/r.divisor), i+1) - start(r::LinSpace) = 1 done(r::LinSpace, i::Int) = length(r) < i -next{T}(r::LinSpace{T}, i::Int) = - (convert(T, ((r.len-i)*r.start + (i-1)*r.stop)/r.divisor), i+1) +function next(r::LinSpace, i::Int) + @_inline_meta + unsafe_getindex(r, i), i+1 +end start(r::StepRange) = oftype(r.start + r.step, r.start) next{T}(r::StepRange{T}, i) = (convert(T,i), i+r.step) @@ -485,6 +419,11 @@ done{T,S}(r::StepRange{T,S}, i) = isempty(r) | (i < min(r.start, r.stop)) | (i > done{T,S}(r::StepRange{T,S}, i::Integer) = isempty(r) | (i == oftype(i, r.stop) + r.step) +# see also twiceprecision.jl +start{T}(r::StepRangeLen{T}) = (unsafe_getindex(r, 1), 1) +next{T}(r::StepRangeLen{T}, s) = s[1], (T(s[1]+r.step), s[2]+1) +done{T}(r::StepRangeLen{T}, s) = s[2] > length(r) + start{T}(r::UnitRange{T}) = oftype(r.start + one(T), r.start) next{T}(r::AbstractUnitRange{T}, i) = (convert(T, i), i + one(T)) done{T}(r::AbstractUnitRange{T}, i) = i == oftype(i, r.stop) + one(T) @@ -529,16 +468,26 @@ function getindex{T}(v::Range{T}, i::Integer) ret end -function getindex{T}(r::FloatRange{T}, i::Integer) +function getindex(r::Union{StepRangeLen,LinSpace}, i::Integer) @_inline_meta @boundscheck checkbounds(r, i) - convert(T, (r.start + (i-1)*r.step)/r.divisor) + unsafe_getindex(r, i) +end + +# This is separate to make it useful even when running with --check-bounds=yes +function unsafe_getindex{T}(r::StepRangeLen{T}, i::Integer) + u = i - r.offset + T(r.ref + u*r.step) end -function getindex{T}(r::LinSpace{T}, i::Integer) +function unsafe_getindex(r::LinSpace, i::Integer) + lerpi.(i-1, r.lendiv, r.start, r.stop) +end + +function lerpi{T}(j::Integer, d::Integer, a::T, b::T) @_inline_meta - @boundscheck checkbounds(r, i) - convert(T, ((r.len-i)*r.start + (i-1)*r.stop)/r.divisor) + t = j/d + T((1-t)*a + t*b) end getindex(r::Range, ::Colon) = copy(r) @@ -571,30 +520,33 @@ function getindex{T<:Integer}(r::StepRange, s::Range{T}) range(st, step(r)*step(s), length(s)) end -function getindex(r::FloatRange, s::OrdinalRange) +function getindex{T<:Integer}(r::StepRangeLen, s::OrdinalRange{T}) @_inline_meta @boundscheck checkbounds(r, s) - FloatRange(r.start + (first(s)-1)*r.step, step(s)*r.step, length(s), r.divisor) + vfirst = unsafe_getindex(r, first(s)) + return StepRangeLen(vfirst, r.step*step(s), length(s)) end -function getindex{T}(r::LinSpace{T}, s::OrdinalRange) +function getindex{T<:Integer}(r::LinSpace, s::OrdinalRange{T}) @_inline_meta @boundscheck checkbounds(r, s) - sl::T = length(s) - ifirst = first(s) - ilast = last(s) - vfirst::T = ((r.len - ifirst) * r.start + (ifirst - 1) * r.stop) / r.divisor - vlast::T = ((r.len - ilast) * r.start + (ilast - 1) * r.stop) / r.divisor - return linspace(vfirst, vlast, sl) + vfirst = unsafe_getindex(r, first(s)) + vlast = unsafe_getindex(r, last(s)) + return linspace(vfirst, vlast, length(s)) end show(io::IO, r::Range) = print(io, repr(first(r)), ':', repr(step(r)), ':', repr(last(r))) show(io::IO, r::UnitRange) = print(io, repr(first(r)), ':', repr(last(r))) show(io::IO, r::OneTo) = print(io, "Base.OneTo(", r.stop, ")") -=={T<:Range}(r::T, s::T) = (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s)) -==(r::OrdinalRange, s::OrdinalRange) = (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s)) -=={T<:LinSpace}(r::T, s::T) = (first(r) == first(s)) & (length(r) == length(s)) & (last(r) == last(s)) +=={T<:Range}(r::T, s::T) = + (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s)) +==(r::OrdinalRange, s::OrdinalRange) = + (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s)) +=={T<:Union{StepRangeLen,LinSpace}}(r::T, s::T) = + (first(r) == first(s)) & (length(r) == length(s)) & (last(r) == last(s)) +=={T}(r::Union{StepRange{T},StepRangeLen{T,T}}, s::Union{StepRange{T},StepRangeLen{T,T}}) = + (first(r) == first(s)) & (last(r) == last(s)) & (step(r) == step(s)) function ==(r::Range, s::Range) lr = length(r) @@ -737,44 +689,43 @@ end ## linear operations on ranges ## -(r::OrdinalRange) = range(-first(r), -step(r), length(r)) --(r::FloatRange) = FloatRange(-r.start, -r.step, r.len, r.divisor) --(r::LinSpace) = LinSpace(-r.start, -r.stop, r.len, r.divisor) +-(r::StepRangeLen) = StepRangeLen(-r.ref, -r.step, length(r), r.offset) +-(r::LinSpace) = LinSpace(-r.start, -r.stop, length(r)) +(x::Real, r::AbstractUnitRange) = range(x + first(r), length(r)) -+(x::Real, r::Range) = (x+first(r)):step(r):(x+last(r)) -#+(x::Real, r::StepRange) = range(x + r.start, r.step, length(r)) -+(x::Real, r::FloatRange) = FloatRange(r.divisor*x + r.start, r.step, r.len, r.divisor) -function +{T}(x::Real, r::LinSpace{T}) - x2 = x * r.divisor / (r.len - 1) - LinSpace(x2 + r.start, x2 + r.stop, r.len, r.divisor) -end -+(r::Range, x::Real) = x + r -#+(r::FloatRange, x::Real) = x + r - --(x::Real, r::Range) = (x-first(r)):-step(r):(x-last(r)) --(x::Real, r::FloatRange) = FloatRange(r.divisor*x - r.start, -r.step, r.len, r.divisor) -function -(x::Real, r::LinSpace) - x2 = x * r.divisor / (r.len - 1) - LinSpace(x2 - r.start, x2 - r.stop, r.len, r.divisor) -end --(r::AbstractUnitRange, x::Real) = range(first(r)-x, length(r)) --(r::StepRange , x::Real) = range(r.start-x, r.step, length(r)) --(r::FloatRange, x::Real) = FloatRange(r.start - r.divisor*x, r.step, r.len, r.divisor) -function -(r::LinSpace, x::Real) - x2 = x * r.divisor / (r.len - 1) - LinSpace(r.start - x2, r.stop - x2, r.len, r.divisor) -end - -*(x::Real, r::OrdinalRange) = range(x*first(r), x*step(r), length(r)) -*(x::Real, r::FloatRange) = FloatRange(x*r.start, x*r.step, r.len, r.divisor) -*(x::Real, r::LinSpace) = LinSpace(x * r.start, x * r.stop, r.len, r.divisor) -*(r::Range, x::Real) = x * r -*(r::FloatRange, x::Real) = x * r -*(r::LinSpace, x::Real) = x * r - -/(r::OrdinalRange, x::Real) = range(first(r)/x, step(r)/x, length(r)) -/(r::FloatRange, x::Real) = FloatRange(r.start/x, r.step/x, r.len, r.divisor) -/(r::LinSpace, x::Real) = LinSpace(r.start / x, r.stop / x, r.len, r.divisor) +# For #18336 we need to prevent promotion of the step type: ++(x::Number, r::AbstractUnitRange) = range(x + first(r), step(r), length(r)) ++(x::Number, r::Range) = (x+first(r)):step(r):(x+last(r)) +function +(x::Number, r::StepRangeLen) + newref = x + r.ref + StepRangeLen{eltype(newref),typeof(newref),typeof(r.step)}(newref, r.step, length(r), r.offset) +end +function +(x::Number, r::LinSpace) + LinSpace(x + r.start, x + r.stop, r.len) +end ++(r::Range, x::Number) = x + r # assumes addition is commutative + +-(x::Number, r::Range) = (x-first(r)):-step(r):(x-last(r)) +-(x::Number, r::StepRangeLen) = +(x, -r) +function -(x::Number, r::LinSpace) + LinSpace(x - r.start, x - r.stop, r.len) +end + +-(r::Range, x::Number) = +(-x, r) + +*(x::Number, r::Range) = range(x*first(r), x*step(r), length(r)) +*(x::Number, r::StepRangeLen) = StepRangeLen(x*r.ref, x*r.step, length(r), r.offset) +*(x::Number, r::LinSpace) = LinSpace(x * r.start, x * r.stop, r.len) +# separate in case of noncommutative multiplication +*(r::Range, x::Number) = range(first(r)*x, step(r)*x, length(r)) +*(r::StepRangeLen, x::Number) = StepRangeLen(r.ref*x, r.step*x, length(r), r.offset) +*(r::LinSpace, x::Number) = LinSpace(r.start * x, r.stop * x, r.len) + +/(r::Range, x::Number) = range(first(r)/x, step(r)/x, length(r)) +/(r::StepRangeLen, x::Number) = StepRangeLen(r.ref/x, r.step/x, length(r), r.offset) +/(r::LinSpace, x::Number) = LinSpace(r.start / x, r.stop / x, r.len) + +/(x::Number, r::Range) = [ x/y for y=r ] promote_rule{T1,T2}(::Type{UnitRange{T1}},::Type{UnitRange{T2}}) = UnitRange{promote_type(T1,T2)} @@ -802,56 +753,37 @@ convert{T1,T2}(::Type{StepRange{T1,T2}}, r::Range) = convert{T}(::Type{StepRange}, r::AbstractUnitRange{T}) = StepRange{T,T}(first(r), step(r), last(r)) -promote_rule{T1,T2}(::Type{FloatRange{T1}},::Type{FloatRange{T2}}) = - FloatRange{promote_type(T1,T2)} -convert{T<:AbstractFloat}(::Type{FloatRange{T}}, r::FloatRange{T}) = r -convert{T<:AbstractFloat}(::Type{FloatRange{T}}, r::FloatRange) = - FloatRange{T}(r.start,r.step,r.len,r.divisor) - -promote_rule{F,OR<:OrdinalRange}(::Type{FloatRange{F}}, ::Type{OR}) = - FloatRange{promote_type(F,eltype(OR))} -convert{T<:AbstractFloat}(::Type{FloatRange{T}}, r::OrdinalRange) = - FloatRange{T}(first(r), step(r), length(r), one(T)) -convert{T}(::Type{FloatRange}, r::OrdinalRange{T}) = - FloatRange{typeof(float(first(r)))}(first(r), step(r), length(r), one(T)) +promote_rule{T1,T2,R1,R2,S1,S2}(::Type{StepRangeLen{T1,R1,S1}},::Type{StepRangeLen{T2,R2,S2}}) = + StepRangeLen{promote_type(T1,T2), promote_type(R1,R2), promote_type(S1,S2)} +convert{T,R,S}(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen{T,R,S}) = r +convert{T,R,S}(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen) = + StepRangeLen{T,R,S}(convert(R, r.ref), convert(S, r.step), length(r), r.offset) +convert{T}(::Type{StepRangeLen{T}}, r::StepRangeLen) = + StepRangeLen(convert(T, r.ref), convert(T, r.step), length(r), r.offset) + +promote_rule{T,R,S,OR<:Range}(::Type{StepRangeLen{T,R,S}}, ::Type{OR}) = + StepRangeLen{promote_type(T,eltype(OR)),promote_type(R,eltype(OR)),promote_type(S,eltype(OR))} +convert{T,R,S}(::Type{StepRangeLen{T,R,S}}, r::Range) = + StepRangeLen{T,R,S}(R(first(r)), S(step(r)), length(r)) +convert{T}(::Type{StepRangeLen{T}}, r::Range) = + StepRangeLen(T(first(r)), T(step(r)), length(r)) +convert(::Type{StepRangeLen}, r::Range) = convert(StepRangeLen{eltype(r)}, r) promote_rule{T1,T2}(::Type{LinSpace{T1}},::Type{LinSpace{T2}}) = LinSpace{promote_type(T1,T2)} -convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::LinSpace{T}) = r -convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::LinSpace) = - LinSpace{T}(r.start, r.stop, r.len, r.divisor) - -promote_rule{F,OR<:OrdinalRange}(::Type{LinSpace{F}}, ::Type{OR}) = - LinSpace{promote_type(F,eltype(OR))} -convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::OrdinalRange) = - linspace(convert(T, first(r)), convert(T, last(r)), convert(T, length(r))) -convert{T}(::Type{LinSpace}, r::OrdinalRange{T}) = - convert(LinSpace{typeof(float(first(r)))}, r) - -# Promote FloatRange to LinSpace -promote_rule{F,OR<:FloatRange}(::Type{LinSpace{F}}, ::Type{OR}) = - LinSpace{promote_type(F,eltype(OR))} -convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::FloatRange) = - linspace(convert(T, first(r)), convert(T, last(r)), convert(T, length(r))) -convert{T<:AbstractFloat}(::Type{LinSpace}, r::FloatRange{T}) = +convert{T}(::Type{LinSpace{T}}, r::LinSpace{T}) = r +convert{T}(::Type{LinSpace{T}}, r::Range) = + LinSpace{T}(first(r), last(r), length(r)) +convert{T}(::Type{LinSpace}, r::Range{T}) = convert(LinSpace{T}, r) +promote_rule{T,OR<:OrdinalRange}(::Type{LinSpace{T}}, ::Type{OR}) = + LinSpace{promote_type(T,eltype(OR))} -# +/- of ranges is defined in operators.jl (to be able to use @eval etc.) - -## non-linear operations on ranges and fallbacks for non-real numbers ## - -+(x::Number, r::Range) = [ x+y for y=r ] -+(r::Range, y::Number) = [ x+y for x=r ] - --(x::Number, r::Range) = [ x-y for y=r ] --(r::Range, y::Number) = [ x-y for x=r ] +promote_rule{L,T,R,S}(::Type{LinSpace{L}}, ::Type{StepRangeLen{T,R,S}}) = + StepRangeLen{promote_type(L,T),promote_type(L,R),promote_type(L,S)} -*(x::Number, r::Range) = [ x*y for y=r ] -*(r::Range, y::Number) = [ x*y for x=r ] - -/(x::Number, r::Range) = [ x/y for y=r ] -/(r::Range, y::Number) = [ x/y for x=r ] +# +/- of ranges is defined in operators.jl (to be able to use @eval etc.) ## concatenation ## @@ -873,8 +805,8 @@ convert{T}(::Type{Array{T,1}}, r::Range{T}) = vcat(r) collect(r::Range) = vcat(r) reverse(r::OrdinalRange) = colon(last(r), -step(r), first(r)) -reverse(r::FloatRange) = FloatRange(r.start + (r.len-1)*r.step, -r.step, r.len, r.divisor) -reverse(r::LinSpace) = LinSpace(r.stop, r.start, r.len, r.divisor) +reverse(r::StepRangeLen) = StepRangeLen(r.ref, -r.step, length(r), length(r)-r.offset+1) +reverse(r::LinSpace) = LinSpace(r.stop, r.start, length(r)) ## sorting ## @@ -896,17 +828,6 @@ function sum{T<:Real}(r::Range{T}) : (step(r) * l) * ((l-1)>>1)) end -function sum(r::FloatRange) - l = length(r) - if iseven(l) - s = r.step * (l-1) * (l>>1) - else - s = (r.step * l) * ((l-1)>>1) - end - return (l * r.start + s)/r.divisor -end - - function mean{T<:Real}(r::Range{T}) isempty(r) && throw(ArgumentError("mean of an empty range is undefined")) (first(r) + last(r)) / 2 diff --git a/base/rational.jl b/base/rational.jl index 82ac9631071f7..05cecc1aa67c0 100644 --- a/base/rational.jl +++ b/base/rational.jl @@ -401,3 +401,7 @@ function ^{T<:Rational}(z::Complex{T}, n::Integer) end iszero(x::Rational) = iszero(numerator(x)) + +function lerpi(j::Integer, d::Integer, a::Rational, b::Rational) + ((d-j)*a)/d + (j*b)/d +end diff --git a/base/sysimg.jl b/base/sysimg.jl index 66c720647737e..5a42220deb2cb 100644 --- a/base/sysimg.jl +++ b/base/sysimg.jl @@ -56,7 +56,9 @@ include("options.jl") # core operations & types include("promotion.jl") include("tuple.jl") +include("traits.jl") include("range.jl") +include("twiceprecision.jl") include("expr.jl") include("error.jl") diff --git a/base/traits.jl b/base/traits.jl new file mode 100644 index 0000000000000..62bc9ecb4d6da --- /dev/null +++ b/base/traits.jl @@ -0,0 +1,22 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +## numeric/object traits +# trait for objects that have an ordering +abstract TypeOrder +immutable HasOrder <: TypeOrder end +immutable Unordered <: TypeOrder end + +(::Type{TypeOrder})(instance) = TypeOrder(typeof(instance)) +(::Type{TypeOrder}){T<:Real}(::Type{T}) = HasOrder() +(::Type{TypeOrder}){T}(::Type{T}) = Unordered() + +# trait for objects that support arithmetic +abstract TypeArithmetic +immutable ArithmeticRounds <: TypeArithmetic end # least significant bits can be lost +immutable ArithmeticOverflows <: TypeArithmetic end # most significant bits can be lost +immutable ArithmeticUnknown <: TypeArithmetic end + +(::Type{TypeArithmetic})(instance) = TypeArithmetic(typeof(instance)) +(::Type{TypeArithmetic}){T<:AbstractFloat}(::Type{T}) = ArithmeticRounds() +(::Type{TypeArithmetic}){T<:Integer}(::Type{T}) = ArithmeticOverflows() +(::Type{TypeArithmetic}){T}(::Type{T}) = ArithmeticUnknown() diff --git a/base/twiceprecision.jl b/base/twiceprecision.jl new file mode 100644 index 0000000000000..80d1917ea599f --- /dev/null +++ b/base/twiceprecision.jl @@ -0,0 +1,508 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +# Twice-precision arithmetic. + +# Necessary for creating nicely-behaved ranges like r = 0.1:0.1:0.3 +# that return r[3] == 0.3. Otherwise, we have roundoff error due to +# 0.1 + 2*0.1 = 0.30000000000000004 + +""" + TwicePrecision{T}(hi::T, lo::T) + TwicePrecision{T}((num, denom)) + +A number with twice the precision of `T`, e.g., quad-precision if `T = +Float64`. `hi` represents the high bits (most significant bits) and +`lo` the low bits (least significant bits). Rational values +`num//denom` can be approximated conveniently using the syntax +`TwicePrecision{T}((num, denom))`. + +When used with `T<:AbstractFloat` to construct an exact +`StepRangeLen`, `ref` should be the range element with smallest +magnitude and `offset` set to the corresponding index. For +efficiency, multiplication of `step` by the index is not performed at +twice precision: `step.hi` should have enough trailing zeros in its +`bits` representation that `(0:len-1)*step.hi` is exact (has no +roundoff error). If `step` has an exact rational representation +`num//denom`, then you can construct `step` using + + step = TwicePrecision{T}((num, denom), nb) + +where `nb` is the number of trailing zero bits of `step.hi`. For +ranges, you can set `nb = ceil(Int, log2(len-1))`. +""" +immutable TwicePrecision{T} + hi::T # most significant bits + lo::T # least significant bits +end + +function (::Type{TwicePrecision{T}}){T,I}(nd::Tuple{I,I}) + n, d = nd + TwicePrecision{T}(n, zero(T)) / d +end + +function (::Type{TwicePrecision{T}}){T,I}(nd::Tuple{I,I}, nb::Integer) + twiceprecision(TwicePrecision{T}(nd), nb) +end + +function twiceprecision{T<:Number}(val::T, nb::Integer) + hi = truncbits(val, nb) + TwicePrecision{T}(hi, val - hi) +end + +function twiceprecision{T<:Number}(val::TwicePrecision{T}, nb::Integer) + hi = truncbits(val.hi, nb) + TwicePrecision{T}(hi, (val.hi - hi) + val.lo) +end + +nbitslen(r::StepRangeLen) = nbitslen(eltype(r), length(r), r.offset) +nbitslen(::Type{Float64}, len, offset) = min(26, nbitslen(len, offset)) +nbitslen(::Type{Float32}, len, offset) = min(12, nbitslen(len, offset)) +nbitslen(::Type{Float16}, len, offset) = min(5, nbitslen(len, offset)) +nbitslen(len, offset) = len < 2 ? 0 : ceil(Int, log2(max(offset-1, len-offset))) + +eltype{T}(::Type{TwicePrecision{T}}) = T + +promote_rule{R,S}(::Type{TwicePrecision{R}}, ::Type{TwicePrecision{S}}) = + TwicePrecision{promote_type(R,S)} +promote_rule{R,S}(::Type{TwicePrecision{R}}, ::Type{S}) = + TwicePrecision{promote_type(R,S)} + +convert{T}(::Type{TwicePrecision{T}}, x::TwicePrecision{T}) = x +convert{T}(::Type{TwicePrecision{T}}, x::TwicePrecision) = + TwicePrecision{T}(convert(T, x.hi), convert(T, x.lo)) + +convert{T<:Number}(::Type{T}, x::TwicePrecision{T}) = convert(T, x.hi + x.lo) +convert{T}(::Type{TwicePrecision{T}}, x::Number) = TwicePrecision{T}(convert(T, x), zero(T)) + +float{T<:AbstractFloat}(x::TwicePrecision{T}) = x +float(x::TwicePrecision) = TwicePrecision(float(x.hi), float(x.lo)) + +big(x::TwicePrecision) = big(x.hi) + big(x.lo) + +-(x::TwicePrecision) = TwicePrecision(-x.hi, -x.lo) + +zero{T}(::Type{TwicePrecision{T}}) = TwicePrecision{T}(0, 0) + +## StepRangeLen + +# If using TwicePrecision numbers, deliberately force user to specify offset +StepRangeLen{T}(ref::TwicePrecision{T}, step::TwicePrecision{T}, len::Integer, offset::Integer) = + StepRangeLen{T,TwicePrecision{T},TwicePrecision{T}}(ref, step, len, offset) + +# Construct range for rational start=start_n/den, step=step_n/den +function floatrange{T}(::Type{T}, start_n::Integer, step_n::Integer, len::Integer, den::Integer) + if len < 2 + return StepRangeLen(TwicePrecision{T}((start_n, den)), + TwicePrecision{T}((step_n, den)), Int(len), 1) + end + # index of smallest-magnitude value + imin = clamp(round(Int, -start_n/step_n+1), 1, Int(len)) + # Compute smallest-magnitude element to 2x precision + ref_n = start_n+(imin-1)*step_n # this shouldn't overflow, so don't check + nb = nbitslen(T, len, imin) + StepRangeLen(TwicePrecision{T}((ref_n, den)), + TwicePrecision{T}((step_n, den), nb), Int(len), imin) +end + +function floatrange(a::AbstractFloat, st::AbstractFloat, len::Real, divisor::AbstractFloat) + T = promote_type(typeof(a), typeof(st), typeof(divisor)) + m = maxintfloat(T) + if abs(a) <= m && abs(st) <= m && abs(divisor) <= m + ia, ist, idivisor = round(Int, a), round(Int, st), round(Int, divisor) + if ia == a && ist == st && idivisor == divisor + # We can return the high-precision range + return floatrange(T, ia, ist, Int(len), idivisor) + end + end + # Fallback (misses the opportunity to set offset different from 1, + # but otherwise this is still high-precision) + StepRangeLen(TwicePrecision{T}((a,divisor)), + TwicePrecision{T}((st,divisor), nbitslen(T, len, 1)), Int(len), 1) +end + +function colon{T<:Union{Float16,Float32,Float64}}(start::T, step::T, stop::T) + step == 0 && throw(ArgumentError("range step cannot be zero")) + len = max(0, floor(Int, (stop-start)/step) + 1) + # Because len might be too small by 1 due to roundoff error, let's + # see if the inputs have exact rational approximations (and if so, + # perform all computations in terms of the rationals) + step_n, step_d = rat(step) + if step_d != 0 && T(step_n/step_d) == step + start_n, start_d = rat(start) + stop_n, stop_d = rat(stop) + if start_d != 0 && stop_d != 0 && + T(start_n/start_d) == start && T(stop_n/stop_d) == stop + den = lcm(start_d, step_d) # use same denominator for start and step + m = maxintfloat(T) + if den != 0 && abs(start*den) <= m && abs(step*den) <= m && # will round succeed? + rem(den, start_d) == 0 && rem(den, step_d) == 0 # check lcm overflow + start_n = round(Int, start*den) + step_n = round(Int, step*den) + len = max(0, div(den*stop_n - stop_d*start_n + step_n*stop_d, step_n*stop_d)) + # Integer ops could overflow, so check that this makes sense + if isbetween(start, start + (len-1)*step, stop + step/2) && + !isbetween(start, start + len*step, stop) + # Return a 2x precision range + return floatrange(T, start_n, step_n, len, den) + end + end + end + end + # Fallback, taking start and step literally + StepRangeLen(TwicePrecision(start, zero(T)), twiceprecision(step, nbitslen(T, len, 1)), len) +end + +function range{T<:Union{Float16,Float32,Float64}}(a::T, st::T, len::Integer) + start_n, start_d = rat(a) + step_n, step_d = rat(st) + if start_d != 0 && step_d != 0 && + T(start_n/start_d) == a && T(step_n/step_d) == st + den = lcm(start_d, step_d) + m = maxintfloat(T) + if abs(den*a) <= m && abs(den*st) <= m && + rem(den, start_d) == 0 && rem(den, step_d) == 0 + start_n = round(Int, den*a) + step_n = round(Int, den*st) + return floatrange(T, start_n, step_n, len, den) + end + end + StepRangeLen(TwicePrecision(a, zero(T)), TwicePrecision(st, zero(T)), len) +end + +step{T,R,S<:TwicePrecision}(r::StepRangeLen{T,R,S}) = convert(eltype(S), r.step) + +start{T,R<:TwicePrecision,S<:TwicePrecision}(r::StepRangeLen{T,R,S}) = 1 +done{T,R<:TwicePrecision,S<:TwicePrecision}(r::StepRangeLen{T,R,S}, i::Int) = length(r) < i +function next{T,R<:TwicePrecision,S<:TwicePrecision}(r::StepRangeLen{T,R,S}, i::Int) + @_inline_meta + unsafe_getindex(r, i), i+1 +end + +# This assumes that r.step has already been split so that (0:len-1)*r.step.hi is exact +function unsafe_getindex{T,R<:TwicePrecision,S<:TwicePrecision}(r::StepRangeLen{T,R,S}, i::Integer) + # Very similar to _getindex_hiprec, but optimized to avoid a 2nd call to add2 + @_inline_meta + u = i - r.offset + shift_hi, shift_lo = u*r.step.hi, u*r.step.lo + x_hi, x_lo = add2(r.ref.hi, shift_hi) + T(x_hi + (x_lo + (shift_lo + r.ref.lo))) +end + +function _getindex_hiprec{T,R<:TwicePrecision,S<:TwicePrecision}(r::StepRangeLen{T,R,S}, i::Integer) + u = i - r.offset + shift_hi, shift_lo = u*r.step.hi, u*r.step.lo + x_hi, x_lo = add2(r.ref.hi, shift_hi) + x_hi, x_lo = add2(x_hi, x_lo + (shift_lo + r.ref.lo)) + TwicePrecision(x_hi, x_lo) +end + +function getindex{T,R<:TwicePrecision,S<:TwicePrecision,I<:Integer}(r::StepRangeLen{T,R,S}, s::OrdinalRange{I}) + @_inline_meta + @boundscheck checkbounds(r, s) + soffset = 1 + round(Int, (r.offset - first(s))/step(s)) + soffset = clamp(soffset, 1, length(s)) + ioffset = start(s) + (soffset-1)*step(s) + if step(s) == 1 || length(s) < 2 + newstep = r.step + else + newstep = twiceprecision(r.step*step(s), nbitslen(T, length(s), soffset)) + end + if ioffset == r.offset + StepRangeLen(r.ref, newstep, length(s), max(1,soffset)) + else + StepRangeLen(r.ref + (ioffset-r.offset)*r.step, newstep, length(s), max(1,soffset)) + end +end + +*{T<:Real,R<:TwicePrecision}(x::Real, r::StepRangeLen{T,R}) = + StepRangeLen(x*r.ref, twiceprecision(x*r.step, nbitslen(r)), length(r), r.offset) +*{T<:Real,R<:TwicePrecision}(r::StepRangeLen{T,R}, x::Real) = x*r +/{T<:Real,R<:TwicePrecision}(r::StepRangeLen{T,R}, x::Real) = + StepRangeLen(r.ref/x, twiceprecision(r.step/x, nbitslen(r)), length(r), r.offset) + +convert{T<:AbstractFloat,R<:TwicePrecision,S<:TwicePrecision}(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen{T,R,S}) = r + +convert{T<:AbstractFloat,R<:TwicePrecision,S<:TwicePrecision}(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen) = + _convertSRL(StepRangeLen{T,R,S}, r) + +convert{T<:Union{Float16,Float32,Float64}}(::Type{StepRangeLen{T}}, r::StepRangeLen) = + _convertSRL(StepRangeLen{T,TwicePrecision{T},TwicePrecision{T}}, r) + +convert{T<:Union{Float16,Float32,Float64}}(::Type{StepRangeLen{T}}, r::Range) = + _convertSRL(StepRangeLen{T,TwicePrecision{T},TwicePrecision{T}}, r) + +function _convertSRL{T,R,S,I<:Integer}(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen{I}) + StepRangeLen{T,R,S}(R(r.ref), S(r.step), length(r), r.offset) +end + +function _convertSRL{T,R,S,I<:Integer}(::Type{StepRangeLen{T,R,S}}, r::Range{I}) + StepRangeLen{T,R,S}(R(first(r)), S(step(r)), length(r)) +end + +function _convertSRL{T,R,S,U}(::Type{StepRangeLen{T,R,S}}, r::Range{U}) + # if start and step have a rational approximation in the old type, + # then we transfer that rational approximation to the new type + f, s = first(r), step(r) + start_n, start_d = rat(f) + step_n, step_d = rat(s) + if start_d != 0 && step_d != 0 && + U(start_n/start_d) == f && U(step_n/step_d) == s + den = lcm(start_d, step_d) + m = maxintfloat(T) + if den != 0 && abs(f*den) <= m && abs(s*den) <= m && + rem(den, start_d) == 0 && rem(den, step_d) == 0 + start_n = round(Int, f*den) + step_n = round(Int, s*den) + return floatrange(T, start_n, step_n, length(r), den) + end + end + __convertSRL(StepRangeLen{T,R,S}, r) +end + +function __convertSRL{T,R,S,U}(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen{U}) + StepRangeLen{T,R,S}(R(r.ref), S(r.step), length(r), r.offset) +end +function __convertSRL{T,R,S,U}(::Type{StepRangeLen{T,R,S}}, r::Range{U}) + StepRangeLen{T,R,S}(R(first(r)), S(step(r)), length(r)) +end + +function sum(r::StepRangeLen) + l = length(r) + # Compute the contribution of step over all indexes. + # Indexes on opposite side of r.offset contribute with opposite sign, + # r.step * (sum(1:np) - sum(1:nn)) + np, nn = l - r.offset, r.offset - 1 # positive, negative + # To prevent overflow in sum(1:n), multiply its factors by the step + sp, sn = sumpair(np), sumpair(nn) + tp = _prod(r.step, sp[1], sp[2]) + tn = _prod(r.step, sn[1], sn[2]) + s_hi, s_lo = add2(tp.hi, -tn.hi) + s_lo += tp.lo - tn.lo + # Add in contributions of ref + ref = r.ref * l + sm_hi, sm_lo = add2(s_hi, ref.hi) + add2(sm_hi, sm_lo + ref.lo)[1] +end + +# sum(1:n) as a product of two integers +sumpair(n::Integer) = iseven(n) ? (n+1, n>>1) : (n, (n+1)>>1) + +function +{T,R<:TwicePrecision}(r1::StepRangeLen{T,R}, r2::StepRangeLen{T,R}) + len = length(r1) + (len == length(r2) || + throw(DimensionMismatch("argument dimensions must match"))) + if r1.offset == r2.offset + imid = r1.offset + ref = r1.ref + r2.ref + else + imid = round(Int, (r1.offset+r2.offset)/2) + ref1mid = _getindex_hiprec(r1, imid) + ref2mid = _getindex_hiprec(r2, imid) + ref = ref1mid + ref2mid + end + step = twiceprecision(r1.step + r2.step, nbitslen(T, len, imid)) + StepRangeLen{T,typeof(ref),typeof(step)}(ref, step, len, imid) +end + +## LinSpace + +# For Float16, Float32, and Float64, linspace returns a StepRangeLen +function linspace{T<:Union{Float16,Float32,Float64}}(start::T, stop::T, len::Integer) + len < 2 && return _linspace1(T, start, stop, len) + # Attempt to find exact rational approximations + start_n, start_d = rat(start) + stop_n, stop_d = rat(stop) + if start_d != 0 && stop_d != 0 + den = lcm(start_d, stop_d) + m = maxintfloat(T) + if den != 0 && abs(den*start) <= m && abs(den*stop) <= m + start_n = round(Int, den*start) + stop_n = round(Int, den*stop) + if T(start_n/den) == start && T(stop_n/den) == stop + return linspace(T, start_n, stop_n, len, den) + end + end + end + _linspace(start, stop, len) +end + +function _linspace{T<:Union{Float16,Float32,Float64}}(start::T, stop::T, len::Integer) + (isfinite(start) && isfinite(stop)) || throw(ArgumentError("start and stop must be finite, got $start and $stop")) + # Find the index that returns the smallest-magnitude element + Δ, Δfac = stop-start, 1 + if !isfinite(Δ) # handle overflow for large endpoints + Δ, Δfac = stop/len - start/len, Int(len) + end + tmin = -(start/Δ)/Δfac # interpolation t such that return value is 0 + imin = round(Int, tmin*(len-1)+1) + if 1 < imin < len + # The smallest-magnitude element is in the interior + t = (imin-1)/(len-1) + ref = T((1-t)*start + t*stop) + step = imin-1 < len-imin ? (ref-start)/(imin-1) : (stop-ref)/(len-imin) + elseif imin <= 1 + imin = 1 + ref = start + step = (Δ/(len-1))*Δfac + else + imin = Int(len) + ref = stop + step = (Δ/(len-1))*Δfac + end + if len == 2 && !isfinite(step) + # For very large endpoints where step overflows, exploit the + # split-representation to handle the overflow + return StepRangeLen(TwicePrecision(start, zero(T)), + TwicePrecision(-start, stop), 2) + end + # 2x calculations to get high precision endpoint matching while also + # preventing overflow in ref_hi+(i-offset)*step_hi + m, k = prevfloat(realmax(T)), max(imin-1, len-imin) + step_hi_pre = clamp(step, max(-(m+ref)/k, (-m+ref)/k), min((m-ref)/k, (m+ref)/k)) + nb = nbitslen(T, len, imin) + step_hi = truncbits(step_hi_pre, nb) + x1_hi, x1_lo = add2((1-imin)*step_hi, ref) + x2_hi, x2_lo = add2((len-imin)*step_hi, ref) + a, b = (start - x1_hi) - x1_lo, (stop - x2_hi) - x2_lo + step_lo = (b - a)/(len - 1) + ref_lo = a - (1 - imin)*step_lo + StepRangeLen(TwicePrecision(ref, ref_lo), TwicePrecision(step_hi, step_lo), Int(len), imin) +end + +# linspace for rational numbers, start = start_n/den, stop = stop_n/den +# Note this returns a StepRangeLen +function linspace{T}(::Type{T}, start_n::Integer, stop_n::Integer, len::Integer, den::Integer) + len < 2 && return _linspace1(T, start_n/den, stop_n/den, len) + start_n == stop_n && return StepRangeLen(TwicePrecision{T}((start_n, den)), zero(TwicePrecision{T}), len) + tmin = -start_n/(Float64(stop_n) - Float64(start_n)) + imin = round(Int, tmin*(len-1)+1) + imin = clamp(imin, 1, Int(len)) + # Compute (1-t)*a and t*b separately in 2x precision (itp = interpolant)... + dent = (den, len-1) # represent products as a tuple to eliminate risk of overflow + start_itp = proddiv(T, (len-imin, start_n), dent) + stop_itp = proddiv(T, (imin-1, stop_n), dent) + # ...and then combine them to make ref + ref = start_itp + stop_itp + # Compute step to 2x precision without risking overflow... + rend = proddiv(T, (stop_n,), dent) + rbeg = proddiv(T, (-start_n,), dent) + step = twiceprecision(rbeg + rend, nbitslen(T, len, imin)) # ...and truncate hi-bits as needed + StepRangeLen(ref, step, Int(len), imin) +end + +# For len < 2 +function _linspace1{T}(::Type{T}, start, stop, len::Integer) + len >= 0 || throw(ArgumentError("linspace($start, $stop, $len): negative length")) + if len <= 1 + len == 1 && (start == stop || throw(ArgumentError("linspace($start, $stop, $len): endpoints differ"))) + # Ensure that first(r)==start and last(r)==stop even for len==0 + return StepRangeLen(TwicePrecision(start, zero(T)), TwicePrecision(start, -stop), len, 1) + end + throw(ArgumentError("should only be called for len < 2, got $len")) +end + +### Numeric utilities + +# Approximate x with a rational representation. Guaranteed to return, +# but not guaranteed to return a precise answer. +# https://en.wikipedia.org/wiki/Continued_fraction#Best_rational_approximations +function rat(x) + y = x + a = d = 1 + b = c = 0 + m = maxintfloat(narrow(typeof(x))) + while abs(y) <= m + f = trunc(Int,y) + y -= f + a, c = f*a + c, a + b, d = f*b + d, b + max(abs(a), abs(b)) <= convert(Int,m) || return c, d + oftype(x,a)/oftype(x,b) == x && break + y = inv(y) + end + return a, b +end + +narrow(::Type{Float64}) = Float32 +narrow(::Type{Float32}) = Float16 +narrow(::Type{Float16}) = Float16 + +function add2{T<:Number}(u::T, v::T) + @_inline_meta + u, v = ifelse(abs(v) > abs(u), (v, u), (u, v)) + w = u + v + w, (u-w) + v +end + +add2(u, v) = _add2(promote(u, v)...) +_add2{T<:Number}(u::T, v::T) = add2(u, v) +_add2(u, v) = error("$u::$(typeof(u)) and $v::$(typeof(v)) cannot be promoted to a common type") + +function +(x::TwicePrecision, y::Number) + s_hi, s_lo = add2(x.hi, y) + TwicePrecision(s_hi, s_lo+x.lo) +end ++(x::Number, y::TwicePrecision) = y+x + +function +{T}(x::TwicePrecision{T}, y::TwicePrecision{T}) + r = x.hi + y.hi + s = abs(x.hi) > abs(y.hi) ? (((x.hi - r) + y.hi) + y.lo) + x.lo : (((y.hi - r) + x.hi) + x.lo) + y.lo + TwicePrecision(r, s) +end ++(x::TwicePrecision, y::TwicePrecision) = _add2(promote(x, y)...) +_add2{T<:TwicePrecision}(x::T, y::T) = x + y +_add2(x::TwicePrecision, y::TwicePrecision) = TwicePrecision(x.hi+y.hi, x.lo+y.lo) + +function *(x::TwicePrecision, v::Integer) + v == 0 && return TwicePrecision(x.hi*v, x.lo*v) + nb = ceil(Int, log2(abs(v))) + u = truncbits(x.hi, nb) + y_hi, y_lo = add2(u*v, ((x.hi-u) + x.lo)*v) + TwicePrecision(y_hi, y_lo) +end + +function _mul2{T<:Union{Float16,Float32,Float64}}(x::TwicePrecision{T}, v::T) + v == 0 && return TwicePrecision(T(0), T(0)) + xhh, xhl = splitprec(x.hi) + vh, vl = splitprec(v) + y_hi, y_lo = add2(xhh*vh, xhh*vl + xhl*vh) + TwicePrecision(y_hi, y_lo + xhl*vl + x.lo*v) +end + +_mul2(x::TwicePrecision, v::Number) = TwicePrecision(x.hi*v, x.lo*v) + +function *{R,S<:Number}(x::TwicePrecision{R}, v::S) + T = promote_type(R, S) + _mul2(convert(TwicePrecision{T}, x), convert(T, v)) +end + +*(v::Number, x::TwicePrecision) = x*v + +function /(x::TwicePrecision, v::Number) + hi = x.hi/v + w = TwicePrecision(hi, zero(hi)) * v + lo = (((x.hi - w.hi) - w.lo) + x.lo)/v + y_hi, y_lo = add2(hi, lo) + TwicePrecision(y_hi, y_lo) +end + +# hi-precision version of prod(num)/prod(den) +# num and den are tuples to avoid risk of overflow +function proddiv(T, num, den) + @_inline_meta + t = TwicePrecision(T(num[1]), zero(T)) + t = _prod(t, tail(num)...) + _divt(t, den...) +end +function _prod(t::TwicePrecision, x, y...) + @_inline_meta + _prod(t * x, y...) +end +_prod(t::TwicePrecision) = t +function _divt(t::TwicePrecision, x, y...) + @_inline_meta + _divt(t / x, y...) +end +_divt(t::TwicePrecision) = t + +isbetween(a, x, b) = a <= x <= b || b <= x <= a diff --git a/doc/src/stdlib/math.md b/doc/src/stdlib/math.md index eab15c7daccb3..dd36ec60d3572 100644 --- a/doc/src/stdlib/math.md +++ b/doc/src/stdlib/math.md @@ -34,6 +34,7 @@ Base.:(>>>) Base.colon Base.range Base.OneTo +Base.StepRangeLen Base.:(==) Base.:(!=) Base.:(!==) diff --git a/test/core.jl b/test/core.jl index cee8edac22a62..69674e5229ed0 100644 --- a/test/core.jl +++ b/test/core.jl @@ -1938,7 +1938,7 @@ end # a method specificity issue c99991{T}(::Type{T},x::T) = 0 -c99991{T}(::Type{UnitRange{T}},x::FloatRange{T}) = 1 +c99991{T}(::Type{UnitRange{T}},x::StepRangeLen{T}) = 1 c99991{T}(::Type{UnitRange{T}},x::Range{T}) = 2 @test c99991(UnitRange{Float64}, 1.0:2.0) == 1 @test c99991(UnitRange{Int}, 1:2) == 2 diff --git a/test/ranges.jl b/test/ranges.jl index 4d8161da0a627..ef7d28a89203b 100644 --- a/test/ranges.jl +++ b/test/ranges.jl @@ -9,12 +9,19 @@ @test length(2:.2:1) == 0 @test length(2.:.2:1.) == 0 +@inferred(colon(10, 1, 0)) +@inferred(colon(1, .2, 2)) +@inferred(colon(1., .2, 2.)) +@inferred(colon(2, -.2, 1)) +@inferred(colon(1, 0)) +@inferred(colon(0.0, -0.5)) + @test length(1:0) == 0 @test length(0.0:-0.5) == 0 @test length(1:2:0) == 0 -L32 = linspace(Int32(1), Int32(4), 4) -L64 = linspace(Int64(1), Int64(4), 4) -@test L32[1] == 1 && L64[1] == 1 +L32 = @inferred(linspace(Int32(1), Int32(4), 4)) +L64 = @inferred(linspace(Int64(1), Int64(4), 4)) +@test @inferred(L32[1]) === 1.0 && @inferred(L64[1]) === 1.0 @test L32[2] == 2 && L64[2] == 2 @test L32[3] == 3 && L64[3] == 3 @test L32[4] == 4 && L64[4] == 4 @@ -31,16 +38,16 @@ r = 5:-1:1 @test length(1.1:1.3:3) == 2 @test length(1:1:1.8) == 1 -@test (1:5)[1:4] == 1:4 +@test @inferred((0.1:0.1:0.3)[2]) === 0.2 +@test @inferred((0.1f0:0.1f0:0.3f0)[2]) === 0.2f0 + +@test @inferred((1:5)[1:4]) === 1:4 +@test @inferred((1.0:5)[1:4]) === 1.0:4 @test (2:6)[1:4] == 2:5 -@test (1:6)[2:5] == 2:5 -@test typeof((1:6)[2:5]) == typeof(2:5) -@test (1:6)[2:2:5] == 2:2:4 -@test typeof((1:6)[2:2:5]) == typeof(2:2:4) -@test (1:2:13)[2:6] == 3:2:11 -@test typeof((1:2:13)[2:6]) == typeof(3:2:11) -@test (1:2:13)[2:3:7] == 3:6:13 -@test typeof((1:2:13)[2:3:7]) == typeof(3:6:13) +@test (1:6)[2:5] === 2:5 +@test (1:6)[2:2:5] === 2:2:4 +@test (1:2:13)[2:6] === 3:2:11 +@test (1:2:13)[2:3:7] === 3:6:13 @test isempty((1:4)[5:4]) @test_throws BoundsError (1:10)[8:-1:-2] @@ -245,7 +252,7 @@ else @test sum(Int64(1):10^9-1) == div(10^9 * (Int64(10^9)-1), 2) end -# Tricky sums of FloatRange #8272 +# Tricky sums of StepRangeLen #8272 @test sum(10000.:-0.0001:0) == 5.00000005e11 @test sum(0:0.001:1) == 500.5 @test sum(0:0.000001:1) == 500000.5 @@ -338,9 +345,6 @@ for T = (Float32, Float64) @test [linspace(-u,-u,1);] == [-u] @test [linspace(-u,u,2);] == [-u,u] @test [linspace(-u,u,3);] == [-u,0,u] - @test [linspace(-u,u,4);] == [-u,0,0,u] - @test [linspace(-u,u,4);][2] === -z - @test [linspace(-u,u,4);][3] === z @test first(linspace(-u,-u,0)) == -u @test last(linspace(-u,-u,0)) == -u @test first(linspace(u,-u,0)) == u @@ -349,18 +353,15 @@ for T = (Float32, Float64) @test [linspace(u,u,1);] == [u] @test [linspace(u,-u,2);] == [u,-u] @test [linspace(u,-u,3);] == [u,0,-u] - @test [linspace(u,-u,4);] == [u,0,0,-u] - @test [linspace(u,-u,4);][2] === z - @test [linspace(u,-u,4);][3] === -z - v = [linspace(-u,u,12);] + v = linspace(-u,u,12) @test length(v) == 12 - @test issorted(v) && unique(v) == [-u,0,0,u] @test [-3u:u:3u;] == [linspace(-3u,3u,7);] == [-3:3;].*u @test [3u:-u:-3u;] == [linspace(3u,-3u,7);] == [3:-1:-3;].*u end # linspace with very large endpoints for T = (Float32, Float64) + largeint = Int(min(maxintfloat(T), typemax(Int))) a = realmax() for i = 1:5 @test [linspace(a,a,1);] == [a] @@ -373,17 +374,15 @@ for T = (Float32, Float64) @test [linspace(a,-b,0);] == [] @test [linspace(a,-b,2);] == [a,-b] @test [linspace(a,-b,3);] == [a,(a-b)/2,-b] - for c = maxintfloat(T)-3:maxintfloat(T) + for c = largeint-3:largeint s = linspace(-a,b,c) @test first(s) == -a @test last(s) == b - c <= typemax(Int) && @test length(s) == c - @test s.len == c + @test length(s) == c s = linspace(a,-b,c) @test first(s) == a @test last(s) == -b - c <= typemax(Int) && @test length(s) == c - @test s.len == c + @test length(s) == c end b = prevfloat(b) end @@ -456,7 +455,7 @@ let r1 = 1.0:0.1:2.0, r2 = 1.0f0:0.2f0:3.0f0, r3 = 1:2:21 @test r1 + r1 == 2*r1 @test r1 + r2 == 2.0:0.3:5.0 @test (r1 + r2) - r2 == r1 - @test r1 + r3 == convert(FloatRange{Float64}, r3) + r1 + @test r1 + r3 == convert(StepRangeLen{Float64}, r3) + r1 @test r3 + r3 == 2 * r3 end @@ -475,7 +474,7 @@ r = linspace(1/3,5/7,6) @test abs(r[end] - 5/7) <= eps(5/7) r = linspace(0.25,0.25,1) @test length(r) == 1 -@test_throws ErrorException linspace(0.25,0.5,1) +@test_throws ArgumentError linspace(0.25,0.5,1) # issue #7426 @test [typemax(Int):1:typemax(Int);] == [typemax(Int)] @@ -486,12 +485,16 @@ r7484 = 0.1:0.1:1 # issue #7387 for r in (0:1, 0.0:1.0) - @test r+im == [r;]+im - @test r-im == [r;]-im - @test r*im == [r;]*im - @test r/im == [r;]/im + @test [r+im;] == [r;]+im + @test [r-im;] == [r;]-im + @test [r*im;] == [r;]*im + @test [r/im;] == [r;]/im end +# Preservation of high precision upon addition +r = (-0.1:0.1:0.3) + ((-0.3:0.1:0.1) + 1e-12) +@test r[3] == 1e-12 + # issue #7709 @test length(map(identity, 0x01:0x05)) == 5 @test length(map(identity, 0x0001:0x0005)) == 5 @@ -554,20 +557,26 @@ end @test_throws ArgumentError StepRange(1.1,1,5.1) @test promote(0f0:inv(3f0):1f0, 0.:2.:5.) === (0:1/3:1, 0.:2.:5.) -@test convert(FloatRange{Float64}, 0:1/3:1) === 0:1/3:1 -@test convert(FloatRange{Float64}, 0f0:inv(3f0):1f0) === 0:1/3:1 + +@test convert(StepRangeLen{Float64}, 0:1/3:1) === 0:1/3:1 +@test convert(StepRangeLen{Float64}, 0f0:inv(3f0):1f0) === 0:1/3:1 @test promote(0:1/3:1, 0:5) === (0:1/3:1, 0.:1.:5.) -@test convert(FloatRange{Float64}, 0:5) === 0.:1.:5. -@test convert(FloatRange{Float64}, 0:1:5) === 0.:1.:5. -@test convert(FloatRange, 0:5) === 0.:1.:5. -@test convert(FloatRange, 0:1:5) === 0.:1.:5. +@test convert(StepRangeLen{Float64}, 0:5) === 0.:1.:5. +@test convert(StepRangeLen{Float64}, 0:1:5) === 0.:1.:5. +@test convert(StepRangeLen, 0:5) == 0:5 +@test convert(StepRangeLen, 0:1:5) == 0:1:5 + +@test convert(LinSpace{Float64}, 0.0:0.1:0.3) === LinSpace{Float64}(0.0, 0.3, 4) +@test convert(LinSpace, 0.0:0.1:0.3) === LinSpace{Float64}(0.0, 0.3, 4) +@test convert(LinSpace, 0:3) === LinSpace{Int}(0, 3, 4) # Issue #11245 let io = IOBuffer() show(io, linspace(1, 2, 3)) str = String(take!(io)) - @test str == "linspace(1.0,2.0,3)" +# @test str == "linspace(1.0,2.0,3)" + @test str == "1.0:0.5:2.0" end # issue 10950 @@ -585,13 +594,16 @@ end replstrmime(x) = sprint((io,x) -> show(IOContext(io, limit=true), MIME("text/plain"), x), x) @test replstrmime(1:4) == "1:4" @test stringmime("text/plain", 1:4) == "1:4" -@test stringmime("text/plain", linspace(1,5,7)) == "7-element LinSpace{Float64}:\n 1.0,1.66667,2.33333,3.0,3.66667,4.33333,5.0" -@test repr(linspace(1,5,7)) == "linspace(1.0,5.0,7)" +@test stringmime("text/plain", linspace(1,5,7)) == "1.0:0.6666666666666666:5.0" +@test stringmime("text/plain", LinSpace{Float64}(1,5,7)) == "7-element LinSpace{Float64}:\n 1.0,1.66667,2.33333,3.0,3.66667,4.33333,5.0" +@test repr(linspace(1,5,7)) == "1.0:0.6666666666666666:5.0" +@test repr(LinSpace{Float64}(1,5,7)) == "linspace(1.0,5.0,7)" @test replstrmime(0:100.) == "0.0:1.0:100.0" # next is to test a very large range, which should be fast because print_range # only examines spacing of the left and right edges of the range, sufficient # to cover the designated screen size. -@test replstrmime(linspace(0,100, 10000)) == "10000-element LinSpace{Float64}:\n 0.0,0.010001,0.020002,0.030003,0.040004,…,99.95,99.96,99.97,99.98,99.99,100.0" +@test replstrmime(linspace(0,100, 10000)) == "0.0:0.010001000100010001:100.0" +@test replstrmime(LinSpace{Float64}(0,100, 10000)) == "10000-element LinSpace{Float64}:\n 0.0,0.010001,0.020002,0.030003,0.040004,…,99.95,99.96,99.97,99.98,99.99,100.0" @test sprint(io -> show(io,UnitRange(1,2))) == "1:2" @test sprint(io -> show(io,StepRange(1,2,5))) == "1:2:5" @@ -600,15 +612,15 @@ replstrmime(x) = sprint((io,x) -> show(IOContext(io, limit=true), MIME("text/pla # Issue 11049 and related @test promote(linspace(0f0, 1f0, 3), linspace(0., 5., 2)) === (linspace(0., 1., 3), linspace(0., 5., 2)) -@test convert(LinSpace{Float64}, linspace(0., 1., 3)) === linspace(0., 1., 3) -@test convert(LinSpace{Float64}, linspace(0f0, 1f0, 3)) === linspace(0., 1., 3) +@test convert(LinSpace{Float64}, linspace(0., 1., 3)) === LinSpace(0., 1., 3) +@test convert(LinSpace{Float64}, linspace(0f0, 1f0, 3)) === LinSpace(0., 1., 3) @test promote(linspace(0., 1., 3), 0:5) === (linspace(0., 1., 3), linspace(0., 5., 6)) -@test convert(LinSpace{Float64}, 0:5) === linspace(0., 5., 6) -@test convert(LinSpace{Float64}, 0:1:5) === linspace(0., 5., 6) -@test convert(LinSpace, 0:5) === linspace(0., 5., 6) -@test convert(LinSpace, 0:1:5) === linspace(0., 5., 6) +@test convert(LinSpace{Float64}, 0:5) === LinSpace(0., 5., 6) +@test convert(LinSpace{Float64}, 0:1:5) === LinSpace(0., 5., 6) +@test convert(LinSpace, 0:5) === LinSpace{Int}(0, 5, 6) +@test convert(LinSpace, 0:1:5) === LinSpace{Int}(0, 5, 6) function test_range_index(r, s) @test typeof(r[s]) == typeof(r) @@ -620,26 +632,26 @@ test_range_index(linspace(1.0, 1.0, 1), 1:1) test_range_index(linspace(1.0, 1.0, 1), 1:0) test_range_index(linspace(1.0, 2.0, 0), 1:0) -function test_linspace_identity{T}(r::LinSpace{T}, mr::LinSpace{T}) +function test_linspace_identity{T}(r::Range{T}, mr) @test -r == mr @test -collect(r) == collect(mr) - @test isa(-r, LinSpace) + @test isa(-r, typeof(r)) @test 1 + r + (-1) == r @test 1 + collect(r) == collect(1 + r) == collect(r + 1) - @test isa(1 + r + (-1), LinSpace) + @test isa(1 + r + (-1), typeof(r)) @test 1 - r - 1 == mr @test 1 - collect(r) == collect(1 - r) == collect(1 + mr) @test collect(r) - 1 == collect(r - 1) == -collect(mr + 1) - @test isa(1 - r - 1, LinSpace) + @test isa(1 - r - 1, typeof(r)) @test 1 * r * 1 == r @test 2 * r * T(0.5) == r - @test isa(1 * r * 1, LinSpace) + @test isa(1 * r * 1, typeof(r)) @test r / 1 == r @test r / 2 * 2 == r @test r / T(0.5) * T(0.5) == r - @test isa(r / 1, LinSpace) + @test isa(r / 1, typeof(r)) @test (2 * collect(r) == collect(r * 2) == collect(2 * r) == collect(r * T(2.0)) == collect(T(2.0) * r) == @@ -786,3 +798,29 @@ io = IOBuffer() show(io, r) str = String(take!(io)) @test str == "Base.OneTo(3)" + +# linspace of other types +r = linspace(0, 3//10, 4) +@test eltype(r) == Rational{Int} +@test r[2] === 1//10 + +a, b = 1.0, nextfloat(1.0) +ba, bb = BigFloat(a), BigFloat(b) +r = linspace(ba, bb, 3) +@test eltype(r) == BigFloat +@test r[1] == a && r[3] == b +@test r[2] == (ba+bb)/2 + +a, b = rand(10), rand(10) +r = linspace(a, b, 5) +@test r[1] == a && r[5] == b +for i = 2:4 + x = ((5-i)//4)*a + ((i-1)//4)*b + @test r[i] == x +end + +# issue #14420 +for r in (linspace(0.10000000000000045, 1), 0.10000000000000045:(1-0.10000000000000045)/49:1) + @test r[1] === 0.10000000000000045 + @test r[end] === 1.0 +end From a12e587f49582a5d9dfdca52e80ee857c53c71ef Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Thu, 26 Jan 2017 13:55:16 -0600 Subject: [PATCH 2/2] Clarify NEWS on StepRangeLen [ci skip] --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index f6dd9408fe99b..0d5a4604e0b14 100644 --- a/NEWS.md +++ b/NEWS.md @@ -150,7 +150,7 @@ This section lists changes that do not have deprecation warnings. `FloatNN`. `LinSpace(start, stop, len)` always returns a `LinSpace`. - + `StepRangeLen(a, step, len)` constructs a single-precision range + + `StepRangeLen(a, step, len)` constructs an ordinary-precision range using the values and types of `a` and `step` as given, whereas `range(a, step, len)` will attempt to match inputs `a::FloatNN` and `step::FloatNN` to rationals and construct a `StepRangeLen`