Add pow, sqrt and root to quaternions

Author: markuswntr

Sources/QuaternionModule/ElementaryFunctions.swift

@@ -163,6 +163,45 @@ extension Quaternion/*: ElementaryFunctions */ {
   public static func tan(_ q: Quaternion) -> Quaternion {
     return sin(q) / cos(q)
   }
+
+  // MARK: - pow-like functions
+  @inlinable
+  public static func pow(_ q: Quaternion, _ p: Quaternion) -> Quaternion {
+    // pow(q, p) = exp(log(q^p)) = exp(p * log(q))
+    return exp(p * log(q))
+  }
+
+  @inlinable
+  public static func pow(_ q: Quaternion, _ n: Int) -> Quaternion {
+    if q.isZero { 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.
+    //
+    // Note that this does not have the same problems that a similar
+    // implementation for a real type would have, because there's no
+    // parity/sign interaction in the complex plane.
+    return exp(log(q).multiplied(by: RealType(n)))
+  }
+
+  @inlinable
+  public static func sqrt(_ q: Quaternion) -> Quaternion<RealType> {
+    if q.isZero { return .zero }
+    // TODO: This is not the fastest implementation available
+    return exp(log(q).divided(by: 2))
+  }
+
+  @inlinable
+  public static func root(_ q: Quaternion, _ n: Int) -> Quaternion {
+    if q.isZero { 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.
+    //
+    // Note that this does not have the same problems that a similar
+    // implementation for a real type would have, because there's no
+    // parity/sign interaction in the complex plane.
+    return exp(log(q).divided(by: RealType(n)))
   }
 }
 

Tests/QuaternionTests/ElementaryFunctionTests.swift

@@ -210,8 +210,7 @@ final class ElementaryFunctionTests: XCTestCase {
   }
 
 
-  func testCosSinIdentity<T: Real & FixedWidthFloatingPoint & SIMDScalar>(_ type: T.Type) {
-    // For randomly-chosen well-scaled finite values, we expect to have cos² + sin² ≈ 1
+  func testCos<T: Real & FixedWidthFloatingPoint & SIMDScalar>(_ type: T.Type) {
     var g = SystemRandomNumberGenerator()
     let values: [Quaternion<T>] = (0..<1000).map { _ in
       Quaternion(
@@ -224,23 +223,57 @@ final class ElementaryFunctionTests: XCTestCase {
     }
     for q in values {
       let c = Quaternion.cos(q)
+
+      // For randomly-chosen well-scaled finite values, we expect to have
+      // cos ≈ (e^(q*||v||)+e^(-q*||v||)) / 2
+      let p = Quaternion(imaginary: q.imaginary / q.imaginary.length)
+      let e = (.exp(p * q) + .exp(-p * q)) / 2
+      XCTAssert(c.isApproximatelyEqual(to: e))
+    }
+  }
+
+  func testSin<T: Real & FixedWidthFloatingPoint & SIMDScalar>(_ type: T.Type) {
+    var g = SystemRandomNumberGenerator()
+    let values: [Quaternion<T>] = (0..<1000).map { _ in
+      Quaternion(
+        real: T.random(in: -2 ... 2, using: &g),
+        imaginary:
+          T.random(in: -2 ... 2, using: &g) / 3,
+          T.random(in: -2 ... 2, using: &g) / 3,
+          T.random(in: -2 ... 2, using: &g) / 3
+      )
+    }
+    for q in values {
       let s = Quaternion.sin(q)
+
+      // For randomly-chosen well-scaled finite values, we expect to have
+      // cos ≈ (e^(q*||v||)+e^(-q*||v||)) / 2
+      let p = Quaternion(imaginary: q.imaginary / q.imaginary.length)
+      let e = (.exp(p * q) - .exp(-p * q)) / (p * 2)
+      XCTAssert(s.isApproximatelyEqual(to: e))
+
+      // For randomly-chosen well-scaled finite values, we expect to have
+      // cos² + sin² ≈ 1
+      let c = Quaternion.cos(q)
       XCTAssert((c*c + s*s).isApproximatelyEqual(to: .one))
     }
   }
+
   func testFloat() {
     testExp(Float32.self)
     testExpMinusOne(Float32.self)
     testCosh(Float32.self)
     testSinh(Float32.self)
-    testCosSinIdentity(Float32.self)
+    testCos(Float32.self)
+    testSin(Float32.self)
   }
 
   func testDouble() {
     testExp(Float64.self)
     testExpMinusOne(Float64.self)
     testCosh(Float64.self)
     testSinh(Float64.self)
-    testCosSinIdentity(Float64.self)
+    testCos(Float64.self)
+    testSin(Float64.self)
   }
 }