Skip to content

Commit 68a86c2

Browse files
markuswntrmarkuswntr
committed
Add cosh, sinh and tanh to quaternions
1 parent f895b11 commit 68a86c2

File tree

2 files changed

+155
-0
lines changed

2 files changed

+155
-0
lines changed

‎Sources/QuaternionModule/ElementaryFunctions.swift

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -76,5 +76,70 @@ extension Quaternion/*: ElementaryFunctions */ {
7676
imaginary: n̂ * .exp(q.real) * .sin(θ)
7777
)
7878
}
79+
80+
// cosh(r + xi + yj + zk) = cosh(r + v)
81+
// cosh(r + v) = cosh(r) cos(||v||) + (v/||v||) sinh(r) sin(||v||).
82+
//
83+
// See cosh on complex numbers for algorithm details.
84+
@inlinable
85+
public static func cosh(_ q: Quaternion) -> Quaternion {
86+
guard q.isFinite else { return q }
87+
// TODO: Replace q.imaginary == .zero with `q.isReal`
88+
let θ = q.imaginary == .zero ? .zero : q.imaginary.length // θ = ||v||
89+
let n̂ = q.imaginary == .zero ? .zero : (q.imaginary / θ) // n̂ = v / ||v||
90+
guard q.real.magnitude < -RealType.log(.ulpOfOne) else {
91+
let rotation = Quaternion(halfAngle: θ, unitAxis: n̂)
92+
let firstScale = RealType.exp(q.real.magnitude/2)
93+
let secondScale = firstScale/2
94+
return rotation.multiplied(by: firstScale).multiplied(by: secondScale)
95+
}
96+
return Quaternion(
97+
real: .cosh(q.real) * .cos(θ),
98+
imaginary: n̂ * .sinh(q.real) * .sin(θ)
99+
)
100+
}
101+
102+
// sinh(r + xi + yj + zk) = sinh(r + v)
103+
// sinh(r + v) = sinh(r) cos(||v||) + (v/||v||) cosh(r) sin(||v||)
104+
//
105+
// See cosh on complex numbers for algorithm details.
106+
@inlinable
107+
public static func sinh(_ q: Quaternion) -> Quaternion {
108+
guard q.isFinite else { return q }
109+
// TODO: Replace q.imaginary == .zero with `q.isReal`
110+
let θ = q.imaginary == .zero ? .zero : q.imaginary.length // θ = ||v||
111+
let n̂ = q.imaginary == .zero ? .zero : (q.imaginary / θ) // n̂ = v / ||v||
112+
guard q.real.magnitude < -RealType.log(.ulpOfOne) else {
113+
let rotation = Quaternion(halfAngle: θ, unitAxis: n̂)
114+
let firstScale = RealType.exp(q.real.magnitude/2)
115+
let secondScale = RealType(signOf: q.real, magnitudeOf: firstScale/2)
116+
return rotation.multiplied(by: firstScale).multiplied(by: secondScale)
117+
}
118+
return Quaternion(
119+
real: .sinh(q.real) * .cos(θ),
120+
imaginary: n̂ * .cosh(q.real) * .sin(θ)
121+
)
122+
}
123+
124+
// tanh(q) = sinh(q) / cosh(q)
125+
//
126+
// See tanh on complex numbers for algorithm details.
127+
@inlinable
128+
public static func tanh(_ q: Quaternion) -> Quaternion {
129+
guard q.isFinite else { return q }
130+
// Note that when |r| is larger than -log(.ulpOfOne),
131+
// sinh(r + v) == ±cosh(r + v), so tanh(r + v) is just ±1.
132+
guard q.real.magnitude < -RealType.log(.ulpOfOne) else {
133+
return Quaternion(
134+
real: RealType(signOf: q.real, magnitudeOf: 1),
135+
imaginary:
136+
RealType(signOf: q.imaginary.x, magnitudeOf: 0),
137+
RealType(signOf: q.imaginary.y, magnitudeOf: 0),
138+
RealType(signOf: q.imaginary.z, magnitudeOf: 0)
139+
)
140+
}
141+
return sinh(q) / cosh(q)
142+
}
143+
79144
}
80145
}

‎Tests/QuaternionTests/ElementaryFunctionTests.swift

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -119,13 +119,103 @@ final class ElementaryFunctionTests: XCTestCase {
119119
}
120120
}
121121

