Merge pull request #9666 from JuliaLang/sk/linspace

linspace: try to "lift" linspace the float ranges are [close #9637]

Author: StefanKarpinski

NEWS.md

@@ -245,6 +245,8 @@ Library improvements
     * `code_llvm` now outputs stripped IR without debug info or other attached metadata.
       Use `code_llvm_raw` for the unstripped output ([#10747]).
 
+    * `linspace` now returns a `LinSpace` object which lazily computes linear interpolation of values between the start and stop values. It "lifts" endpoints which are approximately rational in the same manner as the `colon` operator.
+
 Deprecated or removed
 ---------------------
 
@@ -309,6 +311,8 @@ Deprecated or removed
 
   * the --int-literals compiler option is no longer accepted.
 
+  * Instead of `linrange` use `linspace`.
+
 Julia v0.3.0 Release Notes
 ==========================
 

base/array.jl

@@ -246,28 +246,6 @@ function one{T}(x::AbstractMatrix{T})
     eye(T, m)
 end
 
-linspace(start::Integer, stop::Integer, n::Integer) =
-    linspace(float(start), float(stop), n)
-function linspace(start::Real, stop::Real, n::Integer)
-    (start, stop) = promote(start, stop)
-    T = typeof(start)
-    a = Array(T, Int(n))
-    if n == 1
-        a[1] = start
-        return a
-    end
-    n -= 1
-    S = promote_type(T, Float64)
-    for i=0:n
-        a[i+1] = start*(convert(S, (n-i))/n) + stop*(convert(S, i)/n)
-    end
-    a
-end
-linspace(start::Real, stop::Real) = linspace(start, stop, 100)
-
-logspace(start::Real, stop::Real, n::Integer) = 10.^linspace(start, stop, n)
-logspace(start::Real, stop::Real) = logspace(start, stop, 50)
-
 ## Conversions ##
 
 convert{T,n}(::Type{Array{T}}, x::Array{T,n}) = x

base/deprecated.jl

@@ -522,3 +522,5 @@ export float32_isvalid, float64_isvalid
 @deprecate parseint(s,base)           parse(Int, s, base)
 @deprecate parseint(T::Type, s)       parse(T, s)
 @deprecate parseint(T::Type, s, base) parse(T, s, base)
+
+@deprecate linrange linspace

base/exports.jl

@@ -64,6 +64,7 @@ export
     IO,
     IOBuffer,
     IOStream,
+    LinSpace,
     LocalProcess,
     LowerTriangular,
     MathConst,
@@ -521,7 +522,6 @@ export
     issorted,
     last,
     levicivita,
-    linrange,
     linspace,
     logspace,
     mapslices,

base/range.jl

@@ -163,10 +163,81 @@ range(a::FloatingPoint, st::FloatingPoint, len::Integer) = FloatRange(a,st,len,o
 range(a::Real, st::FloatingPoint, len::Integer) = FloatRange(float(a), st, len, one(st))
 range(a::FloatingPoint, st::Real, len::Integer) = FloatRange(a, float(st), len, one(a))
 
-linrange(a::Real, b::Real, len::Integer) =
-    len >= 2           ? range(a, (b-a)/(len-1), len) :
-    len == 1 && a == b ? range(a, zero((b-a)/(len-1)), 1) :
-                         throw(ArgumentError("invalid range length"))
+## linspace and logspace
+
+immutable LinSpace{T<:FloatingPoint} <: Range{T}
+    start::T
+    stop::T
+    len::T
+    divisor::T
+end
+
+function linspace{T<:FloatingPoint}(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)
+        end
+        return LinSpace(-start, -stop, -one(T), n)
+    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)
+                p = oftype(s,2)^p
+                a /= p; c /= p
+            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)
+        p = oftype(s,2)^p
+        a /= p; c /= p
+    end
+    if a*n/s == start && c*n/s == stop
+        return LinSpace(a, c, len, s)
+    end
+    error("linspace($start, $stop, $len): cannot be constructed")
+end
+function linspace{T<:FloatingPoint}(start::T, stop::T, len::Real)
+    T_len = convert(T, len)
+    T_len == len || throw(InexactError())
+    linspace(start, stop, T_len)
+end
+linspace(start::Real, stop::Real, len::Real=50) =
+    linspace(promote(FloatingPoint(start), FloatingPoint(stop))..., len)
+
+function show(io::IO, r::LinSpace)
+    print(io, "linspace(")
+    show(io, first(r))
+    print(io, ',')
+    show(last(r))
+    print(io, ',')
+    show(length(r))
+    print(io, ')')
+end
+
+logspace(start::Real, stop::Real, n::Integer=50) = 10.^linspace(start, stop, n)
 
 ## interface implementations
 
@@ -178,18 +249,21 @@ size(r::Range) = (length(r),)
 isempty(r::StepRange) =
     (r.start != r.stop) & ((r.step > zero(r.step)) != (r.stop > r.start))
 isempty(r::UnitRange) = r.start > r.stop
-isempty(r::FloatRange) = length(r)==0
+isempty(r::FloatRange) = length(r) == 0
+isempty(r::LinSpace) = length(r) == 0
 
 step(r::StepRange) = r.step
 step(r::UnitRange) = 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)
 
 function length(r::StepRange)
     n = Integer(div(r.stop+r.step - r.start, r.step))
     isempty(r) ? zero(n) : n
 end
 length(r::UnitRange) = Integer(r.stop - r.start + 1)
 length(r::FloatRange) = Integer(r.len)
+length(r::LinSpace) = Integer(r.len + signbit(r.len - 1))
 
 function length{T<:Union(Int,UInt,Int64,UInt64)}(r::StepRange{T})
     isempty(r) && return zero(T)
@@ -220,11 +294,13 @@ let smallint = (Int === Int64 ?
 end
 
 first{T}(r::OrdinalRange{T}) = convert(T, r.start)
-first(r::FloatRange) = r.start/r.divisor
+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)
 
 last{T}(r::StepRange{T}) = r.stop
 last(r::UnitRange) = 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)
 
 minimum(r::UnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : first(r)
 maximum(r::UnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : last(r)
@@ -241,8 +317,14 @@ copy(r::Range) = r
 ## iteration
 
 start(r::FloatRange) = 0
-next{T}(r::FloatRange{T}, i::Int) = (convert(T, (r.start + i*r.step)/r.divisor), i+1)
-done(r::FloatRange, i::Int) = (length(r) <= i)
+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)
 
 # NOTE: For ordinal ranges, we assume start+step might be from a
 # lifted domain (e.g. Int8+Int8 => Int); use that for iterating.
@@ -268,6 +350,10 @@ function getindex{T}(r::FloatRange{T}, i::Integer)
     1 <= i <= length(r) || throw(BoundsError())
     convert(T, (r.start + (i-1)*r.step)/r.divisor)
 end
+function getindex{T}(r::LinSpace{T}, i::Integer)
+    1 <= i <= length(r) || throw(BoundsError())
+    convert(T, ((r.len-i)*r.start + (i-1)*r.stop)/r.divisor)
+end
 
 function check_indexingrange(s, r)
     sl = length(s)

doc/manual/arrays.rst

@@ -89,7 +89,7 @@ Function                                            Description
                                                     distributed random values
 :func:`eye(n) <eye>`                                ``n``-by-``n`` identity matrix
 :func:`eye(m, n) <eye>`                             ``m``-by-``n`` identity matrix
-:func:`linspace(start, stop, n) <linspace>`         vector of ``n`` linearly-spaced elements from ``start`` to ``stop``
+:func:`linspace(start, stop, n) <linspace>`         range of ``n`` linearly spaced elements from ``start`` to ``stop``
 :func:`fill!(A, x) <fill!>`                         fill the array ``A`` with the value ``x``
 :func:`fill(x, dims) <fill>`                        create an array filled with the value ``x``
 =================================================== =====================================================================

doc/manual/noteworthy-differences.rst

@@ -44,8 +44,7 @@ some noteworthy differences that may trip up Julia users accustomed to MATLAB:
   the syntax ``[a b; c d]`` is used to avoid confusion. In Julia v0.4, the
   concatenation syntax ``[x, [y, z]]`` is deprecated in favor of ``[x; [y, z]]``.
 - In Julia, ``a:b`` and ``a:b:c`` construct :obj:`Range` objects. To construct
-  a full vector like in MATLAB, use :func:`collect(a:b) <collect>` or
-  :func:`linspace`.
+  a full vector like in MATLAB, use :func:`collect(a:b) <collect>`.
 - Functions in Julia return values from their last expression or the ``return``
   keyword instead of listing the names of variables to return in the function
   definition (see :ref:`man-return-keyword` for details).

doc/stdlib/arrays.rst

@@ -169,12 +169,11 @@ Constructors
 
 .. function:: linspace(start, stop, n=100)
 
-   Construct a vector of ``n`` linearly-spaced elements from ``start`` to ``stop``.
-   See also: :func:`linrange` that constructs a range object.
+   Construct a range of ``n`` linearly spaced elements from ``start`` to ``stop``.
 
 .. function:: logspace(start, stop, n=50)
 
-   Construct a vector of ``n`` logarithmically-spaced numbers from ``10^start`` to ``10^stop``.
+   Construct a vector of ``n`` logarithmically spaced numbers from ``10^start`` to ``10^stop``.
 
 Mathematical operators and functions
 ------------------------------------

doc/stdlib/math.rst

@@ -189,10 +189,6 @@ Mathematical Operators
 
    Construct a range by length, given a starting value and optional step (defaults to 1).
 
-.. function:: linrange(start, end, length)
-
-   Construct a range by length, given a starting and ending value.
-
 .. _==:
 .. function:: ==(x, y)