Merge pull request #3149 from andrioni/floatrounding

Ability to change rounding mode for all floats (cf. #2976)
JuliaLang · Sep 2, 2013 · c43e84a · c43e84a
2 parents 4467885 + 773b60e
commit c43e84a
Show file tree

Hide file tree

Showing 13 changed files with 258 additions and 24 deletions.
diff --git a/base/.gitignore b/base/.gitignore
@@ -1,5 +1,6 @@
 /pcre_h.jl
 /errno_h.jl
 /build_h.jl
+/fenv_constants.jl
 /file_constants.jl
 /uv_constants.jl
diff --git a/base/Makefile b/base/Makefile
@@ -26,14 +26,17 @@ dirty = $(shell [ -z "$(shell git status --porcelain 2>/dev/null)" ] && echo "fa
 
 
 
-all: pcre_h.jl errno_h.jl build_h.jl.phony file_constants.jl uv_constants.jl
+all: pcre_h.jl errno_h.jl build_h.jl.phony fenv_constants.jl file_constants.jl uv_constants.jl
 
 pcre_h.jl:
 	@$(call PRINT_PERL, $(CPP) -dM $(shell $(PCRE_CONFIG) --prefix)/include/pcre.h | perl -nle '/^\s*#define\s+PCRE_(\w*)\s*\(?($(PCRE_CONST))\)?\s*$$/ and print "const $$1 = uint32($$2)"' | sort > $@)
 
 errno_h.jl:
 	@$(call PRINT_PERL, echo '#include "errno.h"' | $(CPP) -dM - | perl -nle 'print "const $$1 = int32($$2)" if /^#define\s+(E\w+)\s+(\d+)\s*$$/' | sort > $@)
 
+fenv_constants.jl: ../src/fenv_constants.h
+	@$(PRINT_PERL) ${CPP} -P -DJULIA ../src/fenv_constants.h | tail -n 8 | perl -ple 's/\sFE_UN\w+/ 0x10/g; s/\sFE_O\w+/ 0x08/g; s/\sFE_DI\w+/ 0x04/g; s/\sFE_INV\w+/ 0x01/g; s/\sFE_TON\w+/ 0x00/g; s/\sFE_UP\w+/ 0x800/g; s/\sFE_DO\w+/ 0x400/g; s/\sFE_TOW\w+/ 0xc00/g' > $@
+
 file_constants.jl: ../src/file_constants.h
 	@$(call PRINT_PERL, $(CPP) -P -DJULIA ../src/file_constants.h | perl -nle 'print "$$1 0o$$2" if /^(\s*const\s+[A-z_]+\s+=)\s+(0[0-9]*)\s*$$/; print "$$1" if /^\s*(const\s+[A-z_]+\s+=\s+([1-9]|0x)[0-9A-z]*)\s*$$/' > $@)
 
@@ -98,5 +101,6 @@ clean:
 	rm -f pcre_h.jl
 	rm -f errno_h.jl
 	rm -f build_h.jl
+	rm -f fenv_constants.jl
 	rm -f uv_constants.jl
 	rm -f file_constants.jl
diff --git a/base/exports.jl b/base/exports.jl
@@ -82,6 +82,12 @@ export
     Reverse,
     RevString,
     RopeString,
+    RoundFromZero,
+    RoundDown,
+    RoundingMode,
+    RoundNearest,
+    RoundToZero,
+    RoundUp,
     Schur,
     Set,
     SparseMatrixCSC,
@@ -828,6 +834,9 @@ export
     get_bigfloat_rounding,
     set_bigfloat_rounding,
     with_bigfloat_rounding,
+    get_rounding,
+    set_rounding,
+    with_rounding,
 
 # statistics
     cor,

diff --git a/base/mpfr.jl b/base/mpfr.jl
@@ -2,6 +2,7 @@ module MPFR
 
 export
     BigFloat,
+    RoundFromZero,
     get_bigfloat_precision,
     set_bigfloat_precision,
     with_bigfloat_precision,
@@ -19,19 +20,17 @@ import
         gamma, lgamma, digamma, erf, erfc, zeta, log1p, airyai, iceil, ifloor,
         itrunc, eps, signbit, sin, cos, tan, sec, csc, cot, acos, asin, atan,
         cosh, sinh, tanh, sech, csch, coth, acosh, asinh, atanh, atan2,
-        serialize, deserialize, inf, nan, hash
+        serialize, deserialize, inf, nan, hash,
+        RoundingMode, RoundDown, RoundingMode, RoundNearest, RoundToZero, 
+        RoundUp
 
 import Base.Math.lgamma_r
 
 const ROUNDING_MODE = [0]
 const DEFAULT_PRECISION = [256]
 
 # Rounding modes
-const RoundToNearest = 0
-const RoundToZero = 1
-const RoundUp = 2
-const RoundDown = 3
-const RoundAwayZero = 4
+type RoundFromZero <: RoundingMode end
 
 # Basic type and initialization definitions
 
@@ -98,7 +97,7 @@ for to in (Int8, Int16, Int32, Int64)
         function convert(::Type{$to}, x::BigFloat)
             (isinteger(x) && (typemin($to) <= x <= typemax($to))) || throw(InexactError())
             convert($to, ccall((:mpfr_get_si,:libmpfr),
-                               Clong, (Ptr{BigFloat}, Int32), &x, RoundToZero))
+                               Clong, (Ptr{BigFloat}, Int32), &x, 0))
         end
     end
 end
