Need to calculate the complex log of a double precision float

Question

I tried computing the real part of clog(a + i*b) using the following approach

Consider 'x' to be the complex number. x = a + i*b Let z be the complex log of x.

real(x) = 0.5 * log(a^2 + b^2)

This approach gives a huge error in terms of ULP for values between 0.5 and 1.0 especially.

I tried other approaches to avoid squaring of both the real and imaginary parts such as

Let t = b / a; real(x) = log(a) + 0.5 * log1p(t*t)

The error continued to persist with this approach as well. I understand that the error is likely from the squaring of a and b and hence I tried using fma() operations to get the error due to the squaring of 'a' and 'b'

Let a2 = a * a b2 = b * b

err_a2 = fma(a,a, -a2)

err_b2 = fma(b,b,-b2)

I then tried 0.5 * log(((err_a1 + err_b2) + a2) + b2) to get the real value of the complex log of x.

But the result is still inaccurate.

How can I compute log(sqrt(a^2 + b^2)) accurately (error within 2 ULP). I think I need to compute the square root of a^2 + b^2 in higher precision at a higher precision but I am not sure how to proceed from here.

Answer 1

sqrt(a^2 + b^2) is just std::hypot(a,b). With a bit of luck, that's already precise.

Answer 2

... calculate the complex log of a double ...

Code could use double real_part = (double) clog(x).

---

To calculate the real part of a complex log of a double without using clog(x) near |x| == 1.0, consider using log1p()^*1 to form a better precision result.

The core issue is |x| - 1.0 can suffer severe loss of precision and this is the first step in determining log().

0.5 * log(a^2 + b^2) is mathematically like 0.5 * logp1(a^2 + b^2 - 1). When |x| is near 1.0 and |a| > |b|, use 0.5 * logp1((a-1)*(a+1) + b^2). This subtracts the 1.0 from |a| exactly and retains precision with (a-1)*(a+1) + b^2. This differs from subtracting 1.0 from x as calculation of x has already lost important precession.

#include <complex.h>
#include <math.h>
#include <stdio.h>

#define root2 1.4142135623730950488016887242097

double clog_real(double a, double b) {
  double real_x;
  double h = hypot(a, b);
  // |x| near 1.0?
  if (h >= root2 / 2 && h < root2) {
    // Subtract 1 from the larger part
    if (fabs(a) > fabs(b)) {
      real_x = 0.5 * log1p((a - 1) * (a + 1) + b * b);
    } else {
      real_x = 0.5 * log1p((b - 1) * (b + 1) + a * a);
      // or (here and like-wise above IF you have a good fma())
      real_x = 0.5 * log1p(fma(a, a, (b - 1) * (b + 1)));
    }
  } else {
    real_x = log(h);
  }
  return real_x;
}

int main() {
  double a = 0x1.fffffe0000010p-12 * 2;
  double b = 0x1.fffffc0000040p-1;
  printf("%g %g\n", a, b);
  complex double c = a + csqrt(-1) * b;
  printf("%g\n", (double) clog(c));
  printf("%g\n", clog_real(a, b));
}

Output

0.000976562 1
3.57628e-07
3.57628e-07

---

Re: "I tried using fma() ..." --> Some fma() are low quality.

---

^*1 The log1p functions compute the base-e (natural) logarithm of 1 plus the argument.

Answer 3

Have you tried log(sqrt(a*a + b*b)) ? Normally, the square root has the opposite effect of the square, so probably if you try to calculate it you will make the log of a better suited number.

Anyway, to calculate logarithms of numbers close to 1.0, you can probably calculate the derivatives of log(z + 1) for z == 0 and you will get a better approach, as the function is analitic in the circle of radius < 0.5 and so, you will get a good taylor approximation. This approximation is written below (thanks to Wolfram alpha)

log(1+x) ~= x - x^2/2 + x^3/3 - x^4/4 + ... (-1)^(n)*x^(n+1)/(n+1) + O(x^(n+2))

This is a series that converges absolutely on an open circle of radius 1, so to calculate values close to 1 or the logarithm is well suited (indeed, it is what is used in many places).

Of course, if you want to solve this in the complex plane, you need to operate the calculation as complex numbers.