Skip to content

Commit c74d817

Browse files
stephentyronestephentyrone
authored
Merge pull request #208 from markuswntr/two-sum
Add `Augmented.sum(_:_:)` aka "twoSum"
2 parents f956c23 + 05366d7 commit c74d817

File tree

2 files changed

+210
-2
lines changed

2 files changed

+210
-2
lines changed

‎Sources/RealModule/AugmentedArithmetic.swift

Lines changed: 50 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
//
33
// This source file is part of the Swift Numerics open source project
44
//
5-
// Copyright (c) 2020 Apple Inc. and the Swift Numerics project authors
5+
// Copyright (c) 2020-2021 Apple Inc. and the Swift Numerics project authors
66
// Licensed under Apache License v2.0 with Runtime Library Exception
77
//
88
// See https://swift.org/LICENSE.txt for license information
@@ -64,7 +64,7 @@ extension Augmented {
6464
/// never underflow. However, it may not be exactly representable when
6565
/// `a` and `b` differ widely in magnitude.
6666
///
67-
/// This operation is sometimes called "fastTwoSum".
67+
/// This operation is sometimes called ["fastTwoSum"].
6868
///
6969
/// - Parameters:
7070
/// - a: The summand with larger magnitude.
@@ -86,11 +86,59 @@ extension Augmented {
8686
/// -
8787
/// - If `head` is normal, then `abs(tail) < head.ulp`.
8888
/// Assuming IEEE 754 default rounding, `abs(tail) <= head.ulp/2`.
89+
///
90+
/// ["fastTwoSum"]: https://en.wikipedia.org/wiki/2Sum
8991
@_transparent
9092
public static func sum<T:Real>(large a: T, small b: T) -> (head: T, tail: T) {
9193
assert(!(b.magnitude > a.magnitude))
9294
let head = a + b
9395
let tail = a - head + b
9496
return (head, tail)
9597
}
98+
99+
/// The sum `a + b` represented as an implicit sum `head + tail`.
100+
///
101+
/// `head` is the correctly rounded value of `a + b`. `tail` is the
102+
/// error from that computation rounded to the closest representable
103+
/// value.
104+
///
105+
/// Unlike `Augmented.sum(large: a, small: b)`, the magnitude of the summands
106+
/// does not matter and `a.magnitude` might as well be strictly less than
107+
/// `b.magnitude`. However, it is recommended to only use this function over
108+
/// `Augmented.sum(large: a, small: b)` in cases where the ordering of the
109+
/// summands magnitude is unknown at compile time. In cases where either of
110+
/// the summands magnitude is guaranteed to be greater than or equal the
111+
/// magnitude of the other summand, use `Augmented.sum(large: a, small: b)`
112+
/// over this function; as it faster to calculate.
113+
///
114+
/// Unlike `Augmented.product(a, b)`, the rounding error of a sum can
115+
/// never underflow. However, it may not be exactly representable when
116+
/// `a` and `b` differ widely in magnitude.
117+
///
118+
/// This operation is sometimes called ["twoSum"].
119+
///
120+
/// - Parameters:
121+
/// - a: One of the summand
122+
/// - b: The other summand
123+
///
124+
/// Edge Cases:
125+
/// -
126+
/// - `head` is always the IEEE 754 sum `a + b`.
127+
/// - If `head` is not finite, `tail` is unspecified and should not be
128+
/// interpreted as having any meaning (it may be `NaN` or `infinity`).
129+
///
130+
/// Postconditions:
131+
/// -
132+
/// - If `head` is normal, then `abs(tail) < head.ulp`.
133+
/// Assuming IEEE 754 default rounding, `abs(tail) <= head.ulp/2`.
134+
///
135+
/// ["twoSum"]: https://en.wikipedia.org/wiki/2Sum
136+
@_transparent
137+
public static func sum<T: Real>(_ a: T, _ b: T) -> (head: T, tail: T) {
138+
let head = a + b
139+
let x = head - b
140+
let y = head - x
141+
let tail = (a - x) + (b - y)
142+
return (head, tail)
143+
}
96144
}

‎Tests/RealTests/AugmentedArithmeticTests.swift

Lines changed: 160 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,160 @@
1+
//===--- AugmentedArithmeticTests.swift -----------------------*- swift -*-===//
2+
//
3+
// This source file is part of the Swift.org open source project
4+
//
5+
// Copyright (c) 2021 Apple Inc. and the Swift 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+
// See https://swift.org/CONTRIBUTORS.txt for the list of Swift project authors
10+
//
11+
//===----------------------------------------------------------------------===//
12+
13+
import RealModule
14+
import XCTest
15+
import _TestSupport
16+
17+
final class AugmentedArithmeticTests: XCTestCase {
18+
19+
func testTwoSumSpecials<T: Real & FixedWidthFloatingPoint>(_: T.Type) {
20+
// Must be exact and not overflow on outer bounds
21+
var x = T.greatestFiniteMagnitude
22+
var y = -T.greatestFiniteMagnitude
23+
XCTAssertEqual(Augmented.sum(x, y).head, .zero)
24+
XCTAssertEqual(Augmented.sum(x, y).tail, .zero)
25+
XCTAssert(Augmented.sum(x, y).head.isFinite)
26+
XCTAssert(Augmented.sum(x, y).tail.isFinite)
27+
// Must be exact on lower subnormal bounds
28+
x = T.leastNonzeroMagnitude
29+
y = -T.leastNonzeroMagnitude
30+
XCTAssertEqual(Augmented.sum(x, y).head, .zero)
31+
XCTAssertEqual(Augmented.sum(x, y).tail, .zero)
32+
// Must preserve floating point signs for:
33+
// (1) (+0) + (-0) == +0
34+
// (2) (-0) + (+0) == +0
35+
// (3) (-0) + (-0) == -0
36+
x = T(sign: .plus, exponent: 1, significand: 0)
37+
y = T(sign: .minus, exponent: 1, significand: 0)
38+
XCTAssertEqual(Augmented.sum(x, y).head.sign, .plus) // (1)
39+
XCTAssertEqual(Augmented.sum(x, y).tail.sign, .plus)
40+
x = T(sign: .minus, exponent: 1, significand: 0)
41+
y = T(sign: .plus, exponent: 1, significand: 0)
42+
XCTAssertEqual(Augmented.sum(x, y).head.sign, .plus) // (2)
43+
XCTAssertEqual(Augmented.sum(x, y).tail.sign, .plus)
44+
x = T(sign: .minus, exponent: 1, significand: 0)
45+
y = T(sign: .minus, exponent: 1, significand: 0)
46+
XCTAssertEqual(Augmented.sum(x, y).head.sign, .minus) // (3)
47+
XCTAssertEqual(Augmented.sum(x, y).tail.sign, .plus)
48+
// Infinity and NaN are propagated correctly
49+
XCTAssertEqual(Augmented.sum( 0, T.infinity).head, T.infinity)
50+
XCTAssertEqual(Augmented.sum( T.infinity, 0 ).head, T.infinity)
51+
XCTAssertEqual(Augmented.sum( T.infinity, T.infinity).head, T.infinity)
52+
XCTAssertEqual(Augmented.sum( 0, -T.infinity).head, -T.infinity)
53+
XCTAssertEqual(Augmented.sum(-T.infinity, 0 ).head, -T.infinity)
54+
XCTAssertEqual(Augmented.sum(-T.infinity, -T.infinity).head, -T.infinity)
55+
XCTAssert(Augmented.sum( T.infinity, -T.infinity).head.isNaN)
56+
XCTAssert(Augmented.sum(-T.infinity, T.infinity).head.isNaN)
57+
XCTAssert(Augmented.sum( T.infinity, T.nan).head.isNaN)
58+
XCTAssert(Augmented.sum( T.nan, T.infinity).head.isNaN)
59+
XCTAssert(Augmented.sum(-T.infinity, T.nan).head.isNaN)
60+
XCTAssert(Augmented.sum( T.nan, -T.infinity).head.isNaN)
61+
XCTAssert(Augmented.sum( 0, T.nan).head.isNaN)
62+
XCTAssert(Augmented.sum(T.nan, 0).head.isNaN)
63+
XCTAssert(Augmented.sum(T.nan, T.nan).head.isNaN)
64+
}
65+
66+
func testTwoSumRandomValues<T: Real & FixedWidthFloatingPoint>(_: T.Type) {
67+
// For randomly-chosen well-scaled finite values, we expect:
68+
// (1) `head` to be exactly the IEEE 754 sum of `a + b`
69+
// (2) `tail` to be less than or equal `head.ulp/2`
70+
// (3) the result of `twoSum` for unordered input to be exactly equal to
71+
// the result of `fastTwoSum` for ordered input.
72+
var g = SystemRandomNumberGenerator()
73+
let values: [T] = (0 ..< 100).map { _ in
74+
T.random(
75+
in: T.ulpOfOne ..< 1,
76+
using: &g)
77+
}
78+
for a in values {
79+
for b in values {
80+
let twoSum = Augmented.sum(a, b)
81+
XCTAssertEqual(twoSum.head, a + b) // (1)
82+
XCTAssert(twoSum.tail.magnitude <= twoSum.head.ulp/2) // (2)
83+
let x: T = a.magnitude < b.magnitude ? b : a
84+
let y: T = a.magnitude < b.magnitude ? a : b
85+
let fastTwoSum = Augmented.sum(large: x, small: y)
86+
XCTAssertEqual(twoSum.head, fastTwoSum.head) // (3)
87+
XCTAssertEqual(twoSum.tail, fastTwoSum.tail) // (3)
88+
}
89+
}
90+
}
91+
92+
func testTwoSumCancellation<T: Real & FixedWidthFloatingPoint>(_: T.Type) {
93+
// Must be exact for exactly representable values
94+
XCTAssertEqual(Augmented.sum( 0.984375, 1.375).head, 2.359375)
95+
XCTAssertEqual(Augmented.sum( 0.984375, 1.375).tail, 0.0)
96+
XCTAssertEqual(Augmented.sum(-0.984375, 1.375).head, 0.390625)
97+
XCTAssertEqual(Augmented.sum(-0.984375, 1.375).tail, 0.0)
98+
XCTAssertEqual(Augmented.sum( 0.984375,-1.375).head,-0.390625)
99+
XCTAssertEqual(Augmented.sum( 0.984375,-1.375).tail, 0.0)
100+
XCTAssertEqual(Augmented.sum(-0.984375,-1.375).head,-2.359375)
101+
XCTAssertEqual(Augmented.sum(-0.984375,-1.375).tail, 0.0)
102+
XCTAssertEqual(Augmented.sum( 1.375, 0.984375).head, 2.359375)
103+
XCTAssertEqual(Augmented.sum( 1.375, 0.984375).tail, 0.0)
104+
XCTAssertEqual(Augmented.sum( 1.375,-0.984375).head, 0.390625)
105+
XCTAssertEqual(Augmented.sum( 1.375,-0.984375).tail, 0.0)
106+
XCTAssertEqual(Augmented.sum(-1.375, 0.984375).head,-0.390625)
107+
XCTAssertEqual(Augmented.sum(-1.375, 0.984375).tail, 0.0)
108+
XCTAssertEqual(Augmented.sum(-1.375,-0.984375).head,-2.359375)
109+
XCTAssertEqual(Augmented.sum(-1.375,-0.984375).tail, 0.0)
110+
// Must handle cancellation when `b` is not representable in `a` and
111+
// we expect `b` to be lost entirely in the calculation of `a + b`.
112+
var a: T = 1.0
113+
var b: T = .ulpOfOne * .ulpOfOne
114+
var twoSum = Augmented.sum(a, b)
115+
XCTAssertEqual(twoSum.head, a) // a + b = a
116+
XCTAssertEqual(twoSum.tail, b) // Error: b
117+
twoSum = Augmented.sum( a, -b)
118+
XCTAssertEqual(twoSum.head, a)
119+
XCTAssertEqual(twoSum.tail,-b)
120+
twoSum = Augmented.sum(-a, b)
121+
XCTAssertEqual(twoSum.head,-a)
122+
XCTAssertEqual(twoSum.tail, b)
123+
twoSum = Augmented.sum(-a, -b)
124+
XCTAssertEqual(twoSum.head,-a)
125+
XCTAssertEqual(twoSum.tail,-b)
126+
// Must handle cancellation when `b` is only partially representable in `a`.
127+
// We expect the fractional digits of `b` to be cancelled in the following
128+
// example but the fractional digits to be preserved in `tail`.
129+
let exponent = T.Exponent(T.significandBitCount + 1)
130+
a = T(sign: .plus, exponent: exponent, significand: 1.0)
131+
b = 256 + 0.5
132+
twoSum = Augmented.sum( a, b)
133+
XCTAssertEqual(twoSum.head, a + 256)
134+
XCTAssertEqual(twoSum.tail, 0.5)
135+
twoSum = Augmented.sum( a, -b)
136+
XCTAssertEqual(twoSum.head, a - 256)
137+
XCTAssertEqual(twoSum.tail, -0.5)
138+
twoSum = Augmented.sum(-a, b)
139+
XCTAssertEqual(twoSum.head, -a + 256)
140+
XCTAssertEqual(twoSum.tail, 0.5)
141+
twoSum = Augmented.sum(-a, -b)
142+
XCTAssertEqual(twoSum.head, -a - 256)
143+
XCTAssertEqual(twoSum.tail, -0.5)
144+
}
145+
146+
func testTwoSum() {
147+
testTwoSumSpecials(Float32.self)
148+
testTwoSumRandomValues(Float32.self)
149+
testTwoSumCancellation(Float32.self)
150+
testTwoSumSpecials(Float64.self)
151+
testTwoSumRandomValues(Float64.self)
152+
testTwoSumCancellation(Float64.self)
153+
#if (arch(i386) || arch(x86_64)) && !os(Windows) && !os(Android)
154+
testTwoSumSpecials(Float80.self)
155+
testTwoSumRandomValues(Float80.self)
156+
testTwoSumCancellation(Float80.self)
157+
#endif
158+
}
159+
}
160+

0 commit comments

Comments
 (0)