Ranges: use Float64 rather than TwicePrecision for Float16, Float32

Author: timholy

base/float.jl

@@ -885,7 +885,8 @@ 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)
+float(r::StepRangeLen{T}) where {T} =
+    StepRangeLen{typeof(float(T(r.ref)))}(float(r.ref), float(r.step), length(r), r.offset)
 function float(r::LinSpace)
     LinSpace(float(r.start), float(r.stop), length(r))
 end

base/range.jl

@@ -205,6 +205,8 @@ end
 
 StepRangeLen(ref::R, step::S, len::Integer, offset::Integer = 1) where {R,S} =
     StepRangeLen{typeof(ref+0*step),R,S}(ref, step, len, offset)
+StepRangeLen{T}(ref::R, step::S, len::Integer, offset::Integer = 1) where {T,R,S} =
+    StepRangeLen{T,R,S}(ref, step, len, offset)
 
 ## linspace and logspace
 
@@ -370,9 +372,12 @@ julia> step(linspace(2.5,10.9,85))
 """
 step(r::StepRange) = r.step
 step(r::AbstractUnitRange) = 1
-step(r::StepRangeLen) = r.step
+step(r::StepRangeLen{T}) where {T} = T(r.step)
 step(r::LinSpace) = (last(r)-first(r))/r.lendiv
 
+step_hp(r::StepRangeLen) = r.step
+step_hp(r::Range) = step(r)
+
 unsafe_length(r::Range) = length(r)  # generic fallback
 
 function unsafe_length(r::StepRange)
@@ -455,10 +460,9 @@ done(r::StepRange, i) = isempty(r) | (i < min(r.start, r.stop)) | (i > max(r.sta
 done(r::StepRange, i::Integer) =
     isempty(r) | (i == oftype(i, r.stop) + r.step)
 
-# see also twiceprecision.jl
-start(r::StepRangeLen) = (unsafe_getindex(r, 1), 1)
-next(r::StepRangeLen{T}, s) where {T} = s[1], (T(s[1]+r.step), s[2]+1)
-done(r::StepRangeLen, s) = s[2] > length(r)
+start(r::StepRangeLen) = 1
+next(r::StepRangeLen{T}, i) where {T} = unsafe_getindex(r, i), i+1
+done(r::StepRangeLen, i) = i > length(r)
 
 start(r::UnitRange{T}) where {T} = oftype(r.start + oneunit(T), r.start)
 next(r::AbstractUnitRange{T}, i) where {T} = (convert(T, i), i + oneunit(T))
@@ -496,7 +500,7 @@ end
 
 function getindex(v::Range{T}, i::Integer) where T
     @_inline_meta
-    ret = convert(T, first(v) + (i - 1)*step(v))
+    ret = convert(T, first(v) + (i - 1)*step_hp(v))
     ok = ifelse(step(v) > zero(step(v)),
                 (ret <= v.stop) & (ret >= v.start),
                 (ret <= v.start) & (ret >= v.stop))
@@ -516,6 +520,11 @@ function unsafe_getindex(r::StepRangeLen{T}, i::Integer) where T
     T(r.ref + u*r.step)
 end
 
+function _getindex_hiprec(r::StepRangeLen, i::Integer)  # without rounding by T
+    u = i - r.offset
+    r.ref + u*r.step
+end
+
 function unsafe_getindex(r::LinSpace, i::Integer)
     lerpi.(i-1, r.lendiv, r.start, r.stop)
 end
@@ -556,11 +565,14 @@ function getindex(r::StepRange, s::Range{<:Integer})
     range(st, step(r)*step(s), length(s))
 end
 
-function getindex(r::StepRangeLen, s::OrdinalRange{<:Integer})
+function getindex(r::StepRangeLen{T}, s::OrdinalRange{<:Integer}) where {T}
     @_inline_meta
     @boundscheck checkbounds(r, s)
-    vfirst = unsafe_getindex(r, first(s))
-    return StepRangeLen(vfirst, r.step*step(s), length(s))
+    # Find closest approach to offset by s
+    ind = linearindices(s)
+    offset = max(min(1 + round(Int, (r.offset - first(s))/step(s)), last(ind)), first(ind))
+    ref = _getindex_hiprec(r, first(s) + (offset-1)*step(s))
+    return StepRangeLen{T}(ref, r.step*step(s), length(s), offset)
 end
 
 function getindex(r::LinSpace, s::OrdinalRange{<:Integer})
@@ -725,16 +737,17 @@ end
 ## linear operations on ranges ##
 
 -(r::OrdinalRange) = range(-first(r), -step(r), length(r))
--(r::StepRangeLen) = StepRangeLen(-r.ref, -r.step, length(r), r.offset)
+-(r::StepRangeLen{T,R,S}) where {T,R,S} =
+    StepRangeLen{T,R,S}(-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))
 # 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)
+function +(x::Number, r::StepRangeLen{T}) where T
     newref = x + r.ref
-    StepRangeLen{eltype(newref),typeof(newref),typeof(r.step)}(newref, r.step, length(r), r.offset)
+    StepRangeLen{typeof(T(r.ref) + x)}(newref, r.step, length(r), r.offset)
 end
 function +(x::Number, r::LinSpace)
     LinSpace(x + r.start, x + r.stop, r.len)
@@ -750,15 +763,18 @@ 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::StepRangeLen{T}) where {T} =
+    StepRangeLen{typeof(x*T(r.ref))}(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::StepRangeLen{T}, x::Number) where {T} =
+    StepRangeLen{typeof(T(r.ref)*x)}(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::StepRangeLen{T}, x::Number) where {T} =
+    StepRangeLen{typeof(T(r.ref)/x)}(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 ]

base/twiceprecision.jl

@@ -308,23 +308,57 @@ end
 
 ## StepRangeLen
 
-# If using TwicePrecision numbers, deliberately force user to specify offset
-StepRangeLen(ref::TwicePrecision{T}, step::TwicePrecision{T}, len::Integer, offset::Integer) where {T} =
+# Use TwicePrecision only for Float64; use Float64 for T<:Union{Float16,Float32}
+# Ratio-of-integers constructors
+function steprangelen_hp(::Type{Float64}, ref::Tuple{Integer,Integer},
+                         step::Tuple{Integer,Integer}, nb::Integer,
+                         len::Integer, offset::Integer)
+    StepRangeLen(TwicePrecision{Float64}(ref),
+                 TwicePrecision{Float64}(step, nb), Int(len), offset)
+end
+
+function steprangelen_hp(::Type{T}, ref::Tuple{Integer,Integer},
+                         step::Tuple{Integer,Integer}, nb::Integer,
+                         len::Integer, offset::Integer) where {T<:IEEEFloat}
+    StepRangeLen{T}(ref[1]/ref[2], step[1]/step[2], Int(len), offset)
+end
+
+# AbstractFloat constructors (can supply a single number or a 2-tuple
+const F_or_FF = Union{AbstractFloat, Tuple{AbstractFloat,AbstractFloat}}
+asF64(x::AbstractFloat) = Float64(x)
+asF64(x::Tuple{AbstractFloat,AbstractFloat}) = Float64(x[1]) + Float64(x[2])
+
+function steprangelen_hp(::Type{Float64}, ref::F_or_FF,
+                         step::F_or_FF, nb::Integer,
+                         len::Integer, offset::Integer)
+    StepRangeLen(TwicePrecision{Float64}(ref...),
+                 twiceprecision(TwicePrecision{Float64}(step...), nb), Int(len), offset)
+end
+
+function steprangelen_hp(::Type{T}, ref::F_or_FF,
+                         step::F_or_FF, nb::Integer,
+                         len::Integer, offset::Integer) where {T<:IEEEFloat}
+    StepRangeLen{T}(asF64(ref),
+                    asF64(step), Int(len), offset)
+end
+
+
+
+StepRangeLen(ref::TwicePrecision{T}, step::TwicePrecision{T},
+             len::Integer, offset::Integer=1) where {T} =
     StepRangeLen{T,TwicePrecision{T},TwicePrecision{T}}(ref, step, len, offset)
 
 # Construct range for rational start=start_n/den, step=step_n/den
 function floatrange(::Type{T}, start_n::Integer, step_n::Integer, len::Integer, den::Integer) where T
     if len < 2
-        return StepRangeLen(TwicePrecision{T}((start_n, den)),
-                            TwicePrecision{T}((step_n, den)), Int(len), 1)
+        return steprangelen_hp(T, (start_n, den), (step_n, den), 0, 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)
+    steprangelen_hp(T, (ref_n, den), (step_n, den), nb, Int(len), imin)
 end
 
 function floatrange(a::AbstractFloat, st::AbstractFloat, len::Real, divisor::AbstractFloat)
@@ -339,8 +373,7 @@ function floatrange(a::AbstractFloat, st::AbstractFloat, len::Real, divisor::Abs
     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)
+    steprangelen_hp(T, (a,divisor), (st,divisor), nbitslen(T, len, 1), Int(len), 1)
 end
 
 function colon(start::T, step::T, stop::T) where T<:Union{Float16,Float32,Float64}
@@ -381,8 +414,7 @@ function colon(start::T, step::T, stop::T) where T<:Union{Float16,Float32,Float6
         # if we've overshot the end, subtract one:
         len -= (start < stop < stop′) + (start > stop > stop′)
     end
-
-    StepRangeLen(TwicePrecision(start, zero(T)), twiceprecision(step, nbitslen(T, len, 1)), len)
+    steprangelen_hp(T, start, step, 0, len, 1)
 end
 
 function range(a::T, st::T, len::Integer) where T<:Union{Float16,Float32,Float64}
@@ -399,16 +431,7 @@ function range(a::T, st::T, len::Integer) where T<:Union{Float16,Float32,Float64
             return floatrange(T, start_n, step_n, len, den)
         end
     end
-    StepRangeLen(TwicePrecision(a, zero(T)), TwicePrecision(st, zero(T)), len)
-end
-
-step(r::StepRangeLen{T,R,S}) where {T,R,S<:TwicePrecision} = convert(eltype(S), r.step)
-
-start(r::StepRangeLen{<:Any,<:TwicePrecision,<:TwicePrecision}) = 1
-done(r::StepRangeLen{<:Any,<:TwicePrecision,<:TwicePrecision}, i::Int) = length(r) < i
-function next(r::StepRangeLen{<:Any,<:TwicePrecision,<:TwicePrecision}, i::Int)
-    @_inline_meta
-    unsafe_getindex(r, i), i+1
+    steprangelen_hp(T, a, st, 0, len, 1)
 end
 
 # This assumes that r.step has already been split so that (0:len-1)*r.step.hi is exact
@@ -542,7 +565,7 @@ end
 function linspace(start::T, stop::T, len::Integer) where {T<:IEEEFloat}
     len < 2 && return _linspace1(T, start, stop, len)
     if start == stop
-        return StepRangeLen(TwicePrecision(start,zero(T)), TwicePrecision(zero(T),zero(T)), len)
+        return steprangelen_hp(T, start, zero(T), 0, len, 1)
     end
     # Attempt to find exact rational approximations
     start_n, start_d = rat(start)
@@ -587,8 +610,7 @@ function _linspace(start::T, stop::T, len::Integer) where {T<:IEEEFloat}
     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)
+        return steprangelen_hp(T, start, (-start, stop), 0, 2, 1)
     end
     # 2x calculations to get high precision endpoint matching while also
     # preventing overflow in ref_hi+(i-offset)*step_hi
@@ -601,23 +623,22 @@ function _linspace(start::T, stop::T, len::Integer) where {T<:IEEEFloat}
     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)
+    steprangelen_hp(T, (ref, ref_lo), (step_hi, step_lo), 0, Int(len), imin)
 end
 
 # linspace for rational numbers, start = start_n/den, stop = stop_n/den
 # Note this returns a StepRangeLen
 function linspace(::Type{T}, start_n::Integer, stop_n::Integer, len::Integer, den::Integer) where T
     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)
+    start_n == stop_n && return steprangelen_hp(T, (start_n, den), (zero(start_n), den), 0, len)
     tmin = -start_n/(Float64(stop_n) - Float64(start_n))
     imin = round(Int, tmin*(len-1)+1)
     imin = clamp(imin, 1, Int(len))
     ref_num = Int128(len-imin) * start_n + Int128(imin-1) * stop_n
     ref_denom = Int128(len-1) * den
-    ref = TwicePrecision{T}((ref_num, ref_denom))
-    step_full = TwicePrecision{T}((Int128(stop_n) - Int128(start_n), ref_denom))
-    step = twiceprecision(step_full, nbitslen(T, len, imin))
-    StepRangeLen(ref, step, Int(len), imin)
+    ref = (ref_num, ref_denom)
+    step_full = (Int128(stop_n) - Int128(start_n), ref_denom)
+    steprangelen_hp(T, ref, step_full,  nbitslen(T, len, imin), Int(len), imin)
 end
 
 # For len < 2

test/ranges.jl

@@ -521,13 +521,13 @@ end
 for T = (Float32, Float64,),# BigFloat),
     a = -5:25, s = [-5:-1;1:25;], d = 1:25, n = -1:15
     denom = convert(T,d)
-    start = convert(T,a)/denom
-    step  = convert(T,s)/denom
+    strt = convert(T,a)/denom
+    Δ     = convert(T,s)/denom
     stop  = convert(T,(a+(n-1)*s))/denom
     vals  = T[a:s:a+(n-1)*s;]./denom
-    r = start:step:stop
+    r = strt:Δ:stop
     @test [r;] == vals
-    @test [linspace(start, stop, length(r));] == vals
+    @test [linspace(strt, stop, length(r));] == vals
     # issue #7420
     n = length(r)
     @test [r[1:n];] == [r;]
@@ -546,24 +546,24 @@ end
 
 function range_fuzztests(::Type{T}, niter, nrange) where {T}
     for i = 1:niter, n in nrange
-        start, Δ = randn(T), randn(T)
+        strt, Δ = randn(T), randn(T)
         Δ == 0 && continue
-        stop = start + (n-1)*Δ
-        # `n` is not necessarily unique s.t. `start + (n-1)*Δ == stop`
-        # so test that `length(start:Δ:stop)` satisfies this identity
-        # and is the closest value to `(stop-start)/Δ` to do so
+        stop = strt + (n-1)*Δ
+        # `n` is not necessarily unique s.t. `strt + (n-1)*Δ == stop`
+        # so test that `length(strt:Δ:stop)` satisfies this identity
+        # and is the closest value to `(stop-strt)/Δ` to do so
         lo = hi = n
-        while start + (lo-1)*Δ == stop; lo -= 1; end
-        while start + (hi-1)*Δ == stop; hi += 1; end
-        m = clamp(round(Int, (stop-start)/Δ) + 1, lo+1, hi-1)
-        r = start:Δ:stop
+        while strt + (lo-1)*Δ == stop; lo -= 1; end
+        while strt + (hi-1)*Δ == stop; hi += 1; end
+        m = clamp(round(Int, (stop-strt)/Δ) + 1, lo+1, hi-1)
+        r = strt:Δ:stop
         @test m == length(r)
-        @test start == first(r)
+        @test strt == first(r)
         @test Δ == step(r)
         @test_skip stop == last(r)
-        l = linspace(start,stop,n)
+        l = linspace(strt,stop,n)
         @test n == length(l)
-        @test start == first(l)
+        @test strt == first(l)
         @test stop  == last(l)
     end
 end