Skip to content

Commit 4cffd1a

Browse files
JeffreySarnoffJeffreySarnoff
authored
Merge pull request #102 from RalphAS/ras/roundfix1
Fix some rounding and conversion problems
2 parents ccb9337 + 8ec62a2 commit 4cffd1a

3 files changed

Lines changed: 319 additions & 3 deletions

File tree

‎src/Quadmath.jl‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -226,14 +226,15 @@ for Ti in (UInt8, UInt16, UInt32, UInt64, UInt128)
226226
# `Float128(typemax(Ti))+1` is exactly representable
227227
@eval begin
228228
function trunc(::Type{$Ti},x::Float128)
229-
if !signbit(x) && x < Float128(typemax($Ti)) + one(Float128)
229+
if (!signbit(x) && x < Float128(typemax($Ti)) + one(Float128)) || (x == 0)
230230
return unsafe_trunc($Ti,x)
231231
else
232232
throw(InexactError(:trunc, $Ti, x))
233233
end
234234
end
235235
function (::Type{$Ti})(x::Float128)
236-
if (!signbit(x) && x <= Float128(typemax($Ti))) && (round(x, RoundToZero) == x)
236+
if ((!signbit(x) && x <= Float128(typemax($Ti))) && (round(x, RoundToZero) == x)
237+
|| (x == 0))
237238
return unsafe_trunc($Ti,x)
238239
else
239240
throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x))
@@ -303,7 +304,8 @@ end
303304
@assume_effects :foldable round(x::Float128) = Float128(@quad_ccall(libquadmath.rintq(x::Cfloat128)::Cfloat128))
304305
round(x::Float128, r::RoundingMode{:Down}) = floor(x)
305306
round(x::Float128, r::RoundingMode{:Up}) = ceil(x)
306-
round(x::Float128, r::RoundingMode{:ToZero}) = round(x)
307+
round(x::Float128, r::RoundingMode{:Nearest}) = round(x)
308+
round(x::Float128, r::RoundingMode{:ToZero}) = trunc(x)
307309

308310
## two argument
309311
@assume_effects :foldable (^)(x::Float128, y::Float128) =
@@ -498,8 +500,10 @@ function Float128(x::BigFloat)
498500
prec = 113 + (k + 16381)
499501
elseif k == -16381-113 && abs(y) > 0.5
500502
z = reinterpret(Float128, UInt128(1))
503+
return copysign(z,x)
501504
else
502505
z = reinterpret(Float128, UInt128(0))
506+
return copysign(z,x)
503507
end
504508

505509
y = BigFloat(x, precision=prec)

‎test/rounding.jl‎

