From 2b5ec5e43b5dbc6899862e10c83be830fb758eea Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Thu, 26 Mar 2015 14:39:08 +0100 Subject: [PATCH] linspace: "center" computations to handle extreme ranges correctly. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This forces the LinSpace divisor – either n*e or just n – to be in the range [1/2,0), which, in particular, prevents overflow even if the old computation could overflow. Needs more testing of corner cases, but this passes all of the existing tests at least. --- base/range.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/base/range.jl b/base/range.jl index 2224862ac6a5c..a2451ca7c8b0f 100644 --- a/base/range.jl +++ b/base/range.jl @@ -186,20 +186,25 @@ function linspace{T<:FloatingPoint}(start::T, stop::T, len::Int) n = len - 1 a0, b = rat(start) a = convert(T,a0) - if 0 < n && a/convert(T,b) == start + 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) - ne = convert(T,n*e) - if a*n/ne == start && c*n/ne == stop - return LinSpace(a, c, len, ne) + s, p = frexp(convert(T,n*e)) + p = one(p) << p + a /= p; c /= p + if a*n/s == start && c*n/s == stop + return LinSpace(a, c, len, s) end end end - return LinSpace(start, stop, len, max(1,n)) + s, p = frexp(convert(T,n)) + p = one(p) << p + start /= p; stop /= p + return LinSpace(start, stop, len, s) end linspace(start::Real, stop::Real, len::Real=50) = linspace(promote(FloatingPoint(start), FloatingPoint(stop))..., Int(len))