Skip to content

Commit

Permalink
add various divrem combinations to avoid overflow
Browse files Browse the repository at this point in the history
  • Loading branch information
simonbyrne committed Aug 25, 2019
1 parent 1a374a7 commit eb42c70
Show file tree
Hide file tree
Showing 2 changed files with 59 additions and 25 deletions.
76 changes: 51 additions & 25 deletions base/div.jl
Expand Up @@ -113,6 +113,56 @@ julia> divrem(7,3)
```
"""
divrem(x, y) = divrem(x, y, RoundToZero)
function divrem(x::Integer, y::Integer, rnd::typeof(RoundNearest))
(q, r) = divrem(x, y)
if x >= 0
if y >= 0
r >= (y÷2) + (isodd(y) | iseven(q)) ? (q+true, r-y) : (q, r)
else
r >= -(y÷2) + (isodd(y) | iseven(q)) ? (q-true, r+y) : (q, r)
end
else
if y >= 0
r <= -signed(y÷2) - (isodd(y)| iseven(q)) ? (q-true, r+y) : (q, r)
else
r <= (y÷2) - (isodd(y) | iseven(q)) ? (q+true, r-y) : (q, r)
end
end
end
function divrem(x::Integer, y::Integer, rnd:: typeof(RoundNearestTiesAway))
(q, r) = divrem(x, y)
if x >= 0
if y >= 0
r >= (y÷2) + isodd(y) ? (q+true, r-y) : (q, r)
else
r >= -(y÷2) + isodd(y) ? (q-true, r+y) : (q, r)
end
else
if y >= 0
r <= -signed(y÷2) - isodd(y) ? (q-true, r+y) : (q, r)
else
r <= (y÷2) - isodd(y) ? (q+true, r-y) : (q, r)
end
end
end
function divrem(x::Integer, y::Integer, rnd::typeof(RoundNearestTiesUp))
(q, r) = divrem(x, y)
if x >= 0
if y >= 0
r >= (y÷2) + isodd(y) ? (q+true, r-y) : (q, r)
else
r >= -(y÷2) + true ? (q-true, r+y) : (q, r)
end
else
if y >= 0
r <= -signed(y÷2) - true ? (q-true, r+y) : (q, r)
else
r <= (y÷2) - isodd(y) ? (q+true, r-y) : (q, r)
end
end
end


divrem(a, b, r::RoundingMode) = (div(a, b, r), rem(a, b, r))

"""
Expand Down Expand Up @@ -146,31 +196,7 @@ end
function div(x::Integer, y::Integer, rnd::Union{typeof(RoundNearest),
typeof(RoundNearestTiesAway),
typeof(RoundNearestTiesUp)})
(q, r) = divrem(x, y)
# No remainder, we're done
iszero(r) && return q
if isodd(y)
# The divisior is odd - no ties are possible, just
# round to nearest. N.B. 2r == y is impossible because
# y is odd.
return abs(2r) > abs(y) ? q + copysign(one(q), q) : q
end
# y is even, divide y by two and check whether we have a tie.
# If not, as above we just round to the nearest value.
halfy = copysign(y >> 1, r)
c = cmp(r, halfy)
c == -1 && return q
c == 1 && return q + copysign(one(q), q)
# We have a tie (r == y/2). Select tie behavior according to
# rounding mode.
if rnd == RoundNearest
return iseven(q) ? q : q + copysign(one(q), q)
elseif rnd == RoundNearestTiesAway
return q + copysign(one(q), q)
else
@assert rnd == RoundNearestTiesUp
return sign(q) == -1 ? q : q + one(q)
end
divrem(x,y,rnd)[1]
end

# For bootstrapping purposes, we define div for integers directly. Provide the
Expand Down
8 changes: 8 additions & 0 deletions test/int.jl
Expand Up @@ -307,6 +307,14 @@ let i = MyInt26779(1)
end

@testset "rounding division" begin
for x = -100:100
for y = 1:100
for rnd in (RoundNearest, RoundNearestTiesAway, RoundNearestTiesUp)
@test div(x,y,rnd) == round(x/y,rnd)
@test div(x,-y,rnd) == round(x/-y,rnd)
end
end
end
for (a, b, nearest, away, up) in (
(3, 2, 2, 2, 2),
(5, 3, 2, 2, 2),
Expand Down

0 comments on commit eb42c70

Please sign in to comment.