Merge pull request #9344 from simonbyrne/roundeven

Changes behaviour of round to use default rounding mode

Author: JeffBezanson

NEWS.md

@@ -54,6 +54,9 @@ Language changes
       * `Char` + `Int` = `Char`
       * `Char` - `Char` = `Int`
 
+  * `round` rounds to the nearest integer using the default rounding mode,
+    which is ties to even by default ([#8750]).
+
 Compiler improvements
 ---------------------
 

base/exports.jl

@@ -85,6 +85,8 @@ export
     RoundDown,
     RoundingMode,
     RoundNearest,
+    RoundNearestTiesAway,
+    RoundNearestTiesUp,
     RoundToZero,
     RoundUp,
     Set,

base/float.jl

@@ -166,13 +166,17 @@ floor{T<:Integer}(::Type{T}, x::FloatingPoint) = trunc(T,floor(x))
 ceil {T<:Integer}(::Type{T}, x::FloatingPoint) = trunc(T,ceil(x))
 round {T<:Integer}(::Type{T}, x::FloatingPoint) = trunc(T,round(x))
 
+trunc(x::Float64) = box(Float64,trunc_llvm(unbox(Float64,x)))
+trunc(x::Float32) = box(Float32,trunc_llvm(unbox(Float32,x)))
 
-# this is needed very early because it is used by Range and colon
-function round(x::Float64)
-    y = trunc(x)
-    ifelse(x==y,y,trunc(2.0*x-y))
-end
 floor(x::Float64) = box(Float64,floor_llvm(unbox(Float64,x)))
+floor(x::Float32) = box(Float32,floor_llvm(unbox(Float32,x)))
+
+ceil(x::Float64) = box(Float64,ceil_llvm(unbox(Float64,x)))
+ceil(x::Float32) = box(Float32,ceil_llvm(unbox(Float32,x)))
+
+round(x::Float64) = box(Float64,rint_llvm(unbox(Float64,x)))
+round(x::Float32) = box(Float32,rint_llvm(unbox(Float32,x)))
 
 ## floating point promotions ##
 promote_rule(::Type{Float32}, ::Type{Float16}) = Float32
@@ -347,13 +351,6 @@ for Ti in (Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UIn
     end
 end
 
-# adding prevfloat(0.5) will prevent prevfloat(0.5) and odd x with eps(x)=1.0
-# from rounding in the wrong direction in RoundToNearest
-for Tf in (Float32,Float64)
-    @eval function round{T<:Integer}(::Type{T}, x::$Tf)
-        trunc(T,x+copysign($(prevfloat(Tf(0.5))),x))
-    end
-end
 
 
 @eval begin

base/floatfuncs.jl

@@ -43,6 +43,27 @@ end
 @vectorize_1arg Number isinf
 @vectorize_1arg Number isfinite
 
+
+round(x::Real, ::RoundingMode{:ToZero}) = trunc(x)
+round(x::Real, ::RoundingMode{:Up}) = ceil(x)
+round(x::Real, ::RoundingMode{:Down}) = floor(x)
+# C-style round
+function round(x::FloatingPoint, ::RoundingMode{:NearestTiesAway})
+    y = trunc(x)
+    ifelse(x==y,y,trunc(2*x-y))
+end
+# Java-style round
+function round(x::FloatingPoint, ::RoundingMode{:NearestTiesUp})
+    y = floor(x)
+    ifelse(x==y,y,copysign(floor(2*x-y),x))
+end
+round{T<:Integer}(::Type{T}, x::FloatingPoint, r::RoundingMode) = trunc(T,round(x,r))
+
+@vectorize_1arg Real trunc
+@vectorize_1arg Real floor
+@vectorize_1arg Real ceil
+@vectorize_1arg Real round
+
 for f in (:trunc,:floor,:ceil,:round)
     @eval begin
         function ($f){T,R<:Real}(::Type{T}, x::AbstractArray{R,1})
@@ -57,6 +78,25 @@ for f in (:trunc,:floor,:ceil,:round)
     end
 end
 
+function round{R<:Real}(x::AbstractArray{R,1}, r::RoundingMode)
+    [ round(x[i], r) for i = 1:length(x) ]
+end
+function round{R<:Real}(x::AbstractArray{R,2}, r::RoundingMode)
+    [ round(x[i,j], r) for i = 1:size(x,1), j = 1:size(x,2) ]
+end
+function round{R<:Real}(x::AbstractArray{R}, r::RoundingMode)
+    reshape([ round(x[i], r) for i = 1:length(x) ], size(x))
+end
+
+function round{T,R<:Real}(::Type{T}, x::AbstractArray{R,1}, r::RoundingMode)
+    [ round(T, x[i], r) for i = 1:length(x) ]
+end
+function round{T,R<:Real}(::Type{T}, x::AbstractArray{R,2}, r::RoundingMode)
+    [ round(T, x[i,j], r) for i = 1:size(x,1), j = 1:size(x,2) ]
+end
+function round{T,R<:Real}(::Type{T}, x::AbstractArray{R}, r::RoundingMode)
+    reshape([ round(T, x[i], r) for i = 1:length(x) ], size(x))
+end
 
 # adapted from Matlab File Exchange roundsd: http://www.mathworks.com/matlabcentral/fileexchange/26212
 # for round, og is the power of 10 relative to the decimal point

base/math.jl

@@ -9,7 +9,7 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,
        rad2deg, deg2rad,
        log, log2, log10, log1p, exponent, exp, exp2, exp10, expm1,
        cbrt, sqrt, erf, erfc, erfcx, erfi, dawson,
-       ceil, floor, trunc, round, significand,
+       significand,
        lgamma, hypot, gamma, lfact, max, min, minmax, ldexp, frexp,
        clamp, modf, ^, mod2pi,
        airy, airyai, airyprime, airyaiprime, airybi, airybiprime, airyx,
@@ -22,10 +22,10 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,
 
 import Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin,
              acos, atan, asinh, acosh, atanh, sqrt, log2, log10,
-             max, min, minmax, ceil, floor, trunc, round, ^, exp2,
+             max, min, minmax, ^, exp2,
              exp10, expm1, log1p
 
-import Core.Intrinsics: nan_dom_err, ceil_llvm, floor_llvm, trunc_llvm, sqrt_llvm, box, unbox, powi_llvm
+import Core.Intrinsics: nan_dom_err, sqrt_llvm, box, unbox, powi_llvm
 
 # non-type specific math functions
 
@@ -132,30 +132,15 @@ sqrt(x::Float32) = box(Float32,sqrt_llvm(unbox(Float32,x)))
 sqrt(x::Real) = sqrt(float(x))
 @vectorize_1arg Number sqrt
 
-ceil(x::Float64) = box(Float64,ceil_llvm(unbox(Float64,x)))
-ceil(x::Float32) = box(Float32,ceil_llvm(unbox(Float32,x)))
-@vectorize_1arg Real ceil
 
-trunc(x::Float64) = box(Float64,trunc_llvm(unbox(Float64,x)))
-trunc(x::Float32) = box(Float32,trunc_llvm(unbox(Float32,x)))
-@vectorize_1arg Real trunc
-
-for f in (:significand, :rint) # :nearbyint
+for f in (:significand,)
     @eval begin
         ($f)(x::Float64) = ccall(($(string(f)),libm), Float64, (Float64,), x)
         ($f)(x::Float32) = ccall(($(string(f,"f")),libm), Float32, (Float32,), x)
         @vectorize_1arg Real $f
     end
 end
 
-function round(x::Float32)
-    y = trunc(x)
-    ifelse(x==y,y,trunc(2.f0*x-y))
-end
-@vectorize_1arg Real round
-
-floor(x::Float32) = box(Float32,floor_llvm(unbox(Float32,x)))
-@vectorize_1arg Real floor
 
 hypot(x::Real, y::Real) = hypot(promote(float(x), float(y))...)
 function hypot{T<:FloatingPoint}(x::T, y::T)

base/mpfr.jl

@@ -26,7 +26,7 @@ import Base.GMP: ClongMax, CulongMax, CdoubleMax
 
 import Base.Math.lgamma_r
 
-const ROUNDING_MODE = [0]
+const ROUNDING_MODE = Cint[0]
 const DEFAULT_PRECISION = [256]
 
 # Basic type and initialization definitions
@@ -93,29 +93,32 @@ convert(::Type{BigFloat}, x::Real) = BigFloat(x)
 convert(::Type{FloatingPoint}, x::BigInt) = BigFloat(x)
 
 ## BigFloat -> Integer
-function unsafe_cast(::Type{Int64}, x::BigFloat, r::RoundingMode)
+function unsafe_cast(::Type{Int64}, x::BigFloat, ri::Cint)
     ccall((:__gmpfr_mpfr_get_sj,:libmpfr), Cintmax_t,
-          (Ptr{BigFloat}, Int32), &x, to_mpfr(r))
+          (Ptr{BigFloat}, Cint), &x, ri)
 end
-function unsafe_cast(::Type{UInt64}, x::BigFloat, r::RoundingMode)
+function unsafe_cast(::Type{UInt64}, x::BigFloat, ri::Cint)
     ccall((:__gmpfr_mpfr_get_uj,:libmpfr), Cuintmax_t,
-          (Ptr{BigFloat}, Int32), &x, to_mpfr(r))
+          (Ptr{BigFloat}, Cint), &x, ri)
 end
 
-function unsafe_cast{T<:Signed}(::Type{T}, x::BigFloat, r::RoundingMode)
-    unsafe_cast(Int64, x, r) % T
+function unsafe_cast{T<:Signed}(::Type{T}, x::BigFloat, ri::Cint)
+    unsafe_cast(Int64, x, ri) % T
 end
-function unsafe_cast{T<:Unsigned}(::Type{T}, x::BigFloat, r::RoundingMode)
-    unsafe_cast(UInt64, x, r) % T
+function unsafe_cast{T<:Unsigned}(::Type{T}, x::BigFloat, ri::Cint)
+    unsafe_cast(UInt64, x, ri) % T
 end
 
-function unsafe_cast(::Type{BigInt}, x::BigFloat, r::RoundingMode)
+function unsafe_cast(::Type{BigInt}, x::BigFloat, ri::Cint)
     # actually safe, just keep naming consistent
     z = BigInt()
     ccall((:mpfr_get_z, :libmpfr), Int32, (Ptr{BigInt}, Ptr{BigFloat}, Int32),
-          &z, &x, to_mpfr(r))
+          &z, &x, ri)
     z
 end