@@ -108,7 +107,7 @@ for to in (Uint8, Uint16, Uint32, Uint64)
         function convert(::Type{$to}, x::BigFloat)
             (isinteger(x) && (typemin($to) <= x <= typemax($to))) || throw(InexactError())
             convert($to, ccall((:mpfr_get_ui,:libmpfr),
-                               Culong, (Ptr{BigFloat}, Int32), &x, RoundToZero))
+                               Culong, (Ptr{BigFloat}, Int32), &x, 0))
         end
     end
 end
@@ -597,13 +596,26 @@ function set_bigfloat_precision(x::Int)
     DEFAULT_PRECISION[end] = x
 end
 
-get_bigfloat_rounding() = ROUNDING_MODE[end]
-function set_bigfloat_rounding(x::Int)
-    if x < 0 || x > 4
-        throw(DomainError())
+function get_bigfloat_rounding()
+    if ROUNDING_MODE[end] == 0
+        return RoundNearest
+    elseif ROUNDING_MODE[end] == 1
+        return RoundToZero
+    elseif ROUNDING_MODE[end] == 2
+        return RoundUp
+    elseif ROUNDING_MODE[end] == 3
+        return RoundDown
+    elseif ROUNDING_MODE[end] == 4
+        return RoundFromZero
+    else
+        error("Invalid rounding mode")
     end
-    ROUNDING_MODE[end] = x
 end
+set_bigfloat_rounding(::Type{RoundNearest}) = ROUNDING_MODE[end] = 0
+set_bigfloat_rounding(::Type{RoundToZero}) = ROUNDING_MODE[end] = 1
+set_bigfloat_rounding(::Type{RoundUp}) = ROUNDING_MODE[end] = 2
+set_bigfloat_rounding(::Type{RoundDown}) = ROUNDING_MODE[end] = 3
+set_bigfloat_rounding(::Type{RoundFromZero}) = ROUNDING_MODE[end] = 4
 
 function copysign(x::BigFloat, y::BigFloat)
     z = BigFloat()
@@ -635,7 +647,7 @@ end
 
 function itrunc(x::BigFloat)
     z = BigInt()
-    ccall((:mpfr_get_z, :libmpfr), Int32, (Ptr{BigInt}, Ptr{BigFloat}, Int32), &z, &x, RoundToZero)
+    ccall((:mpfr_get_z, :libmpfr), Int32, (Ptr{BigInt}, Ptr{BigFloat}, Int32), &z, &x, 0)
     return z
 end
 
@@ -682,7 +694,7 @@ function with_bigfloat_precision(f::Function, precision::Integer)
     end
 end
 
-function with_bigfloat_rounding(f::Function, rounding::Integer)
+function with_bigfloat_rounding{T<:RoundingMode}(f::Function, rounding::Type{T})
     old_rounding = get_bigfloat_rounding()
     set_bigfloat_rounding(rounding)
     try

