Skip to content

Commit f956c23

Browse files
stephentyronestephentyrone
authored
Merge pull request #203 from stephentyrone/rounding
Right shift and division with rounding.
2 parents af53934 + 24d8f51 commit f956c23

File tree

10 files changed

+1155
-7
lines changed

10 files changed

+1155
-7
lines changed

‎Sources/IntegerUtilities/CMakeLists.txt

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,8 +8,11 @@ See https://swift.org/LICENSE.txt for license information
88
#]]
99

1010
add_library(IntegerUtilities
11+
DivideWithRounding.swift
1112
GCD.swift
12-
Rotate.swift)
13+
Rotate.swift
14+
Rounding.swift
15+
ShiftWithRounding.swift)
1316
set_target_properties(IntegerUtilities PROPERTIES
1417
INTERFACE_INCLUDE_DIRECTORIES ${CMAKE_Swift_MODULE_DIRECTORY})
1518

‎Sources/IntegerUtilities/DivideWithRounding.swift

Lines changed: 314 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,314 @@
1+
//===--- Divide.swift -----------------------------------------*- swift -*-===//
2+
//
3+
// This source file is part of the Swift Numerics open source project
4+
//
5+
// Copyright (c) 2021 Apple Inc. and the Swift Numerics project authors
6+
// Licensed under Apache License v2.0 with Runtime Library Exception
7+
//
8+
// See https://swift.org/LICENSE.txt for license information
9+
//
10+
//===----------------------------------------------------------------------===//
11+
12+
extension BinaryInteger {
13+
/// `self` divided by `other`, rounding the result according to `rule`.
14+
///
15+
/// The default rounding rule is `.down`, which _is not the same_ as the
16+
/// behavior of the `/` operator from the Swift standard library, but is
17+
/// chosen because it generally produces a more useful remainder. In
18+
/// particular, when `b` is positive, the remainder is always positive.
19+
/// To match the behavior of `/`, use the `.towardZero` rounding mode.
20+
///
21+
/// Note that the remainder of division is not always representable in an
22+
/// unsigned type if a rounding rule other than `.down`, `.towardZero`, or
23+
/// `.requireExact` is used. For example:
24+
///
25+
/// let a: UInt = 5
26+
/// let b: UInt = 3
27+
/// let q = a.divided(by: b, rounding: .up) // 2
28+
/// let r = a - b*q // 5 - 3*2 overflows UInt.
29+
///
30+
/// For this reason, there is no `remainder(dividingBy:rounding:)`
31+
/// operation defined on `BinaryInteger`. Signed integers do not have
32+
/// this problem, so it is defined on the `SignedInteger` protocol
33+
/// instead, as is an overload of `divided(by:rounding:)` that returns
34+
/// both quotient and remainder.
35+
@inlinable
36+
public func divided(
37+
by other: Self,
38+
rounding rule: RoundingRule = .down
39+
) -> Self {
40+
// "Normal divsion" rounds toward zero, so we get self = q*other + r
41+
// with |r| < |other| and r matching the sign of self.
42+
let q = self / other
43+
let r = self - q*other
44+
// In every rounding mode, the result is the same when the result is
45+
// exact.
46+
if r == 0 { return q }
47+
// From this point forward, we can assume r != 0.
48+
//
49+
// To get the quotient and remainder rounded as directed by rule, we
50+
// will adjust q and r. Note that the quotient is either q-1 or q (if
51+
// q is negative) or q or q+1 (if q is positive), because q has been
52+
// rounded toward zero.
53+
//
54+
// If we subtract 1 from q, we add other to r to compensate, because:
55+
//
56+
// self = q*other + r
57+
// = (q-1)*other + (r+other)
58+
//
59+
// Similarly, if we add 1 to q, we subtract other from r to compensate.
60+
switch rule {
61+
case .down:
62+
// For rounding down, we want to have r match the sign of other
63+
// rather than self; this means that if the signs of r and other
64+
// disagree, we have to adjust q downward and r to match.
65+
if other.signum() != r.signum() { return q-1 }
66+
return q
67+
case .up:
68+
// For rounding up, we want to have r have the opposite sign of
69+
// other; if not, we adjust q upward and r to match.
70+
if other.signum() == r.signum() { return q+1 }
71+
return q
72+
case .towardZero:
73+
// This is exactly what the `/` operator did for us.
74+
return q
75+
case .toOdd:
76+
// If q is already odd, we're done.
77+
if q._lowWord & 1 == 1 { return q }
78+
// Otherwise, q is even but inexact; it was originally rounded toward
79+
// zero, so rounding away from zero instead will make it odd.
80+
fallthrough
81+
case .awayFromZero:
82+
// To round away from zero, we apply the adjustments for both down
83+
// and up.
84+
if other.signum() != r.signum() { return q-1 }
85+
return q+1
86+
case .toNearestOrAwayFromZero:
87+
// For round to nearest or away, the condition we want to satisfy is
88+
// |r| <= |other/2|, with sign(q) != sign(r) when equality holds.
89+
if r.magnitude < other.magnitude.shifted(rightBy: 1, rounding: .up) {
90+
return q
91+
}
92+
// The (q,r) we have does not satisfy the to nearest or away condition;
93+
// round away from zero to choose the other representative of (q, r).
94+
if other.signum() != r.signum() { return q-1 }
95+
return q+1
96+
case .toNearestOrEven:
97+
// For round to nearest or away, the condition we want to satisfy is
98+
// |r| <= |other/2|, with q even when equality holds.
99+
if r.magnitude > other.magnitude.shifted(rightBy: 1, rounding: .down) ||
100+
2*r.magnitude == other.magnitude && q._lowWord & 1 == 1 {
101+
if (other > 0) != (r > 0) { return q-1 }
102+
return q+1
103+
}
104+
return q
105+
case .stochastically:
106+
var qhi: UInt64
107+
var rhi: UInt64
108+
if other.magnitude <= UInt64.max {
109+
qhi = UInt64(other.magnitude)
110+
rhi = UInt64(r.magnitude)
111+
} else {
112+
// TODO: this is untested currently.
113+
let qmag = other.magnitude
114+
let shift = qmag._msb - 1
115+
qhi = UInt64(truncatingIfNeeded: qmag >> shift)
116+
rhi = UInt64(truncatingIfNeeded: r.magnitude >> shift)
117+
}
118+
let (sum, car) = rhi.addingReportingOverflow(.random(in: 0 ..< qhi))
119+
if car || sum >= qhi {
120+
if (other > 0) != (r > 0) { return q-1 }
121+
return q+1
122+
}
123+
return q
124+
case .requireExact:
125+
preconditionFailure("Division was not exact.")
126+
}
127+
}
128+
129+
// TODO: make this API and make it possible to implement more
130+
// efficiently. Customization point on new/revised integer
131+
// protocol? Shouldn't have to go through .words.
132+
/// The index of the most-significant set bit.
133+
///
134+
/// - Precondition: self is assumed to be non-zero (to be changed
135+
/// if/when this becomes API).
136+
@usableFromInline
137+
internal var _msb: Int {
138+
// a == 0 is never used for division, because this is called
139+
// on the divisor which cannot be zero as a precondition; if
140+
// this becomes API, the behavior for this case will have to
141+
// be defined.
142+
assert(self != 0)
143+
// Because self is non-zero, mswIndex is guaranteed to exist,
144+
// hence force-unwrap is appropriate.
145+
let mswIndex = words.lastIndex { $0 != 0 }!
146+
let mswBits = UInt.bitWidth * words.distance(from: words.startIndex, to: mswIndex)
147+
return mswBits + (UInt.bitWidth - words[mswIndex].leadingZeroBitCount - 1)
148+
}
149+
}
150+
151+
extension SignedInteger {
152+
/// Divides `self` by `other`, rounding the quotient according to `rule`,
153+
/// and returns the remainder.
154+
///
155+
/// The default rounding rule is `.down`, which _is not the same_ as the
156+
/// behavior of the `%` operator from the Swift standard library, but is
157+
/// chosen because it generally produces a more useful remainder. To
158+
/// match the behavior of `%`, use the `.towardZero` rounding mode.
159+
///
160+
/// - Precondition: `other` cannot be zero.
161+
@inlinable
162+
public func remainder(
163+
dividingBy other: Self,
164+
rounding rule: RoundingRule = .down
165+
) -> Self {
166+
// Produce correct remainder for the .min/-1 case, rather than trapping.
167+
if other == -1 { return 0 }
168+
return self.divided(by: other, rounding: rule).remainder
169+
}
170+
171+
/// Divides `self` by `other`, rounding the quotient according to `rule`,
172+
/// and returns both the quotient and remainder.
173+
///
174+
/// The default rounding rule is `.down`, which _is not the same_ as the
175+
/// behavior of the `/` operator from the Swift standard library, but is
176+
/// chosen because it generally produces a more useful remainder. To
177+
/// match the behavior of `/`, use the `.towardZero` rounding mode.
178+
///
179+
/// Because the default rounding mode does not match Swift's standard
180+
/// library, this function is a disfavored overload of `divided(by:)`
181+
/// instead of using the name `quotientAndRemainder(dividingBy:)`, which
182+
/// would shadow the standard library operation and change the behavior
183+
/// of any existing use sites. To call this method, you must explicitly
184+
/// bind the result to a tuple:
185+
///
186+
/// // This calls BinaryInteger's method, which returns only
187+
/// // the quotient.
188+
/// let result = 5.divided(by: 3, rounding: .up) // 2
189+
///
190+
/// // This calls SignedInteger's method, which returns both
191+
/// // the quotient and remainder.
192+
/// let (q, r) = 5.divided(by: 3, rounding: .up) // (q = 2, r = -1)
193+
@inlinable @inline(__always) @_disfavoredOverload
194+
public func divided(
195+
by other: Self,
196+
rounding rule: RoundingRule = .down
197+
) -> (quotient: Self, remainder: Self) {
198+
// "Normal divsion" rounds toward zero, so we get self = q*other + r
199+
// with |r| < |other| and r matching the sign of self.
200+
let q = self / other
201+
let r = self - q*other
202+
// In every rounding mode, the result is the same when the result is
203+
// exact.
204+
if r == 0 { return (q, r) }
205+
// From this point forward, we can assume r != 0.
206+
//
207+
// To get the quotient and remainder rounded as directed by rule, we
208+
// will adjust q and r. Note that the quotient is either q-1 or q (if
209+
// q is negative) or q or q+1 (if q is positive), because q has been
210+
// rounded toward zero.
211+
//
212+
// If we subtract 1 from q, we add other to r to compensate, because:
213+
//
214+
// self = q*other + r
215+
// = (q-1)*other + (r+other)
216+
//
217+
// Similarly, if we add 1 to q, we subtract other from r to compensate.
218+
switch rule {
219+
case .down:
220+
// For rounding down, we want to have r match the sign of other
221+
// rather than self; this means that if the signs of r and other
222+
// disagree, we have to adjust q downward and r to match.
223+
if other.signum() != r.signum() { return (q-1, r+other) }
224+
return (q, r)
225+
case .up:
226+
// For rounding up, we want to have r have the opposite sign of
227+
// other; if not, we adjust q upward and r to match.
228+
if other.signum() == r.signum() { return (q+1, r-other) }
229+
return (q, r)
230+
case .towardZero:
231+
// This is exactly what the `/` operator did for us.
232+
return (q, r)
233+
case .toOdd:
234+
// If q is already odd, we're done.
235+
if q._lowWord & 1 == 1 { return (q, r) }
236+
// Otherwise, q is even but inexact; it was originally rounded toward
237+
// zero, so rounding away from zero instead will make it odd.
238+
fallthrough
239+
case .awayFromZero:
240+
// To round away from zero, we apply the adjustments for both down
241+
// and up.
242+
if other.signum() != r.signum() { return (q-1, r+other) }
243+
return (q+1, r-other)
244+
case .toNearestOrAwayFromZero:
245+
// For round to nearest or away, the condition we want to satisfy is
246+
// |r| <= |other/2|, with sign(q) != sign(r) when equality holds.
247+
if r.magnitude < other.magnitude.shifted(rightBy: 1, rounding: .up) {
248+
return (q, r)
249+
}
250+
// The (q,r) we have does not satisfy the to nearest or away condition;
251+
// round away from zero to choose the other representative of (q, r).
252+
if other.signum() != r.signum() { return (q-1, r+other) }
253+
return (q+1, r-other)
254+
case .toNearestOrEven:
255+
// For round to nearest or away, the condition we want to satisfy is
256+
// |r| <= |other/2|, with q even when equality holds.
257+
if r.magnitude > other.magnitude.shifted(rightBy: 1, rounding: .down) ||
258+
2*r.magnitude == other.magnitude && q._lowWord & 1 == 1 {
259+
if (other > 0) != (r > 0) { return (q-1, r+other) }
260+
return (q+1, r-other)
261+
}
262+
return (q, r)
263+
case .stochastically:
264+
var qhi: UInt64
265+
var rhi: UInt64
266+
if other.magnitude <= UInt64.max {
267+
qhi = UInt64(other.magnitude)
268+
rhi = UInt64(r.magnitude)
269+
} else {
270+
// TODO: this is untested currently.
271+
let qmag = other.magnitude
272+
let shift = qmag._msb - 1
273+
qhi = UInt64(truncatingIfNeeded: qmag >> shift)
274+
rhi = UInt64(truncatingIfNeeded: r.magnitude >> shift)
275+
}
276+
let (sum, car) = rhi.addingReportingOverflow(.random(in: 0 ..< qhi))
277+
if car || sum >= qhi {
278+
if (other > 0) != (r > 0) { return (q-1, r+other) }
279+
return (q+1, r-other)
280+
}
281+
return (q, r)
282+
case .requireExact:
283+
preconditionFailure("Division was not exact.")
284+
}
285+
}
286+
}
287+
288+
/// `a = quotient*b + remainder`, with `remainder >= 0`.
289+
///
290+
/// When `a` and `b` are both positive, `quotient` is `a/b` and `remainder`
291+
/// is `a%b`.
292+
///
293+
/// Rounding the quotient so that the remainder is non-negative is called
294+
/// "Euclidean division". This is not a _rounding rule_, as `quotient`
295+
/// cannot be determined from the unrounded value `a/b`; we need to also
296+
/// know the sign of `a` or `b` or `r` to know which way to round. Because
297+
/// of this, is not present in the `RoundingRule` enum and uses a separate
298+
/// API from the other division operations.
299+
///
300+
/// - Parameters:
301+
/// - a: The dividend
302+
/// - b: The divisor
303+
///
304+
/// - Precondition: `b` must be non-zero, and the quotient `a/b` must be
305+
/// representable. In particular, if `T` is a signed fixed-width integer
306+
/// type, then `euclideanDivision(T.min, -1)` will trap, because `-T.min`
307+
/// is not representable.
308+
///
309+
/// - Returns: `(quotient, remainder)`, with `0 <= remainder < b.magnitude`.
310+
func euclideanDivision<T>(_ a: T, _ b: T) -> (quotient: T, remainder: T)
311+
where T: SignedInteger
312+
{
313+
a.divided(by: b, rounding: a >= 0 ? .towardZero : .awayFromZero)
314+
}