+unsafe_cast(::Type{Int128}, x::BigFloat, ri::Cint) = Int128(unsafe_cast(BigInt,x,ri))
+unsafe_cast(::Type{UInt128}, x::BigFloat, ri::Cint) = UInt128(unsafe_cast(BigInt,x,ri))
+unsafe_cast{T<:Integer}(::Type{T}, x::BigFloat, r::RoundingMode) = unsafe_cast(T,x,to_mpfr(r))
 
 unsafe_trunc{T<:Integer}(::Type{T}, x::BigFloat) = unsafe_cast(T,x,RoundToZero)
 
@@ -132,21 +135,21 @@ function ceil{T<:Union(Signed,Unsigned)}(::Type{T}, x::BigFloat)
     unsafe_cast(T,x,RoundUp)
 end
 
+function round{T<:Union(Signed,Unsigned)}(::Type{T}, x::BigFloat)
+    (typemin(T) <= x <= typemax(T)) || throw(InexactError())
+    unsafe_cast(T,x,ROUNDING_MODE[end])
+end
+
 trunc(::Type{BigInt}, x::BigFloat) = unsafe_cast(BigInt, x, RoundToZero)
 floor(::Type{BigInt}, x::BigFloat) = unsafe_cast(BigInt, x, RoundDown)
 ceil(::Type{BigInt}, x::BigFloat) = unsafe_cast(BigInt, x, RoundUp)
