diff --git a/Sources/QuaternionModule/ImaginaryHelper.swift b/Sources/QuaternionModule/ImaginaryHelper.swift index 64c4bc04..9ffd4cd2 100644 --- a/Sources/QuaternionModule/ImaginaryHelper.swift +++ b/Sources/QuaternionModule/ImaginaryHelper.swift @@ -19,12 +19,6 @@ extension SIMD3 where Scalar: FloatingPoint { SIMD3(repeating: .infinity) } - /// Returns a vector with nan in all lanes - @usableFromInline @inline(__always) - internal static var nan: Self { - SIMD3(repeating: .nan) - } - /// True if all values of this instance are finite @usableFromInline @inline(__always) internal var isFinite: Bool { diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index 12d2ccd6..52e47acb 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -2,7 +2,7 @@ // // This source file is part of the Swift.org open source project // -// Copyright (c) 2019 - 2022 Apple Inc. and the Swift project authors +// 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 @@ -16,7 +16,7 @@ // 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 the comments of the following functions. +// 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 @@ -26,13 +26,12 @@ import RealModule -extension Quaternion/*: ElementaryFunctions */ { - +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` as follows (`let θ = ||v||`): + // 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) @@ -42,7 +41,7 @@ extension Quaternion/*: ElementaryFunctions */ { // 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. + // `exp(r) cos(||v||)` would not be). guard q.isFinite else { return q } let (â, θ) = q.imaginary.unitAxisAndLength let rotation = Quaternion(halfAngle: θ, unitAxis: â) @@ -59,8 +58,8 @@ extension Quaternion/*: ElementaryFunctions */ { @inlinable public static func expMinusOne(_ q: Quaternion) -> Quaternion { - // Mathematically, this operation can be expanded in terms of the `Real` - // operations `exp`, `cos` and `sin` as follows (`let θ = ||v||`): + // 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 @@ -102,7 +101,7 @@ extension Quaternion/*: ElementaryFunctions */ { @inlinable public static func cosh(_ q: Quaternion) -> Quaternion { // Mathematically, this operation can be expanded in terms of - // trigonometric `Real` operations as follows (`let θ = ||v||`): + // trigonometric `Real` operations (`let θ = ||v||`): // // ``` // cosh(q) = (exp(q) + exp(-q)) / 2 @@ -136,7 +135,7 @@ extension Quaternion/*: ElementaryFunctions */ { @inlinable public static func sinh(_ q: Quaternion) -> Quaternion { // Mathematically, this operation can be expanded in terms of - // trigonometric `Real` operations as follows (`let θ = ||v||`): + // trigonometric `Real` operations (`let θ = ||v||`): // // ``` // sinh(q) = (exp(q) - exp(-q)) / 2 @@ -159,7 +158,7 @@ extension Quaternion/*: ElementaryFunctions */ { @inlinable public static func tanh(_ q: Quaternion) -> Quaternion { // Mathematically, this operation can be expanded in terms of - // trigonometric `Real` operations as follows (`let θ = ||v||`): + // quaternionic `sinh` and `cosh` operations: // // ``` // tanh(q) = sinh(q) / cosh(q) @@ -181,7 +180,12 @@ extension Quaternion/*: ElementaryFunctions */ { @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) @@ -189,7 +193,12 @@ extension Quaternion/*: ElementaryFunctions */ { @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) @@ -197,18 +206,289 @@ extension Quaternion/*: ElementaryFunctions */ { @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. - @usableFromInline @inline(__always) - internal var unitAxisAndLength: (Self, Scalar) { + @_alwaysEmitIntoClient + fileprivate var unitAxisAndLength: (Self, Scalar) { if self == .zero { return (SIMD3( Scalar(signOf: x, magnitudeOf: 0), diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index 196be513..ac026a82 100644 --- a/Tests/QuaternionTests/ElementaryFunctionTests.swift +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -2,7 +2,7 @@ // // This source file is part of the Swift Numerics open source project // -// Copyright (c) 2019 - 2021 Apple Inc. and the Swift Numerics project authors +// 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 @@ -21,28 +21,38 @@ final class ElementaryFunctionTests: XCTestCase { func testExp(_ type: T.Type) { // exp(0) = 1 - XCTAssertEqual(1, Quaternion.exp(Quaternion(real: 0, imaginary: 0, 0, 0))) - XCTAssertEqual(1, Quaternion.exp(Quaternion(real:-0, imaginary: 0, 0, 0))) - XCTAssertEqual(1, Quaternion.exp(Quaternion(real:-0, imaginary:-0,-0,-0))) - XCTAssertEqual(1, Quaternion.exp(Quaternion(real: 0, imaginary:-0,-0,-0))) - // In general, exp(Quaternion(r,0,0,0)) should be exp(r), but that breaks - // down when r is infinity or NaN, because we want all non-finite - // quaternions to be semantically a single point at infinity. This is fine - // for most inputs, but exp(Quaternion(-.infinity, 0, 0, 0)) would produce - // 0 if we evaluated it in the usual way. - XCTAssertFalse(Quaternion.exp(Quaternion(real: .infinity, imaginary: .zero)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: 0, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: -.infinity, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: -.infinity, imaginary: .zero)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: -.infinity, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: 0, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: .infinity, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: .nan, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: .infinity, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: .nan, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: -.infinity, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.exp(Quaternion(real: .nan, imaginary: -.infinity)).isFinite) + 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: @@ -79,28 +89,38 @@ final class ElementaryFunctionTests: XCTestCase { func testExpMinusOne(_ type: T.Type) { // expMinusOne(0) = 0 - XCTAssertEqual(0, Quaternion.expMinusOne(Quaternion(real: 0, imaginary: 0, 0, 0))) - XCTAssertEqual(0, Quaternion.expMinusOne(Quaternion(real:-0, imaginary: 0, 0, 0))) - XCTAssertEqual(0, Quaternion.expMinusOne(Quaternion(real:-0, imaginary:-0,-0,-0))) - XCTAssertEqual(0, Quaternion.expMinusOne(Quaternion(real: 0, imaginary:-0,-0,-0))) - // In general, expMinusOne(Quaternion(r,0,0,0)) should be expMinusOne(r), - // but that breaks down when r is infinity or NaN, because we want all non- - // finite Quaternion values to be semantically a single point at infinity. - // This is fine for most inputs, but expMinusOne(Quaternion(-.infinity,0,0,0)) - // would produce 0 if we evaluated it in the usual way. - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .infinity, imaginary: .zero)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: 0, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: -.infinity, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: -.infinity, imaginary: .zero)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: -.infinity, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: 0, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .infinity, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .nan, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .infinity, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .nan, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: -.infinity, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.expMinusOne(Quaternion(real: .nan, imaginary: -.infinity)).isFinite) + 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))) @@ -126,24 +146,38 @@ final class ElementaryFunctionTests: XCTestCase { func testCosh(_ type: T.Type) { // cosh(0) = 1 - XCTAssertEqual(1, Quaternion.cosh(Quaternion(real: 0, imaginary: 0, 0, 0))) - XCTAssertEqual(1, Quaternion.cosh(Quaternion(real:-0, imaginary: 0, 0, 0))) - XCTAssertEqual(1, Quaternion.cosh(Quaternion(real:-0, imaginary:-0,-0,-0))) - XCTAssertEqual(1, Quaternion.cosh(Quaternion(real: 0, imaginary:-0,-0,-0))) + 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: .infinity, imaginary: .zero)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: 0, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: -.infinity, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: -.infinity, imaginary: .zero)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: -.infinity, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: 0, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: .infinity, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: .nan, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: .infinity, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: .nan, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: -.infinity, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.cosh(Quaternion(real: .nan, imaginary: -.infinity)).isFinite) + 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) @@ -162,24 +196,38 @@ final class ElementaryFunctionTests: XCTestCase { func testSinh(_ type: T.Type) { // sinh(0) = 0 - XCTAssertEqual(0, Quaternion.sinh(Quaternion(real: 0, imaginary: 0, 0, 0))) - XCTAssertEqual(0, Quaternion.sinh(Quaternion(real:-0, imaginary: 0, 0, 0))) - XCTAssertEqual(0, Quaternion.sinh(Quaternion(real:-0, imaginary:-0,-0,-0))) - XCTAssertEqual(0, Quaternion.sinh(Quaternion(real: 0, imaginary:-0,-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: .infinity, imaginary: .zero)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: .infinity, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: 0, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: -.infinity, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: -.infinity, imaginary: .zero)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: -.infinity, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: 0, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: .infinity, imaginary: -.infinity)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: .nan, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: .infinity, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: .nan, imaginary: .infinity)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: -.infinity, imaginary: .nan)).isFinite) - XCTAssertFalse(Quaternion.sinh(Quaternion(real: .nan, imaginary: -.infinity)).isFinite) + 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) @@ -234,12 +282,635 @@ final class ElementaryFunctionTests: XCTestCase { } } + // 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() { @@ -248,5 +919,18 @@ final class ElementaryFunctionTests: XCTestCase { 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/TransformationTests.swift b/Tests/QuaternionTests/TransformationTests.swift index 245ab516..c5cf3800 100644 --- a/Tests/QuaternionTests/TransformationTests.swift +++ b/Tests/QuaternionTests/TransformationTests.swift @@ -2,7 +2,7 @@ // // This source file is part of the Swift Numerics open source project // -// Copyright (c) 2020 Apple Inc. and the Swift Numerics project authors +// 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 @@ -367,11 +367,6 @@ final class TransformationTests: XCTestCase { } } -// MARK: - Helper -extension SIMD3 where Scalar: FloatingPoint { - fileprivate var isNaN: Bool { x.isNaN && y.isNaN && z.isNaN } -} - // 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