Intenté calcular la parte real de clog(a + i*b) usando el siguiente enfoque
Considere 'x' como el número complejo. x = a + i*b Sea z el logaritmo complejo de x.
real(x) = 0.5 * log(a^2 + b^2)
Este enfoque da un gran error en términos de ULP para valores entre 0,5 y 1,0 especialmente.
Probé otros enfoques para evitar el cuadrado de las partes real e imaginaria, como
Sea t = b/a; real(x) = log(a) + 0.5 * log1p(t*t)
El error continuó persistiendo con este enfoque también. Entiendo que es probable que el error se deba al cuadrado de a y b y, por lo tanto, intenté usar las operaciones fma() para obtener el error debido al cuadrado de 'a' y 'b'
Sea a2 = a * a b2 = b * b
err_a2 = fma(a,a, -a2)
err_b2 = fma(b,b,-b2)
Luego probé 0.5 * log(((err_a1 + err_b2) + a2) + b2) para obtener el valor real del logaritmo complejo de x.
Pero el resultado sigue siendo inexacto.
¿Cómo puedo calcular log(sqrt(a^2 + b^2)) con precisión (error dentro de 2 ULP). Creo que necesito calcular la raíz cuadrada de a^2 + b^2 con mayor precisión con mayor precisión, pero no estoy seguro de cómo proceder desde aquí.
sqrt(a^2 + b^2) es simplemente std::hypot(a,b) . Con un poco de suerte, eso ya es preciso.
... calcular el logaritmo complejo de un doble ...
El código podría usar double real_part = (double) clog(x) .
Para calcular la parte real de un logaritmo complejo de un doble sin usar clog(x) cerca de |x| == 1.0 , considere usar log1p() *1 para obtener un mejor resultado de precisión.
El problema central es |x| - 1.0 puede sufrir una gran pérdida de precisión y este es el primer paso para determinar log() .
0.5 * log(a^2 + b^2) es matemáticamente como 0.5 * logp1(a^2 + b^2 - 1) . cuando |x| está cerca de 1.0 y |a| > |b| , use 0.5 * logp1((a-1)*(a+1) + b^2) . Esto resta el 1.0 de |a| exactamente y conserva la precisión con (a-1)*(a+1) + b^2 . Esto difiere de restar 1.0 de x ya que el cálculo de x ya ha perdido una precesión importante.
#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)); }Producción
0.000976562 1 3.57628e-07 3.57628e-07 Re: "Intenté usar fma() ..." --> Algunos fma() son de baja calidad .
*1 Las funciones log1p calculan el logaritmo en base e (natural) de 1 más el argumento.
¿Has probado log(sqrt(a*a + b*b)) ? Normalmente, la raíz cuadrada tiene el efecto opuesto al cuadrado, así que probablemente si intentas calcularla, harás el log de un número más adecuado.
De todos modos, para calcular logaritmos de números cercanos a 1.0, probablemente puedas calcular las derivadas de log(z + 1) para z == 0 y obtendrás un mejor enfoque, ya que la función es analítica en el círculo de radio < 0.5 y por lo tanto, obtendrá una buena aproximación de Taylor. Esta aproximación se escribe a continuación (gracias a 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))Esta es una serie que converge absolutamente en un círculo abierto de radio 1, por lo que para calcular valores cercanos a 1 o el logaritmo es muy adecuado (de hecho, es lo que se usa en muchos lugares).
Por supuesto, si desea resolver esto en el plano complejo, debe operar el cálculo como números complejos.