Ranges: use Float64 rather than TwicePrecision for Float16, Float32

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
@@ -457,11 +480,15 @@ convert(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen{T,R,S}) where {T<:AbstractF
 convert(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen) where {T<:AbstractFloat,R<:TwicePrecision,S<:TwicePrecision} =
  _convertSRL(StepRangeLen{T,R,S}, r)
 
+convert(::Type{StepRangeLen{Float64}}, r::StepRangeLen) =
+ _convertSRL(StepRangeLen{Float64,TwicePrecision{Float64},TwicePrecision{Float64}}, r)
 convert(::Type{StepRangeLen{T}}, r::StepRangeLen) where {T<:IEEEFloat} =
- _convertSRL(StepRangeLen{T,TwicePrecision{T},TwicePrecision{T}}, r)
+ _convertSRL(StepRangeLen{T,Float64,Float64}, r)
 
+convert(::Type{StepRangeLen{Float64}}, r::Range) =
+ _convertSRL(StepRangeLen{Float64,TwicePrecision{Float64},TwicePrecision{Float64}}, r)
 convert(::Type{StepRangeLen{T}}, r::Range) where {T<:IEEEFloat} =
- _convertSRL(StepRangeLen{T,TwicePrecision{T},TwicePrecision{T}}, r)
+ _convertSRL(StepRangeLen{T,Float64,Float64}, r)
 
 function _convertSRL(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen{<:Integer}) where {T,R,S}
  StepRangeLen{T,R,S}(R(r.ref), S(r.step), length(r), r.offset)
@@ -542,7 +569,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 +614,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 +627,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
@@ -1078,6 +1078,7 @@ for i = 2:4
  @test r[i] == x
 end
 
+# issue #23178
 r = linspace(Float16(0.1094), Float16(0.9697), 300)
 @test r[1] == Float16(0.1094)
 @test r[end] == Float16(0.9697)
@@ -1148,7 +1149,7 @@ end
  end
 end
 
-# Issue #23300
+# issue #23300
 x = -5:big(1.0):5
 @test map(Float64, x) === -5.0:1.0:5.0
 @test map(Float32, x) === -5.0f0:1.0f0:5.0f0