complex log

[Albert Chan]

(02-07-2024 12:49 AM)Albert Chan Wrote:  .5 < h < 3
min(|x|) = sqrt(.5/2) = 0.5
max(|x|) = sqrt(3/1) ≈ 1.732

We can always scale z to z', with its h = |z'|^2 within this range.
Clog can be coded using only log1p formula, without other branches.

ln(z) = ln(z') + b*ln2      , where z' = z/2^b
ln2 = ln2_hi + ln2_lo

b = binary exponent of z parts = -1074 .. 1023 +1 (*)
if ln2_hi at most 53-11 = 42 bits, then b*ln2_h is exact.

Update: I use constants from https://www.netlib.org/fdlibm/e_log.c
ln2_hi has only 32 bits, which is ok (we have ln2 of 32+53 = 85 bits)
This give more room for future use, say, with fractional b's.

Code:

static const double
ln2_hi  =  6.93147180369123816490e-01,    /* 3fe62e42 fee00000 */
ln2_lo  =  1.90821492927058770002e-10;    /* 3dea39ef 35793c76 */

static double shift_ln2(double x, double b) {    
    double y = b*ln2_hi, z = b*ln2_lo;
    double s = x + y;
    return z - (s-y-x) + s; // Kahan sum x + b*ln2
}

Clog, version 2

Code:

static Complex Clog(Complex z) {      
    double x = Real(z), y = Imag(z);
    double a = atan2(y,x), b;
    if ((x=fabs(x)) < (y=fabs(y))) {b=x; x=y; y=b;}
    if (isinf(b=logb(x))) return CMPLX(b,a);
    x = ldexp(x, -b);
    y = ldexp(y, -b);
    if (x*x+y*y > 2.718281828) {x*=.5; y*=.5; ++b;}
    x -= 1;
    x = log1p(2*x + x*x + y*y)/2;
    return CMPLX(shift_ln2(x,b), a);
}

Which versions are best depends on usage pattern.
It also depends on speed/accuracy of built-in log1p vs log.

(*) b might be +1 higher, to reduce h down to below e
Code doesn't use h=x²+y², but h = [1,8) → [e/4,e)       → x ≈ [.5829, 1.649)

Why e, not 3 Wrote:Let r = ln(|z'|) = ln(h)/2

h = [1, e)      --> r = [0, .5)
h = [e, 8)      --> r = [.5, 1.040)
h = [e/4, 2)   --> r = [-.1931, .3466)

if h < e, then (ulp of r) ≤ (ulp of 1/4)