Skip to content

Commit 40548f2

Browse files
markuswntrmarkuswntr
committed
Fix sign of tan special case and add more comments
1 parent 442f1c6 commit 40548f2

File tree

3 files changed

+178
-143
lines changed

3 files changed

+178
-143
lines changed

‎Sources/ComplexModule/ElementaryFunctions.swift‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -33,10 +33,11 @@
3333
// adapted from Kahan's 1986 paper "Branch Cuts for Complex Elementary
3434
// Functions; or: Much Ado About Nothing's Sign Bit".
3535
//
36-
// As quaternions share the same goals and use adaptations of the elementary
37-
// functions: If you make a modification to either of the following functions,
38-
// you should almost surely make a parallel modification to the same elementary
39-
// function of quaternions (See ElementaryFunctions.swift in QuaternionModule).
36+
// Elementary functions of complex numbers have many similarities with
37+
// elementary functions of quaternions and their definition in terms of real
38+
// operations. Therefore, if you make a modification to one of the following
39+
// functions, you should almost surely make a parallel modification to the same
40+
// elementary function of quaternions.
4041

4142
import RealModule
4243

‎Sources/QuaternionModule/ElementaryFunctions.swift‎

Lines changed: 172 additions & 107 deletions
Original file line numberDiff line numberDiff line change
@@ -10,103 +10,168 @@
1010
//
1111
//===----------------------------------------------------------------------===//
1212

13+
// (r + xi + yj + zk) is a common representation that is often seen for
14+
// quaternions. However, when we want to expand the elementary functions of
15+
// quaternions in terms of real operations it is almost always easier to view
16+
// them as real part (r) and imaginary vector part (v),
17+
// i.e: r + xi + yj + zk = r + v; and so we diverge a little from the
18+
// representation that is used in the documentation in other files and use this
19+
// notation of quaternions in the comments of the following functions.
20+
//
21+
// Quaternionic elementary functions have many similarities with elementary
22+
// functions of complex numbers and their definition in terms of real
23+
// operations. Therefore, if you make a modification to one of the following
24+
// functions, you should almost surely make a parallel modification to the same
25+
// elementary function of complex numbers.
26+
1327
import RealModule
1428

