apple / swift-numerics

Advanced mathematical types and functions for Swift
Apache License 2.0
1.69k stars 147 forks source link

Should complex multiplication use augmented arithmetic? #238

Open NevinBR opened 2 years ago

NevinBR commented 2 years ago

Currently, complex multiplication is implemented with * and +:

https://github.com/apple/swift-numerics/blob/5428505255ad35b8debd2c2c432350b213e4dff5/Sources/ComplexModule/Complex%2BNumeric.swift#L15-L17

However, there are well-known cancellation issues when computing the sum or difference of products in this manner. Here are some references:

Given four numbers (a, b, c, d), we want to compute ab-cd to good accuracy even when the two products are nearly equal. The sources above (among others) describe a way to do so using fused multiply-add instructions, which can be written in Swift as:

// Computes a*b - c*d accurately
func differenceOfProducts<T: FloatingPoint>(_ a: T, _ b: T, _ c: T, _ d: T) -> T {
  let cd = c * d
  let e = cd.addingProduct(-c, d)
  let ab_cd = (-cd).addingProduct(a, b)
  return ab_cd + e
}

This uses one more multiplication than the naive approach, but it is robust against catastrophic cancellation and it mitigates some (but not all) spurious overflows.

The same idea can be applied to the sum of products, and either version can perform both operations by simply negating one of the inputs. Furthermore, there are a variety of options for the name and signature of the function, such as:

func sumOfProducts<T: FloatingPoint>(_ a: T, _ b: T, _ c: T, _ d: T) -> T { ... }

func add<T: FloatingPoint>(product ab: (T, T), plusProduct cd: (T, T)) -> T { ... }

func subtract<T: FloatingPoint>(product ab: (T, T), minusProduct cd: (T, T)) -> T { ... }

func dotProduct<T: FloatingPoint>(_ u: (T, T), _ v: (T, T)) -> T { ... }

func determinant<T: FloatingPoint>(_ M: (T, T, T, T) -> T { ... }

...

Should we add a function along these lines to Numerics, and use it to compute the components of complex multiplication?

If so, what spelling do we prefer, and where should it reside? (Free function / static method on FloatingPoint / namespaced in Augmented / etc.)

stephentyrone commented 1 year ago

Certainly we would never do this for the "normal" multiplication and division; cancellation occurs, but it is harmless--even under cancellation, the result satisfies both componentwise and normwise backwards error bounds, and normwise forwards error bounds, which is all you need for stability of numerical algorithms. This is why no other language uses these formulas.

These are occasionally useful as building blocks, but I would probably tuck them under Augmented for that purpose.