Make LinSpace generic

base/abstractarray.jl

@@ -783,9 +783,7 @@ 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))
 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/abstractarraymath.jl

@@ -100,6 +100,11 @@ imag{T<:Real}(x::AbstractArray{T}) = zero(x)
 
 \(A::Number, B::AbstractArray) = B ./ A
 
+function lerpi(j::Integer, d::Integer, A::AbstractArray, B::AbstractArray)
+ broadcast((a,b) -> lerpi(j, d, a, b), A, B)
+end
+
+
 # index A[:,:,...,i,:,:,...] where "i" is in dimension "d"
 
 """

base/float.jl

@@ -722,9 +722,7 @@ for fn in (:float,:big)
  $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))
+ LinSpace($fn(r.start), $fn(r.stop), length(r))
  end
  end
 end

base/mpfr.jl

@@ -922,4 +922,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

@@ -901,14 +901,8 @@ for f in (:+, :-)
  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},

base/range.jl

@@ -209,68 +209,29 @@ range(a::AbstractFloat, st::Real, len::Integer) = FloatRange(a, float(st), len,
 
 ## linspace and logspace
 
-immutable LinSpace{T<:AbstractFloat} <: Range{T}
+immutable LinSpace{T} <: Range{T}
  start::T
  stop::T
- len::T
- divisor::T
-end
-
-function linspace{T<:AbstractFloat}(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)
- p2 = oftype(s,2)^p
- a /= p2; c /= p2
- end
- if a*n/s == start && c*n/s == stop
- return LinSpace(a, c, len, s)
- end
+ len::Int
+ lendiv::Int
+
+ function LinSpace(start,stop,len)
+ len >= 0 || error("linspace($start, $stop, $len): negative length")
+ if len == 1
+ start == stop || error("linspace($start, $stop, $len): endpoints differ")
+ return new(start, stop, 1, 1)
  end
+ new(start,stop,len,max(len-1,1))
  end
- a, c, s = start, stop, n
- if isinf(a*n) || isinf(c*n)
- s, p = frexp(s)
- p2 = oftype(s,2)^p
- a /= p2; c /= p2
- end
- if a*n/s == start && c*n/s == stop
- return LinSpace(a, c, len, s)
- end
- return LinSpace(start, stop, len, n)
 end
-function linspace{T<:AbstractFloat}(start::T, stop::T, len::Real)
- T_len = convert(T, len)
- T_len == len || throw(InexactError())
- linspace(start, stop, T_len)
+
+function LinSpace(start, stop, len::Integer)
+ T = typeof((stop-start)/len)
+ LinSpace{T}(start, stop, len)
 end
 
 """
- linspace(start::Real, stop::Real, n::Real=50)
+ linspace(start, stop, n=50)
 
 Construct a range of `n` linearly spaced elements from `start` to `stop`.
 
