需要计算双精度浮点数的复数对数

Question

我尝试使用以下方法计算 clog(a + i*b) 的实部

将“x”视为复数。x = a + i*b 令 z 为 x 的复数对数。

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

这种方法在 ULP 方面会产生巨大的误差，特别是对于 0.5 和 1.0 之间的值。

我尝试了其他方法来避免实部和虚部的平方，例如

设 t = b / a；实数(x) = log(a) + 0.5 * log1p(t*t)

使用此方法时错误仍然存​​在。我知道错误可能来自 a 和 b 的平方，因此我尝试使用fma()操作来获取由于 'a' 和 'b' 的平方引起的错误

令 a2 = a * a b2 = b * b

err_a2 = fma(a,a, -a2)

err_b2 = fma(b,b,-b2)

然后我尝试0.5 * log(((err_a1 + err_b2) + a2) + b2)获取 x 的复数对数的实际值。

但结果仍然不准确。

如何才能log(sqrt(a^2 + b^2))准确计算（误差在 2 ULP 以内）。我认为我需要以更高的精度计算更高精度的平方根a^2 + b^2，但我不确定如何从这里开始。

Answer 1

sqrt(a^2 + b^2)只是std::hypot(a,b)。如果运气好的话，这已经很精确了。

Answer 2

...计算双精度数的复数对数...

代码可以使用double real_part = (double) clog(x).

要计算double 的复数对数的实部而不使用clog(x)close |x| == 1.0，请考虑使用log1p()^*1来形成更好的精度结果。

核心问题是|x| - 1.0精度会严重损失，这是确定的第一步log()。

0.5 * log(a^2 + b^2)在数学上就像0.5 * logp1(a^2 + b^2 - 1). 当|x|接近 1.0 且时|a| > |b|，使用0.5 * logp1((a-1)*(a+1) + b^2). |a|这会从精确值中减去 1.0并保留精度(a-1)*(a+1) + b^2。这与减去 1.0 不同，因为x的计算x已经失去了重要的进动。

#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));
}

输出

0.000976562 1
3.57628e-07
3.57628e-07

回复：“我尝试使用fma()...” --> 有些fma()质量低下。

^*1这些log1p函数计算 1 加上参数的以 e（自然）为底的对数。