Merge pull request #18777 from JuliaLang/teh/linspace

Make LinSpace generic

Author: timholy

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}