Merge pull request #181 from stephentyrone/document-augmented-arithmetic

Document the Augmented enum in RealModule in advance of 1.0

Author: stephentyrone

Sources/ComplexModule/ElementaryFunctions.swift

@@ -311,13 +311,13 @@ extension Complex: ElementaryFunctions {
     // as exact head-tail products (u is guaranteed to be well scaled,
     // v may underflow, but if it does it doesn't matter, the u term is
     // all we need).
-    let (a,b) = Augmented.twoProdFMA(u, u)
-    let (c,d) = Augmented.twoProdFMA(v, v)
+    let (a,b) = Augmented.product(u, u)
+    let (c,d) = Augmented.product(v, v)
     // It would be nice if we could simply use a - 1, but unfortunately
     // we don't have a tight enough bound to guarantee that that expression
     // is exact; a may be as small as 1/4, so we could lose a single bit
     // to rounding if we did that.
-    var (s,e) = Augmented.fastTwoSum(-1, a)
+    var (s,e) = Augmented.sum(large: -1, small: a)
     // Now we are ready to assemble the result. If cancellation happens,
     // then |c| > |e| > |b|, |d|, so this assembly order is safe. It's
     // also possible that |c| and |d| are small, but if that happens then
@@ -353,9 +353,9 @@ extension Complex: ElementaryFunctions {
     // So now we need to evaluate (2+x)x + y² accurately. To do this,
     // we employ augmented arithmetic; the key observation here is that
     // cancellation is only a problem when y² ≈ -(2+x)x
-    let xp2 = Augmented.fastTwoSum(2, z.x) // Known that 2 > |x|.
-    let a = Augmented.twoProdFMA(z.x, xp2.head)
-    let y² = Augmented.twoProdFMA(z.y, z.y)
+    let xp2 = Augmented.sum(large: 2, small: z.x) // Known that 2 > |x|.
+    let a = Augmented.product(z.x, xp2.head)
+    let y² = Augmented.product(z.y, z.y)
     let s = (a.head + y².head + a.tail + y².tail).addingProduct(z.x, xp2.tail)
     return Complex(.log(onePlus: s)/2, θ)
   }

Sources/RealModule/AugmentedArithmetic.swift

@@ -9,18 +9,85 @@
 //
 //===----------------------------------------------------------------------===//
 
+/// A namespace for "augmented arithmetic" operations for types conforming to
+/// `Real`.
+///
+/// Augmented arithmetic refers to a family of algorithms that represent
+/// the results of floating-point computations using multiple values such that
+/// either the error is minimized or the result is exact.
 public enum Augmented { }
 
 extension Augmented {
+  /// The product `a * b` represented as an implicit sum `head + tail`.
+  ///
+  /// `head` is the correctly rounded value of `a*b`. If no overflow or
+  /// underflow occurs, `tail` represents the rounding error incurred in
+  /// computing `head`, such that the exact product is the sum of `head`
+  /// and `tail` computed without rounding.
+  ///
+  /// This operation is sometimes called "twoProd" or "twoProduct".
+  ///
+  /// Edge Cases:
+  /// -
+  /// - `head` is always the IEEE 754 product `a * b`.
+  /// - If `head` is not finite, `tail` is unspecified and should not be
+  ///   interpreted as having any meaning (it may be `NaN` or `infinity`).
+  /// - When `head` is close to the underflow boundary, the rounding error
+  ///   may not be representable due to underflow, and `tail` will be rounded.
+  ///   If `head` is very small, `tail` may even be zero, even though the
+  ///   product is not exact.
+  /// - If `head` is zero, `tail` is also a zero with unspecified sign.
+  ///
+  /// Postconditions:
+  /// -
+  /// - If `head` is normal, then `abs(tail) < abs(head.ulp)`.
+  ///   Assuming IEEE 754 default rounding, `abs(tail) <= abs(head.ulp)/2`.
+  /// - If both `head` and `tail` are normal, then `a * b` is exactly
+  ///   equal to `head + tail` when computed as real numbers.
   @_transparent
-  public static func twoProdFMA<T:Real>(_ a: T, _ b: T) -> (head: T, tail: T) {
+  public static func product<T:Real>(_ a: T, _ b: T) -> (head: T, tail: T) {
     let head = a*b
+    // TODO: consider providing an FMA-less implementation for use when
+    // targeting platforms without hardware FMA support. This works everywhere,
+    // falling back on the C math.h fma funcions, but may be slow on older x86.
     let tail = (-head).addingProduct(a, b)
     return (head, tail)
   }
 
+  /// The sum `a + b` represented as an implicit sum `head + tail`.
+  ///
+  /// `head` is the correctly rounded value of `a + b`. `tail` is the
+  /// error from that computation rounded to the closest representable
+  /// value.
+  ///
+  /// Unlike `Augmented.product(a, b)`, the rounding error of a sum can
+  /// never underflow. However, it may not be exactly representable when
+  /// `a` and `b` differ widely in magnitude.
+  ///
+  /// This operation is sometimes called "fastTwoSum".
+  ///
+  /// - Parameters:
+  ///   - a: The summand with larger magnitude.
+  ///   - b: The summand with smaller magnitude.
+  ///
+  /// Preconditions:
+  /// -
+  /// - `large.magnitude` must not be smaller than `small.magnitude`.
+  ///   They may be equal, or one or both may be `NaN`.
+  ///   This precondition is only enforced in debug builds.
+  ///
+  /// Edge Cases:
+  /// -
+  /// - `head` is always the IEEE 754 sum `a + b`.
+  /// - If `head` is not finite, `tail` is unspecified and should not be
+  ///   interpreted as having any meaning (it may be `NaN` or `infinity`).
+  ///
+  /// Postconditions:
+  /// -
+  /// - If `head` is normal, then `abs(tail) < abs(head.ulp)`.
+  ///   Assuming IEEE 754 default rounding, `abs(tail) <= abs(head.ulp)/2`.
   @_transparent
-  public static func fastTwoSum<T:Real>(_ a: T, _ b: T) -> (head: T, tail: T) {
+  public static func sum<T:Real>(large a: T, small b: T) -> (head: T, tail: T) {
     assert(!(b.magnitude > a.magnitude))
     let head = a + b
     let tail = a - head + b