15-
// As the following elementary functions algorithms are adaptations of the
16-
// elementary functions of complex numbers: If you make a modification to either
17-
// of the following functions, you should almost surely make a parallel
18-
// modification to the same elementary function of complex numbers (See
19-
// ElementaryFunctions.swift in ComplexModule).
2029
extension Quaternion/*: ElementaryFunctions */ {
2130

2231
// MARK: - exp-like functions
2332

24-
// exp(r + xi + yj + zk) = exp(r + v) = exp(r) exp(v)
25-
// = exp(r) (cos(||v||) + (v/||v||) sin(||v||))
26-
//
27-
// See exp on complex numbers for algorithm details.
2833
@inlinable
2934
public static func exp(_ q: Quaternion) -> Quaternion {
35+
// Mathematically, this operation can be expanded in terms of the `Real`
36+
// operations `exp`, `cos` and `sin` as follows (`let θ = ||v||`):
37+
//
38+
// ```
39+
// exp(r + v) = exp(r) exp(v)
40+
// = exp(r) (cos(θ) + (v/θ) sin(θ))
41+
// ```
42+
//
43+
// Note that naive evaluation of this expression in floating-point would be
44+
// prone to premature overflow, since `cos` and `sin` both have magnitude
45+
// less than 1 for most inputs (i.e. `exp(r)` may be infinity when
46+
// `exp(r) cos(||v||)` would not be.
3047
guard q.isFinite else { return q }
31-
// For real quaternions we can skip phase and axis calculations
32-
let argument = q.isReal ? .zero : q.imaginary.length
33-
let axis = q.isReal ? .zero : (q.imaginary / argument)
34-
// If real < log(greatestFiniteMagnitude), then exp(q.real) does not overflow.
48+
let θ = q.imaginary.length
49+
let axis = !θ.isZero ? (q.imaginary / θ) : .zero
50+
// If real < log(greatestFiniteMagnitude), then exp(real) does not overflow.
3551
// To protect ourselves against sketchy log or exp implementations in
3652
// an unknown host library, or slight rounding disagreements between
3753
// the two, subtract one from the bound for a little safety margin.
3854
guard q.real < RealType.log(.greatestFiniteMagnitude) - 1 else {
3955
let halfScale = RealType.exp(q.real/2)
40-
let rotation = Quaternion(halfAngle: argument, unitAxis: axis)
56+
let rotation = Quaternion(halfAngle: θ, unitAxis: axis)
4157
return rotation.multiplied(by: halfScale).multiplied(by: halfScale)
4258
}
43-
return Quaternion(halfAngle: argument, unitAxis: axis).multiplied(by: .exp(q.real))
59+
return Quaternion(
60+
halfAngle: θ,
61+
unitAxis: axis
62+
).multiplied(by: .exp(q.real))
4463
}
4564

4665
@inlinable
4766
public static func expMinusOne(_ q: Quaternion) -> Quaternion {
48-
// Note that the imaginary part is just the usual exp(r) sin(argument);
49-
// the only trick is computing the real part, which allows us to borrow
50-
// the derivative of real part for this function from complex numbers.
51-
// See `expMinusOne` in the ComplexModule for implementation details.
67+
// Mathematically, this operation can be expanded in terms of the `Real`
68+
// operations `exp`, `cos` and `sin` as follows (`let θ = ||v||`):
69+
//
70+
// ```
71+
// exp(r + v) - 1 = exp(r) exp(v) - 1
72+
// = exp(r) (cos(θ) + (v/θ) sin(θ)) - 1
73+
// = exp(r) cos(θ) + exp(r) (v/θ) sin(θ) - 1
74+
// = (exp(r) cos(θ) - 1) + exp(r) (v/θ) sin(θ)
75+
// -------- u --------
76+
// ```
77+
//
78+
// Note that the imaginary part is just the usual exp(x) sin(y);
79+
// the only trick is computing the real part ("u"):
80+
//
81+
// ```
82+
// u = exp(r) cos(θ) - 1
83+
// = exp(r) cos(θ) - cos(θ) + cos(θ) - 1
84+
// = (exp(r) - 1) cos(θ) + (cos(θ) - 1)
85+
// = expMinusOne(r) cos(θ) + cosMinusOne(θ)
86+
// ```
87+
//
88+
// See `expMinusOne` on complex numbers for error bounds.
5289
guard q.isFinite else { return q }
53-
let argument = q.isReal ? .zero : q.imaginary.length
54-
let axis = q.isReal ? .zero : (q.imaginary / argument)
90+
let θ = q.imaginary.length
91+
let axis = !θ.isZero ? (q.imaginary / θ) : .zero
5592
// If exp(q) is close to the overflow boundary, we don't need to
5693
// worry about the "MinusOne" part of this function; we're just
57-
// computing exp(q). (Even when q.y is near a multiple of π/2,
58-
// it can't be close enough to overcome the scaling from exp(q.x),
59-
// so the -1 term is _always_ negligable). So we simply handle
60-
// these cases exactly the same as exp(q).
94+
// computing exp(q). (Even when θ is near a multiple of π/2,
95+
// it can't be close enough to overcome the scaling from exp(r),
96+
// so the -1 term is _always_ negligable).
6197
guard q.real < RealType.log(.greatestFiniteMagnitude) - 1 else {
6298
let halfScale = RealType.exp(q.real/2)
63-
let rotation = Quaternion(halfAngle: argument, unitAxis: axis)
99+
let rotation = Quaternion(halfAngle: θ, unitAxis: axis)
64100
return rotation.multiplied(by: halfScale).multiplied(by: halfScale)
65101
}
66102
return Quaternion(
67-
real: RealType._mulAdd(.cos(argument), .expMinusOne(q.real), .cosMinusOne(argument)),
68-
imaginary: axis * .exp(q.real) * .sin(argument)
103+
real: RealType._mulAdd(.cos(θ), .expMinusOne(q.real), .cosMinusOne(θ)),
104+
imaginary: axis * .exp(q.real) * .sin(θ)
69105
)
70106
}
71107

72-
// cosh(r + xi + yj + zk) = cosh(r + v)
73-
// = cosh(r) cos(||v||) + (v/||v||) sinh(r) sin(||v||)
74-
//
75-
// See cosh on complex numbers for algorithm details.
76108
@inlinable
77109
public static func cosh(_ q: Quaternion) -> Quaternion {
110+
// Mathematically, this operation can be expanded in terms of
111+
// trigonometric `Real` operations as follows (`let θ = ||v||`):
112+
//
113+
// ```
114+
// cosh(q) = (exp(q) + exp(-q)) / 2
115+
// = cosh(r) cos(θ) + (v/θ) sinh(r) sin(θ)
116+
// ```
117+
//
118+
// Like exp, cosh is entire, so we do not need to worry about where
119+
// branch cuts fall. Also like exp, cancellation never occurs in the
120+
// evaluation of the naive expression, so all we need to be careful
121+
// about is the behavior near the overflow boundary.
122+
//
123+
// Fortunately, if |x| >= -log(ulpOfOne), cosh(x) and sinh(x) are
124+
// both just exp(|x|)/2, and we already know how to compute that.
125+
//
126+
// This function and sinh should stay in sync; if you make a
127+
// modification here, you should almost surely make a parallel
128+
// modification to sinh below.
78129
guard q.isFinite else { return q }
79-
let argument = q.isReal ? .zero : q.imaginary.length
80-
let axis = q.isReal ? .zero : (q.imaginary / argument)
81-
return cosh(q.real, argument, axis: axis)
130+
let θ = q.imaginary.length
131+
let axis = !θ.isZero ? (q.imaginary / θ) : .zero
132+
guard q.real.magnitude < -RealType.log(.ulpOfOne) else {
133+
let rotation = Quaternion(halfAngle: θ, unitAxis: axis)
134+
let firstScale = RealType.exp(q.real.magnitude/2)
135+
return rotation.multiplied(by: firstScale).multiplied(by: firstScale/2)
136+
}
137+
return Quaternion(
138+
real: .cosh(q.real) * .cos(θ),
139+
imaginary: axis * .sinh(q.real) * .sin(θ)
140+
)
82141
}
83142

84-
// sinh(r + xi + yj + zk) = sinh(r + v)
85-
// = sinh(r) cos(||v||) + (v/||v||) cosh(r) sin(||v||)
86-
//
87-
// See sinh on complex numbers for algorithm details.
88143
@inlinable
89144
public static func sinh(_ q: Quaternion) -> Quaternion {
145+
// Mathematically, this operation can be expanded in terms of
146+
// trigonometric `Real` operations as follows (`let θ = ||v||`):
147+
//
148+
// ```
149+
// sinh(q) = (exp(q) - exp(-q)) / 2
150+
// = sinh(r) cos(θ) + (v/θ) cosh(r) sin(θ)
151+
// ```
90152
guard q.isFinite else { return q }
91-
let argument = q.isReal ? .zero : q.imaginary.length
92-
let axis = q.isReal ? .zero : (q.imaginary / argument)
153+
let θ = q.imaginary.length
154+
let axis = !θ.isZero ? (q.imaginary / θ) : .zero
93155
guard q.real.magnitude < -RealType.log(.ulpOfOne) else {
94-
let rotation = Quaternion(halfAngle: argument, unitAxis: axis)
156+
let rotation = Quaternion(halfAngle: θ, unitAxis: axis)
95157
let firstScale = RealType.exp(q.real.magnitude/2)
96158
let secondScale = RealType(signOf: q.real, magnitudeOf: firstScale/2)
97159
return rotation.multiplied(by: firstScale).multiplied(by: secondScale)
98160
}
99161
return Quaternion(
100-
real: .sinh(q.real) * .cos(argument),
101-
imaginary: axis * .cosh(q.real) * .sin(argument)
162+
real: .sinh(q.real) * .cos(θ),
163+
imaginary: axis * .cosh(q.real) * .sin(θ)
102164
)
103165
}
104166

105-
// tanh(q) = sinh(q) / cosh(q)
106-
//
107-
// See tanh on complex numbers for algorithm details.
108167
@inlinable
109168
public static func tanh(_ q: Quaternion) -> Quaternion {
169+
// Mathematically, this operation can be expanded in terms of
170+
// trigonometric `Real` operations as follows (`let θ = ||v||`):
171+
//
172+
// ```
173+
// tanh(q) = sinh(q) / cosh(q)
174+
// ```
110175
guard q.isFinite else { return q }
111176
// Note that when |r| is larger than -log(.ulpOfOne),
112177
// sinh(r + v) == ±cosh(r + v), so tanh(r + v) is just ±1.
@@ -122,36 +187,75 @@ extension Quaternion/*: ElementaryFunctions */ {
122187
return sinh(q) / cosh(q)
123188
}
124189

125-
// cos(r + xi + yj + zk) = cos(r + v)
126-
// = cos(r) cosh(||v||) - (v/||v||) sin(r) sinh(||v||)
127-
//
128-
// See cosh for algorithm details.
129190
@inlinable
130191
public static func cos(_ q: Quaternion) -> Quaternion {
192+
// Mathematically, this operation can be expanded in terms of
193+
// trigonometric `Real` operations as follows (`let θ = ||v||`):
194+
//
195+
// ```
196+
// cos(r + v) = (exp(q * (v/θ)) + exp(-q * (v/θ))) / 2
197+
// = cos(r) cosh(θ) - (v/θ) sin(r) sinh(θ)
198+
// ```
131199
guard q.isFinite else { return q }
132-
let argument = q.isReal ? .zero : q.imaginary.length
133-
let axis = q.isReal ? .zero : (q.imaginary / argument)
134-
return cosh(-argument, q.real, axis: axis)
200+
let θ = q.imaginary.length
201+
let axis = !θ.isZero ? (q.imaginary / θ) : .zero
202+
guard θ.magnitude < -RealType.log(.ulpOfOne) else {
203+
let rotation = Quaternion(halfAngle: q.real, unitAxis: axis)
204+
let firstScale = RealType.exp(θ.magnitude/2)
205+
let secondScale = firstScale/2
206+
return rotation.multiplied(by: firstScale).multiplied(by: secondScale)
207+
}
208+
return Quaternion(
209+
real: .cosh(θ) * .cos(q.real),
210+
imaginary: -axis * .sinh(θ) * .sin(q.real)
211+
)
135212
}
136213

137-
// sin(r + xi + yj + zk) = sin(r + v)
138-
// = sin(r) cosh(-||v||) - (v/||v||) cos(r) sinh(-||v||)
139-
//
140-
// See sinh for algorithm details.
141214
@inlinable
142215
public static func sin(_ q: Quaternion) -> Quaternion {
216+
// Mathematically, this operation can be expanded in terms of
217+
// trigonometric `Real` operations as follows (`let θ = ||v||`):
218+
//
219+
// ```
220+
// sin(r + v) = -((exp(q * (v/θ)) - exp(-q * (v/θ))) (v/θ * 2)
221+
// = sin(r) cosh(θ) + (v/θ) cos(r) sinh(θ)
222+
// ```
143223
guard q.isFinite else { return q }
144-
let argument = q.isReal ? .zero : q.imaginary.length
145-
let axis = q.isReal ? .zero : (q.imaginary / argument)
146-
let (x, y) = sinh(-argument, q.real)
147-
return Quaternion(real: y, imaginary: axis * -x)
224+
let θ = q.imaginary.length
225+
let axis = !θ.isZero ? (q.imaginary / θ) : .zero
226+
guard θ.magnitude < -RealType.log(.ulpOfOne) else {
227+
let rotation = Quaternion(halfAngle: q.real, unitAxis: axis)
228+
let firstScale = RealType.exp(θ.magnitude/2)
229+
let secondScale = RealType(signOf: θ, magnitudeOf: firstScale/2)
230+
return rotation.multiplied(by: firstScale).multiplied(by: secondScale)
231+
}
232+
return Quaternion(
233+
real: .cosh(θ) * .sin(q.real),
234+
imaginary: axis * .sinh(θ) * .cos(q.real)
235+
)
148236
}
149237

150-
// tan(q) = sin(q) / cos(q)
151-
//
152-
// See tanh for algorithm details.
153238
@inlinable
154239
public static func tan(_ q: Quaternion) -> Quaternion {
240+
// Mathematically, this operation can be expanded in terms of
241+
// trigonometric `Real` operations as follows (`let θ = ||v||`):
242+
//
243+
// ```
244+
// tan(q) = sin(q) / cos(q)
245+
// ```
246+
guard q.isFinite else { return q }
247+
let θ = q.imaginary.length
248+
// Note that when |θ| is larger than -log(.ulpOfOne),
249+
// sin(r + v) == ±cos(r + v), so tan(r + v) is just ±1.
250+
guard θ.magnitude < -RealType.log(.ulpOfOne) else {
251+
return Quaternion(
252+
real: RealType(signOf: q.real, magnitudeOf: 1),
253+
imaginary:
254+
RealType(signOf: q.imaginary.x, magnitudeOf: 0),
255+
RealType(signOf: q.imaginary.y, magnitudeOf: 0),
256+
RealType(signOf: q.imaginary.z, magnitudeOf: 0)
257+
) * Quaternion(RealType(signOf: q.real, magnitudeOf: 1))
258+
}
155259
return sin(q) / cos(q)
156260
}
157261

@@ -162,14 +266,14 @@ extension Quaternion/*: ElementaryFunctions */ {
162266
// the single exceptional value.
163267
guard q.isFinite && !q.isZero else { return .infinity }
164268

165-
let vectorLength = q.imaginary.length
166-
let scale = q.halfAngle / vectorLength
269+
let argument = q.imaginary.length
270+
let axis = q.imaginary / argument
167271

168272
// We deliberatly choose log(length) over the (faster)
169273
// log(lengthSquared) / 2 which is used for complex numbers; as
170274
// the squared length of quaternions is more prone to overflows than the
171275
// squared length of complex numbers.
172-
return Quaternion(real: .log(q.length), imaginary: q.imaginary * scale)
276+
return Quaternion(real: .log(q.length), imaginary: axis * q.halfAngle)
173277
}
174278

175279
// MARK: - pow-like functions
@@ -206,42 +310,3 @@ extension Quaternion/*: ElementaryFunctions */ {
206310
return exp(log(q).divided(by: RealType(n)))
207311
}
208312
}
209-
210-
// MARK: - Hyperbolic trigonometric function helper
211-
extension Quaternion {
212-
213-
// See cosh of complex numbers for algorithm details.
214-
@usableFromInline @_transparent
215-
internal static func cosh(
216-
_ x: RealType,
217-
_ y: RealType,
218-
axis: SIMD3<RealType>
219-
) -> Quaternion {
220-
guard x.magnitude < -RealType.log(.ulpOfOne) else {
221-
let rotation = Quaternion(halfAngle: y, unitAxis: axis)
222-
let firstScale = RealType.exp(x.magnitude/2)
223-
let secondScale = firstScale/2
224-
return rotation.multiplied(by: firstScale).multiplied(by: secondScale)
225-
}
226-
return Quaternion(
227-
real: .cosh(x) * .cos(y),
228-
imaginary: axis * .sinh(x) * .sin(y)
229-
)
230-
}
231-
232-
// See sinh of complex numbers for algorithm details.
233-
@usableFromInline @_transparent
234-
internal static func sinh(
235-
_ x: RealType,
236-
_ y: RealType
237-
) -> (RealType, RealType) {
238-
guard x.magnitude < -RealType.log(.ulpOfOne) else {
239-
var (x, y) = (RealType.cos(y), RealType.sin(y))
240-
let firstScale = RealType.exp(x.magnitude/2)
241-
(x, y) = (x * firstScale, y * firstScale)
242-
let secondScale = RealType(signOf: x, magnitudeOf: firstScale/2)
243-
return (x * secondScale, y * secondScale)
244-
}
245-
return (.sinh(x) * .cos(y), .cosh(x) * .sin(y))
246-
}
247-
}

0 commit comments

Comments
 (0)