‎Sources/IntegerUtilities/README.md

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,9 +8,45 @@ The following API are defined for all integer types conforming to `BinaryInteger
88

99
- The `gcd(_:_:)` free function implements the _Greatest Common Divisor_
1010
operation.
11+
12+
- The `shifted(rightBy:rounding:)` method implements _bitwise shift with
13+
rounding_.
14+
15+
- The `divided(by:rounding:)` method implements division with specified
16+
rounding. (See also `SignedInteger.divided(by:rounding:)`,
17+
`remainder(dividingBy:rounding:)`, and `euclideanDivision(_:_:)` below).
18+
19+
## Utilities defined for `SignedInteger`
20+
21+
The following API are defined for signed integer types:
22+
23+
- The `divided(by:rounding:)` method implementing division with specified
24+
rounding, returning both quotient and remainder. This requires a signed
25+
type because the remainder is not generally representable for unsigned
26+
types. This is a disfavored overload; by default, you will get only the
27+
quotient as the result:
28+
```
29+
let p = 5.divided(by: 3, rounding: .up) // p = 2
30+
let (q, r) = 5.divided(by: 3, rounding: .up) // q = 2, r = -1
31+
```
32+
33+
- The `remainder(dividingBy:rounding:)` method implementing the remainder
34+
operation; the `rounding` argument describes how to round the _quotient_,
35+
which is not returned. (The remainder is always exact, and hence is not
36+
rounded).
37+
38+
- The `euclideanDivision(_:_:)` free function implements _Euclidean division_.
39+
In this operation, the remainder is chosen to always be non-negative. This
40+
does not correspond to any rounding rule on the quotient, which is why it
41+
uses a distinct API.
1142

1243
## Utilities defined for `FixedWidthInteger`
1344

1445
- The `rotated(right:)` and `rotated(left:)` methods implement _bitwise
1546
rotation_ for signed and unsigned integer types. The count parameter may
1647
be any `BinaryInteger` type.
48+
49+
## Types
50+
51+
The `RoundingRule` enum is used with shift, division, and round operations
52+
to specify how to round their results to a representable value.

0 commit comments

Comments
 (0)