diff --git a/base/.gitignore b/base/.gitignore index 77ecc7d193ca7..fd5e9b3ee9f4d 100644 --- 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 index 4bcade1fbad20..0f392ef909faf 100644 --- a/base/Makefile +++ b/base/Makefile @@ -26,7 +26,7 @@ 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 > $@) @@ -34,6 +34,9 @@ pcre_h.jl: 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 index 77d785d2298a7..36f327d81dc0b 100644 --- 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 index bfc28f44e2ab8..30ce5250d9736 100644 --- 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,7 +20,9 @@ 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 @@ -27,11 +30,7 @@ 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 new file mode 100644 index 0000000000000..736a8bbeed06e --- /dev/null +++ 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 \ No newline at end of file diff --git a/base/sysimg.jl b/base/sysimg.jl index 397fb25dde002..82e2a0e015136 100644 --- 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 index 5a666da88ed35..dda4f56c4a2fc 100644 --- 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 `_:: + + + 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 index f7e97f601b9d6..f4e0928811fa1 100644 --- 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 new file mode 100644 index 0000000000000..9c9e96640b597 --- /dev/null +++ b/src/fenv_constants.h @@ -0,0 +1,10 @@ +#include +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 index d61a8d6e1259f..0d07f92cf1c2b 100644 --- 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 index 257d09a7b85d8..18b71cf66df83 100644 --- 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 diff --git a/test/rounding.jl b/test/rounding.jl new file mode 100644 index 0000000000000..0194dd3567b36 --- /dev/null +++ b/test/rounding.jl @@ -0,0 +1,89 @@ +# Small sanity tests to ensure changing the rounding of float functions work +using Base.Test + +## Float64 checks +# a + b returns a number exactly between prevfloat(1.) and 1., so its +# final result depends strongly on the utilized rounding direction. +a = prevfloat(0.5) +b = 0.5 +c = 0x1p-54 +d = prevfloat(1.) + +# Default rounding direction, RoundNearest +@test a + b === 1. +@test - a - b === -1. +@test a - b === -c +@test b - a === c + +# RoundToZero +with_rounding(RoundToZero) do + @test a + b === d + @test - a - b === -d + @test a - b === -c + @test b - a === c +end + +# Sanity check to see if we have returned to RoundNearest +@test a + b === 1. +@test - a - b === -1. +@test a - b == -c +@test b - a == c + +# RoundUp +with_rounding(RoundUp) do + @test a + b === 1. + @test - a - b === -d + @test a - b === -c + @test b - a === c +end + +# RoundDown +with_rounding(RoundDown) do + @test a + b === d + @test - a - b === -1. + @test a - b === -c + @test b - a === c +end + +## Float32 checks + +a32 = prevfloat(0.5f0) +b32 = 0.5f0 +c32 = (1.f0 - prevfloat(1.f0))/2 +d32 = prevfloat(1.0f0) + +# Default rounding direction, RoundNearest +@test a32 + b32 === 1.0f0 +@test - a32 - b32 === -1.0f0 +@test a32 - b32 === -c32 +@test b32 - a32 === c32 + +# RoundToZero +with_rounding(RoundToZero) do + @test a32 + b32 === d32 + @test - a32 - b32 === -d32 + @test a32 - b32 === -c32 + @test b32 - a32 === c32 +end + +# Sanity check to see if we have returned to RoundNearest +@test a32 + b32 === 1.0f0 +@test - a32 - b32 === -1.0f0 +@test a32 - b32 == -c32 +@test b32 - a32 == c32 + +# RoundUp +with_rounding(RoundUp) do + @test a32 + b32 === 1.0f0 + @test - a32 - b32 === -d32 + @test a32 - b32 === -c32 + @test b32 - a32 === c32 +end + +# RoundDown +with_rounding(RoundDown) do + @test a32 + b32 === d32 + @test - a32 - b32 === -1.0f0 + @test a32 - b32 === -c32 + @test b32 - a32 === c32 +end diff --git a/test/runtests.jl b/test/runtests.jl index a2b270081890b..4ef5791992c9a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,7 @@ testnames = ["core", "keywordargs", "numbers", "strings", "unicode", "arpack", "file", "suitesparse", "version", "resolve", "pollfd", "mpfr", "broadcast", "complex", "socket", "floatapprox", "readdlm", "regex", "float16", - "combinatorics", "sysinfo"] + "combinatorics", "sysinfo", "rounding"] tests = ARGS==["all"] ? testnames : ARGS