122+
func testCosh<T: Real & FixedWidthFloatingPoint & SIMDScalar>(_ type: T.Type) {
123+
// cosh(0) = 1
124+
XCTAssertEqual(1, Quaternion<T>.cosh(Quaternion(real: 0, imaginary: 0, 0, 0)))
125+
XCTAssertEqual(1, Quaternion<T>.cosh(Quaternion(real:-0, imaginary: 0, 0, 0)))
126+
XCTAssertEqual(1, Quaternion<T>.cosh(Quaternion(real:-0, imaginary:-0,-0,-0)))
127+
XCTAssertEqual(1, Quaternion<T>.cosh(Quaternion(real: 0, imaginary:-0,-0,-0)))
128+
// cosh is the identity at infinity.
129+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: .infinity, imaginary: .zero)).isFinite)
130+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: .infinity, imaginary: .infinity)).isFinite)
131+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: 0, imaginary: .infinity)).isFinite)
132+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: -.infinity, imaginary: .infinity)).isFinite)
133+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: -.infinity, imaginary: .zero)).isFinite)
134+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: -.infinity, imaginary: -.infinity)).isFinite)
135+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: 0, imaginary: -.infinity)).isFinite)
136+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: .infinity, imaginary: -.infinity)).isFinite)
137+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: .nan, imaginary: .nan)).isFinite)
138+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: .infinity, imaginary: .nan)).isFinite)
139+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: .nan, imaginary: .infinity)).isFinite)
140+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: -.infinity, imaginary: .nan)).isFinite)
141+
XCTAssertFalse(Quaternion<T>.cosh(Quaternion(real: .nan, imaginary: -.infinity)).isFinite)
142+
// Near-overflow test, same as exp() above, but it happens later, because
143+
// for large x, cosh(x + v) ~ exp(x + v)/2.
144+
let x = T.log(.greatestFiniteMagnitude) + T.log(18/8)
145+
let mag = T.greatestFiniteMagnitude/T.sqrt(2) * (9/8)
146+
var huge = Quaternion<T>.cosh(Quaternion(real: x, imaginary: SIMD3(.pi/4, 0, 0)))
147+
XCTAssert(huge.real.isApproximatelyEqual(to: mag))
148+
XCTAssert(huge.imaginary.x.isApproximatelyEqual(to: mag))
149+
XCTAssertEqual(huge.imaginary.y, .zero)
150+
XCTAssertEqual(huge.imaginary.z, .zero)
151+
huge = Quaternion<T>.cosh(Quaternion(real: -x, imaginary: SIMD3(.pi/4, 0, 0)))
152+
XCTAssert(huge.real.isApproximatelyEqual(to: mag))
153+
XCTAssert(huge.imaginary.x.isApproximatelyEqual(to: mag))
154+
XCTAssertEqual(huge.imaginary.y, .zero)
155+
XCTAssertEqual(huge.imaginary.z, .zero)
156+
}
157+
158+
func testSinh<T: Real & FixedWidthFloatingPoint & SIMDScalar>(_ type: T.Type) {
159+
// sinh(0) = 0
160+
XCTAssertEqual(0, Quaternion<T>.sinh(Quaternion(real: 0, imaginary: 0, 0, 0)))
161+
XCTAssertEqual(0, Quaternion<T>.sinh(Quaternion(real:-0, imaginary: 0, 0, 0)))
162+
XCTAssertEqual(0, Quaternion<T>.sinh(Quaternion(real:-0, imaginary:-0,-0,-0)))
163+
XCTAssertEqual(0, Quaternion<T>.sinh(Quaternion(real: 0, imaginary:-0,-0,-0)))
164+
// sinh is the identity at infinity.
165+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: .infinity, imaginary: .zero)).isFinite)
166+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: .infinity, imaginary: .infinity)).isFinite)
167+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: 0, imaginary: .infinity)).isFinite)
168+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: -.infinity, imaginary: .infinity)).isFinite)
169+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: -.infinity, imaginary: .zero)).isFinite)
170+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: -.infinity, imaginary: -.infinity)).isFinite)
171+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: 0, imaginary: -.infinity)).isFinite)
172+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: .infinity, imaginary: -.infinity)).isFinite)
173+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: .nan, imaginary: .nan)).isFinite)
174+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: .infinity, imaginary: .nan)).isFinite)
175+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: .nan, imaginary: .infinity)).isFinite)
176+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: -.infinity, imaginary: .nan)).isFinite)
177+
XCTAssertFalse(Quaternion<T>.sinh(Quaternion(real: .nan, imaginary: -.infinity)).isFinite)
178+
// Near-overflow test, same as exp() above, but it happens later, because
179+
// for large x, sinh(x + v) ~ ±exp(x + v)/2.
180+
let x = T.log(.greatestFiniteMagnitude) + T.log(18/8)
181+
let mag = T.greatestFiniteMagnitude/T.sqrt(2) * (9/8)
182+
var huge = Quaternion<T>.sinh(Quaternion(real: x, imaginary: SIMD3(.pi/4, 0, 0)))
183+
XCTAssert(huge.real.isApproximatelyEqual(to: mag))
184+
XCTAssert(huge.imaginary.x.isApproximatelyEqual(to: mag))
185+
XCTAssertEqual(huge.imaginary.y, .zero)
186+
XCTAssertEqual(huge.imaginary.z, .zero)
187+
huge = Quaternion<T>.sinh(Quaternion(real: -x, imaginary: SIMD3(.pi/4, 0, 0)))
188+
XCTAssert(huge.real.isApproximatelyEqual(to: -mag))
189+
XCTAssert(huge.imaginary.x.isApproximatelyEqual(to: mag))
190+
XCTAssertEqual(huge.imaginary.y, .zero)
191+
XCTAssertEqual(huge.imaginary.z, .zero)
192+
// For randomly-chosen well-scaled finite values, we expect to have
193+
// cosh² - sinh² ≈ 1. Note that this test would break down due to
194+
// catastrophic cancellation as we get further away from the origin.
195+
var g = SystemRandomNumberGenerator()
196+
let values: [Quaternion<T>] = (0..<1000).map { _ in
197+
Quaternion(
198+
real: T.random(in: -2 ... 2, using: &g),
199+
imaginary: SIMD3(repeating: T.random(in: -2 ... 2, using: &g) / 3))
200+
}
201+
for q in values {
202+
let c = Quaternion.cosh(q)
203+
let s = Quaternion.sinh(q)
204+
XCTAssert((c*c - s*s).isApproximatelyEqual(to: .one))
205+
}
206+
}
207+
122208
func testFloat() {
123209
testExp(Float32.self)
124210
testExpMinusOne(Float32.self)
211+
testCosh(Float32.self)
212+
testSinh(Float32.self)
125213
}
126214

127215
func testDouble() {
128216
testExp(Float64.self)
129217
testExpMinusOne(Float64.self)
218+
testCosh(Float64.self)
219+
testSinh(Float64.self)
130220
}
131221
}

0 commit comments

Comments
 (0)