+round(::Type{BigInt}, x::BigFloat) = unsafe_cast(BigInt, x, ROUNDING_MODE[end])
+
 # convert/round/trunc/floor/ceil(Integer, x) should return a BigInt
 trunc(::Type{Integer}, x::BigFloat) = trunc(BigInt, x)
 floor(::Type{Integer}, x::BigFloat) = floor(BigInt, x)
 ceil(::Type{Integer}, x::BigFloat) = ceil(BigInt, x)
-
-for Ti in (Int128,UInt128)
-    @eval begin
-        trunc(::Type{$Ti}, x::BigFloat) = ($Ti)(trunc(BigInt, x))
-        floor(::Type{$Ti}, x::BigFloat) = ($Ti)(floor(BigInt, x))
-        ceil(::Type{$Ti}, x::BigFloat) = ($Ti)(ceil(BigInt, x))
-    end
-end
+round(::Type{Integer}, x::BigFloat) = round(BigInt, x)
 
 convert(::Type{Bool}, x::BigFloat) = (x != 0)
 function convert(::Type{BigInt},x::BigFloat)
@@ -625,11 +628,11 @@ end
 maxintfloat(x::BigFloat) = BigFloat(2)^precision(x)
 maxintfloat(::Type{BigFloat}) = BigFloat(2)^get_bigfloat_precision()
 
