Merge pull request #18777 from JuliaLang/teh/linspace

Make LinSpace generic

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 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`
+ 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

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 ##

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")

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()

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

base/exports.jl

@@ -62,7 +62,7 @@ export
  ExponentialBackOff,
  Factorization,
  FileMonitor,
- FloatRange,
+ StepRangeLen,
  Future,
  Hermitian,
  UniformScaling,

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,26 +786,19 @@ 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
 # TODO: do the definitions below primarily pertaining to integers belong in float.jl?
 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

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

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}