FloatRange: new type for accurate floating-point ranges [fix #2333]

base/range.jl

@@ -35,6 +35,15 @@ immutable Range1{T<:Real} <: Ranges{T}
 end
 Range1{T}(start::T, len::Integer) = Range1{T}(start, len)
 
+immutable FloatRange{T<:FloatingPoint} <: Ranges{T}
+ start::T
+ step::T
+ divisor::T
+ len::T
+end
+FloatRange(a::FloatingPoint, s::FloatingPoint, d::FloatingPoint, l::Real) =
+ FloatRange{promote_type(typeof(a),typeof(s),typeof(d))}(a,s,d,l)
+
 function colon{T<:Integer}(start::T, step::T, stop::T)
  step != 0 || error("step cannot be zero in colon syntax")
  Range{T}(start, step, max(0, 1 + fld(stop-start, step)))
@@ -59,7 +68,7 @@ end
 
 function colon{T<:Real}(start::T, step::T, stop::T)
  step != 0 || error("step cannot be zero in colon syntax")
- if (step<0) != (stop<start)
+ if (step < 0) != (stop < start)
  len = 0
  else
  nf = (stop-start)/step + 1
@@ -104,17 +113,66 @@ end
 colon(start::Real, step::Real, stop::Real) = colon(promote(start, step, stop)...)
 colon(start::Real, stop::Real) = colon(promote(start, stop)...)
 
+# float rationalization helper
+function rat(x)
+ y = x
+ a = d = 1
+ b = c = 0
+ m = typemax(Int) >> 1
+ while max(abs(a),abs(b)) <= m
+ f = itrunc(y)
+ y -= f
+ a, c = f*a + c, a
+ b, d = f*b + d, b
+ (y == 0 || oftype(x,a)/oftype(x,b) == x) && return a, b
+ y = inv(y)
+ end
+ return c, d
+end
+
+# float range "lifting" helper
+function frange{T<:FloatingPoint}(start::T, step::T, stop::T)
+ r = (stop-start)/step
+ n = iround(r)
+ lo = prevfloat((prevfloat(stop)-nextfloat(start))/n)
+ hi = nextfloat((nextfloat(stop)-prevfloat(start))/n)
+ if lo <= step <= hi
+ a, b = rat(start)
+ if convert(T,a)/convert(T,b) == start
+ c, d = rat(step)
+ if convert(T,c)/convert(T,d) == step
+ e = lcm(b,d)
+ a *= div(e,b)
+ c *= div(e,d)
+ if convert(T,a+n*c)/convert(T,e) == stop
+ return convert(T,a), convert(T,c), convert(T,e), convert(T,n+1)
+ end
+ end
+ end
+ end
+ start, step, one(step), floor(r)+1
+end
+
+colon{T<:FloatingPoint}(start::T, step::T, stop::T) =
+ step == 0 ? error("step cannot be zero in colon syntax") :
+ start == stop ? FloatRange{T}(start,step,1,1) :
+ (0 < step) != (start < stop) ? FloatRange{T}(start,step,1,0) :
+ FloatRange{T}(frange(start,step,stop)...)
+
 similar(r::Ranges, T::Type, dims::Dims) = Array(T, dims)
 
-length(r::Ranges) = r.len
-size(r::Ranges) = (r.len,)
+length(r::Ranges) = integer(r.len)
+size(r::Ranges) = (length(r),)
 isempty(r::Ranges) = r.len==0
 first(r::Ranges) = r.start
+first(r::FloatRange) = r.start/r.divisor
 last{T}(r::Range1{T}) = oftype(T, r.start + r.len-1)
 last{T}(r::Range{T}) = oftype(T, r.start + (r.len-1)*r.step)
+last{T}(r::FloatRange{T}) = oftype(T, (r.start + (r.len-1)*r.step)/r.divisor)
 
-step(r::Range)  = r.step
+step(r::Range) = r.step
 step(r::Range1) = one(r.start)
+step(r::FloatRange) = r.step/r.divisor
 
 minimum(r::Range1) = isempty(r) ? error("range must be non-empty") : first(r)
 maximum(r::Range1) = isempty(r) ? error("range must be non-empty") : last(r)
@@ -130,9 +188,13 @@ copy(r::Ranges) = r
 getindex(r::Ranges, i::Real) = getindex(r, to_index(i))
 
 function getindex{T}(r::Ranges{T}, i::Integer)
- if !(1 <= i <= r.len); error(BoundsError); end
+ 1 <= i <= r.len || error(BoundsError)
  oftype(T, r.start + (i-1)*step(r))
 end
+function getindex{T}(r::FloatRange{T}, i::Integer)
+ 1 <= i <= r.len || error(BoundsError)
+ oftype(T, (r.start + (i-1)*r.step)/r.divisor)
+end
 
 function getindex(r::Range1, s::Range1{Int})
  if s.len > 0
@@ -156,18 +218,16 @@ function getindex(r::Ranges, s::Ranges{Int})
  end
 end
 
-function show(io::IO, r::Range)
- if step(r) == 0
- print(io, "Range(",r.start,",",step(r),",",r.len,")")
- else
- print(io, repr(r.start),':',repr(step(r)),':',repr(last(r)))
- end
+function show(io::IO, r::Union(Range,FloatRange))
+ step(r) == 0 ? invoke(show,(IO,Any),io,r) :
+ print(io, repr(first(r)), ':', repr(step(r)), ':', repr(last(r)))
 end
-show(io::IO, r::Range1) = print(io, repr(r.start),':',repr(last(r)))
+show(io::IO, r::Range1) = print(io, repr(first(r)), ':', repr(last(r)))
 
 start(r::Ranges) = 0
-next{T}(r::Range{T},  i) = (oftype(T, r.start + i*step(r)), i+1)
+next{T}(r::Range{T}, i) = (oftype(T, r.start + i*step(r)), i+1)
 next{T}(r::Range1{T}, i) = (oftype(T, r.start + i), i+1)
+next{T}(r::FloatRange{T}, i) = (oftype(T, (r.start + i*r.step)/r.divisor), i+1)
 done(r::Ranges, i) = (length(r) <= i)
 
 # though these look very similar to the above, for some reason LLVM generates
@@ -176,7 +236,7 @@ start{T<:Integer}(r::Range1{T}) = r.start
 next{T<:Integer}(r::Range1{T}, i) = (i, oftype(T, i+1))
 done{T<:Integer}(r::Range1{T}, i) = i==oftype(T, r.start+r.len)
 
-==(r::Ranges, s::Ranges) = (r.start==s.start) & (step(r)==step(s)) & (r.len==s.len)
+==(r::Ranges, s::Ranges) = (first(r)==first(s)) & (step(r)==step(s)) & (length(r)==length(s))
 ==(r::Range1, s::Range1) = (r.start==s.start) & (r.len==s.len)
 
 # TODO: isless?
@@ -389,7 +449,8 @@ function vcat{T}(rs::Ranges{T}...)
  return a
 end
 
-reverse{T<:Real}(r::Ranges{T}) = Range(last(r), -step(r), r.len)
+reverse(r::Ranges) = Range(last(r), -step(r), r.len)
+reverse(r::FloatRange) = FloatRange(last(r), -r.step, r.divisor, r.len)
 
 ## sorting ##
 

base/serialize.jl

@@ -330,9 +330,6 @@ function deserialize(s, ::Type{Module})
  path = deserialize(s)
  m = Main
  for mname in path
- if !isdefined(m,mname)
- warn("Module $mname not defined on process $(myid())") # an error seemingly fails
- end
  m = eval(m,mname)::Module
  end
  m

test/random.jl

@@ -6,12 +6,12 @@
 @test minimum([rand(int32(1):int32(7^7)) for i = 1:100000]) > 0
 @test(typeof(rand(false:true)) == Bool)
 
-for T in (Int8, Uint8, Int16, Uint16, Int32, Uint32, Int64, Uint64, Int128, Uint128, Char, BigInt,
- Float16, Float32, Float64, Rational{Int})
- r = rand(convert(T, 97):convert(T, 122))
+for T in (Int8, Uint8, Int16, Uint16, Int32, Uint32, Int64, Uint64, Int128, Uint128,
+ Char, BigInt, Float32, Float64, Rational{Int})
+ r = rand(convert(T,97):convert(T,122))
  @test typeof(r) == T
  @test 97 <= r <= 122
- r = rand(convert(T, 97):convert(T,2):convert(T, 122),2)[1]
+ r = rand(convert(T,97):convert(T,2):convert(T,122),2)[1]
  @test typeof(r) == T
  @test 97 <= r <= 122
  @test mod(r,2)==1

test/ranges.jl

@@ -159,3 +159,48 @@ else
  @test sum(int64(1:10^9)) == div(10^9 * (int64(10^9)+1), 2)
  @test sum(int64(1:10^9-1)) == div(10^9 * (int64(10^9)-1), 2)
 end
+
+# tricky floating-point ranges
+
+@test 0.1:0.1:0.3 == [1:3]./10
+@test 0.0:0.1:0.3 == [0:3]./10
+@test 0.3:-0.1:-0.1 == [3:-1:-1]./10
+@test 0.1:-0.1:-0.3 == [1:-1:-3]./10
+@test 0.0:0.1:1.0 == [0:10]./10
+@test 0.0:-0.1:1.0 == []
+@test 0.0:0.1:-1.0 == []
+@test 0.0:-0.1:-1.0 == [0:-1:-10]./10
+@test 1.0:1/49:27.0 == [49:1323]./49
+@test 0.0:0.7:2.1 == [0:7:21]./10
+@test 0.0:1.1:3.3 == [0:11:33]./10
+@test 0.1:1.1:3.4 == [1:11:34]./10
+@test 0.0:1.3:3.9 == [0:13:39]./10
+@test 0.1:1.3:4.0 == [1:13:40]./10
+@test 1.1:1.1:3.3 == [11:11:33]./10
+@test 0.3:0.1:1.1 == [3:1:11]./10
+
+@test 0.0:1.0:5.5 == [0:10:55]./10
+@test 0.0:-1.0:0.5 == []
+@test 0.0:1.0:0.5 == [0.0]
+
+@test prevfloat(0.1):0.1:0.3 == [prevfloat(0.1), 0.2, 0.3]
+@test nextfloat(0.1):0.1:0.3 == [nextfloat(0.1), 0.2]
+@test prevfloat(0.0):0.1:0.3 == [prevfloat(0.0), 0.1, 0.2]
+@test nextfloat(0.0):0.1:0.3 == [nextfloat(0.0), 0.1, 0.2]
+@test 0.1:0.1:prevfloat(0.3) == [0.1, 0.2]
+@test 0.1:0.1:nextfloat(0.3) == [0.1, 0.2, nextfloat(0.3)]
+@test 0.0:0.1:prevfloat(0.3) == [0.0, 0.1, 0.2]
+@test 0.0:0.1:nextfloat(0.3) == [0.0, 0.1, 0.2, nextfloat(0.3)]
+@test 0.1:prevfloat(0.1):0.3 == [0.1, 0.2, 0.3]
+@test 0.1:nextfloat(0.1):0.3 == [0.1, 0.2]
+@test 0.0:prevfloat(0.1):0.3 == [0.0, prevfloat(0.1), prevfloat(0.2), 0.3]
+@test 0.0:nextfloat(0.1):0.3 == [0.0, nextfloat(0.1), nextfloat(0.2)]
+
+for T = (Float32, Float64,),# BigFloat),
+ a = -5:25, s = [-5:-1;1:25], d = 1:25, n = -1:15
+ den = convert(T,d)
+ start = convert(T,a)/den
+ step = convert(T,s)/den
+ stop = convert(T,(a+(n-1)*s))/den
+ @test [start:step:stop] == T[a:s:a+(n-1)*s]./den
+end