-to_mpfr(::RoundingMode{:TiesToEven}) = 0
-to_mpfr(::RoundingMode{:TowardZero}) = 1
-to_mpfr(::RoundingMode{:TowardPositive}) = 2
-to_mpfr(::RoundingMode{:TowardNegative}) = 3
-to_mpfr(::RoundingMode{:AwayFromZero}) = 4
+to_mpfr(::RoundingMode{:Nearest}) = Cint(0)
+to_mpfr(::RoundingMode{:ToZero}) = Cint(1)
+to_mpfr(::RoundingMode{:Up}) = Cint(2)
+to_mpfr(::RoundingMode{:Down}) = Cint(3)
+to_mpfr(::RoundingMode{:FromZero}) = Cint(4)
 
 function from_mpfr(c::Integer)
     if c == 0
@@ -687,7 +690,7 @@ function isinteger(x::BigFloat)
     return ccall((:mpfr_integer_p, :libmpfr), Int32, (Ptr{BigFloat},), &x) != 0
 end
 
-for f in (:ceil, :floor, :trunc, :round)
+for f in (:ceil, :floor, :trunc)
     @eval begin
         function ($f)(x::BigFloat)
             z = BigFloat()
@@ -697,6 +700,17 @@ for f in (:ceil, :floor, :trunc, :round)
     end
 end
 
