Add exponential function to quaternions

Author: markuswntr

Sources/QuaternionModule/ElementaryFunctions.swift

@@ -0,0 +1,60 @@
+//===--- ElementaryFunctions.swift ----------------------------*- swift -*-===//
+//
+// This source file is part of the Swift.org open source project
+//
+// Copyright (c) 2019-2021 Apple Inc. and the Swift project authors
+// Licensed under Apache License v2.0 with Runtime Library Exception
+//
+// See https://swift.org/LICENSE.txt for license information
+// See https://swift.org/CONTRIBUTORS.txt for the list of Swift project authors
+//
+//===----------------------------------------------------------------------===//
+
+import RealModule
+
+extension Quaternion/*: ElementaryFunctions */ {
+
+  // MARK: - exp-like functions
+
+  /// The quaternion exponential function e^q whose base `e` is the base of the
+  /// natural logarithm.
+  ///
+  /// Mathematically, this operation can be expanded in terms of the `Real`
+  /// operations `exp`, `cos` and `sin` as follows:
+  /// ```
+  /// exp(r + xi + yj + zk) = exp(r + v) = exp(r) exp(v)
+  ///                       = exp(r) (cos(θ) + (v/θ) sin(θ)) where θ = ||v||
+  /// ```
+  /// Note that naive evaluation of this expression in floating-point would be
+  /// prone to premature overflow, since `cos` and `sin` both have magnitude
+  /// less than 1 for most inputs (i.e. `exp(r)` may be infinity when
+  /// `exp(r) cos(θ)` would not be).
+  public static func exp(_ q: Quaternion<RealType>) -> Quaternion<RealType> {
+    guard q.isFinite else { return q }
+    // Firstly evaluate θ and v/θ where θ = ||v|| (as discussed above)
+    // There are 2 special cases for ||v|| that we need to take care of:
+    // The value of ||v|| may be invalid due to an overflow in `.lengthSquared`.
+    // As the internal `SIMD3.length` helper functions deals with overflow and
+    // underflow of `.lengthSquared`, we can safely ignore this case here.
+    // However, we still have to check for ||v|| = 0 before evaluating v/θ
+    // as it would incorrectly yield a division by zero.
+    let phase = q.imaginary.length
+    let unitAxis = !phase.isZero ? (q.imaginary / phase) : .zero
+    // If real < log(greatestFiniteMagnitude), then exp(q.real) does not overflow.
+    // To protect ourselves against sketchy log or exp implementations in
+    // an unknown host library, or slight rounding disagreements between
+    // the two, subtract one from the bound for a little safety margin.
+    guard q.real < RealType.log(.greatestFiniteMagnitude) - 1 else {
+      let halfScale = RealType.exp(q.real/2)
+      let rotation = Quaternion(
+        halfAngle: phase,
+        unitAxis: unitAxis
+      )
+      return rotation.multiplied(by: halfScale).multiplied(by: halfScale)
+    }
+    return Quaternion(
+      halfAngle: phase,
+      unitAxis: unitAxis
+    ).multiplied(by: .exp(q.real))
+  }
+}

Tests/QuaternionTests/ElementaryFunctionTests.swift

@@ -0,0 +1,85 @@
+//===--- ElementaryFunctionTests.swift ------------------------*- swift -*-===//
+//
+// This source file is part of the Swift Numerics open source project
+//
+// Copyright (c) 2019 - 2021 Apple Inc. and the Swift Numerics project authors
+// Licensed under Apache License v2.0 with Runtime Library Exception
+//
+// See https://swift.org/LICENSE.txt for license information
+//
+//===----------------------------------------------------------------------===//
+
+import XCTest
+import RealModule
+import _TestSupport
+
+@testable import QuaternionModule
+
+final class ElementaryFunctionTests: XCTestCase {
+
+  func testExp<T: Real & FixedWidthFloatingPoint & SIMDScalar>(_ type: T.Type) {
+    // exp(0) = 1
+    XCTAssertEqual(1, Quaternion<T>.exp(Quaternion(real: 0, imaginary: 0, 0, 0)))
+    XCTAssertEqual(1, Quaternion<T>.exp(Quaternion(real:-0, imaginary: 0, 0, 0)))
+    XCTAssertEqual(1, Quaternion<T>.exp(Quaternion(real:-0, imaginary:-0,-0,-0)))
+    XCTAssertEqual(1, Quaternion<T>.exp(Quaternion(real: 0, imaginary:-0,-0,-0)))
+    // In general, exp(Quaternion(r, 0, 0, 0)) should be exp(r), but that breaks
+    // down when r is infinity or NaN, because we want all non-finite
+    // quaternions to be semantically a single point at infinity. This is fine
+    // for most inputs, but exp(Quaternion(-.infinity, 0, 0, 0)) would produce
+    // 0 if we evaluated it in the usual way.
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real:  .infinity, imaginary: .zero)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real:  .infinity, imaginary: .infinity)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real:          0, imaginary: .infinity)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real: -.infinity, imaginary: .infinity)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real: -.infinity, imaginary: .zero)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real: -.infinity, imaginary: -.infinity)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real:          0, imaginary: -.infinity)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real:  .infinity, imaginary: -.infinity)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real:       .nan, imaginary: .nan)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real:  .infinity, imaginary: .nan)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real:       .nan, imaginary: .infinity)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real: -.infinity, imaginary: .nan)).isFinite)
+    XCTAssertFalse(Quaternion<T>.exp(Quaternion(real:       .nan, imaginary: -.infinity)).isFinite)
+    // Find a value of x such that exp(x) just overflows. Then exp((x, π/4))
+    // should not overflow, but will do so if it is not computed carefully.
+    // The correct value is:
+    //
+    //   exp((log(gfm) + log(9/8), π/4) = exp((log(gfm*9/8), π/4))
+    //                                  = gfm*9/8 * (1/sqrt(2), 1/(sqrt(2))
+    let x = T.log(.greatestFiniteMagnitude) + T.log(9/8)
+    let huge = Quaternion<T>.exp(Quaternion(real: x, imaginary: SIMD3(.pi/4, 0, 0)))
+    let mag = T.greatestFiniteMagnitude/T.sqrt(2) * (9/8)
+    XCTAssert(huge.real.isApproximatelyEqual(to: mag))
+    XCTAssert(huge.imaginary.x.isApproximatelyEqual(to: mag))
+    XCTAssertEqual(huge.imaginary.y, .zero)
+    XCTAssertEqual(huge.imaginary.z, .zero)
+    // For randomly-chosen well-scaled finite values, we expect to have the
+    // usual identities:
+    //
+    //   exp(z + w) = exp(z) * exp(w)
+    //   exp(z - w) = exp(z) / exp(w)
+    var g = SystemRandomNumberGenerator()
+    let values: [Quaternion<T>] = (0..<100).map { _ in
+      Quaternion(
+        real: T.random(in: -1 ... 1, using: &g),
+        imaginary: SIMD3(repeating: T.random(in: -.pi ... .pi, using: &g) / 3))
+    }
+    for z in values {
+      for w in values {
+        let p = Quaternion.exp(z) * Quaternion.exp(w)
+        let q = Quaternion.exp(z) / Quaternion.exp(w)
+        XCTAssert(Quaternion.exp(z + w).isApproximatelyEqual(to: p))
+        XCTAssert(Quaternion.exp(z - w).isApproximatelyEqual(to: q))
+      }
+    }
+  }
+
+  func testFloat() {
+    testExp(Float32.self)
+  }
+
+  func testDouble() {
+    testExp(Float64.self)
+  }
+}