Merge pull request #185 from stephentyrone/reciprocal-documentation

Further expand documentation of `reciprocal` with some example use.

Author: stephentyrone

Sources/ComplexModule/Arithmetic.swift

@@ -52,13 +52,31 @@ extension Complex: AdditiveArithmetic {
 // and turn these into operators if/when we have it.
 // (https://github.com/apple/swift-numerics/issues/12)
 extension Complex {
+  /// `self` scaled by `a`.
   @usableFromInline @_transparent
   internal func multiplied(by a: RealType) -> Complex {
+    // This can be viewed in two different ways, which are mathematically
+    // equivalent: either we are computing `self * Complex(a)` (i.e.
+    // converting `a` to be a complex value, and then using the complex
+    // multiplication) or we are using the scalar product of the vector
+    // space structure: `Complex(a*real, a*imaginary)`.
+    //
+    // Although these two interpretations are _mathematically_ equivalent,
+    // they will generate different representations of the point at
+    // infinity in general. For example, suppose `self` is represented by
+    // `(infinity, 0)`. Then `self * Complex(1)` would evaluate as
+    // `(1*infinity - 0*0, 0*infinity + 1*0) = (infinity, nan)`, but
+    // the vector space interpretation produces `(infinity, 0)`. This does
+    // not matter much, because these are two representations of the same
+    // semantic value, but note that one requires four multiplies and two
+    // additions, while the one we use requires only two real multiplications.
     Complex(x*a, y*a)
   }
 
+  /// `self` unscaled by `a`.
   @usableFromInline @_transparent
   internal func divided(by a: RealType) -> Complex {
+    // See implementation notes for `multiplied` above.
     Complex(x/a, y/a)
   }
 }
@@ -147,17 +165,23 @@ extension Complex: AlgebraicField {
   /// isolated division, but if you are dividing many values by a single
   /// denominator, this will often be a significant performance win.
   ///
-  /// Typical use looks like this:
+  /// A typical use case looks something like this:
   /// ```
   /// func divide<T: Real>(data: [Complex<T>], by divisor: Complex<T>) -> [Complex<T>] {
-  ///   // If divisor is well-scaled, use multiply by reciprocal.
+  ///   // If divisor is well-scaled, multiply by reciprocal.
   ///   if let recip = divisor.reciprocal {
   ///     return data.map { $0 * recip }
   ///   }
   ///   // Fallback on using division.
   ///   return data.map { $0 / divisor }
   /// }
   /// ```
+  ///
+  /// Error Bounds:
+  /// -
+  /// Unlike real types, when working with complex types, multiplying by the
+  /// reciprocal instead of dividing cannot change the result. If the
+  /// reciprocal is non-nil, the two computations are always equivalent.
   @inlinable
   public var reciprocal: Complex? {
     let recip = 1/self

Sources/RealModule/AlgebraicField.swift

@@ -69,6 +69,30 @@ public protocol AlgebraicField: SignedNumeric {
   /// Note that `.zero.reciprocal`, somewhat surprisingly, is *not* nil
   /// for `Real` or `Complex` types, because these types have an
   /// `.infinity` value that acts as the reciprocal of `.zero`.
+  ///
+  /// If `b.reciprocal` is non-nil, you may be able to replace division by `b`
+  /// with multiplication by this value. It is not advantageous to do this
+  /// for an isolated division unless it is a compile-time constant visible
+  /// to the compiler, but if you are dividing many values by a single
+  /// denominator, this will often be a significant performance win.
+  ///
+  /// Note that this will slightly perturb results for fields with approximate
+  /// arithmetic, such as real or complex types--using a normal division
+  /// is generally more accurate--but no catastrophic loss of accuracy will
+  /// result. For fields with exact arithmetic, the results are necessarily
+  /// identical.
+  ///
+  /// A typical use case looks something like this:
+  /// ```
+  /// func divide<T: AlgebraicField>(data: [T], by divisor: T) -> [T] {
+  ///   // If divisor is well-scaled, multiply by reciprocal.
+  ///   if let recip = divisor.reciprocal {
+  ///     return data.map { $0 * recip }
+  ///   }
+  ///   // Fallback on using division.
+  ///   return data.map { $0 / divisor }
+  /// }
+  /// ```
   var reciprocal: Self? { get }
 }
 

Sources/RealModule/Real.swift

@@ -103,7 +103,47 @@ extension Real {
   /// if it is representable.
   ///
   /// If `a` if finite and nonzero, and `1/a` overflows or underflows,
-  /// then `a.reciprocal` is `nil`. Otherwise, `a.reciprcoal` is `1/a`.
+  /// then `a.reciprocal` is `nil`. Otherwise, `a.reciprocal` is `1/a`.
+  ///
+  /// If `b.reciprocal` is non-nil, you may be able to replace division by `b`
+  /// with multiplication by this value. It is not advantageous to do this
+  /// for an isolated division unless it is a compile-time constant visible
+  /// to the compiler, but if you are dividing many values by a single
+  /// denominator, this will often be a significant performance win.
+  ///
+  /// A typical use case looks something like this:
+  /// ```
+  /// func divide<T: Real>(data: [T], by divisor: T) -> [T] {
+  ///   // If divisor is well-scaled, multiply by reciprocal.
+  ///   if let recip = divisor.reciprocal {
+  ///     return data.map { $0 * recip }
+  ///   }
+  ///   // Fallback on using division.
+  ///   return data.map { $0 / divisor }
+  /// }
+  /// ```
+  ///
+  /// Error Bounds:
+  /// -
+  /// Multiplying by the reciprocal instead of dividing will slightly
+  /// perturb results. For example `5.0 / 3` is 1.6666666666666667, but
+  /// `5.0 * 3.reciprocal!` is 1.6666666666666665.
+  ///
+  /// The error of a normal division is bounded by half an ulp of the
+  /// result; we can derive a quick error bound for multiplication by
+  /// the real reciprocal (when it exists) as follows (I will use circle
+  /// operators to denote real-number arithmetic, and normal operators
+  /// for floating-point arithmetic):
+  ///
+  ///   a * b.reciprocal! = a * (1/b)
+  ///                     = a * (1 ⊘ b)(1 + δ₁)
+  ///                     = (a ⊘ b)(1 + δ₁)(1 + δ₂)
+  ///                     = (a ⊘ b)(1 + δ₁ + δ₂ + δ₁δ₂)
+  ///
+  /// where 0 < δᵢ <= ulpOfOne/2. This gives a roughly 1-ulp error,
+  /// about twice the error bound we get using division. For most
+  /// purposes this is an acceptable error, but if you need to match
+  /// results obtained using division, you should not use this.
   @inlinable
   public var reciprocal: Self? {
     let recip = 1/self