+function round(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_rint, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Cint), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+function round(x::BigFloat,::RoundingMode{:NearestTiesAway})
+    z = BigFloat()
+    ccall((:mpfr_round, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &x)
+    return z
+end
+
 function isinf(x::BigFloat)
     return ccall((:mpfr_inf_p, :libmpfr), Int32, (Ptr{BigFloat},), &x) != 0
 end

base/printf.jl

@@ -909,7 +909,7 @@ function ini_hex(x::SmallFloatingPoint, n::Int, symbols::Array{UInt8,1})
     else
         s, p = frexp(x)
         sigbits = 4*min(n-1,13)
-        s = 0.25*Base.Math.rint(ldexp(s,1+sigbits))
+        s = 0.25*round(ldexp(s,1+sigbits))
         # ensure last 2 exponent bits either 01 or 10
         u = (reinterpret(UInt64,s) & 0x003f_ffff_ffff_ffff) >> (52-sigbits)
         if n > 14

base/rational.jl

@@ -220,13 +220,25 @@ for op in (:div, :fld, :cld)
     end
 end
 
-trunc{T<:Integer}(::Type{T}, x::Rational) = convert(T,div(x.num,x.den))
-floor{T<:Integer}(::Type{T}, x::Rational) = convert(T,fld(x.num,x.den))
-ceil {T<:Integer}(::Type{T}, x::Rational) = convert(T,cld(x.num,x.den))
-function round{T<:Integer}(::Type{T}, x::Rational)
-    t = trunc(T,x)
-    r = x-t
-    abs(r.num) > (r.den-one(r.den))>>1 ? t + copysign(one(t),x) : t
+trunc{T}(::Type{T}, x::Rational) = convert(T,div(x.num,x.den))
+floor{T}(::Type{T}, x::Rational) = convert(T,fld(x.num,x.den))
+ceil {T}(::Type{T}, x::Rational) = convert(T,cld(x.num,x.den))
+
+function round{T}(::Type{T}, x::Rational, ::RoundingMode{:Nearest})
+    q,r = divrem(x.num,x.den)
+    s = abs(r) < (x.den+one(x.den)+iseven(q))>>1 ? q : q+copysign(one(q),x.num)
+    convert(T,s)
+end
+round{T}(::Type{T}, x::Rational) = round(T,x,RoundNearest)
+function round{T}(::Type{T}, x::Rational, ::RoundingMode{:NearestTiesAway})
+    q,r = divrem(x.num,x.den)
+    s = abs(r) < (x.den+one(x.den))>>1 ? q : q+copysign(one(q),x.num)
+    convert(T,s)
+end
+function round{T}(::Type{T}, x::Rational, ::RoundingMode{:NearestTiesUp})
+    q,r = divrem(x.num,x.den)
+    s = abs(r) < (x.den+one(x.den)+(x.num<0))>>1 ? q : q+copysign(one(q),x.num)
+    convert(T,s)
 end
 
 trunc{T}(x::Rational{T}) = Rational(trunc(T,x))

base/rounding.jl

@@ -3,22 +3,26 @@ include("fenv_constants.jl")
 
 export
     RoundingMode, RoundNearest, RoundToZero, RoundUp, RoundDown, RoundFromZero,
+    RoundNearestTiesAway, RoundNearestTiesUp,
     get_rounding, set_rounding, with_rounding
 
 ## rounding modes ##
 immutable RoundingMode{T} end
 
-const RoundNearest = RoundingMode{:TiesToEven}()
-# const RoundNearestTiesAway = RoundingMode{:TiesToAway}() # currently unsupported
-const RoundToZero = RoundingMode{:TowardZero}()
-const RoundUp = RoundingMode{:TowardPositive}()
-const RoundDown = RoundingMode{:TowardNegative}()
-const RoundFromZero = RoundingMode{:AwayFromZero}() # mpfr only
+const RoundNearest = RoundingMode{:Nearest}()
+const RoundToZero = RoundingMode{:ToZero}()
+const RoundUp = RoundingMode{:Up}()
+const RoundDown = RoundingMode{:Down}()
+const RoundFromZero = RoundingMode{:FromZero}() # mpfr only
+# C-style round behaviour
+const RoundNearestTiesAway = RoundingMode{:NearestTiesAway}()
+# Java-style round behaviour
+const RoundNearestTiesUp = RoundingMode{:NearestTiesUp}()
 
-to_fenv(::RoundingMode{:TiesToEven}) = JL_FE_TONEAREST
-to_fenv(::RoundingMode{:TowardZero}) = JL_FE_TOWARDZERO
-to_fenv(::RoundingMode{:TowardPositive}) = JL_FE_UPWARD
-to_fenv(::RoundingMode{:TowardNegative}) = JL_FE_DOWNWARD
+to_fenv(::RoundingMode{:Nearest}) = JL_FE_TONEAREST
+to_fenv(::RoundingMode{:ToZero}) = JL_FE_TOWARDZERO
+to_fenv(::RoundingMode{:Up}) = JL_FE_UPWARD
+to_fenv(::RoundingMode{:Down}) = JL_FE_DOWNWARD
 
 function from_fenv(r::Integer)
     if r == JL_FE_TONEAREST
@@ -61,16 +65,16 @@ end
 # To avoid ambiguous dispatch with methods in mpfr.jl:
 call{T<:FloatingPoint}(::Type{T},x::Real,r::RoundingMode) = _convert_rounding(T,x,r)
 
-_convert_rounding{T<:FloatingPoint}(::Type{T},x::Real,r::RoundingMode{:TiesToEven}) = convert(T,x)
-function _convert_rounding{T<:FloatingPoint}(::Type{T},x::Real,r::RoundingMode{:TowardNegative})
+_convert_rounding{T<:FloatingPoint}(::Type{T},x::Real,r::RoundingMode{:Nearest}) = convert(T,x)
+function _convert_rounding{T<:FloatingPoint}(::Type{T},x::Real,r::RoundingMode{:Down})
     y = convert(T,x)
     y > x ? prevfloat(y) : y
 end
-function _convert_rounding{T<:FloatingPoint}(::Type{T},x::Real,r::RoundingMode{:TowardPositive})
+function _convert_rounding{T<:FloatingPoint}(::Type{T},x::Real,r::RoundingMode{:Up})
     y = convert(T,x)
     y < x ? nextfloat(y) : y
 end
-function _convert_rounding{T<:FloatingPoint}(::Type{T},x::Real,r::RoundingMode{:TowardZero})
+function _convert_rounding{T<:FloatingPoint}(::Type{T},x::Real,r::RoundingMode{:ToZero})
     y = convert(T,x)
     if x > 0.0
         y > x ? prevfloat(y) : y