complex log

[Albert Chan]

(02-10-2024 09:55 PM)Albert Chan Wrote:  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.

Version 1b with smaller shift (M_LN2*61) avoid most double rounding problem, but not all.
There are still holes where log part and final result ULP size are different.

This version traded hypot() for frexp(), and have no performance penalty.
I still use 732 for rough shift, but any number that get finite (x²+y²) is OK

(x²+y²) = h * 2^n
log(x²+y²)/2 = log(h)/2 + (n/2)*log(2)

With h = [0.5, 1), log(h)/2 = [-.3466, -0). max ulp = (ulp of 1/4), same as version 2

Clog Version 1c

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 = 0x1p732, b = -732;
        if (isinf(h)) {m = 1/m; b = -b;};
        x *= m; y *= m;
        int n;
        h = frexp(x*x+y*y, &n); 
        h = shift_ln2(log(h)/2, b + .5*n);
    }
    return CMPLX(h, a);        
}