Conform Quaternion to Elementary Functions

Package.swift

@@ -20,6 +20,7 @@ let package = Package(
   products: [
     .library(name: "ComplexModule", targets: ["ComplexModule"]),
     .library(name: "Numerics", targets: ["Numerics"]),
+    .library(name: "QuaternionModule", targets: ["QuaternionModule"]),
     .library(name: "RealModule", targets: ["RealModule"]),
   ],
 
@@ -39,9 +40,18 @@ let package = Package(
 
     .target(
       name: "Numerics",
-      dependencies: ["ComplexModule", "IntegerUtilities", "RealModule"],
+      dependencies: [
+        "ComplexModule", "IntegerUtilities",
+        "QuaternionModule", "RealModule"
+      ],
       exclude: excludedFilenames
     ),
+
+    .target(
+      name: "QuaternionModule",
+      dependencies: ["RealModule"],
+      exclude: ["README.md", "Transformation.md"]
+    ),
 
     .target(
       name: "RealModule",
@@ -74,6 +84,11 @@ let package = Package(
       dependencies: ["IntegerUtilities", "_TestSupport"],
       exclude: ["CMakeLists.txt"]
     ),
+
+    .testTarget(
+      name: "QuaternionTests",
+      dependencies: ["_TestSupport"]
+    ),
 
     .testTarget(
       name: "RealTests",

README.md

@@ -106,6 +106,7 @@ Questions about how to use Swift Numerics modules, or issues that are not clearl
 1. [`RealModule`](Sources/RealModule/README.md)
 2. [`ComplexModule`](Sources/ComplexModule/README.md)
 3. [`IntegerUtilities`](Sources/IntegerUtilties/README.md) (on main only, not yet present in a released tag)
+4. [`QuaternionModule`](Sources/QuaternionModule/README.md)
 
 ## Future expansion
 

Sources/ComplexModule/Complex+ElementaryFunctions.swift

@@ -32,6 +32,12 @@
 // Except where derivations are given, the expressions used here are all
 // adapted from Kahan's 1986 paper "Branch Cuts for Complex Elementary
 // Functions; or: Much Ado About Nothing's Sign Bit".
+//
+// Elementary functions of complex numbers have many similarities with
+// elementary functions of quaternions and their definition in terms of real
+// operations. Therefore, if you make a modification to one of the following
+// functions, you should almost surely make a parallel modification to the same
+// elementary function of quaternions.
 
 import RealModule
 

Sources/Numerics/Numerics.swift

@@ -12,4 +12,5 @@
 // A module that re-exports the complete Swift Numerics public API.
 @_exported import ComplexModule
 @_exported import IntegerUtilities
+@_exported import QuaternionModule
 @_exported import RealModule

Sources/QuaternionModule/ImaginaryHelper.swift

@@ -0,0 +1,77 @@
+//===--- ImaginaryHelper.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
+//
+//===----------------------------------------------------------------------===//
+
+// Provides common vector operations on SIMD3 to ease the use of the quaternions
+// imaginary/vector components internally to the module, and in tests.
+extension SIMD3 where Scalar: FloatingPoint {
+  /// Returns a vector with infinity in all lanes
+  @usableFromInline @inline(__always)
+  internal static var infinity: Self {
+    SIMD3(repeating: .infinity)
+  }
+
+  /// True if all values of this instance are finite
+  @usableFromInline @inline(__always)
+  internal var isFinite: Bool {
+    x.isFinite && y.isFinite && z.isFinite
+  }
+
+  /// The ∞-norm of the value (`max(abs(x), abs(y), abs(z))`).
+  @usableFromInline @inline(__always)
+  internal var magnitude: Scalar {
+    Swift.max(x.magnitude, y.magnitude, z.magnitude)
+  }
+
+  /// The 1-norm of the value (`abs(x) + abs(y) + abs(z))`).
+  @usableFromInline @inline(__always)
+  internal var oneNorm: Scalar {
+    x.magnitude + y.magnitude + z.magnitude
+  }
+
+  /// The Euclidean norm (a.k.a. 2-norm, `sqrt(x*x + y*y + z*z)`).
+  @usableFromInline @inline(__always)
+  internal var length: Scalar {
+    let naive = lengthSquared
+    guard naive.isNormal else { return carefulLength }
+    return naive.squareRoot()
+  }
+
+  // Implementation detail of `length`, moving slow path off of the
+  // inline function.
+  @usableFromInline
+  internal var carefulLength: Scalar {
+    guard isFinite else { return .infinity }
+    guard !magnitude.isZero else { return .zero }
+    // Unscale the vector, calculate its length and rescale the result
+    return (self / magnitude).length * magnitude
+  }
+
+  /// Returns the squared length of this instance.
+  @usableFromInline @inline(__always)
+  internal var lengthSquared: Scalar {
+    dot(self)
+  }
+
+  /// Returns the scalar/dot product of this vector with `other`.
+  @usableFromInline @inline(__always)
+  internal func dot(_ other: SIMD3<Scalar>) -> Scalar {
+    (self * other).sum()
+  }
+
+  /// Returns the vector/cross product of this vector with `other`.
+  @usableFromInline @inline(__always)
+  internal func cross(_ other: SIMD3<Scalar>) -> SIMD3<Scalar> {
+    let yzx = SIMD3<Int>(1,2,0)
+    let zxy = SIMD3<Int>(2,0,1)
+    return (self[yzx] * other[zxy]) - (self[zxy] * other[yzx])
+  }
+}

Sources/QuaternionModule/Polar.swift

@@ -0,0 +1,200 @@
+//===--- Polar.swift ------------------------------------------*- swift -*-===//
+//
+// This source file is part of the Swift Numerics open source project
+//
+// Copyright (c) 2019 - 2022 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
+//
+//===----------------------------------------------------------------------===//
+
+// Norms and related quantities defined for Quaternion.
+//
+// The following API are provided by this extension:
+//
+//   var magnitude: RealType     // infinity norm
+//   var length: RealType        // Euclidean norm
+//   var lengthSquared: RealType // Euclidean norm squared
+//
+// For detailed documentation, consult Norms.md or the inline documentation
+// for each operation.
+//
+// Implementation notes:
+//
+// `.magnitude` does not bind the Euclidean norm; it binds the infinity norm
+// instead. There are two reasons for this choice:
+//
+// - It's simply faster to compute in general, because it does not require
+//   a square root.
+//
+// - There exist finite values `q` for which the Euclidean norm is not
+//   representable (consider the quaternion with `r`, `x`, `y` and `z` all
+//   equal to `RealType.greatestFiniteMagnitude`; the Euclidean norm is
+//   `.sqrt(4) * .greatestFiniteMagnitude`, which overflows).
+//
+// The infinity norm is unique among the common vector norms in having
+// the property that every finite vector has a representable finite norm,
+// which makes it the obvious choice to bind `.magnitude`.
+extension Quaternion {
+  /// The Euclidean norm (a.k.a. 2-norm, `sqrt(r*r + x*x + y*y + z*z)`).
+  ///
+  /// This value is highly prone to overflow or underflow.
+  ///
+  /// For most use cases, you can use the cheaper `.magnitude`
+  /// property (which computes the ∞-norm) instead, which always produces
+  /// a representable result.
+  ///
+  /// Edge cases:
+  /// - If a quaternion is not finite, its `.length` is `infinity`.
+  ///
+  /// See also `.magnitude`, `.lengthSquared`, `.polar` and
+  /// `init(length:,phase:,axis:)`.
+  @_transparent
+  public var length: RealType {
+    let naive = lengthSquared
+    guard naive.isNormal else { return carefulLength }
+    return .sqrt(naive)
+  }
+
+  // Internal implementation detail of `length`, moving slow path off
+  // of the inline function.
+  @usableFromInline
+  internal var carefulLength: RealType {
+    guard isFinite else { return .infinity }
+    guard !magnitude.isZero else { return .zero }
+    // Unscale the quaternion, calculate its length and rescale the result
+    return divided(by: magnitude).length * magnitude
+  }
+
+  /// The squared length `(r*r + x*x + y*y + z*z)`.
+  ///
+  /// This value is highly prone to overflow or underflow.
+  /// 
+  /// For many cases, `.magnitude` can be used instead, which is similarly
+  /// cheap to compute and always returns a representable value.
+  ///
+  /// This property is more efficient to compute than `length`.
+  ///
+  /// See also `.magnitude`, `.length`, `.polar` and
+  /// `init(length:,phase:,axis:)`.
+  @_transparent
+  public var lengthSquared: RealType {
+    (components * components).sum()
+  }
+
+  /// The [polar decomposition][wiki].
+  ///
+  /// Returns the length of this quaternion, phase in radians of range *[0, π]*
+  /// and the rotation axis as SIMD3 vector of unit length.
+  ///
+  /// Edge cases:
+  /// - If the quaternion is zero, length is `.zero` and angle and axis
+  /// are `nan`.
+  /// - If the quaternion is non-finite, length is `.infinity` and angle and
+  /// axis are `nan`.
+  /// - For any length other than `.zero` or `.infinity`, if angle is zero, axis
+  /// is `nan`.
+  ///
+  /// See also `.magnitude`, `.length`, `.lengthSquared` and
+  /// `init(length:,phase:,axis:)`.
+  ///
+  /// [wiki]: https://en.wikipedia.org/wiki/Polar_decomposition#Quaternion_polar_decomposition
+  public var polar: (length: RealType, phase: RealType, axis: SIMD3<RealType>) {
+    (length, halfAngle, axis)
+  }
+
+  /// Creates a quaternion specified with [polar coordinates][wiki].
+  ///
+  /// This initializer reads given `length`, `phase` and `axis` values and
+  /// creates a quaternion of equal rotation properties and specified *length*
+  /// using the following equation:
+  ///
+  ///     Q = (cos(phase), axis * sin(phase)) * length
+  ///
+  /// - Note: `axis` must be of unit length, or an assertion failure occurs.
+  ///
+  /// Edge cases:
+  /// - Negative lengths are interpreted as reflecting the point through the
+  ///   origin, i.e.:
+  ///   ```
+  ///   Quaternion(length: -r, phase: θ, axis: axis) == -Quaternion(length: r, phase: θ, axis: axis)
+  ///   ```
+  /// - For any `θ` and any `axis`, even `.infinity` or `.nan`:
+  ///   ```
+  ///   Quaternion(length: .zero, phase: θ, axis: axis) == .zero
+  ///   ```
+  /// - For any `θ` and any `axis`, even `.infinity` or `.nan`:
+  ///   ```
+  ///   Quaternion(length: .infinity, phase: θ, axis: axis) == .infinity
+  ///   ```
+  /// - Otherwise, `θ` must be finite, or a precondition failure occurs.
+  ///
+  /// See also `.magnitude`, `.length`, `.lengthSquared` and `.polar`.
+  ///
+  /// [wiki]: https://en.wikipedia.org/wiki/Polar_decomposition#Quaternion_polar_decomposition
+  @inlinable
+  public init(length: RealType, phase: RealType, axis: SIMD3<RealType>) {
+    guard !length.isZero, length.isFinite else {
+      self = Quaternion(length)
+      return
+    }
+
+    // Length is finite and non-zero, therefore
+    // 1. `phase` must be finite or a precondition failure needs to occur; as
+    //    this is not representable.
+    // 2. `axis` must be of unit length or an assertion failure occurs; while
+    //    "wrong" by definition, it is representable.
+    precondition(
+      phase.isFinite,
+      "Either phase must be finite, or length must be zero or infinite."
+    )
+    assert(
+      // TODO: Replace with `approximateEquality()`
+      abs(.sqrt(axis.lengthSquared)-1) < max(.sqrt(axis.lengthSquared), 1)*RealType.ulpOfOne.squareRoot(),
+      "Given axis must be of unit length."
+    )
+
+    self = Quaternion(halfAngle: phase, unitAxis: axis).multiplied(by: length)
+  }
+}
+
+// MARK: - Operations for working with polar form
+
+extension Quaternion {
+  /// The half rotation angle in radians within *[0, π]* range.
+  ///
+  /// Edge cases:
+  /// - If the quaternion is zero or non-finite, halfAngle is `nan`.
+  @usableFromInline @inline(__always)
+  internal var halfAngle: RealType {
+    guard isFinite else { return .nan }
+    guard imaginary != .zero else {
+      // A zero quaternion does not encode transformation properties.
+      // If imaginary is zero, real must be non-zero or nan is returned.
+      return real.isZero ? .nan : .zero
+    }
+
+    // If lengthSquared computes without over/underflow, everything is fine
+    // and the result is correct. If not, we have to do the computation
+    // carefully and unscale the quaternion first.
+    let lenSq = imaginary.lengthSquared
+    guard lenSq.isNormal else { return divided(by: magnitude).halfAngle }
+    return .atan2(y: .sqrt(lenSq), x: real)
+  }
+
+  /// Creates a new quaternion from given half rotation angle about given
+  /// rotation axis.
+  ///
+  /// The angle-axis values are transformed using the following equation:
+  ///
+  ///     Q = (cos(halfAngle), unitAxis * sin(halfAngle))
+  ///
+  /// - Parameters:
+  ///   - halfAngle: The half rotation angle
+  ///   - unitAxis: The rotation axis of unit length
+  @usableFromInline @inline(__always)
+  internal init(halfAngle: RealType, unitAxis: SIMD3<RealType>) {
+    self.init(real: .cos(halfAngle), imaginary: unitAxis * .sin(halfAngle))
+  }
+}

Sources/QuaternionModule/Quaternion+AdditiveArithmetic.swift

@@ -0,0 +1,41 @@
+//===--- Quaternion+AdditiveArithmetic.swift ------------------*- swift -*-===//
+//
+// This source file is part of the Swift Numerics open source project
+//
+// Copyright (c) 2019 - 2022 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
+//
+//===----------------------------------------------------------------------===//
+
+extension Quaternion: AdditiveArithmetic {
+  /// The additive identity, with real and *all* imaginary parts zero, i.e.:
+  /// `0 + 0i + 0j + 0k`
+  ///
+  /// See also: `one`, `i`, `j`, `k`, `infinity`
+  @_transparent
+  public static var zero: Quaternion {
+    Quaternion(from: SIMD4(repeating: 0))
+  }
+
+  @_transparent
+  public static func + (lhs: Quaternion, rhs: Quaternion) -> Quaternion {
+    Quaternion(from: lhs.components + rhs.components)
+  }
+
+  @_transparent
+  public static func - (lhs: Quaternion, rhs: Quaternion) -> Quaternion {
+    Quaternion(from: lhs.components - rhs.components)
+  }
+
+  @_transparent
+  public static func += (lhs: inout Quaternion, rhs: Quaternion) {
+    lhs = lhs + rhs
+  }
+
+  @_transparent
+  public static func -= (lhs: inout Quaternion, rhs: Quaternion) {
+    lhs = lhs - rhs
+  }
+}