From 3ebee8f5f1c501de6e8703fae6d20fb24aa92e89 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Tue, 26 Apr 2022 15:53:43 +0200 Subject: [PATCH 01/18] Revert "Remove log and pow functions" This reverts commit 3262d4d5c581886f2b76880c9520ec80800598f9. --- .../Quaternion+ElementaryFunctions.swift | 132 +++++++++++++++++- .../ElementaryFunctionTests.swift | 24 ++++ 2 files changed, 155 insertions(+), 1 deletion(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index 12d2ccd6..19cc2802 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -26,7 +26,7 @@ import RealModule -extension Quaternion/*: ElementaryFunctions */ { +extension Quaternion/*: ElementaryFunctions*/ { // MARK: - exp-like functions @inlinable @@ -202,6 +202,136 @@ extension Quaternion/*: ElementaryFunctions */ { 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 } + + let argument = q.imaginary.length + let axis = q.imaginary / argument + + // We deliberatly choose log(length) over the (faster) + // log(lengthSquared) / 2 which is used for complex numbers; as + // the squared length of quaternions is more prone to overflows than the + // squared length of complex numbers. + return Quaternion(real: .log(q.length), imaginary: axis * q.halfAngle) + } + + @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 + // = log1p((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) + } + + // + // 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 as follows: + // + // ``` + // pow(q, p) = exp(log(pow(q, p))) + // = exp(p * log(q)) + // ``` + 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 as follows: + // + // ``` + // pow(q, n) = exp(log(pow(q, 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).multiplied(by: RealType(n))) + } + + @inlinable + public static func sqrt(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of the + // quaternionic `exp` and `log` operations as follows: + // + // ``` + // sqrt(q) = q^(1/2) = exp(log(q^(1/2))) + // = exp(log(q) * (1/2)) + // ``` + guard !q.isZero else { return .zero } + return exp(log(q).divided(by: 2)) + } + + @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 as follows: + // + // ``` + // root(q, n) = exp(log(root(q, 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 { diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index 196be513..d2bf1cea 100644 --- a/Tests/QuaternionTests/ElementaryFunctionTests.swift +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -234,12 +234,34 @@ final class ElementaryFunctionTests: XCTestCase { } } + // MARK: - log-like functions + + func testLog(_ type: T.Type) { + // log(0) = undefined/infinity + XCTAssertFalse(Quaternion.log(Quaternion(real: 0, imaginary: 0, 0, 0)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-0, imaginary: 0, 0, 0)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real:-0, imaginary:-0,-0,-0)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: 0, imaginary:-0,-0,-0)).isFinite) + + 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 q in values { + XCTAssertTrue(q.isApproximatelyEqual(to: .log(.exp(q)))) + } + } + func testFloat() { testExp(Float32.self) testExpMinusOne(Float32.self) testCosh(Float32.self) testSinh(Float32.self) testCosSin(Float32.self) + + testLog(Float32.self) } func testDouble() { @@ -248,5 +270,7 @@ final class ElementaryFunctionTests: XCTestCase { testCosh(Float64.self) testSinh(Float64.self) testCosSin(Float64.self) + + testLog(Float64.self) } } From cee5aaf4a5b4940694a2da503fc40ff28f7a8f61 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Thu, 28 Apr 2022 16:23:43 +0200 Subject: [PATCH 02/18] Add more comments to quaternionic elfns --- .../Quaternion+ElementaryFunctions.swift | 46 +++++++++++++------ 1 file changed, 32 insertions(+), 14 deletions(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index 19cc2802..d0ec6533 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -31,8 +31,8 @@ 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 +42,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 +59,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 +102,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 +136,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 +159,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 +181,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 +194,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,7 +207,12 @@ 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) @@ -241,13 +256,16 @@ extension Quaternion/*: ElementaryFunctions*/ { // 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 - // = log1p((2+r)r + x² + y² + z²)/2 + // 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; + // 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| From 42ad4279ba6b262a76b235b1bd76b0e0bf6edfb3 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Thu, 28 Apr 2022 18:41:50 +0200 Subject: [PATCH 03/18] Specialized quaternionic logarithm --- .../Quaternion+ElementaryFunctions.swift | 53 +++++++++++++++---- 1 file changed, 44 insertions(+), 9 deletions(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index d0ec6533..c4f4b3dd 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -224,15 +224,50 @@ extension Quaternion/*: ElementaryFunctions*/ { // 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 } - - let argument = q.imaginary.length - let axis = q.imaginary / argument - - // We deliberatly choose log(length) over the (faster) - // log(lengthSquared) / 2 which is used for complex numbers; as - // the squared length of quaternions is more prone to overflows than the - // squared length of complex numbers. - return Quaternion(real: .log(q.length), imaginary: axis * q.halfAngle) + // Having eliminated non-finite values and zero, the imaginary part is + // straightforward: + // TODO: There is a potential optimisation hidden here, as length is + // calculated twice (halfAngle, unitAxisAndLength) + let argument = q.halfAngle + let (â, θ) = q.imaginary.unitAxisAndLength + 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 From 2a1fcfadab221ebf0d60b5754c0a8d5f3e9dd49f Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Thu, 28 Apr 2022 18:43:29 +0200 Subject: [PATCH 04/18] Add log1p test cases --- .../ElementaryFunctionTests.swift | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index d2bf1cea..2f1d4539 100644 --- a/Tests/QuaternionTests/ElementaryFunctionTests.swift +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -254,6 +254,51 @@ final class ElementaryFunctionTests: XCTestCase { } } + func testLogOnePlus(_ type: T.Type) { + // log(onePlus: 0) = 0 + XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real: 0, imaginary: 0, 0, 0)).isZero) + XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real:-0, imaginary: 0, 0, 0)).isZero) + XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real:-0, imaginary:-0,-0,-0)).isZero) + XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real: 0, imaginary:-0,-0,-0)).isZero) + // log(onePlus:) is the identity at infinity. + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: .nan, .nan, .nan)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: .nan, .nan, .nan)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: .nan, .nan, .nan)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.infinity, imaginary: .nan, .nan, .nan)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: -.infinity, -.infinity, -.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: -.infinity, -.infinity, -.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: -.infinity, -.infinity, -.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.infinity, imaginary: -.infinity, -.infinity, -.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.ulpOfOne, imaginary: -.infinity, -.infinity, -.infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: .zero, .zero, .zero)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.infinity, imaginary: .zero, .zero, .zero)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: .zero, .zero, .zero)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: .infinity, .infinity, .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: .infinity, .infinity, .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: .infinity, .infinity, .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.infinity, imaginary: .infinity, .infinity, .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.ulpOfOne, imaginary: .infinity, .infinity, .infinity)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: -.ulpOfOne, -.ulpOfOne, -.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.infinity, imaginary: -.ulpOfOne, -.ulpOfOne, -.ulpOfOne)).isFinite) + XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: -.ulpOfOne, -.ulpOfOne, -.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: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g) + ) + } + for q in values { + XCTAssertTrue(q.isApproximatelyEqual(to: .log(onePlus: .expMinusOne(q)))) + } + } + func testFloat() { testExp(Float32.self) testExpMinusOne(Float32.self) @@ -262,6 +307,7 @@ final class ElementaryFunctionTests: XCTestCase { testCosSin(Float32.self) testLog(Float32.self) + testLogOnePlus(Float32.self) } func testDouble() { @@ -272,5 +318,6 @@ final class ElementaryFunctionTests: XCTestCase { testCosSin(Float64.self) testLog(Float64.self) + testLogOnePlus(Float64.self) } } From 6ab47db06b8d86d695fedd2699a1d0f8284ad170 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Thu, 28 Apr 2022 18:44:48 +0200 Subject: [PATCH 05/18] Improve quaternionic logarithm test coverage --- .../QuaternionTests/ElementaryFunctionTests.swift | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index 2f1d4539..e045f835 100644 --- a/Tests/QuaternionTests/ElementaryFunctionTests.swift +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -237,17 +237,24 @@ final class ElementaryFunctionTests: XCTestCase { // MARK: - log-like functions func testLog(_ type: T.Type) { - // log(0) = undefined/infinity + // log(0) = infinity XCTAssertFalse(Quaternion.log(Quaternion(real: 0, imaginary: 0, 0, 0)).isFinite) XCTAssertFalse(Quaternion.log(Quaternion(real:-0, imaginary: 0, 0, 0)).isFinite) XCTAssertFalse(Quaternion.log(Quaternion(real:-0, imaginary:-0,-0,-0)).isFinite) XCTAssertFalse(Quaternion.log(Quaternion(real: 0, imaginary:-0,-0,-0)).isFinite) + XCTAssertFalse(Quaternion.log(Quaternion(real: .nan, imaginary: .nan, .nan, .nan)).isFinite) + // For randomly-chosen well-scaled finite values, we expect to have + // log(exp(q)) ≈ q var g = SystemRandomNumberGenerator() - let values: [Quaternion] = (0..<100).map { _ in + let values: [Quaternion] = (0..<1000).map { _ in Quaternion( - real: T.random(in: -1 ... 1, using: &g), - imaginary: SIMD3(repeating: T.random(in: -.pi ... .pi, using: &g) / 3)) + real: T.random(in: -2 ... 2, using: &g), + imaginary: + T.random(in: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g) + ) } for q in values { XCTAssertTrue(q.isApproximatelyEqual(to: .log(.exp(q)))) From 1953fb3c0d9b27e1a8f9714ccc61ca2868e6a706 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Mon, 2 May 2022 09:50:06 +0200 Subject: [PATCH 06/18] Update quaternionic elfns test cases --- .../QuaternionModule/ImaginaryHelper.swift | 6 + .../ElementaryFunctionTests.swift | 104 +++++++++--------- 2 files changed, 58 insertions(+), 52 deletions(-) diff --git a/Sources/QuaternionModule/ImaginaryHelper.swift b/Sources/QuaternionModule/ImaginaryHelper.swift index 64c4bc04..9033f0b5 100644 --- a/Sources/QuaternionModule/ImaginaryHelper.swift +++ b/Sources/QuaternionModule/ImaginaryHelper.swift @@ -25,6 +25,12 @@ extension SIMD3 where Scalar: FloatingPoint { SIMD3(repeating: .nan) } + /// Returns a vector with .ulpOfOne in all lanes + @usableFromInline @inline(__always) + internal static var ulpOfOne: Self { + SIMD3(repeating: .ulpOfOne) + } + /// True if all values of this instance are finite @usableFromInline @inline(__always) internal var isFinite: Bool { diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index e045f835..88d28c87 100644 --- a/Tests/QuaternionTests/ElementaryFunctionTests.swift +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -21,10 +21,10 @@ 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))) + 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))) // 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 @@ -32,11 +32,11 @@ final class ElementaryFunctionTests: XCTestCase { // 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: .zero, 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: .zero, 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) @@ -79,10 +79,10 @@ 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))) + XCTAssertTrue(Quaternion.expMinusOne(Quaternion(real: .zero, imaginary: .zero)).isZero) + XCTAssertTrue(Quaternion.expMinusOne(Quaternion(real:-.zero, imaginary: .zero)).isZero) + XCTAssertTrue(Quaternion.expMinusOne(Quaternion(real:-.zero, imaginary:-.zero)).isZero) + XCTAssertTrue(Quaternion.expMinusOne(Quaternion(real: .zero, imaginary:-.zero)).isZero) // 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. @@ -90,11 +90,11 @@ final class ElementaryFunctionTests: XCTestCase { // 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: .zero, 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: .zero, 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) @@ -126,18 +126,18 @@ 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: .zero, 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: .zero, 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) @@ -162,18 +162,18 @@ 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))) + XCTAssertTrue(Quaternion.sinh(Quaternion(real: .zero, imaginary: .zero)).isZero) + XCTAssertTrue(Quaternion.sinh(Quaternion(real:-.zero, imaginary: .zero)).isZero) + XCTAssertTrue(Quaternion.sinh(Quaternion(real:-.zero, imaginary:-.zero)).isZero) + XCTAssertTrue(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: .zero, 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: .zero, 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) @@ -238,10 +238,10 @@ final class ElementaryFunctionTests: XCTestCase { func testLog(_ type: T.Type) { // log(0) = infinity - XCTAssertFalse(Quaternion.log(Quaternion(real: 0, imaginary: 0, 0, 0)).isFinite) - XCTAssertFalse(Quaternion.log(Quaternion(real:-0, imaginary: 0, 0, 0)).isFinite) - XCTAssertFalse(Quaternion.log(Quaternion(real:-0, imaginary:-0,-0,-0)).isFinite) - XCTAssertFalse(Quaternion.log(Quaternion(real: 0, imaginary:-0,-0,-0)).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) + XCTAssertFalse(Quaternion.log(Quaternion(real: .zero, imaginary:-.zero)).isFinite) XCTAssertFalse(Quaternion.log(Quaternion(real: .nan, imaginary: .nan, .nan, .nan)).isFinite) // For randomly-chosen well-scaled finite values, we expect to have @@ -263,31 +263,31 @@ final class ElementaryFunctionTests: XCTestCase { func testLogOnePlus(_ type: T.Type) { // log(onePlus: 0) = 0 - XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real: 0, imaginary: 0, 0, 0)).isZero) - XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real:-0, imaginary: 0, 0, 0)).isZero) - XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real:-0, imaginary:-0,-0,-0)).isZero) - XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real: 0, imaginary:-0,-0,-0)).isZero) + XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: .zero)).isZero) + XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real:-.zero, imaginary: .zero)).isZero) + XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real:-.zero, imaginary:-.zero)).isZero) + XCTAssertTrue(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, .nan, .nan)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: .nan, .nan, .nan)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: .nan, .nan, .nan)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.infinity, imaginary: .nan, .nan, .nan)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: -.infinity, -.infinity, -.infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: -.infinity, -.infinity, -.infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: -.infinity, -.infinity, -.infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.infinity, imaginary: -.infinity, -.infinity, -.infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.ulpOfOne, imaginary: -.infinity, -.infinity, -.infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: .zero, .zero, .zero)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.infinity, imaginary: .zero, .zero, .zero)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: .zero, .zero, .zero)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: .infinity, .infinity, .infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: .infinity, .infinity, .infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: .infinity, .infinity, .infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.infinity, imaginary: .infinity, .infinity, .infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.ulpOfOne, imaginary: .infinity, .infinity, .infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .nan, imaginary: -.ulpOfOne, -.ulpOfOne, -.ulpOfOne)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: -.infinity, imaginary: -.ulpOfOne, -.ulpOfOne, -.ulpOfOne)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .infinity, imaginary: -.ulpOfOne, -.ulpOfOne, -.ulpOfOne)).isFinite) + 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: .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: .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: .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 From f54e18bb1ae8676147d4272c2c0ad256d4fa2cba Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Mon, 2 May 2022 10:00:15 +0200 Subject: [PATCH 07/18] Add quaternionic inverse hyperbolic functions --- .../Quaternion+ElementaryFunctions.swift | 48 ++++++-- .../ElementaryFunctionTests.swift | 108 ++++++++++++++++++ 2 files changed, 148 insertions(+), 8 deletions(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index c4f4b3dd..6a2f5b80 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -326,13 +326,45 @@ extension Quaternion/*: ElementaryFunctions*/ { return Quaternion(real: .log(onePlus: s)/2, imaginary: imaginary) } - // - // MARK: - pow-like functions + @inlinable + public static func acosh(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of the + // quaternionic `log` and `sqrt` operations: + // + // ``` + // acosh(q) = log(q + sqrt(q² - 1)) + // ``` + log(q + .sqrt(q*q - .one)) + } + + @inlinable + public static func asinh(_ q: Quaternion) -> Quaternion { + // Mathematically, this operation can be expanded in terms of the + // quaternionic `log` and `sqrt` operations: + // + // ``` + // asinh(q) = log(q + sqrt(q² + 1)) + // ``` + log(q + .sqrt(q*q + .one)) + } + @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 as follows: + // quaternionic `exp` and `log` operations: // // ``` // pow(q, p) = exp(log(pow(q, p))) @@ -344,7 +376,7 @@ extension Quaternion/*: ElementaryFunctions*/ { @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 as follows: + // quaternionic `exp` and `log` operations: // // ``` // pow(q, n) = exp(log(pow(q, n))) @@ -360,7 +392,7 @@ extension Quaternion/*: ElementaryFunctions*/ { @inlinable public static func sqrt(_ q: Quaternion) -> Quaternion { // Mathematically, this operation can be expanded in terms of the - // quaternionic `exp` and `log` operations as follows: + // quaternionic `exp` and `log` operations: // // ``` // sqrt(q) = q^(1/2) = exp(log(q^(1/2))) @@ -373,11 +405,11 @@ extension Quaternion/*: ElementaryFunctions*/ { @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 as follows: + // quaternionic `exp` and `log` operations: // // ``` - // root(q, n) = exp(log(root(q, n))) - // = exp(log(q) / n) + // 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 diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index 88d28c87..7dea6526 100644 --- a/Tests/QuaternionTests/ElementaryFunctionTests.swift +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -306,6 +306,108 @@ final class ElementaryFunctionTests: XCTestCase { } } + func testAcosh(_ type: T.Type) { + // acosh(0) = 0 + XCTAssertTrue(Quaternion.acosh(0).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: .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: .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: .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: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g) + ) + } + for q in values { + XCTAssertTrue(q.isApproximatelyEqual(to: .cosh(.acosh(q)))) + } + } + + func testAsinh(_ type: T.Type) { + // asinh(0) = 0 + XCTAssertTrue(Quaternion.asinh(0).isZero) + // 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: .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: .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: .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: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g) + ) + } + for q in values { + XCTAssertTrue(q.isApproximatelyEqual(to: .sinh(.asinh(q)))) + } + } + + func testAtanh(_ type: T.Type) { + // For randomly-chosen well-scaled finite values, we expect to have + // atanh(tanh(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: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g) + ) + } + for q in values { + XCTAssertTrue(q.isApproximatelyEqual(to: .atanh(.tanh(q)))) + } + } + func testFloat() { testExp(Float32.self) testExpMinusOne(Float32.self) @@ -315,6 +417,9 @@ final class ElementaryFunctionTests: XCTestCase { testLog(Float32.self) testLogOnePlus(Float32.self) + testAcosh(Float32.self) + testAsinh(Float32.self) + testAtanh(Float32.self) } func testDouble() { @@ -326,5 +431,8 @@ final class ElementaryFunctionTests: XCTestCase { testLog(Float64.self) testLogOnePlus(Float64.self) + testAcosh(Float64.self) + testAsinh(Float64.self) + testAtanh(Float64.self) } } From c6f285353959f7d096d1f485310db69a5e39de83 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Mon, 2 May 2022 14:16:47 +0200 Subject: [PATCH 08/18] Add more test cases and cleanup code formatting --- .../ElementaryFunctionTests.swift | 377 ++++++++++++------ 1 file changed, 249 insertions(+), 128 deletions(-) diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index 7dea6526..5586a4cf 100644 --- a/Tests/QuaternionTests/ElementaryFunctionTests.swift +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -25,24 +25,34 @@ final class ElementaryFunctionTests: XCTestCase { 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))) - // 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: .zero, 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: .zero, 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) + // 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 - XCTAssertTrue(Quaternion.expMinusOne(Quaternion(real: .zero, imaginary: .zero)).isZero) - XCTAssertTrue(Quaternion.expMinusOne(Quaternion(real:-.zero, imaginary: .zero)).isZero) - XCTAssertTrue(Quaternion.expMinusOne(Quaternion(real:-.zero, imaginary:-.zero)).isZero) - XCTAssertTrue(Quaternion.expMinusOne(Quaternion(real: .zero, imaginary:-.zero)).isZero) - // 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: .zero, 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: .zero, 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))) @@ -131,19 +151,33 @@ final class ElementaryFunctionTests: XCTestCase { 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: .zero, 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: .zero, 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 - XCTAssertTrue(Quaternion.sinh(Quaternion(real: .zero, imaginary: .zero)).isZero) - XCTAssertTrue(Quaternion.sinh(Quaternion(real:-.zero, imaginary: .zero)).isZero) - XCTAssertTrue(Quaternion.sinh(Quaternion(real:-.zero, imaginary:-.zero)).isZero) - XCTAssertTrue(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) + 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: .zero, 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: .zero, 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) @@ -242,7 +290,34 @@ final class ElementaryFunctionTests: XCTestCase { 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: .nan, imaginary: .nan, .nan, .nan)).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 @@ -257,35 +332,42 @@ final class ElementaryFunctionTests: XCTestCase { ) } for q in values { - XCTAssertTrue(q.isApproximatelyEqual(to: .log(.exp(q)))) + XCTAssert(q.isApproximatelyEqual(to: .log(.exp(q)))) } } func testLogOnePlus(_ type: T.Type) { // log(onePlus: 0) = 0 - XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary: .zero)).isZero) - XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real:-.zero, imaginary: .zero)).isZero) - XCTAssertTrue(Quaternion.log(onePlus: Quaternion(real:-.zero, imaginary:-.zero)).isZero) - XCTAssertTrue(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) + 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: .nan, imaginary:-.infinity)).isFinite) - XCTAssertFalse(Quaternion.log(onePlus: Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + 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: .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: .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) @@ -302,35 +384,41 @@ final class ElementaryFunctionTests: XCTestCase { ) } for q in values { - XCTAssertTrue(q.isApproximatelyEqual(to: .log(onePlus: .expMinusOne(q)))) + XCTAssert(q.isApproximatelyEqual(to: .log(onePlus: .expMinusOne(q)))) } } func testAcosh(_ type: T.Type) { - // acosh(0) = 0 - XCTAssertTrue(Quaternion.acosh(0).isZero) + // 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: .nan, imaginary:-.infinity)).isFinite) - XCTAssertFalse(Quaternion.acosh(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + 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: .nan, imaginary:-.ulpOfOne)).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() @@ -344,35 +432,54 @@ final class ElementaryFunctionTests: XCTestCase { ) } for q in values { - XCTAssertTrue(q.isApproximatelyEqual(to: .cosh(.acosh(q)))) + 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)) + XCTAssertEqual(Quaternion.asin(1).imaginary, SIMD3(repeating: 0)) // asinh(0) = 0 - XCTAssertTrue(Quaternion.asinh(0).isZero) + XCTAssert(Quaternion.asinh(0).isZero) + // asinh(-1) = -π/2 + XCTAssert(Quaternion.asin(-1).real.isApproximatelyEqual(to: -.pi/2)) + XCTAssertEqual(Quaternion.asin(-1).imaginary, SIMD3(repeating: 0)) // 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: .nan, imaginary:-.infinity)).isFinite) - XCTAssertFalse(Quaternion.asinh(Quaternion(real: .zero, imaginary:-.infinity)).isFinite) + 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: .nan, imaginary:-.ulpOfOne)).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() @@ -386,13 +493,20 @@ final class ElementaryFunctionTests: XCTestCase { ) } for q in values { - XCTAssertTrue(q.isApproximatelyEqual(to: .sinh(.asinh(q)))) + 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 - // atanh(tanh(z)) ≈ z + // tanh(atanh(q)) ≈ q var g = SystemRandomNumberGenerator() let values: [Quaternion] = (0..<1000).map { _ in Quaternion( @@ -404,7 +518,14 @@ final class ElementaryFunctionTests: XCTestCase { ) } for q in values { - XCTAssertTrue(q.isApproximatelyEqual(to: .atanh(.tanh(q)))) + 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() + } } } From 55bfe0e8d45503c21fe8b150cc51dba5912493f9 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Mon, 2 May 2022 14:21:31 +0200 Subject: [PATCH 09/18] Add quaternionic acos, asin and atan Adds quaternionic inverse cosine, inverse sine and inverse tangent. Specialized inverse hyperbolic cosine and inverse hyperbolic sine based on these functions. --- .../Quaternion+ElementaryFunctions.swift | 85 +++++++++---- .../ElementaryFunctionTests.swift | 119 +++++++++++++++++- 2 files changed, 177 insertions(+), 27 deletions(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index 6a2f5b80..84dce5bb 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -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` (`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) @@ -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` (`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 @@ -326,32 +325,67 @@ extension Quaternion/*: ElementaryFunctions*/ { 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 `log` and `sqrt` operations: + // Mathematically, this operation can be expanded in terms of + // the quaternionic `acos` operation (`let θ = ||v||`): // // ``` - // acosh(q) = log(q + sqrt(q² - 1)) + // acosh(q) = (v/θ) * acos(q) // ``` - log(q + .sqrt(q*q - .one)) + 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 `log` and `sqrt` operations: + // Mathematically, this operation can be expanded in terms of + // the quaternionic `asin` operation (`let θ = ||v||`): // // ``` - // asinh(q) = log(q + sqrt(q² + 1)) + // sin(q) = -(v/θ) * asin(q * (v/θ))) // ``` - log(q + .sqrt(q*q + .one)) + 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: + // Mathematically, this operation can be expanded in terms of + // the quaternionic `log` operation: // // ``` // atanh(q) = (log(1 + q) - log(1 - q))/2 @@ -363,8 +397,8 @@ extension Quaternion/*: ElementaryFunctions*/ { // 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: + // Mathematically, this operation can be expanded in terms of + // the quaternionic `exp` and `log` operations: // // ``` // pow(q, p) = exp(log(pow(q, p))) @@ -375,8 +409,8 @@ extension Quaternion/*: ElementaryFunctions*/ { @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: + // Mathematically, this operation can be expanded in terms of + // the quaternionic `exp` and `log` operations: // // ``` // pow(q, n) = exp(log(pow(q, n))) @@ -391,8 +425,8 @@ extension Quaternion/*: ElementaryFunctions*/ { @inlinable public static func sqrt(_ q: Quaternion) -> Quaternion { - // Mathematically, this operation can be expanded in terms of the - // quaternionic `exp` and `log` operations: + // Mathematically, this operation can be expanded in terms of + // the quaternionic `exp` and `log` operations: // // ``` // sqrt(q) = q^(1/2) = exp(log(q^(1/2))) @@ -404,8 +438,8 @@ extension Quaternion/*: ElementaryFunctions*/ { @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: + // 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))) @@ -420,7 +454,6 @@ extension Quaternion/*: ElementaryFunctions*/ { } extension SIMD3 where Scalar: FloatingPoint { - /// Returns the normalized axis and the length of this vector. @usableFromInline @inline(__always) internal var unitAxisAndLength: (Self, Scalar) { diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index 5586a4cf..039d0235 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 @@ -388,8 +388,121 @@ final class ElementaryFunctionTests: XCTestCase { } } + 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: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/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: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/2, using: &g), + T.random(in: -.pi/2 ... .pi/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 + XCTAssertEqual(Quaternion.acosh(1).imaginary, .zero) XCTAssert(Quaternion.acosh(1).isZero) // acosh is the identity at infinity. XCTAssertFalse(Quaternion.acosh(Quaternion(real: .nan, imaginary: .nan)).isFinite) @@ -538,6 +651,8 @@ final class ElementaryFunctionTests: XCTestCase { testLog(Float32.self) testLogOnePlus(Float32.self) + testAcos(Float32.self) + testAsin(Float32.self) testAcosh(Float32.self) testAsinh(Float32.self) testAtanh(Float32.self) @@ -552,6 +667,8 @@ final class ElementaryFunctionTests: XCTestCase { testLog(Float64.self) testLogOnePlus(Float64.self) + testAcos(Float64.self) + testAsin(Float64.self) testAcosh(Float64.self) testAsinh(Float64.self) testAtanh(Float64.self) From 4f3ef185f5aa649265c3874edb01d9eaec5d4f71 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Mon, 2 May 2022 17:22:21 +0200 Subject: [PATCH 10/18] Add quaternionic pow, sqrt and root --- .../Quaternion+ElementaryFunctions.swift | 3 +- .../ElementaryFunctionTests.swift | 256 +++++++++++++++++- 2 files changed, 256 insertions(+), 3 deletions(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index 84dce5bb..aa7d50c9 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -404,7 +404,8 @@ extension Quaternion: ElementaryFunctions { // pow(q, p) = exp(log(pow(q, p))) // = exp(p * log(q)) // ``` - exp(p * log(q)) + guard !q.isZero else { return .zero } + return exp(p * log(q)) } @inlinable diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index 039d0235..6274ff35 100644 --- a/Tests/QuaternionTests/ElementaryFunctionTests.swift +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -318,7 +318,6 @@ final class ElementaryFunctionTests: XCTestCase { 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() @@ -370,7 +369,6 @@ final class ElementaryFunctionTests: XCTestCase { 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() @@ -642,6 +640,250 @@ final class ElementaryFunctionTests: XCTestCase { } } + // MARK: - pow-like functions + + func testPowQuaternion(_ type: T.Type) { + // pow(0, 0) = 0 + let zero: Quaternion = .zero + XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary:-.zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero).isZero) + // 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) = 0 + let zero: Int = .zero + XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary:-.zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero).isZero) + XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero).isZero) + // pow(0, x) = 0 for x > 0 + let n: Int = 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) @@ -656,6 +898,11 @@ final class ElementaryFunctionTests: XCTestCase { testAcosh(Float32.self) testAsinh(Float32.self) testAtanh(Float32.self) + + testPowQuaternion(Float32.self) + testPowInt(Float32.self) + testSqrt(Float32.self) + testRoot(Float32.self) } func testDouble() { @@ -672,5 +919,10 @@ final class ElementaryFunctionTests: XCTestCase { testAcosh(Float64.self) testAsinh(Float64.self) testAtanh(Float64.self) + + testPowQuaternion(Float64.self) + testPowInt(Float64.self) + testSqrt(Float64.self) + testRoot(Float64.self) } } From 6dff4a68958c28751b329a72ad7bb227d2d54457 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Mon, 2 May 2022 17:29:40 +0200 Subject: [PATCH 11/18] Seperate test helper from quaternion module --- .../QuaternionModule/ImaginaryHelper.swift | 12 -------- .../QuaternionTests/ImaginaryTestHelper.swift | 28 +++++++++++++++++++ .../QuaternionTests/TransformationTests.swift | 7 +---- 3 files changed, 29 insertions(+), 18 deletions(-) create mode 100644 Tests/QuaternionTests/ImaginaryTestHelper.swift diff --git a/Sources/QuaternionModule/ImaginaryHelper.swift b/Sources/QuaternionModule/ImaginaryHelper.swift index 9033f0b5..9ffd4cd2 100644 --- a/Sources/QuaternionModule/ImaginaryHelper.swift +++ b/Sources/QuaternionModule/ImaginaryHelper.swift @@ -19,18 +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) - } - - /// Returns a vector with .ulpOfOne in all lanes - @usableFromInline @inline(__always) - internal static var ulpOfOne: Self { - SIMD3(repeating: .ulpOfOne) - } - /// True if all values of this instance are finite @usableFromInline @inline(__always) internal var isFinite: Bool { diff --git a/Tests/QuaternionTests/ImaginaryTestHelper.swift b/Tests/QuaternionTests/ImaginaryTestHelper.swift new file mode 100644 index 00000000..9fa60a46 --- /dev/null +++ b/Tests/QuaternionTests/ImaginaryTestHelper.swift @@ -0,0 +1,28 @@ +//===--- 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 +// +//===----------------------------------------------------------------------===// + +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 From 514418bcc0d0237f49d3deed88e1f5abd272223c Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Tue, 3 May 2022 08:45:59 +0200 Subject: [PATCH 12/18] Hide implementation detail of unit axis and argument --- Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index aa7d50c9..55d55764 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -456,8 +456,8 @@ extension Quaternion: ElementaryFunctions { 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), From b546eb892b8048105c1545cf4e67e85cdfaf75a9 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Tue, 3 May 2022 08:46:22 +0200 Subject: [PATCH 13/18] Improve test cases and coverage --- .../ElementaryFunctionTests.swift | 51 +++++++++---------- 1 file changed, 25 insertions(+), 26 deletions(-) diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index 6274ff35..9caa94c0 100644 --- a/Tests/QuaternionTests/ElementaryFunctionTests.swift +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -325,13 +325,13 @@ final class ElementaryFunctionTests: XCTestCase { Quaternion( real: T.random(in: -2 ... 2, using: &g), imaginary: - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g) + 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: .log(.exp(q)))) + XCTAssert(q.isApproximatelyEqual(to: .exp(.log(q)))) } } @@ -376,13 +376,13 @@ final class ElementaryFunctionTests: XCTestCase { Quaternion( real: T.random(in: -2 ... 2, using: &g), imaginary: - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g) + 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: .log(onePlus: .expMinusOne(q)))) + XCTAssert(q.isApproximatelyEqual(to: .expMinusOne(.log(onePlus: q)))) } } @@ -430,9 +430,9 @@ final class ElementaryFunctionTests: XCTestCase { Quaternion( real: T.random(in: -2 ... 2, using: &g), imaginary: - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g) + 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 { @@ -486,9 +486,9 @@ final class ElementaryFunctionTests: XCTestCase { Quaternion( real: T.random(in: -2 ... 2, using: &g), imaginary: - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g) + 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 { @@ -500,7 +500,6 @@ final class ElementaryFunctionTests: XCTestCase { func testAcosh(_ type: T.Type) { // acosh(1) = 0 - XCTAssertEqual(Quaternion.acosh(1).imaginary, .zero) XCTAssert(Quaternion.acosh(1).isZero) // acosh is the identity at infinity. XCTAssertFalse(Quaternion.acosh(Quaternion(real: .nan, imaginary: .nan)).isFinite) @@ -537,9 +536,9 @@ final class ElementaryFunctionTests: XCTestCase { Quaternion( real: T.random(in: -2 ... 2, using: &g), imaginary: - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g) + 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 { @@ -557,12 +556,12 @@ final class ElementaryFunctionTests: XCTestCase { func testAsinh(_ type: T.Type) { // asinh(1) = π/2 XCTAssert(Quaternion.asin(1).real.isApproximatelyEqual(to: .pi/2)) - XCTAssertEqual(Quaternion.asin(1).imaginary, SIMD3(repeating: 0)) + 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)) - XCTAssertEqual(Quaternion.asin(-1).imaginary, SIMD3(repeating: 0)) + 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) @@ -598,9 +597,9 @@ final class ElementaryFunctionTests: XCTestCase { Quaternion( real: T.random(in: -2 ... 2, using: &g), imaginary: - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g) + 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 { @@ -623,9 +622,9 @@ final class ElementaryFunctionTests: XCTestCase { Quaternion( real: T.random(in: -2 ... 2, using: &g), imaginary: - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g), - T.random(in: -.pi/2 ... .pi/2, using: &g) + 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 { From cf6188523729d8af67c3bb4b6e1af83b0865ec31 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Tue, 3 May 2022 15:54:13 +0200 Subject: [PATCH 14/18] Specialised quaternionic sqrt --- .../Quaternion+ElementaryFunctions.swift | 48 +++++++++++++++---- 1 file changed, 39 insertions(+), 9 deletions(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index 55d55764..973f7783 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -426,15 +426,45 @@ extension Quaternion: ElementaryFunctions { @inlinable public static func sqrt(_ q: Quaternion) -> Quaternion { - // Mathematically, this operation can be expanded in terms of - // the quaternionic `exp` and `log` operations: - // - // ``` - // sqrt(q) = q^(1/2) = exp(log(q^(1/2))) - // = exp(log(q) * (1/2)) - // ``` - guard !q.isZero else { return .zero } - return exp(log(q).divided(by: 2)) + 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 From 4b13da31a4ac03170b7e63cefc67f03c58d19fd1 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Tue, 3 May 2022 16:08:14 +0200 Subject: [PATCH 15/18] Improve argument calculation in quaternionic log --- Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index 973f7783..5b3fb190 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -225,10 +225,8 @@ extension Quaternion: ElementaryFunctions { guard q.isFinite && !q.isZero else { return .infinity } // Having eliminated non-finite values and zero, the imaginary part is // straightforward: - // TODO: There is a potential optimisation hidden here, as length is - // calculated twice (halfAngle, unitAxisAndLength) - let argument = q.halfAngle 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 From 442d4a75dd023eb59dc30bdcfe867f081b5d6280 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Tue, 3 May 2022 16:12:28 +0200 Subject: [PATCH 16/18] Fix license header --- Tests/QuaternionTests/ImaginaryTestHelper.swift | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Tests/QuaternionTests/ImaginaryTestHelper.swift b/Tests/QuaternionTests/ImaginaryTestHelper.swift index 9fa60a46..0e51b059 100644 --- a/Tests/QuaternionTests/ImaginaryTestHelper.swift +++ b/Tests/QuaternionTests/ImaginaryTestHelper.swift @@ -1,8 +1,8 @@ -//===--- ImaginaryHelper.swift --------------------------------*- swift -*-===// +//===--- ImaginaryTestHelper.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 +// 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 From 4af024b2333023f3e4e6a8090f33d52d694c7999 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Thu, 16 Jun 2022 09:29:30 +0200 Subject: [PATCH 17/18] Remove superfluous type definition of quaternion sqrt --- Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index 5b3fb190..a20b4a04 100644 --- a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift +++ b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift @@ -423,7 +423,7 @@ extension Quaternion: ElementaryFunctions { } @inlinable - public static func sqrt(_ q: Quaternion) -> Quaternion { + 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||`): From f975305391c6c1607b7fc35437b0615776721121 Mon Sep 17 00:00:00 2001 From: Markus Winter Date: Wed, 25 Jan 2023 20:51:35 +0100 Subject: [PATCH 18/18] Improvement to Quaternion.pow edge cases --- .../Quaternion+ElementaryFunctions.swift | 15 ++++--- .../ElementaryFunctionTests.swift | 43 +++++++++++-------- 2 files changed, 35 insertions(+), 23 deletions(-) diff --git a/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift b/Sources/QuaternionModule/Quaternion+ElementaryFunctions.swift index a20b4a04..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 @@ -252,7 +252,10 @@ extension Quaternion: ElementaryFunctions { // 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) + 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 @@ -402,7 +405,7 @@ extension Quaternion: ElementaryFunctions { // pow(q, p) = exp(log(pow(q, p))) // = exp(p * log(q)) // ``` - guard !q.isZero else { return .zero } + if q.isZero { return p.real > 0 ? zero : infinity } return exp(p * log(q)) } @@ -415,7 +418,7 @@ extension Quaternion: ElementaryFunctions { // pow(q, n) = exp(log(pow(q, n))) // = exp(log(q) * n) // ``` - guard !q.isZero else { return .zero } + 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. @@ -431,8 +434,8 @@ extension Quaternion: ElementaryFunctions { // 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. + // 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) diff --git a/Tests/QuaternionTests/ElementaryFunctionTests.swift b/Tests/QuaternionTests/ElementaryFunctionTests.swift index 9caa94c0..ac026a82 100644 --- a/Tests/QuaternionTests/ElementaryFunctionTests.swift +++ b/Tests/QuaternionTests/ElementaryFunctionTests.swift @@ -642,15 +642,15 @@ final class ElementaryFunctionTests: XCTestCase { // MARK: - pow-like functions func testPowQuaternion(_ type: T.Type) { - // pow(0, 0) = 0 + // pow(0, 0) = inf let zero: Quaternion = .zero - XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary:-.zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero).isZero) + 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) @@ -713,17 +713,26 @@ final class ElementaryFunctionTests: XCTestCase { } func testPowInt(_ type: T.Type) { - // pow(0, 0) = 0 + // pow(0, 0) = 1 let zero: Int = .zero - XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary: .zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary: .zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real:-.zero, imaginary:-.zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero).isZero) - XCTAssert(Quaternion.pow(Quaternion(real: .zero, imaginary:-.zero), zero).isZero) + 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 - let n: Int = 2 + 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)