complex log

[Albert Chan]

Here is updated OP version 1, fixed its subnormal accuracy issue.

Instead of LN2 shift for overflowed h, (DBL_MANT_DIG*M_LN2) for subnormal, I combined them.

relative error = (exact - approx) / exact

relerr(ln(2)) ≈ 3.346E-17
relerr(ln(2)*53) ≈ 1.835E-17
relerr(ln(2)*61) ≈ -6.781E-19

Last constant has 59 bits accuracy (~ 18 digits), the most accurate for our domain.

float(ln(2)*61) = 0x1.52417db067f38p+5 = 42.28197801415666390...

Clog Version 1b

Code:

static Complex Clog(Complex z) {
    double x = Real(z), y = Imag(z);
    double h = x*x + y*y;
    double a = atan2(y,x);
    if (0.5 < h && h < 3.0) {
        if ((h=fabs(x)) < (y=fabs(y))) {h=y; y=x;}
        h -= 1; /* x^2-1 = 2*(x-1) + (x-1)^2 */
        h = log1p(2*h + h*h + y*y)/2;
    } else if (isnormal(h)) {
        h = log(h)/2;
    } else {
        double m = 0x1p61, b = -M_LN2*61;
        if (isinf(h)) {m = 1/m; b = -b;};
        h = log(hypot(m*x,m*y)) + b;
    }
    return CMPLX(h, a);        
}

(02-08-2024 09:04 PM)Albert Chan Wrote:  lua> t = 0x1p-1074
lua> I(3*t, 5*t) :log() -- version 2
(-742.6768916590731+1.0303768265243125*I)

We get exactly the same answer for version 1b

Update: previously, I used shift of (M_LN2*732) to avoid hypot() function.
With ULP of log part different than final result, answer is always double rounded.
Double rounding is acceptable if log part ULP is much smaller, but it isn't.

Also, for big |log(x)|, error caused by hypot() inside log may be ignored.

log(x*(1-ε)) ≈ log(x) + (-xε) / x = log(x) - ε

if relerr(x) is ε, then relerr(log(x)) is ε/log(x), much smaller than ε