Skip to content

Commit 2128216

Browse files
stephentyronestephentyrone
authored
Merge pull request #152 from stephentyrone/complex-trig
First pass over hyperbolics and trig functions for Complex.
2 parents f49d1af + 0f02c65 commit 2128216

File tree

2 files changed

+241
-75
lines changed

2 files changed

+241
-75
lines changed

‎Sources/ComplexModule/ElementaryFunctions.swift

Lines changed: 206 additions & 75 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,10 @@
2222
// to do that for all of these functions off the top of my head, and
2323
// I don't think that other libraries have tried to do so in general,
2424
// so this is a research project. We should not sacrifice 1-4 for it.
25+
// Note that multiplication and division don't even provide good
26+
// componentwise relative accuracy, so it's _totally OK_ to not get
27+
// it for these functions too. But: it's a dynamite long-term research
28+
// project.
2529
// 6. Give the best performance we can. We should care about performance,
2630
// but with lower precedence than the other considerations.
2731

@@ -31,104 +35,170 @@ import RealModule
3135
extension Complex /*: ElementaryFunctions */ {
3236

3337
// MARK: - exp-like functions
34-
/// Checks if x is bounded away overflowing exp(x).
35-
///
36-
/// This is a conservative (imprecise) check; if it returns `true`, `exp(x)` is definitely safe, but
37-
/// it will return `false` even in some cases where `exp(x)` would not overflow.
38-
@usableFromInline @inline(__always)
39-
internal static func expIsSafe(_ x: RealType) -> Bool {
40-
// If x < log(greatestFiniteMagnitude), then exp(x) does not overflow.
41-
// To protect ourselves against sketchy log or exp implementations in
42-
// an unknown host library, we round down to the nearest integer to get
43-
// some margin of safety.
44-
return x < RealType.log(.greatestFiniteMagnitude).rounded(.down)
45-
}
4638

47-
/// Computes exp(z) with extra care near the overflow boundary.
39+
/// The complex exponential function e^z whose base `e` is the base of the natural logarithm.
4840
///
49-
/// When x = z.real is large, exp(x) may overflow even when exp(z) is finite,
50-
/// because exp(z) = exp(x) * (cos(y) + i sin(y)), and max(cos(y),sin(y)) may
51-
/// be as small as 1/sqrt(2).
52-
///
53-
/// - Parameter z: a complex number with large real part.
54-
@usableFromInline
55-
internal static func expNearOverflow(_ z: Complex) -> Complex {
56-
let xm1 = z.x - 1
57-
let y = z.y
58-
let r = Complex(.cos(y), .sin(y)).multiplied(by: .exp(1))
59-
return r.multiplied(by: .exp(xm1))
60-
}
61-
62-
// exp(x + iy) = exp(x)(cos(y) + i sin(y))
41+
/// Mathematically, this operation can be expanded in terms of the `Real` operations `exp`,
42+
/// `cos` and `sin` as follows:
43+
/// ```
44+
/// exp(x + iy) = exp(x) exp(iy)
45+
/// = exp(x) cos(y) + i exp(x) sin(y)
46+
/// ```
47+
/// Note that naive evaluation of this expression in floating-point would be prone to premature
48+
/// overflow, since `cos` and `sin` both have magnitude less than 1 for most inputs (i.e.
49+
/// `exp(x)` may be infinity when `exp(x) cos(y)` would not be.
6350
@inlinable
6451
public static func exp(_ z: Complex) -> Complex {
65-
// Naively we would let exp(-∞,0) fall out as 0, matching the real
66-
// type behavior, but that breaks the single-point-at-infinity
67-
// semantics, because exp has an essential singularity at infinity.
68-
guard z.x != -.infinity else { return .infinity }
69-
guard expIsSafe(z.x) else { return expNearOverflow(z) }
52+
guard z.isFinite else { return z }
53+
// If x < log(greatestFiniteMagnitude), then exp(x) does not overflow.
54+
// To protect ourselves against sketchy log or exp implementations in
55+
// an unknown host library, or slight rounding disagreements between
56+
// the two, subtract one from the bound for a little safety margin.
57+
guard z.x < RealType.log(.greatestFiniteMagnitude) - 1 else {
58+
let halfScale = RealType.exp(z.x/2)
59+
let phase = Complex(RealType.cos(z.y), RealType.sin(z.y))
60+
return phase.multiplied(by: halfScale).multiplied(by: halfScale)
61+
}
7062
return Complex(.cos(z.y), .sin(z.y)).multiplied(by: .exp(z.x))
7163
}
7264

73-
// exp(x + iy) - 1 = (exp(x) cos(y) - 1) + i exp(x) sin(y)
74-
// -------- u --------
75-
// Note that the imaginary part is just the usual exp(x) sin(y);
76-
// the only trick is computing the real part ("u"):
77-
//
78-
// u = exp(x) cos(y) - 1
79-
// = exp(x) cos(y) - cos(y) + cos(y) - 1
80-
// = (exp(x) - 1) cos(y) + (cos(y) - 1)
81-
// = expMinusOne(x) cos(y) + cosMinusOne(y)
82-
//
83-
// Note: most implementations of expm1 for complex (e.g. Julia's)
84-
// factor the real part as follows instead:
85-
//
86-
// exp(x) cosMinuxOne(y) + expMinusOne(y)
87-
//
88-
// This expression gives good accuracy close to zero, but suffers from
89-
// catastrophic cancellation when z.x is large and z.y is near an odd
90-
// multiple of π/2. This is _OK_ (the componentwise error is bad, but
91-
// the error in a complex norm is acceptable), but we can do better by
92-
// factoring on cosine instead of exp.
93-
//
94-
// The other implementation that is sometimes seen, 2*exp(z/2)*sinh(z/2),
95-
// has the same weaknesses.
96-
//
97-
// The approach used here achieves good componentwise worst-case error
98-
// (7e-5 for Float) as well as normwise error (2.9e-7) in structured
99-
// and randomized tests. The alternative factorization achieves
100-
// comparable normwise error (3.9e-7), but dramatically worse
101-
// componentwise errors, e.g. Complex(18, -3π/2) produces (4.0, 6.57e7)
102-
// while the reference result would be (-0.22, 6.57e7).
10365
@inlinable
10466
public static func expMinusOne(_ z: Complex) -> Complex {
105-
// Naively we would let exp(-∞,0) fall out as 0, matching the real
106-
// type behavior, but that breaks the single-point-at-infinity
107-
// semantics, because exp has an essential singularity at infinity.
108-
guard z.x != -.infinity else { return .infinity }
67+
// exp(x + iy) - 1 = (exp(x) cos(y) - 1) + i exp(x) sin(y)
68+
// -------- u --------
69+
// Note that the imaginary part is just the usual exp(x) sin(y);
70+
// the only trick is computing the real part ("u"):
71+
//
72+
// u = exp(x) cos(y) - 1
73+
// = exp(x) cos(y) - cos(y) + cos(y) - 1
74+
// = (exp(x) - 1) cos(y) + (cos(y) - 1)
75+
// = expMinusOne(x) cos(y) + cosMinusOne(y)
76+
//
77+
// Note: most implementations of expm1 for complex (e.g. Julia's)
78+
// factor the real part as follows instead:
79+
//
80+
// exp(x) cosMinuxOne(y) + expMinusOne(x)
81+
//
82+
// The other implementation that is sometimes seen is:
83+
//
84+
// expMinusOne(z) = 2*exp(z/2)*sinh(z/2)
85+
//
86+
// All three of these implementations provide good relative error
87+
// bounds _in the complex norm_, but the cosineMinusOne-based
88+
// implementation has the best _componentwise_ error characteristics,
89+
// which is why we use it here:
90+
//
91+
// Implementation | Real | Imaginary |
92+
// ---------------+--------------------+----------------+
93+
// Ours | Hybrid bound | Relative bound |
94+
// Standard | No bound | Relative bound |
95+
// Half Angle | Hybrid bound | Hybrid bound |
96+
//
97+
// FUTURE WORK: devise an algorithm that achieves good _relative_ error
98+
// in the real component as well. Doing this efficiently is a research
99+
// project--exp(x) cos(y) - 1 can be very nearly zero along a curve in
100+
// the complex plane, not only at zero. Evaluating it accurately
101+
// _without_ depending on arbitrary-precision exp and cos is an
102+
// interesting challenge.
103+
guard z.isFinite else { return z }
109104
// If exp(z) is close to the overflow boundary, we don't need to
110-
// worry about the m1 part; we're just computing exp(z). (Even when
111-
// z.y is near a multiple of π/2, it can't be close enough to
112-
// overcome the scaling from exp(z.x), so the -1 term is _always_
113-
// negligable).
114-
guard expIsSafe(z.x) else { return expNearOverflow(z) }
105+
// worry about the "MinusOne" part of this function; we're just
106+
// computing exp(z). (Even when z.y is near a multiple of π/2,
107+
// it can't be close enough to overcome the scaling from exp(z.x),
108+
// so the -1 term is _always_ negligable). So we simply handle
109+
// these cases exactly the same as exp(z).
110+
guard z.x < RealType.log(.greatestFiniteMagnitude) - 1 else {
111+
let halfScale = RealType.exp(z.x/2)
112+
let phase = Complex(RealType.cos(z.y), RealType.sin(z.y))
113+
return phase.multiplied(by: halfScale).multiplied(by: halfScale)
114+
}
115115
// Special cases out of the way, evaluate as discussed above.
116116
return Complex(
117117
RealType._mulAdd(.cos(z.y), .expMinusOne(z.x), .cosMinusOne(z.y)),
118118
.exp(z.x) * .sin(z.y)
119119
)
120120
}
121121

122+
// cosh(x + iy) = cosh(x) cos(y) + i sinh(x) sin(y).
123+
//
124+
// Like exp, cosh is entire, so we do not need to worry about where
125+
// branch cuts fall. Also like exp, cancellation never occurs in the
126+
// evaluation of the naive expression, so all we need to be careful
127+
// about is the behavior near the overflow boundary.
128+
//
129+
// Fortunately, if |x| >= -log(ulpOfOne), cosh(x) and sinh(x) are
130+
// both just exp(|x|)/2, and we already know how to compute that.
131+
//
132+
// This function and sinh should stay in sync; if you make a
133+
// modification here, you should almost surely make a parallel
134+
// modification to sinh below.
135+
@inlinable @inline(__always)
122136
public static func cosh(_ z: Complex) -> Complex {
123-
fatalError()
137+
guard z.isFinite else { return z }
138+
guard z.x.magnitude < -RealType.log(.ulpOfOne) else {
139+
let phase = Complex(RealType.cos(z.y), RealType.sin(z.y))
140+
let firstScale = RealType.exp(z.x.magnitude/2)
141+
let secondScale = firstScale/2
142+
return phase.multiplied(by: firstScale).multiplied(by: secondScale)
143+
}
144+
// Future optimization opportunity: expm1 is faster than cosh/sinh
145+
// on most platforms, and division is now commonly pipelined, so we
146+
// might replace the check above with a much more conservative one,
147+
// and then evaluate cosh(x) and sinh(x) as
148+
//
149+
// cosh(x) = 1 + 0.5*expm1(x)*expm1(x) / (1 + expm1(x))
150+
// sinh(x) = expm1(x) + 0.5*expm1(x) / (1 + expm1(x))
151+
//
152+
// This won't be a _big_ win except on platforms with a crappy sinh
153+
// and cosh, and for those we should probably just provide our own
154+
// implementations of _those_, so for now let's keep it simple and
155+
// obviously correct.
156+
return Complex(
157+
RealType.cosh(z.x) * RealType.cos(z.y),
158+
RealType.sinh(z.x) * RealType.sin(z.y)
159+
)
124160
}
125161

162+
// sinh(x + iy) = sinh(x) cos(y) + i cosh(x) sinh(y)
163+
//
164+
// See cosh above for algorithm details.
165+
@inlinable @inline(__always)
126166
public static func sinh(_ z: Complex) -> Complex {
127-
fatalError()
167+
guard z.isFinite else { return z }
168+
guard z.x.magnitude < -RealType.log(.ulpOfOne) else {
169+
let phase = Complex(RealType.cos(z.y), RealType.sin(z.y))
170+
let firstScale = RealType.exp(z.x.magnitude/2)
171+
let secondScale = RealType(signOf: z.x, magnitudeOf: firstScale/2)
172+
return phase.multiplied(by: firstScale).multiplied(by: secondScale)
173+
}
174+
return Complex(
175+
RealType.sinh(z.x) * RealType.cos(z.y),
176+
RealType.cosh(z.x) * RealType.sin(z.y)
177+
)
128178
}
129179

180+
// tanh(z) = sinh(z) / cosh(z)
181+
@inlinable
130182
public static func tanh(_ z: Complex) -> Complex {
131-
fatalError()
183+
guard z.isFinite else { return z }
184+
// Note that when |x| is larger than -log(.ulpOfOne),
185+
// sinh(x + iy) == ±cosh(x + iy), so tanh(x + iy) is just ±1.
186+
guard z.x.magnitude < -RealType.log(.ulpOfOne) else {
187+
return Complex(
188+
RealType(signOf: z.x, magnitudeOf: 1),
189+
RealType(signOf: z.y, magnitudeOf: 0)
190+
)
191+
}
192+
// Now we have z in a vertical strip where exp(x) is reasonable,
193+
// and y is finite, so we can simply evaluate sinh(z) and cosh(z).
194+
//
195+
// TODO: Kahan uses a different expression for evaluation here; it
196+
// isn't strictly necessary for numerics reasons--it's to avoid
197+
// doing the complex division, but it probably provides better
198+
// componentwise error bounds, and is likely more efficient (because
199+
// it avoids the complex division, which is painful even when well-
200+
// scaled). This suffices to get us up and running.
201+
return sinh(z) / cosh(z)
132202
}
133203

134204
// cos(z) = cosh(iz)
@@ -149,6 +219,7 @@ extension Complex /*: ElementaryFunctions */ {
149219
}
150220

151221
// MARK: - log-like functions
222+
@inlinable
152223
public static func log(_ z: Complex) -> Complex {
153224
// If z is zero or infinite, the phase is undefined, so the result is
154225
// the single exceptional value.
@@ -163,8 +234,68 @@ extension Complex /*: ElementaryFunctions */ {
163234
return Complex(.log(z.magnitude) + .log(w.lengthSquared)/2, θ)
164235
}
165236

237+
@inlinable
166238
public static func log(onePlus z: Complex) -> Complex {
167-
fatalError()
239+
// Nevin proposed the idea for this implementation on the Swift forums:
240+
// https://forums.swift.org/t/elementaryfunctions-compliance-for-complex/37903/3
241+
//
242+
// Here's a quick explainer on why it works: in exact arithmetic,
243+
//
244+
// log(1+z) = (log |1+z|, atan2(y, 1+x))
245+
//
246+
// where x and y are the real and imaginary parts of z, respectively.
247+
//
248+
// The first thing to note is that the expression for the imaginary
249+
// part works fine as is. If cancellation occurs (because x ≈ -1),
250+
// then 1+x is exact, and so we have good componentwise relative
251+
// accuracy. Otherwise, x is bounded away from -1 and 1+x has good
252+
// relative accuracy, and therefore so does atan2(y, 1+x).
253+
//
254+
// So the real part is the hard part (no surprise, just like expPlusOne).
255+
// Nevin's clever idea is simply to take advantage of the expansion:
256+
//
257+
// Re(log 1+z) = (log 1+z + Conj(log 1+z))/2
258+
//
259+
// Log commutes with conjugation, so this becomes:
260+
//
261+
// Re(log 1+z) = (log 1+z + log 1+z̅)/2
262+
// = log((1+z)(1+z̅)/2
263+
// = log(1+z+z̅+zz̅)/2
264+
//
265+
// This behaves well close to zero, because the z+z̅ term dominates
266+
// and is computed exactly. Away from zero, cancellation occurs near
267+
// the circle x(x+2) + y^2 = 0, but everywhere along this curve we
268+
// have |Im(log 1+z)| >= π/2, so the relative error in the complex
269+
// norm is well-controlled. We can take advantage of FMA to further
270+
// reduce the cancellation error and recover a good error bound.
271+
//
272+
// The other common implementation choice for log1p is Kahan's trick:
273+
//
274+
// w := 1+z
275+
// return z/(w-1) * log(w)
276+
//
277+
// But this actually doesn't do as well as Nevin's approach does,
278+
// and requires a complex division, which we want to avoid when we
279+
// can do so.
280+
var a = 2*z.x
281+
// We want to add the larger term first (contra usual guidance for
282+
// floating-point error optimization), because we're optimizing for
283+
// the catastrophic cancellation case; when that happens adding the
284+
// larger term via FMA is always exact. When cancellation doesn't
285+
// happen, the simple relative error bound carries through the
286+
// rest of the computation.
287+
let large = max(z.x.magnitude, z.y.magnitude)
288+
let small = min(z.x.magnitude, z.y.magnitude)
289+
a.addProduct(large, large)
290+
a.addProduct(small, small)
291+
// If r2 overflowed, then |z| ≫ 1, and so log(1+z) = log(z).
292+
guard a.isFinite else { return log(z) }
293+
// Unlike log(z), we do not need to worry about what happens if a
294+
// underflows.
295+
return Complex(
296+
RealType.log(onePlus: a)/2,
297+
RealType.atan2(y: z.y, x: 1+z.x)
298+
)
168299
}
169300

170301
public static func acos(_ z: Complex) -> Complex {

‎Tests/ComplexTests/ElementaryFunctionTests.swift

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -112,20 +112,55 @@ final class ElementaryFunctionTests: XCTestCase {
112112
}
113113
}
114114

115+
func testCosh<T: Real & FixedWidthFloatingPoint>(_ type: T.Type) {
116+
// cosh(0) = 1
117+
XCTAssertEqual(1, Complex<T>.cosh(Complex( 0, 0)))
118+
XCTAssertEqual(1, Complex<T>.cosh(Complex(-0, 0)))
119+
XCTAssertEqual(1, Complex<T>.cosh(Complex(-0,-0)))
120+
XCTAssertEqual(1, Complex<T>.cosh(Complex( 0,-0)))
121+
// cosh is the identity at infinity.
122+
XCTAssertFalse(Complex<T>.cosh(Complex( .infinity, 0)).isFinite)
123+
XCTAssertFalse(Complex<T>.cosh(Complex( .infinity, .infinity)).isFinite)
124+
XCTAssertFalse(Complex<T>.cosh(Complex( 0, .infinity)).isFinite)
125+
XCTAssertFalse(Complex<T>.cosh(Complex(-.infinity, .infinity)).isFinite)
126+
XCTAssertFalse(Complex<T>.cosh(Complex(-.infinity, 0)).isFinite)
127+
XCTAssertFalse(Complex<T>.cosh(Complex(-.infinity,-.infinity)).isFinite)
128+
XCTAssertFalse(Complex<T>.cosh(Complex( 0,-.infinity)).isFinite)
129+
XCTAssertFalse(Complex<T>.cosh(Complex( .infinity,-.infinity)).isFinite)
130+
XCTAssertFalse(Complex<T>.cosh(Complex( .nan, .nan)).isFinite)
131+
XCTAssertFalse(Complex<T>.cosh(Complex( .infinity, .nan)).isFinite)
132+
XCTAssertFalse(Complex<T>.cosh(Complex( .nan, .infinity)).isFinite)
133+
XCTAssertFalse(Complex<T>.cosh(Complex(-.infinity, .nan)).isFinite)
134+
XCTAssertFalse(Complex<T>.cosh(Complex( .nan,-.infinity)).isFinite)
135+
// Near-overflow test, same as exp() above, but it happens later, because
136+
// for large x, cosh(x + iy) ~ exp(x + iy)/2.
137+
let x = T.log(.greatestFiniteMagnitude) + T.log(18/8)
138+
let mag = T.greatestFiniteMagnitude/T.sqrt(2) * (9/8)
139+
var huge = Complex<T>.cosh(Complex(x, .pi/4))
140+
XCTAssert(huge.real.isApproximatelyEqual(to: mag))
141+
XCTAssert(huge.imaginary.isApproximatelyEqual(to: mag))
142+
huge = Complex<T>.cosh(Complex(-x, .pi/4))
143+
XCTAssert(huge.real.isApproximatelyEqual(to: mag))
144+
XCTAssert(huge.imaginary.isApproximatelyEqual(to: mag))
145+
}
146+
115147
func testFloat() {
116148
testExp(Float.self)
117149
testExpMinusOne(Float.self)
150+
testCosh(Float.self)
118151
}
119152

120153
func testDouble() {
121154
testExp(Double.self)
122155
testExpMinusOne(Double.self)
156+
testCosh(Double.self)
123157
}
124158

125159
#if (arch(i386) || arch(x86_64)) && !os(Windows) && !os(Android)
126160
func testFloat80() {
127161
testExp(Float80.self)
128162
testExpMinusOne(Float80.self)
163+
testCosh(Float80.self)
129164
}
130165
#endif
131166
}

0 commit comments

Comments
 (0)