@@ -280,8 +241,7 @@ julia> linspace(1.3,2.9,9)
  1.3,1.5,1.7,1.9,2.1,2.3,2.5,2.7,2.9
 ```
 """
-linspace(start::Real, stop::Real, len::Real=50) =
- linspace(promote(AbstractFloat(start), AbstractFloat(stop))..., len)
+linspace(start, stop, len::Real=50) = LinSpace(start, stop, Int(len))
 
 function show(io::IO, r::LinSpace)
  print(io, "linspace(")
@@ -398,7 +358,7 @@ julia> step(linspace(2.5,10.9,85))
 step(r::StepRange) = r.step
 step(r::AbstractUnitRange) = 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)
+step(r::LinSpace) = (last(r)-first(r))/r.lendiv
 
 unsafe_length(r::Range) = length(r) # generic fallback
 
@@ -412,7 +372,7 @@ unsafe_length(r::OneTo) = r.stop
 length(r::AbstractUnitRange) = unsafe_length(r)
 length(r::OneTo) = unsafe_length(r)
 length(r::FloatRange) = Integer(r.len)
-length(r::LinSpace) = Integer(r.len + signbit(r.len - 1))
+length(r::LinSpace) = r.len
 
 function length{T<:Union{Int,UInt,Int64,UInt64}}(r::StepRange{T})
  isempty(r) && return zero(T)
@@ -450,11 +410,11 @@ end
 first{T}(r::OrdinalRange{T}) = convert(T, r.start)
 first{T}(r::OneTo{T}) = one(T)
 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)
+first(r::LinSpace) = r.start
 
 last{T}(r::OrdinalRange{T}) = convert(T, 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)
+last(r::LinSpace) = r.stop
 
 minimum(r::AbstractUnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : first(r)
 maximum(r::AbstractUnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : last(r)
@@ -477,8 +437,10 @@ next{T}(r::FloatRange{T}, i::Int) =
 
 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)
+function next(r::LinSpace, i::Int)
+ @_inline_meta
+ unsafe_getindex(r, i), i+1
+end
 
 start(r::StepRange) = oftype(r.start + r.step, r.start)
 next{T}(r::StepRange{T}, i) = (convert(T,i), i+r.step)
@@ -536,10 +498,24 @@ function getindex{T}(r::FloatRange{T}, i::Integer)
  convert(T, (r.start + (i-1)*r.step)/r.divisor)
 end
 
-function getindex{T}(r::LinSpace{T}, i::Integer)
+function getindex(r::LinSpace, i::Integer)
  @_inline_meta
  @boundscheck checkbounds(r, i)
- convert(T, ((r.len-i)*r.start + (i-1)*r.stop)/r.divisor)
+ unsafe_getindex(r, i)
+end
+
+# This is separate to make it useful even when running with --check-bounds=yes
+function unsafe_getindex(r::LinSpace, i::Integer)
+ d = r.lendiv
+ j, a, b = ifelse(2i >= length(r), (i-1, r.start, r.stop), (length(r)-i, r.stop, r.start))
+ lerpi(j, d, a, b)
+end
+
+# High-precision interpolation. Accurate for t ∈ [0.5,1], so that 1-t is exact.
+@inline function lerpi{T}(j::Integer, d::Integer, a::T, b::T)
+ t = j/d
+ # computes (1-t)*a + t*b
+ T(fma(t, b, fma(-t, a, a)))
 end
 
 getindex(r::Range, ::Colon) = copy(r)
@@ -581,11 +557,11 @@ end
 function getindex{T}(r::LinSpace{T}, s::OrdinalRange)
  @_inline_meta
  @boundscheck checkbounds(r, s)
- sl::T = length(s)
+ sl = length(s)
  ifirst = first(s)
  ilast = last(s)
- vfirst::T = ((r.len - ifirst) * r.start + (ifirst - 1) * r.stop) / r.divisor
- vlast::T = ((r.len - ilast) * r.start + (ilast - 1) * r.stop) / r.divisor
+ vfirst = unsafe_getindex(r, ifirst)
+ vlast = unsafe_getindex(r, ilast)
  return linspace(vfirst, vlast, sl)
 end
 
@@ -739,43 +715,58 @@ end
 
 -(r::OrdinalRange) = range(-first(r), -step(r), length(r))
 -(r::FloatRange) = FloatRange(-r.start, -r.step, r.len, r.divisor)
--(r::LinSpace)  = LinSpace(-r.start, -r.stop, r.len, r.divisor)
+-(r::LinSpace) = LinSpace(-r.start, -r.stop, r.len)
 
 .+(x::Real, r::AbstractUnitRange) = range(x + first(r), length(r))
 .+(x::Real, r::Range) = (x+first(r)):step(r):(x+last(r))
 #.+(x::Real, r::StepRange) = range(x + r.start, r.step, length(r))
 .+(x::Real, r::FloatRange) = FloatRange(r.divisor*x + r.start, r.step, r.len, r.divisor)
-function .+{T}(x::Real, r::LinSpace{T})
- x2 = x * r.divisor / (r.len - 1)
- LinSpace(x2 + r.start, x2 + r.stop, r.len, r.divisor)
+function .+(x::Real, r::LinSpace)
+ LinSpace(x + r.start, x + r.stop, r.len)
+end
+function .+(x::Number, r::LinSpace)
+ LinSpace(x + r.start, x + r.stop, r.len)
+end
+function .+{T}(x::Ref{T}, r::LinSpace{T})
+ LinSpace(x + r.start, x + r.stop, r.len)
 end
 .+(r::Range, x::Real) = x + r
 #.+(r::FloatRange, x::Real) = x + r
 
 .-(x::Real, r::Range) = (x-first(r)):-step(r):(x-last(r))
 .-(x::Real, r::FloatRange) = FloatRange(r.divisor*x - r.start, -r.step, r.len, r.divisor)
 function .-(x::Real, r::LinSpace)
- x2 = x * r.divisor / (r.len - 1)
- LinSpace(x2 - r.start, x2 - r.stop, r.len, r.divisor)
+ LinSpace(x - r.start, x - r.stop, r.len)
+end
+function .-(x::Number, r::LinSpace)
+ LinSpace(x - r.start, x - r.stop, r.len)
+end
+function .-{T}(x::Ref{T}, r::LinSpace{T})
+ LinSpace(x - r.start, x - r.stop, r.len)
 end
 .-(r::AbstractUnitRange, x::Real) = range(first(r)-x, length(r))
 .-(r::StepRange , x::Real) = range(r.start-x, r.step, length(r))
 .-(r::FloatRange, x::Real) = FloatRange(r.start - r.divisor*x, r.step, r.len, r.divisor)
 function .-(r::LinSpace, x::Real)
- x2 = x * r.divisor / (r.len - 1)
- LinSpace(r.start - x2, r.stop - x2, r.len, r.divisor)
+ LinSpace(r.start - x, r.stop - x, r.len)
+end
+function .-(r::LinSpace, x::Number)
+ LinSpace(r.start - x, r.stop - x, r.len)
+end
+function .-{T}(r::LinSpace{T}, x::Ref{T})
+ LinSpace(r.start - x, r.stop - x, r.len)
 end
 
 .*(x::Real, r::OrdinalRange) = range(x*first(r), x*step(r), length(r))
 .*(x::Real, r::FloatRange) = FloatRange(x*r.start, x*r.step, r.len, r.divisor)
-.*(x::Real, r::LinSpace) = LinSpace(x * r.start, x * r.stop, r.len, r.divisor)
+.*(x::Real, r::LinSpace) = LinSpace(x * r.start, x * r.stop, r.len)
 .*(r::Range, x::Real) = x .* r
 .*(r::FloatRange, x::Real) = x .* r
 .*(r::LinSpace, x::Real) = x .* r
 
 ./(r::OrdinalRange, x::Real) = range(first(r)/x, step(r)/x, length(r))
 ./(r::FloatRange, x::Real) = FloatRange(r.start/x, r.step/x, r.len, r.divisor)
-./(r::LinSpace, x::Real) = LinSpace(r.start / x, r.stop / x, r.len, r.divisor)
+./(r::LinSpace, x::Real) = LinSpace(r.start / x, r.stop / x, r.len)
 
 promote_rule{T1,T2}(::Type{UnitRange{T1}},::Type{UnitRange{T2}}) =
  UnitRange{promote_type(T1,T2)}
@@ -820,20 +811,20 @@ promote_rule{T1,T2}(::Type{LinSpace{T1}},::Type{LinSpace{T2}}) =
  LinSpace{promote_type(T1,T2)}
 convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::LinSpace{T}) = r
 convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::LinSpace) =
- LinSpace{T}(r.start, r.stop, r.len, r.divisor)
+ LinSpace{T}(r.start, r.stop, r.len)
 
 promote_rule{F,OR<:OrdinalRange}(::Type{LinSpace{F}}, ::Type{OR}) =
  LinSpace{promote_type(F,eltype(OR))}
 convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::OrdinalRange) =
- linspace(convert(T, first(r)), convert(T, last(r)), convert(T, length(r)))
+ linspace(convert(T, first(r)), convert(T, last(r)), length(r))
 convert{T}(::Type{LinSpace}, r::OrdinalRange{T}) =
  convert(LinSpace{typeof(float(first(r)))}, r)
 
 # Promote FloatRange to LinSpace
 promote_rule{F,OR<:FloatRange}(::Type{LinSpace{F}}, ::Type{OR}) =
  LinSpace{promote_type(F,eltype(OR))}
 convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::FloatRange) =
- linspace(convert(T, first(r)), convert(T, last(r)), convert(T, length(r)))
+ linspace(convert(T, first(r)), convert(T, last(r)), length(r))
 convert{T<:AbstractFloat}(::Type{LinSpace}, r::FloatRange{T}) =
  convert(LinSpace{T}, r)
 
@@ -878,7 +869,7 @@ collect(r::Range) = vcat(r)
 
 reverse(r::OrdinalRange) = colon(last(r), -step(r), first(r))
 reverse(r::FloatRange) = FloatRange(r.start + (r.len-1)*r.step, -r.step, r.len, r.divisor)
-reverse(r::LinSpace) = LinSpace(r.stop, r.start, r.len, r.divisor)
+reverse(r::LinSpace) = LinSpace(r.stop, r.start, length(r))
 
 ## sorting ##
 

base/rational.jl

@@ -385,3 +385,7 @@ end
 function ^{T<:Rational}(z::Complex{T}, n::Integer)
  n >= 0 ? power_by_squaring(z,n) : power_by_squaring(inv(z),-n)
 end
+
+@inline function lerpi(j::Integer, d::Integer, a::Rational, b::Rational)
+ ((d-j)*a)/d + (j*b)/d
+end