complex log

[Albert Chan]

(04-10-2023 06:47 PM)Albert Chan Wrote:  I considered Kahan's clog(z), but Python's cmath.log code is simpler.

Code:

function complex.log(z)
    local x, y = I.real(z), I.imag(z)
    local h = hypot(x, y)
    if 0.71 <= h and h <= 1.73 then
        if abs(x) < abs(y) then x,y = y,x end
        h = log1p((x+1)*(x-1) + y*y)/2
    else
        h = log(h)
    end
    return I.new(h, I.arg(z))
end

I found a better way to calculate (x^2-1) = (x-1)*2 + (x-1)^2
If x is slightly bigger than 1, RHS is sum of 2 positive terms, no cancellation errors.
If x is slightly below 1, statistically, errors is much smaller than doing (x+1)*(x-1)

How to calculate 1-x·x in floating-point
To SSE, or not to SSE, that is the question

I might as well fix when (x,y) both finite, yet hypot(x,y) overflowed.
Also, optimized a bit, not to use hypot() until really needed.

Clog code for lua lcomplex.c

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 if (isinf(h)) {
        h = log(hypot(x*.5,y*.5)) + M_LN2;
    } else
        h = log(hypot(x,y));
    return CMPLX(h, a);
}