diff --git a/base/rounding.jl b/base/rounding.jl
@@ -0,0 +1,45 @@
+module Rounding
+include("fenv_constants.jl")
+
+export
+    RoundingMode, RoundNearest, RoundToZero, RoundUp, RoundDown,
+    get_rounding, set_rounding, with_rounding
+
+## rounding modes ##
+abstract RoundingMode
+type RoundNearest <: RoundingMode end
+type RoundToZero <: RoundingMode end
+type RoundUp <: RoundingMode end
+type RoundDown <: RoundingMode end
+
+set_rounding(::Type{RoundNearest}) = ccall(:fesetround, Cint, (Cint, ), JL_FE_TONEAREST)
+set_rounding(::Type{RoundToZero}) = ccall(:fesetround, Cint, (Cint, ), JL_FE_TOWARDZERO)
+set_rounding(::Type{RoundUp}) = ccall(:fesetround, Cint, (Cint, ), JL_FE_UPWARD)
+set_rounding(::Type{RoundDown}) = ccall(:fesetround, Cint, (Cint, ), JL_FE_DOWNWARD)
+
+function get_rounding()
+    r = ccall(:fegetround, Cint, ())
+    if r == JL_FE_TONEAREST
+        return RoundNearest
+    elseif r == JL_FE_DOWNWARD
+        return RoundDown
+    elseif r == JL_FE_UPWARD
+        return RoundUp
+    elseif r == JL_FE_TOWARDZERO
+        return RoundToZero
+    else
+        error()
+    end
+end
+
+function with_rounding{T<:RoundingMode}(f::Function, rounding::Type{T})
+    old_rounding = get_rounding()
+    set_rounding(rounding)
+    try
+        return f()
+    finally
+        set_rounding(old_rounding)
+    end
+end
+
+end #module
diff --git a/base/sysimg.jl b/base/sysimg.jl
@@ -182,6 +182,10 @@ include("fftw.jl")
 include("dsp.jl")
 importall .DSP
 
+# rounding utilities
+include("rounding.jl")
+importall .Rounding
+
 # BigInts and BigFloats
 include("gmp.jl")
 importall .GMP

diff --git a/doc/manual/integers-and-floating-point-numbers.rst b/doc/manual/integers-and-floating-point-numbers.rst
@@ -447,6 +447,26 @@ argument respectively: ::
 This example highlights the general principle that the adjacent representable
 floating-point numbers also have adjacent binary integer representations.
 
+Rounding modes
+~~~~~~~~~~~~~~
+
+If a number doesn't have an exact floating-point representation, it must be
+rounded to an appropriate representable value, however, if wanted, the manner
+in which this rounding is done can be changed according to the rounding modes
+presented in the `IEEE 754 standard <http://en.wikipedia.org/wiki/IEEE_754-2008>`_::
+    
+
+    julia> 1.1 + 0.1
+    1.2000000000000002
+
+    julia> with_rounding(RoundDown) do
+           1.1 + 0.1
+           end
+    1.2
+
+The default mode used is always ``RoundNearest``, which rounds to the nearest
+representable value, with ties rounded towards the nearest value with an even
+least significant bit.
 
 Background and References
 ~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -503,10 +523,10 @@ Once created, they participate in arithmetic with all other numeric types thanks
     123456789012345678901234567891
 
     julia> BigFloat("1.23456789012345678901")
-    1.23456789012345678901
+    1.234567890123456789010000000000000000000000000000000000000000000000000000000004e+00 with 256 bits of precision
 
     julia> BigFloat(2.0^66) / 3
-    24595658764946068821.3
+    2.459565876494606882133333333333333333333333333333333333333333333333333333333344e+19 with 256 bits of precision
 
     julia> factorial(BigInt(40))
     815915283247897734345611269596115894272000000000
@@ -531,6 +551,27 @@ However, type promotion between the primitive types above and
     
     julia> typeof(y)
     BigInt
+
+The default precision (in number of bits of the significand) and rounding
+mode of `BigFloat` operations can be changed, and all further calculations 
+will take these changes in account::
+
+    julia> with_bigfloat_rounding(RoundUp) do
+           BigFloat(1) + BigFloat("0.1")
+           end
+    1.100000000000000000000000000000000000000000000000000000000000000000000000000003e+00 with 256 bits of precision
+
+    julia> with_bigfloat_rounding(RoundDown) do
+           BigFloat(1) + BigFloat("0.1")
+           end
+    1.099999999999999999999999999999999999999999999999999999999999999999999999999986e+00 with 256 bits of precision
+
+    julia> with_bigfloat_precision(40) do
+           BigFloat(1) + BigFloat("0.1")
+           end
+    1.0999999999985e+00 with 40 bits of precision
+
+
    
 .. _man-numeric-literal-coefficients:
 

diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst
@@ -2981,6 +2981,25 @@ Numbers
    Create an arbitrary precision floating point number. ``x`` may be an ``Integer``, a ``Float64``, a ``String`` or a ``BigInt``. The
    usual mathematical operators are defined for this type, and results are promoted to a ``BigFloat``.
 
+.. function:: get_rounding()
+
+   Get the current floating point rounding mode. Valid modes are ``RoundNearest``, ``RoundToZero``, ``RoundUp`` and ``RoundDown``.
+
+.. function:: set_rounding(mode)
+
+   Set the floating point rounding mode. See ``get_rounding`` for available modes
+
+.. function:: with_rounding(f::Function,mode)
+
+   Change the floating point rounding mode for the duration of ``f``. It is logically equivalent to::
+
+       old = get_rounding()
+       set_rounding(mode)
+       f()
+       set_rounding(old)
+
+   See ``get_rounding`` for available rounding modes.
+
 Integers
 ~~~~~~~~
 
@@ -3069,7 +3088,7 @@ The `BigFloat` type implements arbitrary-precision floating-point aritmetic usin
 
 .. function:: get_bigfloat_rounding()
 
-   Get the current BigFloat rounding mode. Valid modes are ``RoundToNearest``, ``RoundToZero``, ``RoundUp``, ``RoundDown``, ``RoundAwayZero``
+   Get the current BigFloat rounding mode. Valid modes are ``RoundNearest``, ``RoundToZero``, ``RoundUp``, ``RoundDown``, ``RoundFromZero``
 
 .. function:: set_bigfloat_rounding(mode)
 

diff --git a/src/fenv_constants.h b/src/fenv_constants.h
@@ -0,0 +1,10 @@
+#include <fenv.h>
+const JL_FE_INEXACT = FE_INEXACT
+const JL_FE_UNDERFLOW = FE_UNDERFLOW
+const JL_FE_OVERFLOW = FE_OVERFLOW
+const JL_FE_DIVBYZERO = FE_DIVBYZERO
+const JL_FE_INVALID = FE_INVALID
+const JL_FE_TONEAREST = FE_TONEAREST
+const JL_FE_UPWARD = FE_UPWARD
+const JL_FE_DOWNWARD = FE_DOWNWARD
+const JL_FE_TOWARDZERO = FE_TOWARDZERO
diff --git a/test/Makefile b/test/Makefile
@@ -6,7 +6,7 @@ TESTS = all core keywordargs numbers strings unicode collections hashing	\
 	math functional bigint sorting statistics spawn parallel arpack file	\
 	git pkg resolve suitesparse complex version pollfd mpfr	broadcast       \
 	socket floatapprox priorityqueue readdlm regex float16 combinatorics    \
-	sysinfo
+	sysinfo rounding
 
 $(TESTS) ::
 	@$(call PRINT_JULIA, $(call spawn,$(JULIA_EXECUTABLE)) ./runtests.jl $@)

diff --git a/test/mpfr.jl b/test/mpfr.jl
@@ -82,14 +82,14 @@ z = BigFloat(30)
 # rounding modes
 with_bigfloat_precision(4) do
     # default mode is round to nearest
-    down, up =  with_bigfloat_rounding(MPFR.RoundToNearest) do
+    down, up =  with_bigfloat_rounding(RoundNearest) do
         BigFloat("0.0938"), BigFloat("0.102")
     end
-    with_bigfloat_rounding(MPFR.RoundDown) do
+    with_bigfloat_rounding(RoundDown) do
         @test BigFloat(0.1) == down
         @test BigFloat(0.1) != up
     end
-    with_bigfloat_rounding(MPFR.RoundUp) do
+    with_bigfloat_rounding(RoundUp) do
         @test BigFloat(0.1) != down
         @test BigFloat(0.1) == up
     end