diff --git a/Package.swift b/Package.swift index e8c31c55..ce64c6ed 100644 --- a/Package.swift +++ b/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", diff --git a/README.md b/README.md index 843b7707..15ae3064 100644 --- a/README.md +++ b/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 diff --git a/Sources/ComplexModule/Complex+ElementaryFunctions.swift b/Sources/ComplexModule/Complex+ElementaryFunctions.swift index efa685b4..a0225ff5 100644 --- a/Sources/ComplexModule/Complex+ElementaryFunctions.swift +++ b/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 diff --git a/Sources/Numerics/Numerics.swift b/Sources/Numerics/Numerics.swift index 2aa99036..15dd9b27 100644 --- a/Sources/Numerics/Numerics.swift +++ b/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 diff --git a/Sources/QuaternionModule/ImaginaryHelper.swift b/Sources/QuaternionModule/ImaginaryHelper.swift new file mode 100644 index 00000000..9ffd4cd2 --- /dev/null +++ b/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 { + (self * other).sum() + } + + /// Returns the vector/cross product of this vector with `other`. + @usableFromInline @inline(__always) + internal func cross(_ other: SIMD3) -> SIMD3 { + let yzx = SIMD3(1,2,0) + let zxy = SIMD3(2,0,1) + return (self[yzx] * other[zxy]) - (self[zxy] * other[yzx]) + } +} diff --git a/Sources/QuaternionModule/Polar.swift b/Sources/QuaternionModule/Polar.swift new file mode 100644 index 00000000..9d05bd17 --- /dev/null +++ b/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) { + (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) { + 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) { + self.init(real: .cos(halfAngle), imaginary: unitAxis * .sin(halfAngle)) + } +} diff --git a/Sources/QuaternionModule/Quaternion+AdditiveArithmetic.swift b/Sources/QuaternionModule/Quaternion+AdditiveArithmetic.swift new file mode 100644 index 00000000..f8e45b24 --- /dev/null +++ b/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 + } +} diff --git a/Sources/QuaternionModule/Quaternion+AlgebraicField.swift b/Sources/QuaternionModule/Quaternion+AlgebraicField.swift new file mode 100644 index 00000000..66c633ba --- /dev/null +++ b/Sources/QuaternionModule/Quaternion+AlgebraicField.swift @@ -0,0 +1,120 @@ +//===--- Quaternion+AlgebraicField.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 +// +//===----------------------------------------------------------------------===// + +import RealModule + +extension Quaternion: AlgebraicField { + /// The multiplicative identity, with real part one and *all* imaginary parts + /// zero, i.e.: `1 + 0i + 0j + 0k` + /// + /// See also: `zero`, `i`, `j`, `k`, `infinity` + @_transparent + public static var one: Quaternion { + Quaternion(from: SIMD4(0,0,0,1)) + } + + /// The [conjugate][conj] of this value. + /// + /// [conj]: https://en.wikipedia.org/wiki/Quaternion#Conjugation,_the_norm,_and_reciprocal + @_transparent + public var conjugate: Quaternion { + Quaternion(from: components * [-1, -1, -1, 1]) + } + + @_transparent + public static func / (lhs: Quaternion, rhs: Quaternion) -> Quaternion { + // Try the naive expression lhs/rhs = lhs*conj(rhs) / |rhs|^2; if we can compute + // this without over/underflow, everything is fine and the result is + // correct. If not, we have to rescale and do the computation carefully. + let lengthSquared = rhs.lengthSquared + guard lengthSquared.isNormal else { return rescaledDivide(lhs, rhs) } + return lhs * (rhs.conjugate.divided(by: lengthSquared)) + } + + @_transparent + public static func /= (lhs: inout Quaternion, rhs: Quaternion) { + lhs = lhs / rhs + } + + @usableFromInline @_alwaysEmitIntoClient @inline(never) + internal static func rescaledDivide(_ lhs: Quaternion, _ rhs: Quaternion) -> Quaternion { + if rhs.isZero { return .infinity } + if lhs.isZero || !rhs.isFinite { return .zero } + // TODO: detect when RealType is Float and just promote to Double, then + // use the naive algorithm. + let lhsScale = lhs.magnitude + let rhsScale = rhs.magnitude + let lhsNorm = lhs.divided(by: lhsScale) + let rhsNorm = rhs.divided(by: rhsScale) + let r = (lhsNorm * rhsNorm.conjugate).divided(by: rhsNorm.lengthSquared) + // At this point, the result is (r * lhsScale)/rhsScale computed without + // undue overflow or underflow. We know that r is close to unity, so + // the question is simply what order in which to do this computation + // to avoid spurious overflow or underflow. There are three options + // to choose from: + // + // - r * (lhsScale / rhsScale) + // - (r * lhsScale) / rhsScale + // - (r / rhsScale) * lhsScale + // + // The simplest case is when lhsScale / rhsScale is normal: + if (lhsScale / rhsScale).isNormal { + return r.multiplied(by: lhsScale / rhsScale) + } + // Otherwise, we need to compute either rNorm * lhsScale or rNorm / rhsScale + // first. Choose the first if the first scaling behaves well, otherwise + // choose the other one. + if (r.magnitude * lhsScale).isNormal { + return r.multiplied(by: lhsScale).divided(by: rhsScale) + } + return r.divided(by: rhsScale).multiplied(by: lhsScale) + } + + /// A normalized quaternion with the same direction and phase as this value. + /// + /// If such a value cannot be produced, `nil` is returned. + @inlinable + public var normalized: Quaternion? { + if length.isNormal { + return divided(by: length) + } + if isZero || !isFinite { + return nil + } + return divided(by: magnitude).normalized + } + + /// The reciprocal of this value, if it can be computed without undue overflow or underflow. + /// + /// If z.reciprocal is non-nil, you can safely replace division by z with multiplication by this + /// value. It is not advantageous to do this for an 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: + /// ``` + /// func divide(data: [Quaternion], by divisor: Quaternion) -> [Quaternion] { + /// // If divisor is well-scaled, use multiply by reciprocal. + /// if let recip = divisor.reciprocal { + /// return data.map { $0 * recip } + /// } + /// // Fallback on using division. + /// return data.map { $0 / divisor } + /// } + /// ``` + @inlinable + public var reciprocal: Quaternion? { + let recip = Quaternion(1)/self + if recip.isNormal || isZero || !isFinite { + return recip + } + return nil + } +} diff --git a/Sources/QuaternionModule/Quaternion+Codable.swift b/Sources/QuaternionModule/Quaternion+Codable.swift new file mode 100644 index 00000000..5c43de41 --- /dev/null +++ b/Sources/QuaternionModule/Quaternion+Codable.swift @@ -0,0 +1,25 @@ +//===--- Quaternion+Codable.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 +// +//===----------------------------------------------------------------------===// + +import RealModule + +// FloatingPoint does not refine Codable, so this is a conditional conformance. +extension Quaternion: Decodable where RealType: Decodable { + public init(from decoder: Decoder) throws { + try self.init(from: SIMD4(from: decoder)) + } +} + +extension Quaternion: Encodable where RealType: Encodable { + public func encode(to encoder: Encoder) throws { + try components.encode(to: encoder) + } +} diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift new file mode 100644 index 00000000..52e47acb --- /dev/null +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -0,0 +1,501 @@ +//===--- Quaternion+ElementaryFunctions.swift -----------------*- swift -*-===// +// +// This source file is part of the Swift.org open source project +// +// Copyright (c) 2019 - 2023 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 +// +//===----------------------------------------------------------------------===// + +// (r + xi + yj + zk) is a common representation that is often seen for +// quaternions. However, when we want to expand the elementary functions of +// quaternions in terms of real operations it is almost always easier to view +// them as real part (r) and imaginary vector part (v), +// i.e: r + xi + yj + zk = r + v; and so we diverge a little from the +// representation that is used in the documentation in other files and use this +// notation of quaternions in (internal) comments of the following functions. +// +// Quaternionic elementary functions have many similarities with elementary +// functions of complex numbers 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 complex numbers. + +import RealModule + +extension Quaternion: ElementaryFunctions { + // MARK: - exp-like functions + @inlinable + public static func exp(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // the `Real` operations `exp`, `cos` and `sin` (`let θ = ||v||`): + // + // ``` + // exp(r + v) = exp(r) exp(v) + // = exp(r) (cos(θ) + (v/θ) sin(θ)) + // ``` + // + // 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(||v||)` would not be). + guard q.isFinite else { return q } + let (â, θ) = q.imaginary.unitAxisAndLength + let rotation = Quaternion(halfAngle: θ, unitAxis: â) + // If real < log(greatestFiniteMagnitude), then exp(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) + return rotation.multiplied(by: halfScale).multiplied(by: halfScale) + } + return rotation.multiplied(by: .exp(q.real)) + } + + @inlinable + public static func expMinusOne(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // the `Real` operations `exp`, `cos` and `sin` (`let θ = ||v||`): + // + // ``` + // exp(r + v) - 1 = exp(r) exp(v) - 1 + // = exp(r) (cos(θ) + (v/θ) sin(θ)) - 1 + // = exp(r) cos(θ) + exp(r) (v/θ) sin(θ) - 1 + // = (exp(r) cos(θ) - 1) + exp(r) (v/θ) sin(θ) + // -------- u -------- + // ``` + // + // Note that the imaginary part is just the usual exp(x) sin(y); + // the only trick is computing the real part ("u"): + // + // ``` + // u = exp(r) cos(θ) - 1 + // = exp(r) cos(θ) - cos(θ) + cos(θ) - 1 + // = (exp(r) - 1) cos(θ) + (cos(θ) - 1) + // = expMinusOne(r) cos(θ) + cosMinusOne(θ) + // ``` + // + // See `expMinusOne` on complex numbers for error bounds. + guard q.isFinite else { return q } + let (â, θ) = q.imaginary.unitAxisAndLength + // If exp(q) is close to the overflow boundary, we don't need to + // worry about the "MinusOne" part of this function; we're just + // computing exp(q). (Even when θ is near a multiple of π/2, + // it can't be close enough to overcome the scaling from exp(r), + // so the -1 term is _always_ negligable). + guard q.real < RealType.log(.greatestFiniteMagnitude) - 1 else { + let halfScale = RealType.exp(q.real/2) + let rotation = Quaternion(halfAngle: θ, unitAxis: â) + return rotation.multiplied(by: halfScale).multiplied(by: halfScale) + } + return Quaternion( + real: RealType._mulAdd(.cos(θ), .expMinusOne(q.real), .cosMinusOne(θ)), + imaginary: â * .exp(q.real) * .sin(θ) + ) + } + + @inlinable + public static func cosh(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // trigonometric `Real` operations (`let θ = ||v||`): + // + // ``` + // cosh(q) = (exp(q) + exp(-q)) / 2 + // = cosh(r) cos(θ) + (v/θ) sinh(r) sin(θ) + // ``` + // + // Like exp, cosh is entire, so we do not need to worry about where + // branch cuts fall. Also like exp, cancellation never occurs in the + // evaluation of the naive expression, so all we need to be careful + // about is the behavior near the overflow boundary. + // + // Fortunately, if |r| >= -log(ulpOfOne), cosh(r) and sinh(r) are + // both just exp(|r|)/2, and we already know how to compute that. + // + // This function and sinh should stay in sync; if you make a + // modification here, you should almost surely make a parallel + // modification to sinh below. + guard q.isFinite else { return q } + let (â, θ) = q.imaginary.unitAxisAndLength + guard q.real.magnitude < -RealType.log(.ulpOfOne) else { + let rotation = Quaternion(halfAngle: θ, unitAxis: â) + let firstScale = RealType.exp(q.real.magnitude/2) + return rotation.multiplied(by: firstScale).multiplied(by: firstScale/2) + } + return Quaternion( + real: .cosh(q.real) * .cos(θ), + imaginary: â * .sinh(q.real) * .sin(θ) + ) + } + + @inlinable + public static func sinh(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // trigonometric `Real` operations (`let θ = ||v||`): + // + // ``` + // sinh(q) = (exp(q) - exp(-q)) / 2 + // = sinh(r) cos(θ) + (v/θ) cosh(r) sin(θ) + // ``` + guard q.isFinite else { return q } + let (â, θ) = q.imaginary.unitAxisAndLength + guard q.real.magnitude < -RealType.log(.ulpOfOne) else { + let rotation = Quaternion(halfAngle: θ, unitAxis: â) + let firstScale = RealType.exp(q.real.magnitude/2) + let secondScale = RealType(signOf: q.real, magnitudeOf: firstScale/2) + return rotation.multiplied(by: firstScale).multiplied(by: secondScale) + } + return Quaternion( + real: .sinh(q.real) * .cos(θ), + imaginary: â * .cosh(q.real) * .sin(θ) + ) + } + + @inlinable + public static func tanh(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // quaternionic `sinh` and `cosh` operations: + // + // ``` + // tanh(q) = sinh(q) / cosh(q) + // ``` + guard q.isFinite else { return q } + // Note that when |r| is larger than -log(.ulpOfOne), + // sinh(r + v) == ±cosh(r + v), so tanh(r + v) is just ±1. + guard q.real.magnitude < -RealType.log(.ulpOfOne) else { + return Quaternion( + real: RealType(signOf: q.real, magnitudeOf: 1), + imaginary: + RealType(signOf: q.components.x, magnitudeOf: 0), + RealType(signOf: q.components.y, magnitudeOf: 0), + RealType(signOf: q.components.z, magnitudeOf: 0) + ) + } + return sinh(q) / cosh(q) + } + + @inlinable + public static func cos(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // quaternionic `cosh` operations (`let θ = ||v||`): + // + // ``` + // cos(q) = cosh(q * (v/θ))) + // ``` + let (â,_) = q.imaginary.unitAxisAndLength + let p = Quaternion(imaginary: â) + return cosh(q * p) + } + + @inlinable + public static func sin(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // quaternionic `sinh` operations (`let θ = ||v||`): + // + // ``` + // sin(q) = -(v/θ) * sinh(q * (v/θ))) + // ``` + let (â,_) = q.imaginary.unitAxisAndLength + let p = Quaternion(imaginary: â) + return -p * sinh(q * p) + } + + @inlinable + public static func tan(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // quaternionic `tanh` operations (`let θ = ||v||`): + // + // ``` + // tan(q) = -(v/θ) * tanh(q * (v/θ))) + // ``` + let (â,_) = q.imaginary.unitAxisAndLength + let p = Quaternion(imaginary: â) + return -p * tanh(q * p) + } + + // MARK: - log-like functions + @inlinable + public static func log(_ q: Quaternion) -> Quaternion { + // If q is zero or infinite, the phase is undefined, so the result is + // the single exceptional value. + guard q.isFinite && !q.isZero else { return .infinity } + // Having eliminated non-finite values and zero, the imaginary part is + // straightforward: + let (â, θ) = q.imaginary.unitAxisAndLength + let argument = RealType.atan2(y: θ, x: q.real) + let imaginary = â * argument + // The real part of the result is trickier and we employ the same approach + // as we did for the complex numbers logarithm to improve the relative error + // bounds (`Complex.log`). There you may also find a lot more details to + // the following approach. + // + // To handle very large arguments without overflow, _rescale the problem. + // This is done by finding whichever part has greater magnitude, and + // dividing through by it. + let u = max(q.real.magnitude, θ) + let v = min(q.real.magnitude, θ) + // Now expand out log |w|: + // + // log |w| = log(u² + v²)/2 + // = log(u + log(onePlus: (u/v)²))/2 + // + // This handles overflow well, because log(u) is finite for every finite u, + // and we have 0 ≤ v/u ≤ 1. Unfortunately, it does not handle all points + // close to the unit circle so well, as cancellation might occur. + // + // We are not trying for sub-ulp accuracy, just a good relative error + // bound, so for our purposes it suffices to have log u dominate the + // result: + if u >= 1 || u >= RealType._mulAdd(u,u,v*v) { + let r = v / u + return Quaternion( + real: .log(u) + .log(onePlus: r*r)/2, + imaginary: imaginary + ) + } + // Here we're in the tricky case; cancellation is likely to occur. + // Instead of the factorization used above, we will want to evaluate + // log(onePlus: u² + v² - 1)/2. This all boils down to accurately + // evaluating u² + v² - 1. + let (a,b) = Augmented.product(u, u) + let (c,d) = Augmented.product(v, v) + 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. + s = (s + c) + e + b + d + return Quaternion(real: .log(onePlus: s)/2, imaginary: imaginary) + } + + @inlinable + public static func log(onePlus q: Quaternion) -> Quaternion { + // If either |r| or ||v||₁ is bounded away from the origin, we don't need + // any extra precision, and can just literally compute log(1+z). Note + // that this includes part of the sphere |1+q| = 1 where log(onePlus:) + // vanishes (where r <= -0.5), but on this portion of the sphere 1+r + // is always exact by Sterbenz' lemma, so as long as log( ) produces + // a good result, log(1+q) will too. + guard 2*q.real.magnitude < 1 && q.imaginary.oneNorm < 1 else { + return log(.one + q) + } + // q is in (±0.5, ±1), so we need to evaluate more carefully. + // The imaginary part is straightforward: + let argument = (.one + q).halfAngle + let (â,_) = q.imaginary.unitAxisAndLength + let imaginary = â * argument + // For the real part, we _could_ use the same approach that we do for + // log( ), but we'd need an extra-precise (1+r)², which can potentially + // be quite painful to calculate. Instead, we can use an approach that + // NevinBR suggested on the Swift forums for complex numbers: + // + // Re(log(1+q)) = (log(1+q) + log(1+q̅)) / 2 + // = log((1+q)(1+q̅)) / 2 + // = log(1 + q + q̅ + qq̅) / 2 + // = log(1 + 2r + r² + v²)) / 2 + // = log(1 + (2+r)r + v²)) / 2 + // = log(1 + (2+r)r + x² + y² + z²)) / 2 + // = log(onePlus: (2+r)r + x² + y² + z²) / 2 + // + // So now we need to evaluate (2+r)r + x² + y² + z² accurately. + // To do this, we employ augmented arithmetic + // (2+r)r + x² + y² + z² + // --↓-- + let rp2 = Augmented.sum(large: 2, small: q.real) // Known that 2 > |r| + var (head, δ) = Augmented.product(q.real, rp2.head) + var tail = δ + // head + x² + y² + z² + // ----↓---- + let x² = Augmented.product(q.imaginary.x, q.imaginary.x) + (head, δ) = Augmented.sum(head, x².head) + tail += (δ + x².tail) + // head + y² + z² + // ----↓---- + let y² = Augmented.product(q.imaginary.y, q.imaginary.y) + (head, δ) = Augmented.sum(head, y².head) + tail += (δ + y².tail) + // head + z² + // ----↓---- + let z² = Augmented.product(q.imaginary.z, q.imaginary.z) + (head, δ) = Augmented.sum(head, z².head) + tail += (δ + z².tail) + + let s = (head + tail).addingProduct(q.real, rp2.tail) + return Quaternion(real: .log(onePlus: s)/2, imaginary: imaginary) + } + + @inlinable + public static func acos(_ q: Quaternion) -> Quaternion { + let (â, θ) = (sqrt(1+q).conjugate * sqrt(1-q)).imaginary.unitAxisAndLength + return Quaternion( + real: 2*RealType.atan2(y: sqrt(1-q).real, x: sqrt(1+q).real), + imaginary: â * RealType.asinh(θ) + ) + } + + @inlinable + public static func asin(_ q: Quaternion) -> Quaternion { + let (â, θ) = (sqrt(1-q).conjugate * sqrt(1+q)).imaginary.unitAxisAndLength + return Quaternion( + real: RealType.atan2(y: q.real, x: (sqrt(1-q) * sqrt(1+q)).real), + imaginary: â * RealType.asinh(θ) + ) + } + + @inlinable + public static func atan(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // the quaternionic `atanh` operation (`let θ = ||v||`): + // + // ``` + // atan(q) = -(v/θ) * atanh(q * (v/θ)) + // ``` + let (â, _) = q.imaginary.unitAxisAndLength + let p = Quaternion(imaginary: â) + return -p * .atanh(q * p) + } + + @inlinable + public static func acosh(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // the quaternionic `acos` operation (`let θ = ||v||`): + // + // ``` + // acosh(q) = (v/θ) * acos(q) + // ``` + let (â,_) = q.imaginary.unitAxisAndLength + let p = Quaternion(imaginary: â) + return p * acos(q) + } + + @inlinable + public static func asinh(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // the quaternionic `asin` operation (`let θ = ||v||`): + // + // ``` + // sin(q) = -(v/θ) * asin(q * (v/θ))) + // ``` + let (â,_) = q.imaginary.unitAxisAndLength + let p = Quaternion(imaginary: â) + return -p * .asin(q * p) + } + + @inlinable + public static func atanh(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // the quaternionic `log` operation: + // + // ``` + // atanh(q) = (log(1 + q) - log(1 - q))/2 + // = (log(onePlus: q) - log(onePlus: -q))/2 + // ``` + (log(onePlus: q) - log(onePlus:-q))/2 + } + + // MARK: - pow-like functions + @inlinable + public static func pow(_ q: Quaternion, _ p: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // the quaternionic `exp` and `log` operations: + // + // ``` + // pow(q, p) = exp(log(pow(q, p))) + // = exp(p * log(q)) + // ``` + if q.isZero { return p.real > 0 ? zero : infinity } + return exp(p * log(q)) + } + + @inlinable + public static func pow(_ q: Quaternion, _ n: Int) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // the quaternionic `exp` and `log` operations: + // + // ``` + // pow(q, n) = exp(log(pow(q, n))) + // = exp(log(q) * n) + // ``` + if q.isZero { return n < 0 ? infinity : n == 0 ? one : zero } + // TODO: this implementation is not quite correct, because n may be + // rounded in conversion to RealType. This only effects very extreme + // cases, so we'll leave it alone for now. + return exp(log(q).multiplied(by: RealType(n))) + } + + @inlinable + public static func sqrt(_ q: Quaternion) -> Quaternion { + let lengthSquared = q.lengthSquared + if lengthSquared.isNormal { + // If |q|^2 doesn't overflow, then define s and t by (`let θ = ||v||`): + // + // s = sqrt((|q|+|r|) / 2) + // t = θ/2s + // + // If r is positive, the result is just w = (s, (v/θ) * t). If r is + // negative, the result is (|t|, (v/θ) * copysign(s, θ)) instead. + let (â, θ) = q.imaginary.unitAxisAndLength + let norm: RealType = .sqrt(lengthSquared) + let s: RealType = .sqrt((norm + q.real.magnitude) / 2) + let t: RealType = θ / (2*s) + if q.real.sign == .plus { + return Quaternion( + real: s, + imaginary: â * t) + } else { + return Quaternion( + real: t.magnitude, + imaginary: â * RealType(signOf: θ, magnitudeOf: s) + ) + } + } + // Handle edge cases: + guard !q.isZero else { + return Quaternion( + real: 0, + imaginary: + RealType(signOf: q.components.x, magnitudeOf: 0), + RealType(signOf: q.components.y, magnitudeOf: 0), + RealType(signOf: q.components.z, magnitudeOf: 0) + ) + } + guard q.isFinite else { return q } + // q is finite but badly-scaled. Rescale and replay by factoring out + // the larger of r and v. + let scale = q.magnitude + return Quaternion.sqrt(q.divided(by: scale)).multiplied(by: .sqrt(scale)) + } + + @inlinable + public static func root(_ q: Quaternion, _ n: Int) -> Quaternion { + // Mathematically, this operation can be expanded in terms of + // the quaternionic `exp` and `log` operations: + // + // ``` + // root(q, n) = q^(1/n) = exp(log(q^(1/n))) + // = exp(log(q) / n) + // ``` + guard !q.isZero else { return .zero } + // TODO: this implementation is not quite correct, because n may be + // rounded in conversion to RealType. This only effects very extreme + // cases, so we'll leave it alone for now. + return exp(log(q).divided(by: RealType(n))) + } +} + +extension SIMD3 where Scalar: FloatingPoint { + /// Returns the normalized axis and the length of this vector. + @_alwaysEmitIntoClient + fileprivate var unitAxisAndLength: (Self, Scalar) { + if self == .zero { + return (SIMD3( + Scalar(signOf: x, magnitudeOf: 0), + Scalar(signOf: y, magnitudeOf: 0), + Scalar(signOf: z, magnitudeOf: 0) + ), .zero) + } + return (self/length, length) + } +} diff --git a/Sources/QuaternionModule/Quaternion+Hashable.swift b/Sources/QuaternionModule/Quaternion+Hashable.swift new file mode 100644 index 00000000..7e7a2270 --- /dev/null +++ b/Sources/QuaternionModule/Quaternion+Hashable.swift @@ -0,0 +1,72 @@ +//===--- Quaternion+Hashable.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: Hashable { + /// Returns a Boolean value indicating whether two values are equal. + /// + /// - Important: + /// Quaternions are frequently used to represent 3D transformations. It's + /// important to be aware that, when used this way, any quaternion and its + /// negation represent the same transformation, but they do not compare + /// equal using `==` because they are not the same quaternion. You can + /// compare quaternions as 3D transformations using `equals(as3DTransform:)`. + @_transparent + public static func == (lhs: Quaternion, rhs: Quaternion) -> Bool { + // Identify all numbers with either component non-finite as a single "point at infinity". + guard lhs.isFinite || rhs.isFinite else { return true } + // For finite numbers, equality is defined componentwise. Cases where + // only one of lhs or rhs is infinite fall through to here as well, but this + // expression correctly returns false for them so we don't need to handle + // them explicitly. + return lhs.components == rhs.components + } + + /// Returns a Boolean value indicating whether the 3D transformation of the + /// two quaternions are equal. + /// + /// Use this method to test for equality of the 3D transformation properties + /// of quaternions; where for any quaternion `q`, its negation represent the + /// same 3D transformation; i.e. `q.equals(as3DTransform: q)` as well as + /// `q.equals(as3DTransform: -q)` are both `true`. + /// + /// - Parameter other: The value to compare. + /// - Returns: True if the 3D transformation of this quaternion equals `other`. + @_transparent + public func equals(as3DTransform other: Quaternion) -> Bool { + // Identify all numbers with either component non-finite as a single "point at infinity". + guard isFinite || other.isFinite else { return true } + // For finite numbers, equality is defined componentwise. Cases where only + // one of lhs or rhs is infinite fall through to here as well, but this + // expression correctly returns false for them so we don't need to handle + // them explicitly. + return components == other.components || components == -other.components + } + + @_transparent + public func hash(into hasher: inout Hasher) { + // There are two equivalence classes to which we owe special attention: + // All zeros should hash to the same value, regardless of sign, and all + // non-finite numbers should hash to the same value, regardless of + // representation. The correct behavior for zero falls out for free from + // the hash behavior of floating-point, but we need to use a + // representative member for any non-finite values. + // For any normal values we use the "canonical transform" representation, + // where real is always non-negative. This allows people who are using + // quaternions as rotations to get the expected semantics out of collections + // (while unfortunately producing some collisions for people who are not, + // but not in too catastrophic of a fashion). + if isFinite { + canonicalizedTransform.components.hash(into: &hasher) + } else { + hasher.combine(RealType.infinity) + } + } +} diff --git a/Sources/QuaternionModule/Quaternion+IntegerLiteral.swift b/Sources/QuaternionModule/Quaternion+IntegerLiteral.swift new file mode 100644 index 00000000..cafc3256 --- /dev/null +++ b/Sources/QuaternionModule/Quaternion+IntegerLiteral.swift @@ -0,0 +1,19 @@ +//===--- Quaternion+IntegerLiteral.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: ExpressibleByIntegerLiteral { + public typealias IntegerLiteralType = RealType.IntegerLiteralType + + @inlinable + public init(integerLiteral value: IntegerLiteralType) { + self.init(RealType(integerLiteral: value)) + } +} diff --git a/Sources/QuaternionModule/Quaternion+Numeric.swift b/Sources/QuaternionModule/Quaternion+Numeric.swift new file mode 100644 index 00000000..b08fec2a --- /dev/null +++ b/Sources/QuaternionModule/Quaternion+Numeric.swift @@ -0,0 +1,66 @@ +//===--- Quaternion+Numeric.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: Numeric { + @_transparent + public static func * (lhs: Quaternion, rhs: Quaternion) -> Quaternion { + + let rhsX = SIMD4(+rhs.components.w, +rhs.components.z, -rhs.components.y, +rhs.components.x) + let rhsY = SIMD4(-rhs.components.z, +rhs.components.w, +rhs.components.x, +rhs.components.y) + let rhsZ = SIMD4(+rhs.components.y, -rhs.components.x, +rhs.components.w, +rhs.components.z) + let rhsR = SIMD4(-rhs.components.x, -rhs.components.y, -rhs.components.z, +rhs.components.w) + + let x = (lhs.components * rhsX).sum() + let y = (lhs.components * rhsY).sum() + let z = (lhs.components * rhsZ).sum() + let r = (lhs.components * rhsR).sum() + + return Quaternion(from: SIMD4(x,y,z,r)) + } + + @_transparent + public static func *= (lhs: inout Quaternion, rhs: Quaternion) { + lhs = lhs * rhs + } + + /// The quaternion with specified real part and zero imaginary part. + /// + /// Equivalent to `Quaternion(RealType(real))`. + @inlinable + public init(_ real: Other) { + self.init(RealType(real)) + } + + /// The quaternion with specified real part and zero imaginary part, + /// if it can be constructed without rounding. + @inlinable + public init?(exactly real: Other) { + guard let real = RealType(exactly: real) else { return nil } + self.init(real) + } + + /// The ∞-norm of the value (`max(abs(r), abs(x), abs(y), abs(z))`). + /// + /// If you need the Euclidean norm (a.k.a. 2-norm) use the `length` or `lengthSquared` + /// properties instead. + /// + /// Edge cases: + /// - If `q` is not finite, `q.magnitude` is `.infinity`. + /// - If `q` is zero, `q.magnitude` is `0`. + /// - Otherwise, `q.magnitude` is finite and non-zero. + /// + /// See also `.length` and `.lengthSquared` + @_transparent + public var magnitude: RealType { + guard isFinite else { return .infinity } + return max(abs(components.max()), abs(components.min())) + } +} diff --git a/Sources/QuaternionModule/Quaternion+StringConvertible.swift b/Sources/QuaternionModule/Quaternion+StringConvertible.swift new file mode 100644 index 00000000..fbfa89f4 --- /dev/null +++ b/Sources/QuaternionModule/Quaternion+StringConvertible.swift @@ -0,0 +1,28 @@ +//===--- Quaternion+StringConvertible.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: CustomStringConvertible { + public var description: String { + guard isFinite else { return "inf" } + return "(\(components.w), \(components.x), \(components.y), \(components.z))" + } +} + +extension Quaternion: CustomDebugStringConvertible { + public var debugDescription: String { + let x = String(reflecting: components.x) + let y = String(reflecting: components.y) + let z = String(reflecting: components.z) + let r = String(reflecting: components.w) + return "Quaternion<\(RealType.self)>(\(r), \(x), \(y), \(z))" + } +} + diff --git a/Sources/QuaternionModule/Quaternion.swift b/Sources/QuaternionModule/Quaternion.swift new file mode 100644 index 00000000..a74affa5 --- /dev/null +++ b/Sources/QuaternionModule/Quaternion.swift @@ -0,0 +1,302 @@ +//===--- Quaternion.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 +// +//===----------------------------------------------------------------------===// + +import RealModule + +/// A quaternion represented by a real (or scalar) part and three imaginary (or vector) parts. +/// +/// TODO: introductory text on quaternions +/// +/// Implementation notes: +/// - +/// This type does not provide heterogeneous real/quaternion arithmetic, +/// not even the natural vector-space operations like real * quaternion. +/// There are two reasons for this choice: first, Swift broadly avoids +/// mixed-type arithmetic when the operation can be adequately expressed +/// by a conversion and homogeneous arithmetic. Second, with the current +/// typechecker rules, it would lead to undesirable ambiguity in common +/// expressions (see README.md for more details). +/// +/// `.magnitude` does not return the Euclidean norm; it uses the "infinity +/// norm" (`max(|r|,|x|,|y|,|z|)`) instead. There are two reasons for this +/// choice: first, it's simply faster to compute on most hardware. Second, +/// there exist values for which the Euclidean norm cannot be represented. +/// Using the infinity norm avoids this problem entirely without significant +/// downsides. You can access the Euclidean norm using the `length` property. +/// See `Complex` type of the swift-numerics package for additional details. +/// +/// Quaternions are frequently used to represent 3D transformations. It's +/// important to be aware that, when used this way, any quaternion and its +/// negation represent the same transformation, but they do not compare equal +/// using `==` because they are not the same quaternion. +/// You can compare quaternions as 3D transformations using `equals(as3DTransform:)`. +public struct Quaternion where RealType: Real & SIMDScalar { + + /// The components of the 4-dimensional vector space of the quaternion. + /// + /// Components are stored within a 4-dimensional SIMD vector with the + /// scalar part last, i.e. in the form of: + /// + /// xi + yj + zk + r // SIMD(x,y,z,r) + @usableFromInline @inline(__always) + internal var components: SIMD4 + + /// A quaternion constructed from given 4-dimensional SIMD vector. + /// + /// Creates a new quaternion by reading the values of the SIMD vector + /// as components of a quaternion with the scalar part last, i.e. in the form of: + /// + /// xi + yj + zk + r // SIMD(x,y,z,r) + @usableFromInline @inline(__always) + internal init(from components: SIMD4) { + self.components = components + } +} + +// MARK: - Basic Property +extension Quaternion { + /// The real part of this quaternion. + /// + /// If `q` is not finite, `q.real` is `.nan`. + public var real: RealType { + @_transparent + get { isFinite ? components[3] : .nan } + + @_transparent + set { components[3] = newValue } + } + + /// The imaginary part of this quaternion. + /// + /// If `q` is not finite, `q.imaginary` is `.nan` in all lanes. + public var imaginary: SIMD3 { + @_transparent + get { isFinite ? components[SIMD3(0,1,2)] : SIMD3(repeating: .nan) } + + @_transparent + set { + components = SIMD4(newValue, components[3]) + } + } + + /// The quaternion with the imaginary unit **i** one, i.e. `0 + i + 0j + 0k`. + /// + /// See also `.zero`, `.one`, `.j`, `.k` and `.infinity`. + @_transparent + public static var i: Quaternion { + Quaternion(imaginary: SIMD3(1,0,0)) + } + + /// The quaternion with the imaginary unit **j** one, i.e. `0 + 0i + j + 0k`. + /// + /// See also `.zero`, `.one`, `.i`, `.k` and `.infinity`. + @_transparent + public static var j: Quaternion { + Quaternion(imaginary: SIMD3(0,1,0)) + } + + /// The quaternion with the imaginary unit **k** one, i.e. `0 + 0i + 0j + k`. + /// + /// See also `.zero`, `.one`, `.i`, `.j` and `.infinity`. + @_transparent + public static var k: Quaternion { + Quaternion(imaginary: SIMD3(0,0,1)) + } + + /// The point at infinity. + /// + /// See also `.zero`, `.one`, `.i`, `.j` and `.k`. + @_transparent + public static var infinity: Quaternion { + Quaternion(.infinity) + } + + /// True if this value is finite. + /// + /// A quaternion is finite if neither component is an infinity or nan. + /// + /// See also `.isNormal`, `.isSubnormal`, `.isZero`, `.isReal`, `.isPure`. + @_transparent + public var isFinite: Bool { + return components.x.isFinite + && components.y.isFinite + && components.z.isFinite + && components.w.isFinite + } + + /// True if this value is normal. + /// + /// A quaternion is normal if it is finite and *any* of its real or imaginary components + /// are normal. A floating-point number representing one of the components is normal + /// if its exponent allows a full-precision representation. + /// + /// See also `.isFinite`, `.isSubnormal`, `.isZero`, `.isReal`, `.isPure`. + @_transparent + public var isNormal: Bool { + return isFinite && ( + components.x.isNormal || + components.y.isNormal || + components.z.isNormal || + components.w.isNormal + ) + } + + /// True if this value is subnormal. + /// + /// A quaternion is subnormal if it is finite, not normal, and not zero. When the result of a + /// computation is subnormal, underflow has occurred and the result generally does not have full + /// precision. + /// + /// See also `.isFinite`, `.isNormal`, `.isZero`, `.isReal`, `.isPure`. + @_transparent + public var isSubnormal: Bool { + isFinite && !isNormal && !isZero + } + + /// True if this value is zero. + /// + /// A quaternion is zero if the real and *all* imaginary components are zero. + /// + /// See also `.isFinite`, `.isNormal`, `.isSubnormal`, `.isReal`, `.isPure`. + @_transparent + public var isZero: Bool { + components == .zero + } + + /// True if this quaternion is real. + /// + /// A quaternion is real if *all* imaginary components are zero. + /// + /// See also `.isFinite`, `.isNormal`, `.isSubnormal`, `.isZero`, `.isPure`. + @_transparent + public var isReal: Bool { + imaginary == .zero + } + + /// True if this quaternion is pure. + /// + /// A quaternion is pure if the real component is zero. + /// + /// See also `.isFinite`, `.isNormal`, `.isSubnormal`, `.isZero`, `.isReal`. + @_transparent + public var isPure: Bool { + real.isZero + } + + /// A "canonical" representation of the value. + /// + /// For normal quaternion instances with a RealType conforming to + /// BinaryFloatingPoint (the common case), the result is simply this value + /// unmodified. For zeros, the result has the representation (+0, +0, +0, +0). + /// For infinite values, the result has the representation (+inf, +0, +0, +0). + /// + /// If the RealType admits non-canonical representations, the x, y, z and r + /// components are canonicalized in the result. + /// + /// This is mainly useful for interoperation with other languages, where + /// you may want to reduce each equivalence class to a single representative + /// before passing across language boundaries, but it may also be useful + /// for some serialization tasks. It's also a useful implementation detail for + /// some primitive operations. + /// + /// See also `.canonicalizedTransform`. + @_transparent + public var canonicalized: Self { + guard !isZero else { return .zero } + guard isFinite else { return .infinity } + return self.multiplied(by: 1) + } + + /// A "canonical transformation" representation of the value. + /// + /// For normal quaternion instances with a RealType conforming to + /// BinaryFloatingPoint (the common case) and a non-negative real component, + /// the result is simply this value unmodified. For instances with a negative + /// real component, the result is this quaternion negated -(r, x, y, z); so + /// the real component is always positive. + /// For zeros, the result has the representation (+0, +0, +0, +0). For + /// infinite values, the result has the representation (+inf, +0, +0, +0). + /// + /// If the RealType admits non-canonical representations, the x, y, z and r + /// components are canonicalized in the result. + /// + /// See also `.canonicalized`. + @_transparent + public var canonicalizedTransform: Self { + let canonical = canonicalized + if canonical.real.sign == .plus { return canonical } + return -canonical + } +} + +// MARK: - Additional Initializers +extension Quaternion { + /// The quaternion with specified real part and zero imaginary part. + /// + /// Equivalent to `Quaternion(real: real, imaginary: SIMD3(repeating: 0))`. + @inlinable + public init(_ real: RealType) { + self.init(real: real, imaginary: SIMD3(repeating: 0)) + } + + /// The quaternion with specified imaginary part and zero real part. + /// + /// Equivalent to `Quaternion(real: 0, imaginary: imaginary)`. + @inlinable + public init(imaginary: SIMD3) { + self.init(real: 0, imaginary: imaginary) + } + + /// The quaternion with specified imaginary part and zero real part. + /// + /// Equivalent to `Quaternion(real: 0, imaginary: imaginary)`. + @inlinable + public init(imaginary x: RealType, _ y: RealType, _ z: RealType) { + self.init(imaginary: SIMD3(x, y, z)) + } + + /// The quaternion with specified real part and imaginary parts. + @inlinable + public init(real: RealType, imaginary: SIMD3) { + self.init(from: SIMD4(imaginary, real)) + } + + /// The quaternion with specified real part and imaginary parts. + @inlinable + public init(real: RealType, imaginary x: RealType, _ y: RealType, _ z: RealType) { + self.init(real: real, imaginary: SIMD3(x, y, z)) + } +} + +extension Quaternion where RealType: BinaryFloatingPoint { + /// `other` rounded to the nearest representable value of this type. + @inlinable + public init(_ other: Quaternion) { + self.init(from: SIMD4( + RealType(other.components.x), + RealType(other.components.y), + RealType(other.components.z), + RealType(other.components.w) + )) + } + + /// `other`, if it can be represented exactly in this type; otherwise `nil`. + @inlinable + public init?(exactly other: Quaternion) { + guard + let x = RealType(exactly: other.components.x), + let y = RealType(exactly: other.components.y), + let z = RealType(exactly: other.components.z), + let r = RealType(exactly: other.components.w) + else { return nil } + self.init(from: SIMD4(x, y, z, r)) + } +} diff --git a/Sources/QuaternionModule/README.md b/Sources/QuaternionModule/README.md new file mode 100644 index 00000000..cebccc8f --- /dev/null +++ b/Sources/QuaternionModule/README.md @@ -0,0 +1,30 @@ +# Quaternion + +This module provides a `Quaternion` type generic over an underlying `RealType`: + +```swift +1> import QuaternionModule +2> let q = Quaternion(1, (1,1,1)) // q = 1 + i + j + k +``` + +The usual arithmetic operators are provided for Quaternions, many useful properties, plus conformances to the +obvious usual protocols: `Equatable`, `Hashable`, `Codable` (if the underlying `RealType` is), and `AlgebraicField` +(hence also `AdditiveArithmetic` and `SignedNumeric`). + +### Dependencies: +- `RealModule`. + +### The magnitude property +The `Numeric` protocol requires a `.magnitude` property, but (deliberately) does not fully specify the semantics. +The most obvious choice for `Quaternion` would be to use the Euclidean norm (aka the "2-norm", given by `sqrt(real*real + i*i + k*k + j*j)`). +However, in practice there are good reasons to use something else instead: + +- The 2-norm requires special care to avoid spurious overflow/underflow, but the naive expressions for the 1-norm ("taxicab norm") or ∞-norm ("sup norm") are always correct. +- Even when care is used, near the overflow boundary the 2-norm and the 1-norm are not representable. + As an example, consider `q = Quaternion(big, (big, big, big))`, where `big` is `Double.greatestFiniteMagnitude`. The 1-norm and 2-norm of `q` both overflow (the 1-norm would be `4*big`, and the 2-norm would be `sqrt(4)*big`, neither of which are representable as `Double`), but the ∞-norm is always equal to either `real`, `i`, `j` or `k`, so it is guaranteed to be representable. +Because of this, the ∞-norm is the obvious alternative; it gives the nicest API surface. +- If we consider the magnitude of more exotic types, like operators, the 1-norm and ∞-norm are significantly easier to compute than the 2-norm (O(n) vs. "no closed form expression, but O(n^3) iterative methods"), so it is nice to establish a precedent of `.magnitude` binding one of these cheaper-to-compute norms. +- The ∞-norm is heavily used in other computational libraries; for example, it is used by the `izamax` and `icamax` functions in BLAS. + +The 2-norm still needs to be available, of course, because sometimes you need it. +This functionality is accessed via the `.length` and `.lengthSquared` properties. diff --git a/Sources/QuaternionModule/Scale.swift b/Sources/QuaternionModule/Scale.swift new file mode 100644 index 00000000..25c15346 --- /dev/null +++ b/Sources/QuaternionModule/Scale.swift @@ -0,0 +1,39 @@ +//===--- Scale.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 +// +//===----------------------------------------------------------------------===// + +// Policy: deliberately not using the * and / operators for these at the +// moment, because then there's an ambiguity in expressions like 2*p; is +// that Quaternion(2) * p or is it RealType(2) * p? This is especially +// problematic in type inference: suppose we have: +// +// let a: RealType = 1 +// let b = 2*a +// +// what is the type of b? If we don't have a type context, it's ambiguous. +// If we have a Quaternion type context, then b will be inferred to have type +// Quaternion! Obviously, that doesn't help anyone. +// +// TODO: figure out if there's some way to avoid these surprising results +// and turn these into operators if/when we have it. +// (https://github.com/apple/swift-numerics/issues/12) +extension Quaternion { + /// `self` scaled by `scalar`. + @usableFromInline @_transparent + internal func multiplied(by scalar: RealType) -> Quaternion { + Quaternion(from: components * scalar) + } + + /// `self` unscaled by `scalar`. + @usableFromInline @_transparent + internal func divided(by scalar: RealType) -> Quaternion { + Quaternion(from: components / scalar) + } +} diff --git a/Sources/QuaternionModule/Transformation.md b/Sources/QuaternionModule/Transformation.md new file mode 100644 index 00000000..1dfa2cb0 --- /dev/null +++ b/Sources/QuaternionModule/Transformation.md @@ -0,0 +1,25 @@ +# Transformation + +In computer science, quaternions are frequently used to represent three-dimensional rotations; as quaternions have some dvantages over other representations][advantages]. + +`Transformation.swift` encapsulates an API to interact with the three-dimensional transformation properties of quaternions. It provides conversions to and from other rotation representations, namely [*Angle-Axis*][angle_axis_wiki], [*Rotation Vector*][rotation_vector_wiki] and [*Polar decomposition*][polar_wiki], as well as it provides methods to directly transform arbitrary vectors by quaternions. + +## Policies + +- zero and non-finite quaternions have indeterminate transformation properties and can not be converted to other representations. Thus, + + - The `angle` of `.zero` or `.infinity` is `RealType.nan`. + - The `axis` of `.zero` or `.infinity` is `RealType.nan` in all lanes. + - The `rotationVector` of `.zero` or `.infinity` is `RealType.nan` in all lanes. + - The polar `phase` of `.zero` or `.infinity` is `RealType.nan` + +- Quaternions with an `angle` of `.zero` have an indeterminate rotation axis. Thus, + + - the `axis` of `angle == .zero` is `RealType.nan` in all lanes. + - the `rotationVector` of `angle == .zero` is `RealType.nan` in all lanes. + + +[advantages]: https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Advantages_of_quaternions +[angle_axis_wiki]: https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Recovering_the_axis-angle_representation +[polar_wiki]: https://en.wikipedia.org/wiki/Polar_decomposition#Quaternion_polar_decomposition +[rotation_vector_wiki]: https://en.wikipedia.org/wiki/Axis–angle_representation#Rotation_vector diff --git a/Sources/QuaternionModule/Transformation.swift b/Sources/QuaternionModule/Transformation.swift new file mode 100644 index 00000000..0c64e7a2 --- /dev/null +++ b/Sources/QuaternionModule/Transformation.swift @@ -0,0 +1,269 @@ +//===--- Transformation.swift ---------------------------------*- swift -*-===// +// +// This source file is part of the Swift Numerics open source project +// +// Copyright (c) 2020 - 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 { + /// The [rotation angle][wiki] of the Angle-Axis representation. + /// + /// Returns the rotation angle about the rotation *axis* in radians + /// within *[0, 2π]* range. + /// + /// Edge cases: + /// - If the quaternion is zero or non-finite, angle is `nan`. + /// + /// See also `.axis`, `.angleAxis`, `.rotationVector`, + /// `init(length:angle:axis:)` and `init(rotation:)`. + /// + /// [wiki]: https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Recovering_the_axis-angle_representation + @inlinable + public var angle: RealType { + 2 * halfAngle + } + + /// The [rotation axis][wiki] of the Angle-Axis representation. + /// + /// Returns the *(x,y,z)* rotation axis encoded in the quaternion + /// as SIMD3 vector of unit length. + /// + /// Edge cases: + /// - If the quaternion is zero or non-finite, axis is `nan` in all lanes. + /// - If the rotation angle is zero, axis is `nan` in all lanes. + /// + /// See also `.angle`, `.angleAxis`, `.rotationVector`, + /// `init(length:angle:axis:)` and `init(rotation:)`. + /// + /// [wiki]: https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Recovering_the_axis-angle_representation + @inlinable + public var axis: SIMD3 { + guard isFinite, imaginary != .zero else { return SIMD3(repeating: .nan) } + + // 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).axis } + return imaginary / .sqrt(lenSq) + } + + /// The [Angle-Axis][wiki] representation. + /// + /// Returns the length of the quaternion, the rotation angle in radians + /// within *[0, 2π]* and the rotation axis as SIMD3 vector of unit length. + /// + /// Edge cases: + /// - If the quaternion is zero or non-finite, angle and axis are `nan`. + /// - If the angle is zero, axis is `nan` in all lanes. + /// + /// See also `.angle`, `.axis`, `.rotationVector`, `init(length:angle:axis:)` + /// and `init(rotation:)`. + /// + /// [wiki]: https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Recovering_the_axis-angle_representation + public var angleAxis: (length: RealType, angle: RealType, axis: SIMD3) { + (length, angle, axis) + } + + /// The [rotation vector][rotvector]. + /// + /// A rotation vector is a vector of same direction as the rotation axis, + /// whose length is the rotation angle of an Angle-Axis representation. It + /// is effectively the product of multiplying the rotation `axis` by the + /// rotation `angle`. Rotation vectors are often called "scaled axis" — this + /// is a different name for the same concept. + /// + /// Edge cases: + /// - If the quaternion is zero or non-finite, the rotation vector is `nan` + /// in all lanes. + /// - If the rotation angle is zero, the rotation vector is `nan` + /// in all lanes. + /// + /// See also `.angle`, `.axis`, `.angleAxis`, `init(length:angle:axis:)` + /// and `init(rotation:)`. + /// + /// [rotvector]: https://en.wikipedia.org/wiki/Axis–angle_representation#Rotation_vector + @_transparent + public var rotationVector: SIMD3 { + axis * angle + } + + /// Creates a unit quaternion specified with [Angle-Axis][wiki] values. + /// + /// Angle-Axis is a representation of a three-dimensional rotation using two + /// different quantities: an angle describing the magnitude of rotation, and + /// a vector of unit length indicating the axis direction to rotate along. + /// The optional length parameter scales the quaternion after the conversion. + /// + /// This initializer reads given `length`, `angle` and `axis` values and + /// creates a quaternion of equal rotation properties and of specified length + /// using the following equation: + /// + /// Q = (cos(angle/2), axis * sin(angle/2)) * length + /// + /// If `length` is not specified, it defaults to *1*; and the final + /// quaternion is of unit 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, angle: θ, axis: axis) == -Quaternion(length: r, angle: θ, axis: axis) + /// ``` + /// - For any `θ` and any `axis`, even `.infinity` or `.nan`: + /// ``` + /// Quaternion(length: .zero, angle: θ, axis: axis) == .zero + /// ``` + /// - For any `θ` and any `axis`, even `.infinity` or `.nan`: + /// ``` + /// Quaternion(length: .infinity, angle: θ, axis: axis) == .infinity + /// ``` + /// - Otherwise, `θ` must be finite, or a precondition failure occurs. + /// + /// See also `.angle`, `.axis`, `.angleAxis`, `.rotationVector` + /// and `init(rotation:)`. + /// + /// - Parameter length: The length of the quaternion. Defaults to `1`. + /// - Parameter angle: The rotation angle about the rotation axis in radians. + /// - Parameter axis: The rotation axis. Must be of unit length. + /// + /// [wiki]: https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Recovering_the_axis-angle_representation + @inlinable + public init(length: RealType = 1, angle: RealType, axis: SIMD3) { + guard !length.isZero, length.isFinite else { + self = Quaternion(length) + return + } + + // Length is finite and non-zero, therefore + // 1. `angle` 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( + angle.isFinite, + "Either angle 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: angle/2, unitAxis: axis).multiplied(by: length) + } + + /// Creates a unit quaternion specified with given [rotation vector][wiki]. + /// + /// A rotation vector is a vector of same direction as the rotation axis, + /// whose length is the rotation angle of an Angle-Axis representation. It + /// is effectively the product of multiplying the rotation `axis` by the + /// rotation `angle`. + /// + /// This initializer reads the angle and axis values of given rotation vector + /// and creates a quaternion of equal rotation properties using the following + /// equation: + /// + /// Q = (cos(angle/2), axis * sin(angle/2)) + /// + /// Rotation vectors are sometimes referred to as *scaled axis* — this is a + /// different name for the same concept. + /// + /// The final quaternion is of unit length. + /// + /// Edge cases: + /// - If `vector` is `.zero`, the quaternion is `.zero`: + /// ``` + /// Quaternion(rotation: .zero) == .zero + /// ``` + /// - If `vector` is `.infinity` or `-.infinity`, the quaternion is `.infinity`: + /// ``` + /// Quaternion(rotation: -.infinity) == .infinity + /// ``` + /// + /// See also `.angle`, `.axis`, `.angleAxis`, `.rotationVector` + /// and `init(length:angle:axis:)`. + /// + /// - Parameter vector: The rotation vector. + /// + /// [wiki]: https://en.wikipedia.org/wiki/Axis–angle_representation#Rotation_vector + @inlinable + public init(rotation vector: SIMD3) { + let angle: RealType = .sqrt(vector.lengthSquared) + if !angle.isZero, angle.isFinite { + self = Quaternion(halfAngle: angle/2, unitAxis: vector/angle) + } else { + self = Quaternion(angle) + } + } + + /// Transforms a vector by this quaternion. + /// + /// Quaternions are frequently used to represent three-dimensional + /// transformations, and thus are used to transform vectors in + /// three-dimensional space. The transformation of an arbitrary vector + /// by a quaternion is known as an action. + /// + /// The canonical way of transforming an arbitrary three-dimensional vector + /// `v` by a quaternion `q` is given by the following [formula][wiki] + /// + /// p' = qpq⁻¹ + /// + /// where `p` is a *pure* quaternion (`real == .zero`) with imaginary part + /// equal to vector `v`, and where `p'` is another pure quaternion with + /// imaginary part equal to the transformed vector `v'`. The implementation + /// uses this formular but boils down to a simpler and faster implementation + /// as `p` is known to be pure and `q` is assumed to have unit length – which + /// allows for simplification. + /// + /// - Note: This method assumes this quaternion is of unit length. + /// + /// Edge cases: + /// - For any quaternion `q`, even `.zero` or `.infinity`, if `vector` is + /// `.infinity` or `-.infinity` in any of the lanes or all, the returning + /// vector is `.infinity` in all lanes: + /// ``` + /// SIMD3(-.infinity,0,0) * q == SIMD3(.infinity,.infinity,.infinity) + /// ``` + /// - For any quaternion `q`, even `.zero` or `.infinity`, if `vector` is + /// `.zero`, the returning vector is also `.zero`. + /// ``` + /// SIMD3(0,0,0) * q == .zero + /// ``` + /// + /// - Parameter vector: A vector to rotate by this quaternion + /// - Returns: The vector rotated by this quaternion + /// + /// [wiki]: https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Using_quaternion_as_rotations + @inlinable + public func act(on vector: SIMD3) -> SIMD3 { + guard vector.isFinite else { return .infinity } + guard vector != .zero else { return .zero } + + // The following expression have been split up so the type-checker + // can resolve them in a reasonable time. + let p1 = vector * (real*real - imaginary.lengthSquared) + let p2 = imaginary * imaginary.dot(vector) + let p3 = imaginary.cross(vector) * real + let rotatedVector = p1 + (p2 + p3) * 2 + + // If the rotation computes without over/underflow, everything is fine + // and the result is correct. If not, we have to do the computation + // carefully and first unscale the vector, rotate it again and then + // rescale the vector + if + (rotatedVector.x.isNormal || rotatedVector.x.isZero) && + (rotatedVector.y.isNormal || rotatedVector.y.isZero) && + (rotatedVector.z.isNormal || rotatedVector.z.isZero) + { + return rotatedVector + } + let scale = max(abs(vector.max()), abs(vector.min())) + return act(on: vector/scale) * scale + } +} diff --git a/Tests/QuaternionTests/ArithmeticTests.swift b/Tests/QuaternionTests/ArithmeticTests.swift new file mode 100644 index 00000000..37779c9e --- /dev/null +++ b/Tests/QuaternionTests/ArithmeticTests.swift @@ -0,0 +1,67 @@ +//===--- ArithmeticTests.swift --------------------------------*- swift -*-===// +// +// This source file is part of the Swift Numerics open source project +// +// Copyright (c) 2019 - 2020 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 + +@testable import QuaternionModule + +final class ArithmeticTests: XCTestCase { + + func testMultiplication(_ type: T.Type) { + for value: T in [-3, -2, -1, +1, +2, +3] { + let q = Quaternion(real: value, imaginary: value, value, value) + XCTAssertEqual(q * .one, q) + XCTAssertEqual(q * 1, q) + XCTAssertEqual(1 * q, q) + } + } + + func testMultiplication() { + testMultiplication(Float32.self) + testMultiplication(Float64.self) + } + + func testDivision(_ type: T.Type) { + for value: T in [-3, -2, -1, +1, +2, +3] { + let q = Quaternion(real: value, imaginary: value, value, value) + XCTAssertEqual(q/q, .one) + XCTAssertEqual(0/q, .zero) + + for q2: Quaternion in [-3, -2, -1, 0, +1, +2, +3] { + XCTAssertEqual(q2/q, q2 * q.reciprocal!) + } + + for q2: T in [-3, -2, -1, +1, +2, +3] { + XCTAssertEqual(q.divided(by: q2), q.multiplied(by: 1/q2)) + } + } + } + + func testDivision() { + testDivision(Float32.self) + testDivision(Float64.self) + } + + func testDivisionByZero(_ type: T.Type) { + XCTAssertFalse((Quaternion(real: 0, imaginary: 0, 0, 0) / Quaternion(real: 0, imaginary: 0, 0, 0)).isFinite) + XCTAssertFalse((Quaternion(real: 1, imaginary: 1, 1, 1) / Quaternion(real: 0, imaginary: 0, 0, 0)).isFinite) + XCTAssertFalse((Quaternion.infinity / Quaternion(real: 0, imaginary: 0, 0, 0)).isFinite) + XCTAssertFalse((Quaternion.i / Quaternion(real: 0, imaginary: 0, 0, 0)).isFinite) + XCTAssertFalse((Quaternion.j / Quaternion(real: 0, imaginary: 0, 0, 0)).isFinite) + XCTAssertFalse((Quaternion.k / Quaternion(real: 0, imaginary: 0, 0, 0)).isFinite) + } + + func testDivisionByZero() { + testDivisionByZero(Float32.self) + testDivisionByZero(Float64.self) + } +} diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift new file mode 100644 index 00000000..ac026a82 --- /dev/null +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -0,0 +1,936 @@ +//===--- ElementaryFunctionTests.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 +// +//===----------------------------------------------------------------------===// + +import XCTest +import RealModule +import _TestSupport + +@testable import QuaternionModule + +final class ElementaryFunctionTests: XCTestCase { + + // MARK: - exp-like functions + + func testExp(_ type: T.Type) { + // exp(0) = 1 + XCTAssertEqual(1, Quaternion.exp(Quaternion(real: .zero, imaginary: .zero))) + XCTAssertEqual(1, Quaternion.exp(Quaternion(real:-.zero, imaginary: .zero))) + XCTAssertEqual(1, Quaternion.exp(Quaternion(real:-.zero, imaginary:-.zero))) + XCTAssertEqual(1, Quaternion.exp(Quaternion(real: .zero, imaginary:-.zero))) + // exp is the identity at infinity. + XCTAssertFalse(Quaternion.exp(Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.exp(Quaternion(real: .infinity, imaginary:-.ulpOfOne)).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.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] = (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 testExpMinusOne(_ type: T.Type) { + // expMinusOne(0) = 0 + XCTAssert(Quaternion.expMinusOne(Quaternion(real: .zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.expMinusOne(Quaternion(real:-.zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.expMinusOne(Quaternion(real:-.zero, imaginary:-.zero)).isZero) + XCTAssert(Quaternion.expMinusOne(Quaternion(real: .zero, imaginary:-.zero)).isZero) + // expMinusOne is the identity at infinity + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .infinity, imaginary:-.ulpOfOne)).isFinite) + // Near-overflow test, same as exp() above. + let x = T.log(.greatestFiniteMagnitude) + T.log(9/8) + let huge = Quaternion.expMinusOne(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 small values, expMinusOne should be approximately the identity. + var g = SystemRandomNumberGenerator() + let small = T.ulpOfOne + for _ in 0 ..< 100 { + let q = Quaternion( + real: T.random(in: -small ... small, using: &g), + imaginary: + T.random(in: -small ... small, using: &g), + T.random(in: -small ... small, using: &g), + T.random(in: -small ... small, using: &g) + ) + XCTAssert(q.isApproximatelyEqual(to: Quaternion.expMinusOne(q), relativeTolerance: 16 * .ulpOfOne)) + } + } + + func testCosh(_ type: T.Type) { + // cosh(0) = 1 + XCTAssertEqual(1, Quaternion.cosh(Quaternion(real: .zero, imaginary: .zero))) + XCTAssertEqual(1, Quaternion.cosh(Quaternion(real:-.zero, imaginary: .zero))) + XCTAssertEqual(1, Quaternion.cosh(Quaternion(real:-.zero, imaginary:-.zero))) + XCTAssertEqual(1, Quaternion.cosh(Quaternion(real: .zero, imaginary:-.zero))) + // cosh is the identity at infinity. + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.cosh(Quaternion(real: .infinity, imaginary:-.ulpOfOne)).isFinite) + // Near-overflow test, same as exp() above, but it happens later, because + // for large x, cosh(x + v) ~ exp(x + v)/2. + let x = T.log(.greatestFiniteMagnitude) + T.log(18/8) + let mag = T.greatestFiniteMagnitude/T.sqrt(2) * (9/8) + var huge = Quaternion.cosh(Quaternion(real: x, imaginary: SIMD3(.pi/4, 0, 0))) + XCTAssert(huge.real.isApproximatelyEqual(to: mag)) + XCTAssert(huge.imaginary.x.isApproximatelyEqual(to: mag)) + XCTAssertEqual(huge.imaginary.y, .zero) + XCTAssertEqual(huge.imaginary.z, .zero) + huge = Quaternion.cosh(Quaternion(real: -x, imaginary: SIMD3(.pi/4, 0, 0))) + XCTAssert(huge.real.isApproximatelyEqual(to: mag)) + XCTAssert(huge.imaginary.x.isApproximatelyEqual(to: mag)) + XCTAssertEqual(huge.imaginary.y, .zero) + XCTAssertEqual(huge.imaginary.z, .zero) + } + + func testSinh(_ type: T.Type) { + // sinh(0) = 0 + XCTAssert(Quaternion.sinh(Quaternion(real: .zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.sinh(Quaternion(real:-.zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.sinh(Quaternion(real:-.zero, imaginary:-.zero)).isZero) + XCTAssert(Quaternion.sinh(Quaternion(real: .zero, imaginary:-.zero)).isZero) + // sinh is the identity at infinity. + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.sinh(Quaternion(real: .infinity, imaginary:-.ulpOfOne)).isFinite) + // Near-overflow test, same as exp() above, but it happens later, because + // for large x, sinh(x + v) ~ ±exp(x + v)/2. + let x = T.log(.greatestFiniteMagnitude) + T.log(18/8) + let mag = T.greatestFiniteMagnitude/T.sqrt(2) * (9/8) + var huge = Quaternion.sinh(Quaternion(real: x, imaginary: SIMD3(.pi/4, 0, 0))) + XCTAssert(huge.real.isApproximatelyEqual(to: mag)) + XCTAssert(huge.imaginary.x.isApproximatelyEqual(to: mag)) + XCTAssertEqual(huge.imaginary.y, .zero) + XCTAssertEqual(huge.imaginary.z, .zero) + huge = Quaternion.sinh(Quaternion(real: -x, imaginary: SIMD3(.pi/4, 0, 0))) + 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 + // cosh² - sinh² ≈ 1. Note that this test would break down due to + // catastrophic cancellation as we get further away from the origin. + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: -2 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + let c = Quaternion.cosh(q) + let s = Quaternion.sinh(q) + XCTAssert((c*c - s*s).isApproximatelyEqual(to: .one)) + } + } + + func testCosSin(_ type: T.Type) { + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: -2 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + // For randomly-chosen well-scaled finite values, we expect to have + // cos² + sin² ≈ 1 + let s = Quaternion.sin(q) + let c = Quaternion.cos(q) + XCTAssert((c*c + s*s).isApproximatelyEqual(to: .one)) + } + } + + // MARK: - log-like functions + + func testLog(_ type: T.Type) { + // log(0) = infinity + XCTAssertFalse(Quaternion.log(Quaternion(real: .zero, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.zero, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.zero, imaginary:-.zero)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .zero, imaginary:-.zero)).isFinite) + // log is the identity at infinity + XCTAssertFalse(Quaternion.log(Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .infinity, imaginary:-.ulpOfOne)).isFinite) + // For randomly-chosen well-scaled finite values, we expect to have + // log(exp(q)) ≈ q + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: -2 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + XCTAssert(q.isApproximatelyEqual(to: .exp(.log(q)))) + } + } + + func testLogOnePlus(_ type: T.Type) { + // log(onePlus: 0) = 0 + XCTAssert(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.log(onePlus: Quaternion(real:-.zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.log(onePlus: Quaternion(real:-.zero, imaginary:-.zero)).isZero) + XCTAssert(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary:-.zero)).isZero) + // log(onePlus:) is the identity at infinity. + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary:-.ulpOfOne)).isFinite) + // For randomly-chosen well-scaled finite values, we expect to have + // log(onePlus: expMinusOne(q)) ≈ q + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: -2 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + XCTAssert(q.isApproximatelyEqual(to: .expMinusOne(.log(onePlus: q)))) + } + } + + func testAcos(_ type: T.Type) { + // acos(1) = 0 + XCTAssert(Quaternion.acos(1).isZero) + // acos(0) = π/2 + XCTAssert(Quaternion.acos(0).real.isApproximatelyEqual(to: .pi/2)) + XCTAssertEqual(Quaternion.acos(0).imaginary, .zero) + // acos(-1) = π + XCTAssert(Quaternion.acos(-1).real.isApproximatelyEqual(to: .pi)) + XCTAssertEqual(Quaternion.acos(-1).imaginary, .zero) + // acos is the identity at infinity. + XCTAssertFalse(Quaternion.acos(Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.acos(Quaternion(real: .infinity, imaginary:-.ulpOfOne)).isFinite) + // For randomly-chosen well-scaled finite values, we expect to have + // cos(acos(q)) ≈ q and acos(q) ≈ π - acos(-q) + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: -2 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + let p = Quaternion.acos(q) + XCTAssert(Quaternion.cos(p).isApproximatelyEqual(to: q)) + XCTAssert(p.isApproximatelyEqual(to: Quaternion(.pi) - .acos(-q))) + } + } + + func testAsin(_ type: T.Type) { + // asin(1) = π/2 + XCTAssert(Quaternion.asin(1).real.isApproximatelyEqual(to: .pi/2)) + XCTAssertEqual(Quaternion.asin(1).imaginary, .zero) + // asin(0) = 0 + XCTAssert(Quaternion.asin(0).isZero) + // asin(-1) = -π/2 + XCTAssert(Quaternion.asin(-1).real.isApproximatelyEqual(to: -.pi/2)) + XCTAssertEqual(Quaternion.asin(-1).imaginary, .zero) + // asin is the identity at infinity. + XCTAssertFalse(Quaternion.asin(Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.asin(Quaternion(real: .infinity, imaginary:-.ulpOfOne)).isFinite) + // For randomly-chosen well-scaled finite values, we expect to have + // sin(asin(q)) ≈ q and asin(q) ≈ -asin(-q) + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: -2 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + let p = Quaternion.asin(q) + XCTAssert(Quaternion.sin(p).isApproximatelyEqual(to: q)) + XCTAssert(p.isApproximatelyEqual(to: -.asin(-q))) + } + } + + func testAcosh(_ type: T.Type) { + // acosh(1) = 0 + XCTAssert(Quaternion.acosh(1).isZero) + // acosh is the identity at infinity. + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.acosh(Quaternion(real: .infinity, imaginary:-.ulpOfOne)).isFinite) + // For randomly-chosen well-scaled finite values, we expect to have + // cosh(acosh(q)) ≈ q + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: -2 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + let r = Quaternion.acosh(q) + let s = Quaternion.cosh(r) + if !q.isApproximatelyEqual(to: s) { + print("cosh(acosh()) was not close to identity at q = \(q).") + print("acosh(\(q)) = \(r).") + print("cosh(\(r)) = \(s).") + XCTFail() + } + } + } + + func testAsinh(_ type: T.Type) { + // asinh(1) = π/2 + XCTAssert(Quaternion.asin(1).real.isApproximatelyEqual(to: .pi/2)) + XCTAssert(Quaternion.asin(1).isReal) + // asinh(0) = 0 + XCTAssert(Quaternion.asinh(0).isZero) + // asinh(-1) = -π/2 + XCTAssert(Quaternion.asin(-1).real.isApproximatelyEqual(to: -.pi/2)) + XCTAssert(Quaternion.asin(-1).isReal) + // asinh is the identity at infinity. + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.asinh(Quaternion(real: .infinity, imaginary:-.ulpOfOne)).isFinite) + // For randomly-chosen well-scaled finite values, we expect to have + // sinh(asinh(z)) ≈ z + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: -2 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + let r = Quaternion.asinh(q) + let s = Quaternion.sinh(r) + if !q.isApproximatelyEqual(to: s) { + print("sinh(asinh()) was not close to identity at q = \(q).") + print("asinh(\(q)) = \(r).") + print("sinh(\(r)) = \(s).") + XCTFail() + } + } + } + + func testAtanh(_ type: T.Type) { + // For randomly-chosen well-scaled finite values, we expect to have + // tanh(atanh(q)) ≈ q + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: -2 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + let r = Quaternion.atanh(q) + let s = Quaternion.tanh(r) + if !q.isApproximatelyEqual(to: s) { + print("tanh(atanh()) was not close to identity at q = \(q).") + print("atanh(\(q)) = \(r).") + print("tanh(\(r)) = \(s).") + XCTFail() + } + } + } + + // MARK: - pow-like functions + + func testPowQuaternion(_ type: T.Type) { + // pow(0, 0) = inf + let zero: Quaternion = .zero + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.zero, imaginary:-.zero), zero).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero).isFinite) + // pow(0, x) = 0 for x > 0 + let n: Quaternion = 2 + XCTAssertTrue(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real:-.zero, imaginary:-.zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), n).isZero) + // pow is the identity at infinity. + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .ulpOfOne, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.ulpOfOne, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .ulpOfOne, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.ulpOfOne, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary: .zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary: .zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary: .zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .ulpOfOne, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.ulpOfOne, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary: .ulpOfOne), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary: .ulpOfOne), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary: .ulpOfOne), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary:-.ulpOfOne), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary:-.ulpOfOne), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary:-.ulpOfOne), n).isFinite) + // For randomly-chosen well-scaled finite values, we expect to have + // pow(q, 0) = 1 and pow(q, 1) ≈ q, as well as + // pow(sqrt(q), 2) ≈ q and pow(root(q, n), n) ≈ q for n > 0 + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<100).map { _ in + Quaternion( + real: T.random(in: 0 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + XCTAssertEqual(Quaternion.pow(q, zero), .one) + XCTAssert(q.isApproximatelyEqual(to: .pow(q, .one))) + XCTAssert(q.isApproximatelyEqual(to: .pow(.sqrt(q), Quaternion(2)))) + for n in 1 ... 10 { + let p = Quaternion(n) + XCTAssert(q.isApproximatelyEqual(to: .pow(.root(q, n), p))) + } + } + } + + func testPowInt(_ type: T.Type) { + // pow(0, 0) = 1 + let zero: Int = .zero + XCTAssertEqual(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero), .one) + XCTAssertEqual(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero), .one) + XCTAssertEqual(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero), .one) + XCTAssertEqual(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero), .one) + XCTAssertEqual(Quaternion.pow(Quaternion(real:-.zero, imaginary:-.zero), zero), .one) + XCTAssertEqual(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero), .one) + XCTAssertEqual(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero), .one) + // pow(0, x) = inf for x < 0 + var n: Int = -1 + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.zero, imaginary:-.zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), n).isFinite) + // pow(0, x) = 0 for x > 0 + n = 2 + XCTAssertTrue(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real:-.zero, imaginary:-.zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), n).isZero) + XCTAssertTrue(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), n).isZero) + // pow is the identity at infinity. + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .ulpOfOne, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.ulpOfOne, imaginary: .nan), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .ulpOfOne, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.ulpOfOne, imaginary:-.infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary: .zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary: .zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary: .zero), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .zero, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .ulpOfOne, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.ulpOfOne, imaginary: .infinity), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary: .ulpOfOne), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary: .ulpOfOne), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary: .ulpOfOne), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .nan, imaginary:-.ulpOfOne), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real:-.infinity, imaginary:-.ulpOfOne), n).isFinite) + XCTAssertFalse(Quaternion.pow(Quaternion(real: .infinity, imaginary:-.ulpOfOne), n).isFinite) + // For randomly-chosen well-scaled finite values, we expect to have + // pow(q, 0) = 1 and pow(q, 1) ≈ q, as well as + // pow(sqrt(q), 2) ≈ q and pow(root(q, n), n) ≈ q for n > 0 + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: 0 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + XCTAssertEqual(Quaternion.pow(q, zero), .one) + XCTAssert(q.isApproximatelyEqual(to: .pow(q, 1))) + XCTAssert(q.isApproximatelyEqual(to: .pow(.sqrt(q), 2))) + for n in 1 ... 10 { + XCTAssert(q.isApproximatelyEqual(to: .pow(.root(q, n), n))) + } + } + } + + func testSqrt(_ type: T.Type) { + // sqrt(0) = 0 + XCTAssert(Quaternion.sqrt(Quaternion(real: .zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real: .zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real:-.zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real:-.zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real:-.zero, imaginary:-.zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real: .zero, imaginary:-.zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real: .zero, imaginary:-.zero)).isZero) + // sqrt is the identity at infinity. + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .nan, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .zero, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real:-.infinity, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real:-.ulpOfOne, imaginary: .nan)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .nan, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real:-.infinity, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real:-.ulpOfOne, imaginary:-.infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .nan, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real:-.infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .infinity, imaginary: .zero)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .nan, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .zero, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real:-.infinity, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real:-.ulpOfOne, imaginary: .infinity)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .nan, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real:-.infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .infinity, imaginary: .ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .nan, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real:-.infinity, imaginary:-.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.sqrt(Quaternion(real: .infinity, imaginary:-.ulpOfOne)).isFinite) + } + + func testRoot(_ type: T.Type) { + // root(0, 0) = 0 + XCTAssert(Quaternion.sqrt(Quaternion(real: .zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real: .zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real:-.zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real:-.zero, imaginary: .zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real:-.zero, imaginary:-.zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real: .zero, imaginary:-.zero)).isZero) + XCTAssert(Quaternion.sqrt(Quaternion(real: .zero, imaginary:-.zero)).isZero) + // root(x, 0) = undefined + XCTAssertFalse(Quaternion.root(Quaternion(real: .ulpOfOne, imaginary: .ulpOfOne), .zero).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .ulpOfOne, imaginary: .ulpOfOne), .zero).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.ulpOfOne, imaginary: .ulpOfOne), .zero).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.ulpOfOne, imaginary: .ulpOfOne), .zero).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.ulpOfOne, imaginary:-.ulpOfOne), .zero).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .ulpOfOne, imaginary:-.ulpOfOne), .zero).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .ulpOfOne, imaginary:-.ulpOfOne), .zero).isFinite) + // root is the identity at infinity. + XCTAssertFalse(Quaternion.root(Quaternion(real: .nan, imaginary: .nan), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .zero, imaginary: .nan), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .infinity, imaginary: .nan), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.infinity, imaginary: .nan), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .ulpOfOne, imaginary: .nan), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.ulpOfOne, imaginary: .nan), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .nan, imaginary:-.infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .zero, imaginary:-.infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .infinity, imaginary:-.infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.infinity, imaginary:-.infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .ulpOfOne, imaginary:-.infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.ulpOfOne, imaginary:-.infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .nan, imaginary: .zero), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.infinity, imaginary: .zero), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .infinity, imaginary: .zero), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .nan, imaginary: .infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .zero, imaginary: .infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .infinity, imaginary: .infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.infinity, imaginary: .infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .ulpOfOne, imaginary: .infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.ulpOfOne, imaginary: .infinity), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .nan, imaginary: .ulpOfOne), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.infinity, imaginary: .ulpOfOne), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .infinity, imaginary: .ulpOfOne), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .nan, imaginary:-.ulpOfOne), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real:-.infinity, imaginary:-.ulpOfOne), 2).isFinite) + XCTAssertFalse(Quaternion.root(Quaternion(real: .infinity, imaginary:-.ulpOfOne), 2).isFinite) + // For randomly-chosen well-scaled finite values, we expect to have + // root(q, 1) ≈ q + var g = SystemRandomNumberGenerator() + let values: [Quaternion] = (0..<1000).map { _ in + Quaternion( + real: T.random(in: 0 ... 2, using: &g), + imaginary: + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g), + T.random(in: -2 ... 2, using: &g) + ) + } + for q in values { + XCTAssert(q.isApproximatelyEqual(to: .root(q, 1))) + } + } + + func testFloat() { + testExp(Float32.self) + testExpMinusOne(Float32.self) + testCosh(Float32.self) + testSinh(Float32.self) + testCosSin(Float32.self) + + testLog(Float32.self) + testLogOnePlus(Float32.self) + testAcos(Float32.self) + testAsin(Float32.self) + testAcosh(Float32.self) + testAsinh(Float32.self) + testAtanh(Float32.self) + + testPowQuaternion(Float32.self) + testPowInt(Float32.self) + testSqrt(Float32.self) + testRoot(Float32.self) + } + + func testDouble() { + testExp(Float64.self) + testExpMinusOne(Float64.self) + testCosh(Float64.self) + testSinh(Float64.self) + testCosSin(Float64.self) + + testLog(Float64.self) + testLogOnePlus(Float64.self) + testAcos(Float64.self) + testAsin(Float64.self) + testAcosh(Float64.self) + testAsinh(Float64.self) + testAtanh(Float64.self) + + testPowQuaternion(Float64.self) + testPowInt(Float64.self) + testSqrt(Float64.self) + testRoot(Float64.self) + } +} diff --git a/Tests/QuaternionTests/ImaginaryTestHelper.swift b/Tests/QuaternionTests/ImaginaryTestHelper.swift new file mode 100644 index 00000000..0e51b059 --- /dev/null +++ b/Tests/QuaternionTests/ImaginaryTestHelper.swift @@ -0,0 +1,28 @@ +//===--- ImaginaryTestHelper.swift ----------------------------*- swift -*-===// +// +// This source file is part of the Swift.org open source project +// +// Copyright (c) 2019-2022 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 +// +//===----------------------------------------------------------------------===// + +extension SIMD3 where Scalar: FloatingPoint { + /// Returns a vector with .ulpOfOne in all lanes + static var ulpOfOne: Self { + Self(repeating: .ulpOfOne) + } + + /// Returns a vector with nan in all lanes + static var nan: Self { + SIMD3(repeating: .nan) + } + + /// Returns true if all lanes are NaN + var isNaN: Bool { + x.isNaN && y.isNaN && z.isNaN + } +} diff --git a/Tests/QuaternionTests/PropertyTests.swift b/Tests/QuaternionTests/PropertyTests.swift new file mode 100644 index 00000000..4536fcea --- /dev/null +++ b/Tests/QuaternionTests/PropertyTests.swift @@ -0,0 +1,236 @@ +//===--- PropertyTests.swift ----------------------------------*- swift -*-===// +// +// This source file is part of the Swift Numerics open source project +// +// Copyright (c) 2019 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 + +@testable import QuaternionModule + +final class PropertyTests: XCTestCase { + + func testProperties(_ type: T.Type) { + // The real and imaginary parts of a non-finite value should be nan. + XCTAssertTrue(Quaternion.infinity.real.isNaN) + XCTAssertTrue(Quaternion.infinity.imaginary.x.isNaN) + XCTAssertTrue(Quaternion.infinity.imaginary.y.isNaN) + XCTAssertTrue(Quaternion.infinity.imaginary.z.isNaN) + XCTAssertTrue(Quaternion(real: .infinity, imaginary: .nan, .nan, .nan).real.isNaN) + XCTAssertTrue(Quaternion(real: .nan, imaginary: 0, 0, 0).imaginary.x.isNaN) + XCTAssertTrue(Quaternion(real: .nan, imaginary: 0, 0, 0).imaginary.y.isNaN) + XCTAssertTrue(Quaternion(real: .nan, imaginary: 0, 0, 0).imaginary.z.isNaN) + // The length of a non-finite value should be infinity. + XCTAssertEqual(Quaternion.infinity.length, .infinity) + XCTAssertEqual(Quaternion(real: .infinity, imaginary: .nan, .nan, .nan).length, .infinity) + XCTAssertEqual(Quaternion(real: .nan, imaginary: 0, 0, 0).length, .infinity) + // The length of a zero value should be zero. + XCTAssertEqual(Quaternion.zero.length, .zero) + XCTAssertEqual(Quaternion(real: .zero, imaginary: -.zero).length, .zero) + XCTAssertEqual(Quaternion(real: -.zero, imaginary: -.zero).length, .zero) + // Check for overflow and underflow when calculating the length + XCTAssertEqual(Quaternion(real: .greatestFiniteMagnitude, imaginary: 0, 0, 0).length, .greatestFiniteMagnitude) + XCTAssertEqual(Quaternion(real: 0, imaginary: -.greatestFiniteMagnitude, 0, 0).length, .greatestFiniteMagnitude) + XCTAssertEqual(Quaternion(real: 0, imaginary: 0, .leastNormalMagnitude, 0).length, .leastNormalMagnitude) + XCTAssertEqual(Quaternion(real: 0, imaginary: 0, 0, -.leastNormalMagnitude).length, .leastNormalMagnitude) + // The properties of pure and real + XCTAssertTrue(Quaternion.zero.isPure) // zero quaternion is both, pure... + XCTAssertTrue(Quaternion.zero.isReal) // and real + XCTAssertFalse(Quaternion(1).isPure) + XCTAssertTrue(Quaternion(1).isReal) + XCTAssertTrue(Quaternion(real: .zero, imaginary: 1, 0, 0).isPure) + XCTAssertFalse(Quaternion(real: .zero, imaginary: 1, 0, 0).isReal) + XCTAssertFalse(Quaternion(from: SIMD4(repeating: 1)).isPure) + XCTAssertFalse(Quaternion(from: SIMD4(repeating: 1)).isReal) + } + + func testProperties() { + testProperties(Float32.self) + testProperties(Float64.self) + } + + func testEquatableHashable(_ type: T.Type) { + // Validate that all zeros compare and hash equal, and all non-finites + // do too. + let zeros = [ + Quaternion(real: .zero, imaginary: .zero, .zero, .zero), + Quaternion(real: .zero, imaginary: -.zero, .zero, .zero), + Quaternion(real: .zero, imaginary: .zero, -.zero, .zero), + Quaternion(real: .zero, imaginary: .zero, .zero, -.zero), + Quaternion(real: .zero, imaginary: -.zero, -.zero, .zero), + Quaternion(real: .zero, imaginary: -.zero, .zero, -.zero), + Quaternion(real: .zero, imaginary: .zero, -.zero, -.zero), + Quaternion(real: .zero, imaginary: -.zero, -.zero, -.zero), + + Quaternion(real: -.zero, imaginary: .zero, .zero, .zero), + Quaternion(real: -.zero, imaginary: -.zero, .zero, .zero), + Quaternion(real: -.zero, imaginary: .zero, -.zero, .zero), + Quaternion(real: -.zero, imaginary: .zero, .zero, -.zero), + Quaternion(real: -.zero, imaginary: -.zero, -.zero, .zero), + Quaternion(real: -.zero, imaginary: -.zero, .zero, -.zero), + Quaternion(real: -.zero, imaginary: .zero, -.zero, -.zero), + Quaternion(real: -.zero, imaginary: -.zero, -.zero, -.zero) + ] + for z in zeros[1...] { + XCTAssertEqual(zeros[0], z) + XCTAssertEqual(zeros[0].hashValue, z.hashValue) + } + let infs = [ + Quaternion(real: .nan, imaginary: .nan, .nan, .nan), + Quaternion(real: -.infinity, imaginary: .nan, .nan, .nan), + Quaternion(real: -.ulpOfOne, imaginary: .nan, .nan, .nan), + Quaternion(real: .zero, imaginary: .nan, .nan, .nan), + Quaternion(real: .pi, imaginary: .nan, .nan, .nan), + Quaternion(real: .infinity, imaginary: .nan, .nan, .nan), + Quaternion(real: .nan, imaginary: -.infinity, -.infinity, -.infinity), + Quaternion(real: -.infinity, imaginary: -.infinity, -.infinity, -.infinity), + Quaternion(real: -.ulpOfOne, imaginary: -.infinity, -.infinity, -.infinity), + Quaternion(real: .zero, imaginary: -.infinity, -.infinity, -.infinity), + Quaternion(real: .pi, imaginary: -.infinity, -.infinity, -.infinity), + Quaternion(real: .infinity, imaginary: -.infinity, -.infinity, -.infinity), + Quaternion(real: .nan, imaginary: -.ulpOfOne, -.ulpOfOne, -.ulpOfOne), + Quaternion(real: -.infinity, imaginary: -.ulpOfOne, -.ulpOfOne, -.ulpOfOne), + Quaternion(real: .infinity, imaginary: -.ulpOfOne, -.ulpOfOne, -.ulpOfOne), + Quaternion(real: .nan, imaginary: .zero, .zero, .zero), + Quaternion(real: -.infinity, imaginary: .zero, .zero, .zero), + Quaternion(real: .infinity, imaginary: .zero, .zero, .zero), + Quaternion(real: .nan, imaginary: .pi, .pi, .pi), + Quaternion(real: -.infinity, imaginary: .pi, .pi, .pi), + Quaternion(real: .infinity, imaginary: .pi, .pi, .pi), + Quaternion(real: .nan, imaginary: .infinity, .infinity, .infinity), + Quaternion(real: -.infinity, imaginary: .infinity, .infinity, .infinity), + Quaternion(real: -.ulpOfOne, imaginary: .infinity, .infinity, .infinity), + Quaternion(real: .zero, imaginary: .infinity, .infinity, .infinity), + Quaternion(real: .pi, imaginary: .infinity, .infinity, .infinity), + Quaternion(real: .infinity, imaginary: .infinity, .infinity, .infinity), + ] + for i in infs[1...] { + XCTAssertEqual(infs[0], i) + XCTAssertEqual(infs[0].hashValue, i.hashValue) + } + // Validate that all *normal* values hash their absolute components, so + // that rotations in *R³* of `q` and `-q` will hash to same value. + let pairs: [(lhs: Quaternion, rhs: Quaternion)] = [ + ( + Quaternion(real: -.pi, imaginary: .pi, .pi, .pi), + Quaternion(real: .pi, imaginary: -.pi, -.pi, -.pi) + ), ( + Quaternion(real: .pi, imaginary: -.pi, .pi, .pi), + Quaternion(real: -.pi, imaginary: .pi, -.pi, -.pi) + ), ( + Quaternion(real: .pi, imaginary: .pi, -.pi, .pi), + Quaternion(real: -.pi, imaginary: -.pi, .pi, -.pi) + ), ( + Quaternion(real: .pi, imaginary: .pi, .pi, -.pi), + Quaternion(real: -.pi, imaginary: -.pi, -.pi, .pi) + ), ( + Quaternion(real: -.pi, imaginary: -.pi, .pi, .pi), + Quaternion(real: .pi, imaginary: .pi, -.pi, -.pi) + ), ( + Quaternion(real: .pi, imaginary: -.pi, -.pi, .pi), + Quaternion(real: -.pi, imaginary: .pi, .pi, -.pi) + ), ( + Quaternion(real: .pi, imaginary: .pi, -.pi, -.pi), + Quaternion(real: -.pi, imaginary: -.pi, .pi, .pi) + ), ( + Quaternion(real: .pi, imaginary: .pi, .pi, .pi), + Quaternion(real: -.pi, imaginary: -.pi, -.pi, -.pi) + ) + ] + for pair in pairs { + XCTAssertEqual(pair.lhs.hashValue, pair.rhs.hashValue) + } + } + + func testEquatableHashable() { + testEquatableHashable(Float32.self) + testEquatableHashable(Float64.self) + } + + func testTransformationEquals(_ type: T.Type) { + let rotations: [(lhs: Quaternion, rhs: Quaternion)] = [ + ( + Quaternion(real: -.pi, imaginary: -.pi, -.pi, -.pi), + Quaternion(real: .pi, imaginary: .pi, .pi, .pi) + ), ( + Quaternion(real: .ulpOfOne, imaginary: .ulpOfOne, .ulpOfOne, .ulpOfOne), + Quaternion(real: -.ulpOfOne, imaginary: -.ulpOfOne, -.ulpOfOne, -.ulpOfOne) + ), ( + Quaternion(real: .pi, imaginary: -.pi, .pi, -.pi), + Quaternion(real: -.pi, imaginary: .pi, -.pi, .pi) + ), ( + Quaternion(real: -.ulpOfOne, imaginary: -.ulpOfOne, .ulpOfOne, .ulpOfOne), + Quaternion(real: .ulpOfOne, imaginary: .ulpOfOne, -.ulpOfOne, -.ulpOfOne) + ), + + // Zero and infinity must have equal rotations too + ( + Quaternion.zero, + -Quaternion.zero + ), ( + -Quaternion.infinity, + Quaternion.infinity + ), + ] + for (lhs, rhs) in rotations { + XCTAssertTrue(lhs.equals(as3DTransform: rhs)) + } + + let signDifferentAxis: [(lhs: Quaternion, rhs: Quaternion)] = [ + ( + Quaternion(real: -.pi, imaginary: -.pi, -.pi, -.pi), + Quaternion(real: -.pi, imaginary: .pi, .pi, .pi) + ), ( + Quaternion(real: -.ulpOfOne, imaginary: .ulpOfOne, .ulpOfOne, .ulpOfOne), + Quaternion(real: -.ulpOfOne, imaginary: -.ulpOfOne, -.ulpOfOne, -.ulpOfOne) + ), ( + Quaternion(real: -.pi, imaginary: -.pi, .pi, -.pi), + Quaternion(real: -.pi, imaginary: .pi, -.pi, .pi) + ), ( + Quaternion(real: -.ulpOfOne, imaginary: -.ulpOfOne, .ulpOfOne, .ulpOfOne), + Quaternion(real: -.ulpOfOne, imaginary: .ulpOfOne, -.ulpOfOne, -.ulpOfOne) + ) + ] + for (lhs, rhs) in signDifferentAxis { + XCTAssertFalse(lhs.equals(as3DTransform: rhs)) + } + } + + func testTransformationEquals() { + testTransformationEquals(Float32.self) + testTransformationEquals(Float64.self) + } + + func testCodable(_ type: T.Type) throws { + let encoder = JSONEncoder() + encoder.nonConformingFloatEncodingStrategy = .convertToString( + positiveInfinity: "inf", + negativeInfinity: "-inf", + nan: "nan" + ) + + let decoder = JSONDecoder() + decoder.nonConformingFloatDecodingStrategy = .convertFromString( + positiveInfinity: "inf", + negativeInfinity: "-inf", + nan: "nan" + ) + + for expected: Quaternion in [.zero, .one, .i, .infinity] { + let data = try encoder.encode(expected) + let actual = try decoder.decode(Quaternion.self, from: data) + XCTAssertEqual(actual, expected) + } + } + + func testCodable() throws { + try testCodable(Float32.self) + try testCodable(Float64.self) + } +} diff --git a/Tests/QuaternionTests/TransformationTests.swift b/Tests/QuaternionTests/TransformationTests.swift new file mode 100644 index 00000000..c5cf3800 --- /dev/null +++ b/Tests/QuaternionTests/TransformationTests.swift @@ -0,0 +1,374 @@ +//===--- TransformationTests.swift ----------------------------*- swift -*-===// +// +// This source file is part of the Swift Numerics open source project +// +// Copyright (c) 2020 - 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 +// +//===----------------------------------------------------------------------===// + +import XCTest +import RealModule + +@testable import QuaternionModule + +final class TransformationTests: XCTestCase { + + // MARK: Angle/Axis + + func testAngleAxis(_ type: T.Type) { + let xAxis = SIMD3(1,0,0) + // Positive angle, positive axis + XCTAssertEqual(Quaternion(angle: .pi, axis: xAxis).angle, .pi) + XCTAssertEqual(Quaternion(angle: .pi, axis: xAxis).axis, xAxis) + // Negative angle, positive axis + XCTAssertEqual(Quaternion(angle: -.pi, axis: xAxis).angle, .pi) + XCTAssertEqual(Quaternion(angle: -.pi, axis: xAxis).axis, -xAxis) + // Positive angle, negative axis + XCTAssertEqual(Quaternion(angle: .pi, axis: -xAxis).angle, .pi) + XCTAssertEqual(Quaternion(angle: .pi, axis: -xAxis).axis, -xAxis) + // Negative angle, negative axis + XCTAssertEqual(Quaternion(angle: -.pi, axis: -xAxis).angle, .pi) + XCTAssertEqual(Quaternion(angle: -.pi, axis: -xAxis).axis, xAxis) + } + + func testAngleAxis() { + testAngleAxis(Float32.self) + testAngleAxis(Float64.self) + } + + func testAngleAxisMultipleOfPi(_ type: T.Type) { + let xAxis = SIMD3(1,0,0) + // 2π + let pi2 = Quaternion(angle: .pi * 2, axis: xAxis) + XCTAssertEqual(pi2.angle, .pi * 2) + XCTAssertEqual(pi2.axis, xAxis) + // 3π - axis inverted + let pi3 = Quaternion(angle: .pi * 3, axis: xAxis) + XCTAssertTrue(closeEnough(pi3.angle, .pi, ulps: 1)) + XCTAssertEqual(pi3.axis, -xAxis) + // 5π - axis restored + let pi5 = Quaternion(angle: .pi * 5, axis: xAxis) + XCTAssertTrue(closeEnough(pi5.angle, .pi, ulps: 5)) + XCTAssertEqual(pi5.axis, xAxis) + } + + func testAngleAxisMultipleOfPi() { + testAngleAxisMultipleOfPi(Float32.self) + testAngleAxisMultipleOfPi(Float64.self) + } + + func testAngleAxisEdgeCases(_ type: T.Type) { + // Zero/Zero + XCTAssertTrue(Quaternion(length: .zero, angle: .zero, axis: .zero).axis.isNaN) + XCTAssertTrue(Quaternion(length: .zero, angle: .zero, axis: .zero).angle.isNaN) + XCTAssertEqual(Quaternion(length: .zero, angle: .zero, axis: .zero), .zero) + // Inf/Zero + XCTAssertTrue(Quaternion(length: .zero, angle: .infinity, axis: .infinity).axis.isNaN) + XCTAssertTrue(Quaternion(length: .zero, angle: .infinity, axis: .infinity).angle.isNaN) + XCTAssertEqual(Quaternion(length: .zero, angle: .infinity, axis: .infinity), .zero) + // -Inf/Zero + XCTAssertTrue(Quaternion(length: .zero, angle: -.infinity, axis: -.infinity).axis.isNaN) + XCTAssertTrue(Quaternion(length: .zero, angle: -.infinity, axis: -.infinity).angle.isNaN) + XCTAssertEqual(Quaternion(length: .zero, angle: -.infinity, axis: -.infinity), .zero) + // NaN/Zero + XCTAssertTrue(Quaternion(length: .zero, angle: .nan, axis: .nan).axis.isNaN) + XCTAssertTrue(Quaternion(length: .zero, angle: .nan, axis: .nan).angle.isNaN) + XCTAssertEqual(Quaternion(length: .zero, angle: .nan, axis: .nan), .zero) + // Zero/Inf + XCTAssertTrue(Quaternion(length: .infinity, angle: .zero, axis: .zero).axis.isNaN) + XCTAssertTrue(Quaternion(length: .infinity, angle: .zero, axis: .zero).angle.isNaN) + XCTAssertEqual(Quaternion(length: .infinity, angle: .zero, axis: .zero), .infinity) + // Inf/Inf + XCTAssertTrue(Quaternion(length: .infinity, angle: .infinity, axis: .infinity).axis.isNaN) + XCTAssertTrue(Quaternion(length: .infinity, angle: .infinity, axis: .infinity).angle.isNaN) + XCTAssertEqual(Quaternion(length: .infinity, angle: .infinity, axis: .infinity), .infinity) + // -Inf/Inf + XCTAssertTrue(Quaternion(length: .infinity, angle: -.infinity, axis: -.infinity).axis.isNaN) + XCTAssertTrue(Quaternion(length: .infinity, angle: -.infinity, axis: -.infinity).angle.isNaN) + XCTAssertEqual(Quaternion(length: .infinity, angle: -.infinity, axis: -.infinity), .infinity) + // NaN/Inf + XCTAssertTrue(Quaternion(length: .infinity, angle: .nan, axis: .nan).axis.isNaN) + XCTAssertTrue(Quaternion(length: .infinity, angle: .nan, axis: .nan).angle.isNaN) + XCTAssertEqual(Quaternion(length: .infinity, angle: .nan, axis: .nan), .infinity) + // Zero/-Inf + XCTAssertTrue(Quaternion(length: -.infinity, angle: .zero, axis: .zero).axis.isNaN) + XCTAssertTrue(Quaternion(length: -.infinity, angle: .zero, axis: .zero).angle.isNaN) + XCTAssertEqual(Quaternion(length: -.infinity, angle: .zero, axis: .zero), .infinity) + // Inf/-Inf + XCTAssertTrue(Quaternion(length: -.infinity, angle: .infinity, axis: .infinity).axis.isNaN) + XCTAssertTrue(Quaternion(length: -.infinity, angle: .infinity, axis: .infinity).angle.isNaN) + XCTAssertEqual(Quaternion(length: -.infinity, angle: .infinity, axis: .infinity), .infinity) + // -Inf/-Inf + XCTAssertTrue(Quaternion(length: -.infinity, angle: -.infinity, axis: -.infinity).axis.isNaN) + XCTAssertTrue(Quaternion(length: -.infinity, angle: -.infinity, axis: -.infinity).angle.isNaN) + XCTAssertEqual(Quaternion(length: -.infinity, angle: -.infinity, axis: -.infinity), .infinity) + // NaN/-Inf + XCTAssertTrue(Quaternion(length: -.infinity, angle: .nan, axis: .nan).axis.isNaN) + XCTAssertTrue(Quaternion(length: -.infinity, angle: .nan, axis: .nan).angle.isNaN) + XCTAssertEqual(Quaternion(length: -.infinity, angle: .nan, axis: .nan), .infinity) + } + + func testAngleAxisEdgeCases() { + testAngleAxisEdgeCases(Float32.self) + testAngleAxisEdgeCases(Float64.self) + } + + func testHalfAngleAndAxisOverflow(_ type: T.Type) { + let unscaled = Quaternion(real: 1, imaginary: .one) + let scaled = Quaternion( + real: .greatestFiniteMagnitude, + imaginary: SIMD3(repeating: .greatestFiniteMagnitude) + ) + XCTAssertEqual(scaled.angle, unscaled.angle) + XCTAssertEqual(scaled.axis, unscaled.axis) + } + + func testHalfAngleAndAxisOverflow() { + testHalfAngleAndAxisOverflow(Float32.self) + testHalfAngleAndAxisOverflow(Float64.self) + } + + func testHalfAngleAndAxisUnderflow(_ type: T.Type) { + let unscaled = Quaternion(real: 1, imaginary: .one) + let scaled = Quaternion( + real: .leastNormalMagnitude, + imaginary: SIMD3(repeating: .leastNormalMagnitude) + ) + XCTAssertEqual(scaled.angle, unscaled.angle) + XCTAssertEqual(scaled.axis, unscaled.axis) + } + + func testHalfAngleAndAxisUnderflow() { + testHalfAngleAndAxisUnderflow(Float32.self) + testHalfAngleAndAxisUnderflow(Float64.self) + } + + // MARK: Rotation Vector + + func testRotationVector(_ type: T.Type) { + let vector = SIMD3(0,-1,0) * .pi + XCTAssertEqual(Quaternion(rotation: vector).rotationVector.x, .zero) + XCTAssertEqual(Quaternion(rotation: vector).rotationVector.y, -.pi) + XCTAssertEqual(Quaternion(rotation: vector).rotationVector.z, .zero) + + XCTAssertEqual(Quaternion(rotation: vector).axis, SIMD3(0,-1,0)) + XCTAssertEqual(Quaternion(rotation: vector).angle, .pi) + } + + func testRotationVector() { + testRotationVector(Float32.self) + testRotationVector(Float64.self) + } + + func testRotationVectorEdgeCases(_ type: T.Type) { + XCTAssertEqual(Quaternion(rotation: .zero), .zero) + XCTAssertEqual(Quaternion(rotation: .infinity), .infinity) + XCTAssertEqual(Quaternion(rotation: -.infinity), .infinity) + XCTAssertTrue(Quaternion(rotation: .nan).real.isNaN) + XCTAssertTrue(Quaternion(rotation: .nan).imaginary.isNaN) + } + + func testRotationVectorEdgeCases() { + testRotationVectorEdgeCases(Float32.self) + testRotationVectorEdgeCases(Float64.self) + } + + // MARK: Polar Decomposition + + func testPolarDecomposition(_ type: T.Type) { + let axis = SIMD3(0,-1,0) + + let q = Quaternion(length: 5, phase: .pi, axis: axis) + XCTAssertEqual(q.axis, axis) + XCTAssertEqual(q.angle, .pi * 2) + + XCTAssertEqual(q.polar.length, 5) + XCTAssertEqual(q.polar.phase, .pi) + XCTAssertEqual(q.polar.axis, axis) + } + + func testPolarDecomposition() { + testPolarDecomposition(Float32.self) + testPolarDecomposition(Float64.self) + } + + func testPolarDecompositionEdgeCases(_ type: T.Type) { + XCTAssertEqual(Quaternion(length: .zero, phase: .zero , axis: .zero ), .zero) + XCTAssertEqual(Quaternion(length: .zero, phase: .infinity, axis: .infinity), .zero) + XCTAssertEqual(Quaternion(length: .zero, phase:-.infinity, axis: -.infinity), .zero) + XCTAssertEqual(Quaternion(length: .zero, phase: .nan , axis: .nan ), .zero) + XCTAssertEqual(Quaternion(length: .infinity, phase: .zero , axis: .zero ), .infinity) + XCTAssertEqual(Quaternion(length: .infinity, phase: .infinity, axis: .infinity), .infinity) + XCTAssertEqual(Quaternion(length: .infinity, phase:-.infinity, axis: -.infinity), .infinity) + XCTAssertEqual(Quaternion(length: .infinity, phase: .nan , axis: .infinity), .infinity) + XCTAssertEqual(Quaternion(length:-.infinity, phase: .zero , axis: .zero ), .infinity) + XCTAssertEqual(Quaternion(length:-.infinity, phase: .infinity, axis: .infinity), .infinity) + XCTAssertEqual(Quaternion(length:-.infinity, phase:-.infinity, axis: -.infinity), .infinity) + XCTAssertEqual(Quaternion(length:-.infinity, phase: .nan , axis: .infinity), .infinity) + } + + func testPolarDecompositionEdgeCases() { + testPolarDecompositionEdgeCases(Float32.self) + testPolarDecompositionEdgeCases(Float64.self) + } + + // MARK: Act on Vector + + func testActOnVector(_ type: T.Type) { + let vector = SIMD3(1,1,1) + let xAxis = SIMD3(1,0,0) + + let singleAxis = Quaternion(angle: .pi/2, axis: SIMD3(0, 1, 0)) + XCTAssertTrue(singleAxis.act(on: xAxis).x.isApproximatelyEqual( + to: .zero, absoluteTolerance: T.ulpOfOne.squareRoot() + )) + XCTAssertTrue(singleAxis.act(on: xAxis).y.isApproximatelyEqual( + to: .zero, absoluteTolerance: T.ulpOfOne.squareRoot() + )) + XCTAssertEqual(singleAxis.act(on: xAxis).z, -1) + + let piHalf = Quaternion(angle: .pi/2, axis: xAxis) + XCTAssertTrue(closeEnough(piHalf.act(on: vector).x, 1, ulps: 0)) + XCTAssertTrue(closeEnough(piHalf.act(on: vector).y, -1, ulps: 1)) + XCTAssertTrue(closeEnough(piHalf.act(on: vector).z, 1, ulps: 1)) + + let pi = Quaternion(angle: .pi, axis: xAxis) + XCTAssertTrue(closeEnough(pi.act(on: vector).x, 1, ulps: 0)) + XCTAssertTrue(closeEnough(pi.act(on: vector).y, -1, ulps: 2)) + XCTAssertTrue(closeEnough(pi.act(on: vector).z, -1, ulps: 2)) + + let twoPi = Quaternion(angle: .pi * 2, axis: xAxis) + XCTAssertTrue(closeEnough(twoPi.act(on: vector).x, 1, ulps: 0)) + XCTAssertTrue(closeEnough(twoPi.act(on: vector).y, 1, ulps: 3)) + XCTAssertTrue(closeEnough(twoPi.act(on: vector).z, 1, ulps: 3)) + } + + func testActOnVector() { + testActOnVector(Float32.self) + testActOnVector(Float64.self) + } + + func testActOnVectorRandom(_ type: T.Type) + where T: Real, T: BinaryFloatingPoint, T: SIMDScalar, + T.Exponent: FixedWidthInteger, T.RawSignificand: FixedWidthInteger + { + // Generate random angles, axis and vector to test rotation properties + // - angle are selected from range -π to π + // - axis values are selected from -1 to 1; axis length is unity + // - vector values are selected from 10 to 10000 + let inputs: [(angle: T, axis: SIMD3, vector: SIMD3)] = (0..<100).map { _ in + let angle = T.random(in: -.pi ... .pi) + var axis = SIMD3.random(in: -1 ... 1) + axis /= .sqrt((axis * axis).sum()) // Normalize + var vector = SIMD3.random(in: -1 ... 1) + vector /= .sqrt((vector * vector).sum()) // Normalize + vector *= T.random(in: 10 ... 10000) // Scale + return (angle, axis, vector) + } + + for (angle, axis, vector) in inputs { + let q = Quaternion(angle: angle, axis: axis) + // The following equation in the form of v' = qvq⁻¹ is the mathmatical + // definition for how a quaternion rotates a vector (by promoting it to + // a quaternion) and goes "the full and long way" to calculate the + // rotation of vector by a quaternion. The result is used to test the + // rotation properties of "act(on:)" + let vrot = (q // q + * Quaternion(imaginary: vector) // v (pure quaternion) + * q.conjugate // q⁻¹ (as q is of unit length, q⁻¹ == q*) + ).imaginary // the result is a pure quaternion with v' == imaginary + + XCTAssertTrue(q.act(on: vector).x.isFinite) + XCTAssertTrue(q.act(on: vector).y.isFinite) + XCTAssertTrue(q.act(on: vector).z.isFinite) + // Test for sign equality on the components to see if the vector rotated + // to the correct quadrant and if the vector is of equal length, instead + // of testing component equality – as they are hard to compare with + // proper tolerance + XCTAssertEqual(q.act(on: vector).x.sign, vrot.x.sign) + XCTAssertEqual(q.act(on: vector).y.sign, vrot.y.sign) + XCTAssertEqual(q.act(on: vector).z.sign, vrot.z.sign) + XCTAssertTrue(closeEnough(q.act(on: vector).lengthSquared, vrot.lengthSquared, ulps: 16)) + } + } + + func testActOnVectorRandom() { + testActOnVectorRandom(Float32.self) + testActOnVectorRandom(Float64.self) + } + + func testActOnVectorEdgeCase(_ type: T.Type) { + + /// Test for zero, infinity + let q = Quaternion(angle: .pi, axis: SIMD3(1,0,0)) + XCTAssertEqual(q.act(on: .zero), .zero) + XCTAssertEqual(q.act(on: -.zero), .zero) + XCTAssertEqual(q.act(on: .nan ), .infinity) + XCTAssertEqual(q.act(on: .infinity), .infinity) + XCTAssertEqual(q.act(on: -.infinity), .infinity) + } + + func testActOnVectorEdgeCase() { + testActOnVectorEdgeCase(Float32.self) + testActOnVectorEdgeCase(Float64.self) + } + + func testActOnVectorOverflow(_ type: T.Type) { + // Create a vector (somewhat) close to greatestFiniteMagnitude on all lanes. + // We can not use greatestFiniteMagnitude here to test the careful rotation + // path, as we lose some precision in the process and it will overflow after + // rescaling the vector. + let scalar = T( + sign: .plus, exponent: T.greatestFiniteMagnitude.exponent, + significand: 1.999999 + ) + let closeToBounds = SIMD3(repeating: scalar) + + // An axis perpendicular to the vector, so all lanes change equally + let axis = SIMD3(1,0,-1) / .sqrt(2) + // Perform a 180° rotation on all components + let pi = Quaternion(angle: .pi, axis: axis).act(on: closeToBounds) + // Must be finite after the rotation + XCTAssertTrue(pi.x.isFinite) + XCTAssertTrue(pi.y.isFinite) + XCTAssertTrue(pi.z.isFinite) + XCTAssertTrue(closeEnough(pi.x, -scalar, ulps: 4)) + XCTAssertTrue(closeEnough(pi.y, -scalar, ulps: 4)) + XCTAssertTrue(closeEnough(pi.z, -scalar, ulps: 4)) + } + + func testActOnVectorOverflow() { + testActOnVectorOverflow(Float32.self) + testActOnVectorOverflow(Float64.self) + } + + func testActOnVectorUnderflow(_ type: T.Type) { + let scalar = T.leastNormalMagnitude + let closeToZero = SIMD3(repeating: scalar) + // An axis perpendicular to the vector, so all lanes change equally + let axis = SIMD3(1,0,-1) / .sqrt(2) + // Perform a 180° rotation on all components + let pi = Quaternion(angle: .pi, axis: axis).act(on: closeToZero) + // Must be finite after the rotation + XCTAssertTrue(!pi.x.isZero) + XCTAssertTrue(!pi.y.isZero) + XCTAssertTrue(!pi.z.isZero) + XCTAssertTrue(closeEnough(pi.x, -scalar, ulps: 2)) + XCTAssertTrue(closeEnough(pi.y, -scalar, ulps: 2)) + XCTAssertTrue(closeEnough(pi.z, -scalar, ulps: 2)) + } + + func testActOnVectorUnderflow() { + testActOnVectorUnderflow(Float32.self) + testActOnVectorUnderflow(Float64.self) + } +} + +// TODO: replace with approximately equals +func closeEnough(_ a: T, _ b: T, ulps allowed: T) -> Bool { + let scale = max(a.magnitude, b.magnitude, T.leastNormalMagnitude).ulp + return (a - b).magnitude <= allowed * scale +}