diff --git a/src/Quadmath.jl b/src/Quadmath.jl index efcbe07..b3af96a 100644 --- a/src/Quadmath.jl +++ b/src/Quadmath.jl @@ -226,14 +226,15 @@ for Ti in (UInt8, UInt16, UInt32, UInt64, UInt128) # `Float128(typemax(Ti))+1` is exactly representable @eval begin function trunc(::Type{$Ti},x::Float128) - if !signbit(x) && x < Float128(typemax($Ti)) + one(Float128) + if (!signbit(x) && x < Float128(typemax($Ti)) + one(Float128)) || (x == 0) return unsafe_trunc($Ti,x) else throw(InexactError(:trunc, $Ti, x)) end end function (::Type{$Ti})(x::Float128) - if (!signbit(x) && x <= Float128(typemax($Ti))) && (round(x, RoundToZero) == x) + if ((!signbit(x) && x <= Float128(typemax($Ti))) && (round(x, RoundToZero) == x) + || (x == 0)) return unsafe_trunc($Ti,x) else throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x)) @@ -303,7 +304,8 @@ end @assume_effects :foldable round(x::Float128) = Float128(@quad_ccall(libquadmath.rintq(x::Cfloat128)::Cfloat128)) round(x::Float128, r::RoundingMode{:Down}) = floor(x) round(x::Float128, r::RoundingMode{:Up}) = ceil(x) -round(x::Float128, r::RoundingMode{:ToZero}) = round(x) +round(x::Float128, r::RoundingMode{:Nearest}) = round(x) +round(x::Float128, r::RoundingMode{:ToZero}) = trunc(x) ## two argument @assume_effects :foldable (^)(x::Float128, y::Float128) = @@ -498,8 +500,10 @@ function Float128(x::BigFloat) prec = 113 + (k + 16381) elseif k == -16381-113 && abs(y) > 0.5 z = reinterpret(Float128, UInt128(1)) + return copysign(z,x) else z = reinterpret(Float128, UInt128(0)) + return copysign(z,x) end y = BigFloat(x, precision=prec) diff --git a/test/rounding.jl b/test/rounding.jl new file mode 100644 index 0000000..7015bcc --- /dev/null +++ b/test/rounding.jl @@ -0,0 +1,310 @@ +# Derived from Julia Base test suite. License is MIT: https://julialang.org/license + +# Small sanity tests to ensure changing the rounding of float functions work +using Base.MathConstants + +q1 = one(Float128) +q0 = zero(Float128) + +@testset "Arithmetic rounding checks" begin + # a + b returns a number exactly between prevfloat(1.) and 1., so its + # final result depends strongly on the utilized rounding direction. + a = prevfloat(Float128(0.5)) + b = Float128(0.5) + c = eps(Float128)/4 + d = prevfloat(Float128(1)) + + @testset "Default rounding direction, RoundNearest" begin + @test a + b === q1 + @test - a - b === -q1 + @test a - b === -c + @test b - a === c + end +end + +@testset "convert with rounding" begin + for v = [sqrt(Float128(2)),-1/Float128(3),nextfloat(q1),prevfloat(q1),nextfloat(-q1), + prevfloat(-q1),nextfloat(q0),prevfloat(q0)] + + pn = Float64(v,RoundNearest) + @test pn == convert(Float64,v) + + pz = Float64(v,RoundToZero) + @test abs(pz) <= abs(v) < nextfloat(abs(pz)) + @test signbit(pz) == signbit(v) + + pd = Float64(v,RoundDown) + @test pd <= v < nextfloat(pd) + + pu = Float64(v,RoundUp) + @test prevfloat(pu) < v <= pu + + @test pn == pd || pn == pu + @test v > 0 ? pz == pd : pz == pu + @test pu - pd == eps(pz) + end + + # numbers which are not represented exactly in Float128 + for T in [Float128,] + for v in [sqrt(big(2.0)),-big(1.0)/big(3.0),nextfloat(big(1.0)), + prevfloat(big(1.0)), nextfloat(big(0.0)),prevfloat(big(0.0)), + pi,ℯ,eulergamma,catalan,golden, + typemax(Int128),typemax(UInt128), + ] + pn = T(v,RoundNearest) + @test pn == convert(T,BigFloat(v)) + pz = T(v,RoundToZero) + @test pz == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundToZero) + pd = T(v,RoundDown) + if v == prevfloat(big(0.0)) + @test_broken pd == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundDown) + else + @test pd == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundDown) + end + pu = T(v,RoundUp) + if v == nextfloat(big(0.0)) + @test_broken pu == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundUp) + else + @test pu == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundUp) + end + @test pn == pd || pn == pu + @test v > 0 ? pz == pd : pz == pu + @test isinf(pu) || pu - pd == eps(pz) + end + end +end + +@testset "rounding properties" for Tf in [Float128,] + # these should hold for all u, but we just test the smallest and largest + # of each binade + + for i in exponent(floatmin(Tf)):exponent(floatmax(Tf)) + for u in [ldexp(Tf(1.0), i), -ldexp(Tf(1.0), i), + ldexp(prevfloat(Tf(2.0)), i), -ldexp(prevfloat(Tf(2.0)), i)] + + r = round(u, RoundNearest) + if isfinite(u) + @test isfinite(r) + @test isinteger(r) + @test abs(r-u) < 0.5 || abs(r-u) == 0.5 && isinteger(r/2) + @test signbit(u) == signbit(r) + else + @test u === r + end + + r = round(u, RoundNearestTiesAway) + if isfinite(u) + @test isfinite(r) + @test isinteger(r) + @test abs(r-u) < 0.5 || (r-u) == copysign(0.5,u) + @test signbit(u) == signbit(r) + else + @test u === r + end + + r = round(u, RoundNearestTiesUp) + if isfinite(u) + @test isfinite(r) + @test isinteger(r) + @test -0.5 < r-u <= 0.5 + @test signbit(u) == signbit(r) + else + @test u === r + end + + r = round(u, RoundFromZero) + if isfinite(u) + @test isfinite(r) + @test isinteger(r) + @test signbit(u) ? (r == floor(u)) : (r == ceil(u)) + @test signbit(u) == signbit(r) + else + @test u === r + end + end + end +end + +@testset "rounding difficult values" begin + for x = Int128(2)^113-10:Int128(2)^113+10 + y = Float128(x) + i = trunc(Int128,y) + @test Int128(trunc(y)) == i + @test Int128(round(y)) == i + @test Int128(floor(y)) == i + @test Int128(ceil(y)) == i + + @test round(Int128,y) == i + @test floor(Int128,y) == i + @test ceil(Int128,y) == i + end + + # rounding vectors + let ≈(x,y) = x==y && typeof(x)==typeof(y) + for t in [Float128,] + # try different vector lengths + for n in [0,3,255,256] + r = (1:n) .- div(n,2) + y = t[x/4 for x in r] + @test trunc.(y) ≈ t[div(i,4) for i in r] + @test floor.(y) ≈ t[i>>2 for i in r] + @test ceil.(y) ≈ t[(i+3)>>2 for i in r] + @test round.(y) ≈ t[(i+1+isodd(i>>2))>>2 for i in r] + @test broadcast(x -> round(x, RoundNearestTiesAway), y) ≈ t[(i+1+(i>=0))>>2 for i in r] + @test broadcast(x -> round(x, RoundNearestTiesUp), y) ≈ t[(i+2)>>2 for i in r] + @test broadcast(x -> round(x, RoundFromZero), y) ≈ t[(i+3*(i>=0))>>2 for i in r] + end + end + end + + @test_throws InexactError round(Int,Float128(Inf)) + @test_throws InexactError round(Int,Float128(NaN)) + @test round(Int,Float128(2.5)) == 2 + @test round(Int,Float128(1.5)) == 2 + @test round(Int,Float128(-2.5)) == -2 + @test round(Int,Float128(-1.5)) == -2 + @test round(Int,Float128(2.5),RoundNearestTiesAway) == 3 + @test round(Int,Float128(1.5),RoundNearestTiesAway) == 2 + @test round(Int,Float128(2.5),RoundNearestTiesUp) == 3 + @test round(Int,Float128(1.5),RoundNearestTiesUp) == 2 + @test round(Int,Float128(-2.5),RoundNearestTiesAway) == -3 + @test round(Int,Float128(-1.5),RoundNearestTiesAway) == -2 + @test round(Int,Float128(-2.5),RoundNearestTiesUp) == -2 + @test round(Int,Float128(-1.5),RoundNearestTiesUp) == -1 + @test round(Int,Float128(-1.9)) == -2 + @test round(Int,nextfloat(q1),RoundFromZero) == 2 + @test round(Int,-nextfloat(q1),RoundFromZero) == -2 + @test round(Int,prevfloat(q1),RoundFromZero) == 1 + @test round(Int,-prevfloat(q1),RoundFromZero) == -1 + #= + FIXME: want analogs + @test_throws InexactError round(Int64, 9.223372036854776e18) + @test round(Int64, 9.223372036854775e18) == 9223372036854774784 + @test_throws InexactError round(Int64, -9.223372036854778e18) + @test round(Int64, -9.223372036854776e18) == typemin(Int64) + @test_throws InexactError round(UInt64, 1.8446744073709552e19) + @test round(UInt64, 1.844674407370955e19) == 0xfffffffffffff800 + =# + for Ti in [Int,UInt] + for Tf in [Float128,] + # Note: these threw InexactError + @test round(Ti,Tf(-0.0)) == 0 + @test round(Ti,Tf(-0.0),RoundNearestTiesAway) == 0 + @test round(Ti,Tf(-0.0),RoundNearestTiesUp) == 0 + + @test round(Ti, Tf(0.5)) == 0 + @test round(Ti, Tf(0.5), RoundNearestTiesAway) == 1 + @test round(Ti, Tf(0.5), RoundNearestTiesUp) == 1 + + @test round(Ti, prevfloat(Tf(0.5))) == 0 + @test round(Ti, prevfloat(Tf(0.5)), RoundNearestTiesAway) == 0 + @test round(Ti, prevfloat(Tf(0.5)), RoundNearestTiesUp) == 0 + + @test round(Ti, nextfloat(Tf(0.5))) == 1 + @test round(Ti, nextfloat(Tf(0.5)), RoundNearestTiesAway) == 1 + @test round(Ti, nextfloat(Tf(0.5)), RoundNearestTiesUp) == 1 + # Note: these threw InexactError + @test round(Ti, Tf(-0.5)) == 0 + @test round(Ti, Tf(-0.5), RoundNearestTiesUp) == 0 + + @test round(Ti, nextfloat(Tf(-0.5))) == 0 + @test round(Ti, nextfloat(Tf(-0.5)), RoundNearestTiesAway) == 0 + @test round(Ti, nextfloat(Tf(-0.5)), RoundNearestTiesUp) == 0 + + if Ti <: Signed + @test round(Ti, Tf(-0.5), RoundNearestTiesAway) == -1 + @test round(Ti, prevfloat(Tf(-0.5))) == -1 + @test round(Ti, prevfloat(Tf(-0.5)), RoundNearestTiesAway) == -1 + @test round(Ti, prevfloat(Tf(-0.5)), RoundNearestTiesUp) == -1 + else + @test_throws InexactError round(Ti, Tf(-0.5), RoundNearestTiesAway) + @test_throws InexactError round(Ti, prevfloat(Tf(-0.5))) + @test_throws InexactError round(Ti, prevfloat(Tf(-0.5)), RoundNearestTiesAway) + @test_throws InexactError round(Ti, prevfloat(Tf(-0.5)), RoundNearestTiesUp) + end + end + end + + # numbers that can't be rounded by trunc(x+0.5) + @test round(Int128, Float128(2.0)^112 + 1) == (UInt128(1) << 112) + 1 +end + +# custom rounding and significant-digit ops +@testset "rounding to digits relative to the decimal point" begin + @test round(Float128(pi)) ≈ 3. + @test round(Float128(pi), base=10) ≈ 3. + @test round(Float128(pi), digits=0) ≈ 3. + @test round(Float128(pi), digits=1) ≈ 3.1 + @test round(Float128(pi), digits=3, base=2) ≈ 3.125 + @test round(Float128(pi), sigdigits=1) ≈ 3. + @test round(Float128(pi), sigdigits=3) ≈ 3.14 + @test round(Float128(pi), sigdigits=4, base=2) ≈ 3.25 + @test round(10*Float128(pi), digits=-1) ≈ 30. + @test round(Float128(.1), digits=0) == 0. + @test round(Float128(-.1), digits=0) == -0. + @test isnan(round(Float128(NaN), digits=2)) + @test isinf(round(Float128(Inf), digits=2)) + @test isinf(round(Float128(-Inf), digits=2)) +end +@testset "round vs trunc vs floor vs ceil" begin + @test round(Float128(123.456), digits=1) ≈ 123.5 + @test round(-Float128(123.456), digits=1) ≈ -123.5 + @test trunc(Float128(123.456), digits=1) ≈ 123.4 + @test trunc(-Float128(123.456), digits=1) ≈ -123.4 + @test ceil(Float128(123.456), digits=1) ≈ 123.5 + @test ceil(-Float128(123.456), digits=1) ≈ -123.4 + @test floor(Float128(123.456), digits=1) ≈ 123.4 + @test floor(-Float128(123.456), digits=1) ≈ -123.5 +end +@testset "rounding with too much (or too few) precision" begin + for x in (12345.6789, 0, -12345.6789) + y = Float128(x) + @test y == trunc(x, digits=1000) + @test y == round(x, digits=1000) + @test y == floor(x, digits=1000) + @test y == ceil(x, digits=1000) + end + let x = Float128(12345.6789) + @test 0.0 == trunc(x, digits=-1000) + @test 0.0 == round(x, digits=-1000) + @test 0.0 == floor(x, digits=-1000) + @test Inf128 == ceil(x, digits=-10000) + end + let x = Float128(-12345.6789) + @test -0.0 == trunc(x, digits=-1000) + @test -0.0 == round(x, digits=-1000) + @test -Inf128 == floor(x, digits=-10000) + @test -0.0 == ceil(x, digits=-1000) + end + let x = q0 + @test 0.0 == trunc(x, digits=-1000) + @test 0.0 == round(x, digits=-1000) + @test 0.0 == floor(x, digits=-1000) + @test 0.0 == ceil(x, digits=-1000) + end +end +@testset "rounding in other bases" begin + @test round(Float128(pi), digits = 2, base = 2) ≈ 3.25 + @test round(Float128(pi), digits = 3, base = 2) ≈ 3.125 + @test round(Float128(pi), digits = 3, base = 5) ≈ 3.144 +end + +if VERSION >= v"1.11" +@testset "rounding floats with specified return type #50778" begin + @test round(Float128, Float128(1.2)) === q1 +end + +# The libmpfr shipped with Julia is built without float128 support, so we can't +# easily do comparisons to conversion and rounding by library calls. +# If this situation changes, consider adapting more tests from Base. + +@testset "floor(<:AbstractFloat, large_number) (#52355)" begin + for i in [-BigInt(floatmax(Float128)), -BigInt(floatmax(Float128))*100, + BigInt(floatmax(Float128)), BigInt(floatmax(Float128))*100] + f = ceil(Float128, i) + @test f >= i + @test isinteger(f) || isinf(f) + @test prevfloat(f) < i + end +end +end diff --git a/test/runtests.jl b/test/runtests.jl index cc85c1b..49c81ac 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -273,5 +273,7 @@ include("specfun.jl") include("printf.jl") +@testset verbose=true "rounding" include("rounding.jl") + using Aqua Aqua.test_all(Quadmath)