Need to calculate the complex log of a double precision float-CodePudding

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.

CodePudding user response：

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

CodePudding user response：

... 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 thus retains precision in (a-1)*(a 1) b^2.

#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)
      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.