@@ -224,78 +224,131 @@ extension Complex /*: ElementaryFunctions */ {
224224 // If z is zero or infinite, the phase is undefined, so the result is
225225 // the single exceptional value.
226226 guard z. isFinite && !z. isZero else { return . infinity }
227- // Otherwise, try computing lengthSquared; if the result is normal,
228- // we can just take its log to get the real part of the result.
229- let r2 = z . lengthSquared
227+ // Having eliminated non-finite values and zero, the imaginary part is
228+ // easy; it's just the phase, which is always computed with good
229+ // relative accuracy via atan2.
230230 let θ = z. phase
231- if r2. isNormal { return Complex ( . log( r2) / 2 , θ) }
232- // z is finite, but z.lengthSquared is not normal. Rescale and recompute.
233- let w = z. divided ( by: z. magnitude)
234- return Complex ( . log( z. magnitude) + . log( w. lengthSquared) / 2 , θ)
235- }
236-
237- @inlinable
238- public static func log( onePlus z: Complex ) -> Complex {
239- // Nevin proposed the idea for this implementation on the Swift forums:
240- // https://forums.swift.org/t/elementaryfunctions-compliance-for-complex/37903/3
231+ // The real part of the result is trickier. In exact arithmetic, the
232+ // real part is just log |z|--many implementations of complex functions
233+ // simply use this expression as is. However, there are two problems
234+ // lurking here:
241235 //
242- // Here's a quick explainer on why it works: in exact arithmetic,
236+ // - length can overflow even when log(z) is finite.
243237 //
244- // log(1+z) = (log |1+z|, atan2(y, 1+x))
238+ // - when length is close to 1, catastrophic cancellation is hidden
239+ // in this expression. Consider, e.g. z = 1 + δi for small δ.
245240 //
246- // where x and y are the real and imaginary parts of z, respectively.
241+ // Because δ ≪ 1, |z| rounds to 1, and so log |z| produces zero.
242+ // We can expand using Taylor series to see that the result should
243+ // be:
247244 //
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).
245+ // log |z| = log √(1 + δ²)
246+ // = log(1 + δ²)/2
247+ // = δ²/2 + O(δ⁴)
253248 //
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:
249+ // So naively using log |z| results in a total loss of relative
250+ // accuracy for this case. Note that this is _not_ constrained near
251+ // a single point; it occurs everywhere close to the circle |z| = 1.
256252 //
257- // Re(log 1+z) = (log 1+z + Conj(log 1+z))/2
253+ // Note that this case still _does_ deliver a result with acceptable
254+ // relative accuracy in the complex norm, because
258255 //
259- // Log commutes with conjugation, so this becomes:
256+ // Im(log z) ≈ δ ≫ δ²/2 ≈ Re(log z).
260257 //
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
258+ // There are a number of ways to try to tackle this problem. I'll begin
259+ // with a simple one that solves the first issue, and _sometimes_ the
260+ // second, then analyze when it doesn't work for the second case.
264261 //
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.
262+ // To handle very large arguments without overflow, the standard
263+ // approach is to _rescale_ the problem. We can do this by finding
264+ // whichever of x and y has greater magnitude, and dividing through
265+ // by it. You can think of this as changing coordinates by reflections
266+ // so that we get a new value w = u + iv with |w| = |z| (and hence
267+ // Re(log w) = Re(log z), and 0 ≤ u, 0 ≤ v ≤ u.
268+ let u = max ( z. x. magnitude, z. y. magnitude)
269+ let v = min ( z. x. magnitude, z. y. magnitude)
270+ // Now expand out log |w|:
271271 //
272- // The other common implementation choice for log1p is Kahan's trick:
272+ // log |w| = log(u² + v²)/2
273+ // = log u + log(onePlus: (v/u)²)/2
273274 //
274- // w := 1+z
275- // return z/(w-1) * log(w)
275+ // This looks promising! It handles overflow well, because log(u) is
276+ // finite for every finite u, and we have 0 ≤ v/u ≤ 1, so the second
277+ // term is bounded by 0 ≤ log(1 + (v/u)²)/2 ≤ (log 2)/2. It also
278+ // handles the example I gave above well: we have u = 1, v = δ, and
276279 //
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- )
280+ // log(1) + log(onePlus: δ²)/2 = 0 + δ²/2
281+ //
282+ // as expected.
283+ //
284+ // Unfortunately, it does not handle all points close to the unit
285+ // circle so well; it's easy to see why if we look at the two terms
286+ // that contribute to the result. Cancellation occurs when the result
287+ // is close to zero and the terms have opposing signs. By construction,
288+ // the second term is always positive, so the easiest observation is
289+ // that cancellation is only a problem for u < 1 (because otherwise
290+ // log u is also positive, and there can be no cancellation).
291+ //
292+ // We are not trying for sub-ulp accuracy, just a good relative error
293+ // bound, so for our purposes it suffices to have log u dominate the
294+ // result:
295+ if u >= 1 || u >= RealType . _mulAdd ( u, u, v*v) {
296+ let r = v / u
297+ return Complex ( . log( u) + . log( onePlus: r*r) / 2 , θ)
298+ }
299+ // Here we're in the tricky case; cancellation is likely to occur.
300+ // Instead of the factorization used above, we will want to evaluate
301+ // log(onePlus: u² + v² - 1)/2. This all boils down to accurately
302+ // evaluating u² + v² - 1. To begin, calculate both squared terms
303+ // as exact head-tail products (u is guaranteed to be well scaled,
304+ // v may underflow, but if it does it doesn't matter, the u term is
305+ // all we need).
306+ let ( a, b) = Augmented . twoProdFMA ( u, u)
307+ let ( c, d) = Augmented . twoProdFMA ( v, v)
308+ // It would be nice if we could simply use a - 1, but unfortunately
309+ // we don't have a tight enough bound to guarantee that that expression
310+ // is exact; a may be as small as 1/4, so we could lose a single bit
311+ // to rounding if we did that.
312+ var ( s, e) = Augmented . fastTwoSum ( - 1 , a)
313+ // Now we are ready to assemble the result. If cancellation happens,
314+ // then |c| > |e| > |b|, |d|, so this assembly order is safe. It's
315+ // also possible that |c| and |d| are small, but if that happens then
316+ // there is no significant cancellation, and the exact assembly doesn't
317+ // matter.
318+ s = ( s + c) + e + b + d
319+ return Complex ( . log( onePlus: s) / 2 , θ)
320+ }
321+
322+ @inlinable
323+ public static func log( onePlus z: Complex ) -> Complex {
324+ // If either |x| or |y| is bounded away from the origin, we don't need
325+ // any extra precision, and can just literally compute log(1+z). Note
326+ // that this includes part of the circle |1+z| = 1 where log(onePlus:)
327+ // vanishes (where x <= -0.5), but on this portion of the circle 1+x
328+ // is always exact by Sterbenz' lemma, so as long as log( ) produces
329+ // a good result, log(1+z) will too.
330+ guard 2 * z. x. magnitude < 1 && z. y. magnitude < 1 else { return log ( 1 + z) }
331+ // z is in (±0.5, ±1), so we need to evaluate more carefully.
332+ // The imaginary part is straightforward:
333+ let θ = z. phase
334+ // For the real part, we _could_ use the same approach that we do for
335+ // log( ), but we'd need an extra-precise (1+x)², which can potentially
336+ // be quite painful to calculate. Instead, we can use an approach that
337+ // NevinBR suggested on the Swift forums:
338+ //
339+ // Re(log 1+z) = (log 1+z + log 1+z̅)/2
340+ // = log((1+z)(1+z̅)/2
341+ // = log(1+z+z̅+zz̅)/2
342+ // = log((2+x)x + y²)/2
343+ //
344+ // So now we need to evaluate (2+x)x + y² accurately. To do this,
345+ // we employ augmented arithmetic; the key observation here is that
346+ // cancellation is only a problem when y² ≈ -(2+x)x
347+ let xp2 = Augmented . fastTwoSum ( 2 , z. x) // Known that 2 > |x|.
348+ let a = Augmented . twoProdFMA ( z. x, xp2. head)
349+ let y ² = Augmented . twoProdFMA ( z. y, z. y)
350+ let s = ( a. head + y ². head + a . tail + y ². tail) . addingProduct ( z. x, xp2. tail)
351+ return Complex ( . log( onePlus: s) / 2 , θ)
299352 }
300353
301354 public static func acos( _ z: Complex ) -> Complex {
0 commit comments