Lines changed: 310 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,310 @@
1+
# Derived from Julia Base test suite. License is MIT: https://julialang.org/license
2+
3+
# Small sanity tests to ensure changing the rounding of float functions work
4+
using Base.MathConstants
5+
6+
q1 = one(Float128)
7+
q0 = zero(Float128)
8+
9+
@testset "Arithmetic rounding checks" begin
10+
# a + b returns a number exactly between prevfloat(1.) and 1., so its
11+
# final result depends strongly on the utilized rounding direction.
12+
a = prevfloat(Float128(0.5))
13+
b = Float128(0.5)
14+
c = eps(Float128)/4
15+
d = prevfloat(Float128(1))
16+
17+
@testset "Default rounding direction, RoundNearest" begin
18+
@test a + b === q1
19+
@test - a - b === -q1
20+
@test a - b === -c
21+
@test b - a === c
22+
end
23+
end
24+
25+
@testset "convert with rounding" begin
26+
for v = [sqrt(Float128(2)),-1/Float128(3),nextfloat(q1),prevfloat(q1),nextfloat(-q1),
27+
prevfloat(-q1),nextfloat(q0),prevfloat(q0)]
28+
29+
pn = Float64(v,RoundNearest)
30+
@test pn == convert(Float64,v)
31+
32+
pz = Float64(v,RoundToZero)
33+
@test abs(pz) <= abs(v) < nextfloat(abs(pz))
34+
@test signbit(pz) == signbit(v)
35+
36+
pd = Float64(v,RoundDown)
37+
@test pd <= v < nextfloat(pd)
38+
39+
pu = Float64(v,RoundUp)
40+
@test prevfloat(pu) < v <= pu
41+
42+
@test pn == pd || pn == pu
43+
@test v > 0 ? pz == pd : pz == pu
44+
@test pu - pd == eps(pz)
45+
end
46+
47+
# numbers which are not represented exactly in Float128
48+
for T in [Float128,]
49+
for v in [sqrt(big(2.0)),-big(1.0)/big(3.0),nextfloat(big(1.0)),
50+
prevfloat(big(1.0)), nextfloat(big(0.0)),prevfloat(big(0.0)),
51+
pi,ℯ,eulergamma,catalan,golden,
52+
typemax(Int128),typemax(UInt128),
53+
]
54+
pn = T(v,RoundNearest)
55+
@test pn == convert(T,BigFloat(v))
56+
pz = T(v,RoundToZero)
57+
@test pz == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundToZero)
58+
pd = T(v,RoundDown)
59+
if v == prevfloat(big(0.0))
60+
@test_broken pd == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundDown)
61+
else
62+
@test pd == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundDown)
63+
end
64+
pu = T(v,RoundUp)
65+
if v == nextfloat(big(0.0))
66+
@test_broken pu == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundUp)
67+
else
68+
@test pu == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundUp)
69+
end
70+
@test pn == pd || pn == pu
71+
@test v > 0 ? pz == pd : pz == pu
72+
@test isinf(pu) || pu - pd == eps(pz)
73+
end
74+
end
75+
end
76+
77+
@testset "rounding properties" for Tf in [Float128,]
78+
# these should hold for all u, but we just test the smallest and largest
79+
# of each binade
80+
81+
for i in exponent(floatmin(Tf)):exponent(floatmax(Tf))
82+
for u in [ldexp(Tf(1.0), i), -ldexp(Tf(1.0), i),
83+
ldexp(prevfloat(Tf(2.0)), i), -ldexp(prevfloat(Tf(2.0)), i)]
84+
85+
r = round(u, RoundNearest)
86+
if isfinite(u)
87+
@test isfinite(r)
88+
@test isinteger(r)
89+
@test abs(r-u) < 0.5 || abs(r-u) == 0.5 && isinteger(r/2)
90+
@test signbit(u) == signbit(r)
91+
else
92+
@test u === r
93+
end
94+
95+
r = round(u, RoundNearestTiesAway)
96+
if isfinite(u)
97+
@test isfinite(r)
98+
@test isinteger(r)
99+
@test abs(r-u) < 0.5 || (r-u) == copysign(0.5,u)
100+
@test signbit(u) == signbit(r)
101+
else
102+
@test u === r
103+
end
104+
105+
r = round(u, RoundNearestTiesUp)
106+
if isfinite(u)
107+
@test isfinite(r)
108+
@test isinteger(r)
109+
@test -0.5 < r-u <= 0.5
110+
@test signbit(u) == signbit(r)
111+
else
112+
@test u === r
113+
end
114+
115+
r = round(u, RoundFromZero)
116+
if isfinite(u)
117+
@test isfinite(r)
118+
@test isinteger(r)
119+
@test signbit(u) ? (r == floor(u)) : (r == ceil(u))
120+
@test signbit(u) == signbit(r)
121+
else
122+
@test u === r
123+
end
124+
end
125+
end
126+
end
127+
128+
@testset "rounding difficult values" begin
129+
for x = Int128(2)^113-10:Int128(2)^113+10
130+
y = Float128(x)
131+
i = trunc(Int128,y)
132+
@test Int128(trunc(y)) == i
133+
@test Int128(round(y)) == i
134+
@test Int128(floor(y)) == i
135+
@test Int128(ceil(y)) == i
136+
137+
@test round(Int128,y) == i
138+
@test floor(Int128,y) == i
139+
@test ceil(Int128,y) == i
140+
end
141+
142+
# rounding vectors
143+
let ≈(x,y) = x==y && typeof(x)==typeof(y)
144+
for t in [Float128,]
145+
# try different vector lengths
146+
for n in [0,3,255,256]
147+
r = (1:n) .- div(n,2)
148+
y = t[x/4 for x in r]
149+
@test trunc.(y) ≈ t[div(i,4) for i in r]
150+
@test floor.(y) ≈ t[i>>2 for i in r]
151+
@test ceil.(y) ≈ t[(i+3)>>2 for i in r]
152+
@test round.(y) ≈ t[(i+1+isodd(i>>2))>>2 for i in r]
153+
@test broadcast(x -> round(x, RoundNearestTiesAway), y) ≈ t[(i+1+(i>=0))>>2 for i in r]
154+
@test broadcast(x -> round(x, RoundNearestTiesUp), y) ≈ t[(i+2)>>2 for i in r]
155+
@test broadcast(x -> round(x, RoundFromZero), y) ≈ t[(i+3*(i>=0))>>2 for i in r]
156+
end
157+
end
158+
end
159+
160+
@test_throws InexactError round(Int,Float128(Inf))
161+
@test_throws InexactError round(Int,Float128(NaN))
162+
@test round(Int,Float128(2.5)) == 2
163+
@test round(Int,Float128(1.5)) == 2
164+
@test round(Int,Float128(-2.5)) == -2
165+
@test round(Int,Float128(-1.5)) == -2
166+
@test round(Int,Float128(2.5),RoundNearestTiesAway) == 3
167+
@test round(Int,Float128(1.5),RoundNearestTiesAway) == 2
168+
@test round(Int,Float128(2.5),RoundNearestTiesUp) == 3
169+
@test round(Int,Float128(1.5),RoundNearestTiesUp) == 2
170+
@test round(Int,Float128(-2.5),RoundNearestTiesAway) == -3
171+
@test round(Int,Float128(-1.5),RoundNearestTiesAway) == -2
172+
@test round(Int,Float128(-2.5),RoundNearestTiesUp) == -2
173+
@test round(Int,Float128(-1.5),RoundNearestTiesUp) == -1
174+
@test round(Int,Float128(-1.9)) == -2
175+
@test round(Int,nextfloat(q1),RoundFromZero) == 2
176+
@test round(Int,-nextfloat(q1),RoundFromZero) == -2
177+
@test round(Int,prevfloat(q1),RoundFromZero) == 1
178+
@test round(Int,-prevfloat(q1),RoundFromZero) == -1
179+
#=
180+
FIXME: want analogs
181+
@test_throws InexactError round(Int64, 9.223372036854776e18)
182+
@test round(Int64, 9.223372036854775e18) == 9223372036854774784
183+
@test_throws InexactError round(Int64, -9.223372036854778e18)
184+
@test round(Int64, -9.223372036854776e18) == typemin(Int64)
185+
@test_throws InexactError round(UInt64, 1.8446744073709552e19)
186+
@test round(UInt64, 1.844674407370955e19) == 0xfffffffffffff800
187+
=#
188+
for Ti in [Int,UInt]
189+
for Tf in [Float128,]
190+
# Note: these threw InexactError
191+
@test round(Ti,Tf(-0.0)) == 0
192+
@test round(Ti,Tf(-0.0),RoundNearestTiesAway) == 0
193+
@test round(Ti,Tf(-0.0),RoundNearestTiesUp) == 0
194+
195+
@test round(Ti, Tf(0.5)) == 0
196+
@test round(Ti, Tf(0.5), RoundNearestTiesAway) == 1
197+
@test round(Ti, Tf(0.5), RoundNearestTiesUp) == 1
198+
199+
@test round(Ti, prevfloat(Tf(0.5))) == 0
200+
@test round(Ti, prevfloat(Tf(0.5)), RoundNearestTiesAway) == 0
201+
@test round(Ti, prevfloat(Tf(0.5)), RoundNearestTiesUp) == 0
202+
203+
@test round(Ti, nextfloat(Tf(0.5))) == 1
204+
@test round(Ti, nextfloat(Tf(0.5)), RoundNearestTiesAway) == 1
205+
@test round(Ti, nextfloat(Tf(0.5)), RoundNearestTiesUp) == 1
206+
# Note: these threw InexactError
207+
@test round(Ti, Tf(-0.5)) == 0
208+
@test round(Ti, Tf(-0.5), RoundNearestTiesUp) == 0
209+
210+
@test round(Ti, nextfloat(Tf(-0.5))) == 0
211+
@test round(Ti, nextfloat(Tf(-0.5)), RoundNearestTiesAway) == 0
212+
@test round(Ti, nextfloat(Tf(-0.5)), RoundNearestTiesUp) == 0
213+
214+
if Ti <: Signed
215+
@test round(Ti, Tf(-0.5), RoundNearestTiesAway) == -1
216+
@test round(Ti, prevfloat(Tf(-0.5))) == -1
217+
@test round(Ti, prevfloat(Tf(-0.5)), RoundNearestTiesAway) == -1
218+
@test round(Ti, prevfloat(Tf(-0.5)), RoundNearestTiesUp) == -1
219+
else
220+
@test_throws InexactError round(Ti, Tf(-0.5), RoundNearestTiesAway)
221+
@test_throws InexactError round(Ti, prevfloat(Tf(-0.5)))
222+
@test_throws InexactError round(Ti, prevfloat(Tf(-0.5)), RoundNearestTiesAway)
223+
@test_throws InexactError round(Ti, prevfloat(Tf(-0.5)), RoundNearestTiesUp)
224+
end
225+
end
226+
end
227+
228+
# numbers that can't be rounded by trunc(x+0.5)
229+
@test round(Int128, Float128(2.0)^112 + 1) == (UInt128(1) << 112) + 1
230+
end
231+
232+
# custom rounding and significant-digit ops
233+
@testset "rounding to digits relative to the decimal point" begin
234+
@test round(Float128(pi)) ≈ 3.
235+
@test round(Float128(pi), base=10) ≈ 3.
236+
@test round(Float128(pi), digits=0) ≈ 3.
237+
@test round(Float128(pi), digits=1) ≈ 3.1
238+
@test round(Float128(pi), digits=3, base=2) ≈ 3.125
239+
@test round(Float128(pi), sigdigits=1) ≈ 3.
240+
@test round(Float128(pi), sigdigits=3) ≈ 3.14
241+
@test round(Float128(pi), sigdigits=4, base=2) ≈ 3.25
242+
@test round(10*Float128(pi), digits=-1) ≈ 30.
243+
@test round(Float128(.1), digits=0) == 0.
244+
@test round(Float128(-.1), digits=0) == -0.
245+
@test isnan(round(Float128(NaN), digits=2))
246+
@test isinf(round(Float128(Inf), digits=2))
247+
@test isinf(round(Float128(-Inf), digits=2))
248+
end
249+
@testset "round vs trunc vs floor vs ceil" begin
250+
@test round(Float128(123.456), digits=1) ≈ 123.5
251+
@test round(-Float128(123.456), digits=1) ≈ -123.5
252+
@test trunc(Float128(123.456), digits=1) ≈ 123.4
253+
@test trunc(-Float128(123.456), digits=1) ≈ -123.4
254+
@test ceil(Float128(123.456), digits=1) ≈ 123.5
255+
@test ceil(-Float128(123.456), digits=1) ≈ -123.4
256+
@test floor(Float128(123.456), digits=1) ≈ 123.4
257+
@test floor(-Float128(123.456), digits=1) ≈ -123.5
258+
end
259+
@testset "rounding with too much (or too few) precision" begin
260+
for x in (12345.6789, 0, -12345.6789)
261+
y = Float128(x)
262+
@test y == trunc(x, digits=1000)
263+
@test y == round(x, digits=1000)
264+
@test y == floor(x, digits=1000)
265+
@test y == ceil(x, digits=1000)
266+
end
267+
let x = Float128(12345.6789)
268+
@test 0.0 == trunc(x, digits=-1000)
269+
@test 0.0 == round(x, digits=-1000)
270+
@test 0.0 == floor(x, digits=-1000)
271+
@test Inf128 == ceil(x, digits=-10000)
272+
end
273+
let x = Float128(-12345.6789)
274+
@test -0.0 == trunc(x, digits=-1000)
275+
@test -0.0 == round(x, digits=-1000)
276+
@test -Inf128 == floor(x, digits=-10000)
277+
@test -0.0 == ceil(x, digits=-1000)
278+
end
279+
let x = q0
280+
@test 0.0 == trunc(x, digits=-1000)
281+
@test 0.0 == round(x, digits=-1000)
282+
@test 0.0 == floor(x, digits=-1000)
283+
@test 0.0 == ceil(x, digits=-1000)
284+
end
285+
end
286+
@testset "rounding in other bases" begin
287+
@test round(Float128(pi), digits = 2, base = 2) ≈ 3.25
288+
@test round(Float128(pi), digits = 3, base = 2) ≈ 3.125
289+
@test round(Float128(pi), digits = 3, base = 5) ≈ 3.144
290+
end
291+
292+
if VERSION >= v"1.11"
293+
@testset "rounding floats with specified return type #50778" begin
294+
@test round(Float128, Float128(1.2)) === q1
295+
end
296+
297+
# The libmpfr shipped with Julia is built without float128 support, so we can't
298+
# easily do comparisons to conversion and rounding by library calls.
299+
# If this situation changes, consider adapting more tests from Base.
300+
301+
@testset "floor(<:AbstractFloat, large_number) (#52355)" begin
302+
for i in [-BigInt(floatmax(Float128)), -BigInt(floatmax(Float128))*100,
303+
BigInt(floatmax(Float128)), BigInt(floatmax(Float128))*100]
304+
f = ceil(Float128, i)
305+
@test f >= i
306+
@test isinteger(f) || isinf(f)
307+
@test prevfloat(f) < i
308+
end
309+
end
310+
end

‎test/runtests.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -273,5 +273,7 @@ include("specfun.jl")
273273

274274
include("printf.jl")
275275

276+
@testset verbose=true "rounding" include("rounding.jl")
277+
276278
using Aqua
277279
Aqua.test_all(Quadmath)

0 commit comments

Comments
 (0)