Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Document the Augmented enum in RealModule in advance of 1.0 #181

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions Sources/ComplexModule/ElementaryFunctions.swift
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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, θ)
}
Expand Down
71 changes: 69 additions & 2 deletions Sources/RealModule/AugmentedArithmetic.swift
Original file line number Diff line number Diff line change
Expand Up @@ -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:
/// -

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is a blank bullet here and also in some lists below. Is it intentional?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this produces the desired formatting in Xcode:
Screen Shot 2021-08-19 at 2 21 14 PM

/// - `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`.
Comment on lines +43 to +44
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are the “abs” on the right-hand sides necessary? I was under the impression that ulp was guaranteed to be non-negative for normal values. If so, the same thing applies to the doc comment for sum(large:small:).

Copy link
Member Author

@stephentyrone stephentyrone Aug 19, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's right (I'll roll this in with another upcoming PR)

/// - 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
Expand Down