better complex tan()

[achan001]

core_command6.cc use identity: tan(A+B) = (sin(2A) + sin(2B)) / (cos(2A) + cos(2B))

tan(x+iy) = (sin(2x) + i sinh(2y)) / (cos(2x) + cosh(2y))

With cos() range of -1 to 1, denominator might suffer cancellation errors. Example, for tan(1.57079632679+1.57079632679E-10i)

Mathematica (34 sig. diigts) = 198259844.2684457351073074341883681 + 6360017400.061800948953381036777052i

Free42: (lost 20 sig. digits !) = 198259844.2684462500753790941269250 + 6360017400.061817468717674802797319i

Better formula is tan(A+B) = (sin(A)cos(A) + cos(B)sin(B)) / (cos(A)cos(A) - sin(B)sin(B))

tan(x+iy) = (sin(x)cos(x) + i cosh(y)sinh(y)) / (cos(x)^2 + sinh(y)^2)

Denominator is now safe from cancellation errors. Divide num, den by cos(x)^2: (with finite-precision, cos(x) ≠ 0, thus will not hit with divide-by-zero issue)

tan(x+iy) = (tan(x) + i cosh(y)*sinh(y)/cos(x)^2) / (1 + sinh(y)^2/cos(x)^2)

1/cos(x)^2 = 1 + tan(x)^2 cosh(y) = sqrt(1 + sinh(y)^2)

Simplify with t = tan(x), S = sinh(y), a = S + St^2, b = 1 + Sa

tan(x+iy) = t/b + i sqrt(1+SS) a/b

If b overflowed, it implied cosh(y) ≈ |sinh(y)|, imaginery part = sign(y) = ±1

Redo example with new formula, using Free42 Decimal:

x + iy = 1.57079632679+1.57079632679E-10i t = tan(x) = 204222536562.2478227241769513521222 S = sinh(y) = 1.570796326790000000006459640975002e-10 a = S+Stt = 6551295804822.614830404305805658496 b = 1+Sa = 1030.075138593010014316438138894843

t/b + i sqrt(1+SS) * a/b = 198259844.2684457351073074341883681 + 6360017400.061800948953381036777053i

[achan001]

Simplify with t = tan(x), S = sinh(y), a = S + St^2, b = 1 + Sa

tan(x+i_y) = t/b + i sqrt(1+SS) * a/b

If b overflowed, it implied cosh(y) ≈ |sinh(y)|, imaginery part = sign(y) = ±1

With finite precision, and SS ≫ 1, such that 1+SS → S*S

cosh(y) = sqrt(1+SS) → sqrt(SS) → abs(S)

This assumed binary math, and correctly rounded square root.

For, decimal system, this might not be true. A safer test is to test finiteness of numerator too.

If both finite, imag_part = sqrt(1+SS)a / b , else imag_part = sign(a) = sign(y)

[achan001]

tan(x+iy) = (tan(x) + i cosh(y)*sinh(y)/cos(x)^2) / (1 + sinh(y)^2/cos(x)^2)

Proof:, Let t = tan(A), T = tan(B)

tan(A+B) = (t+T) / (1 - t*T)
= (t+T) * (1+t*T) / (1 - t^2*T^2) = (t*(1+T^2) + T*(1+t^2)) / ((1+T^2) - T^2*(1+t^2)) = (t + T/(1+T^2)*(1+t^2)) / (1 - T^2/(1+T^2)*(1+t^2))

(1+t^2) = 1/cos(A)^2 T/(1+T^2) = sin(B)*cos(B) T^2/(1+T^2) = sin(B)*sin(B)

Let A=x, B=i*y, we derived the quoted identity.

[thomasokken]

Cool, thanks! I'll take a closer look at this and update complex TAN and TANH with these improved algorithms. I should be able to do this sometime next week; stay tuned...

[thomasokken]

I'm back after dealing with my move for the last few weeks.

Albert, I see you've presented several possible improvements for complex tan. Which one of those would you recommend I use? Accuracy and robustness are more important than speed...

[achan001]

I would pick the example from first post, tan(x+i*y) using t=tan(x) and S=sinh(y)

It has the advantage of tan(x+0.0i) returning tan(x)+0.0i, without rounding errors.

[thomasokken]

If b overflows, what does that mean for the real part?

[achan001]

If b overflows, what does that mean for the real part?

a = S + Stt b = 1 + Sa = 1 + SS * (1+tt)

t = tan(x), or its inverse, with correct angle reduction code, should never be huge. If b overflowed, it is mostly due to huge SS, we can drop the +1:

re(tan(x+yi) ≈ t/(Sa) = 1/(a/t S)

Note that denominator (a/t S) = (t+1/t) SS, and abs(t+1/t) ≥ 2 If denominator still overflow, under-flowed 0.0 real part is correct.

[achan001]

If denominator still overflow, under-flowed 0.0 real part is correct.

Floating point system may have wider range for small numbers. Overflowed denominator might still be able to produce tiny sub-normal reciprocals.
Instead of flipping t to the denominator, we can scale down exponents.

To be safe, we over-scaled down, for both numerator and denominator. For IEEE double, let k = 2^-500. For Intel decimal128, k = 10^-3000

If b overflowed, re(tan(x+yi) ≈ t/(Sa) = (kt) / (kSa)

[thomasokken]

How about this? I'm trying to keep the error checking simple, this is close to the previous version in that regard, but does use the formula that avoids cancellation in the denominator...

static int mappable_tan_c(phloat xre, phloat xim, phloat *yre, phloat *yim) {
    /* NOTE: DEG/RAD/GRAD mode does not apply here. */

    if (xim == 0) {
        *yim = 0;
        return math_tan(xre, yre, true);
    } else if (xre == 0) {
        *yre = 0;
        *yim = tanh(xim);
        return ERR_NONE;
    }
    if (p_isnan(xre) || p_isnan(xim)) {
        *yre = NAN_PHLOAT;
        *yim = NAN_PHLOAT;
        return ERR_NONE;
    }

    phloat sinxre, cosxre;
    p_sincos(xre, &sinxre, &cosxre);
    phloat sinhxim = sinh(xim);
    phloat coshxim = cosh(xim);
    phloat d = cosxre * cosxre + sinhxim * sinhxim;
    if (d == 0) {
        if (flags.f.range_error_ignore) {
            *yre = POS_HUGE_PHLOAT;
            *yim = POS_HUGE_PHLOAT;
            return ERR_NONE;
        } else
            return ERR_OUT_OF_RANGE;
    }
    if (p_isinf(d) != 0) {
        *yre = 0;
        *yim = xim < 0 ? -1 : 1;
        return ERR_NONE;
    }
    *yre = sinxre * cosxre / d;
    *yim = sinhxim * coshxim / d;
    int inf;
    if ((inf = p_isinf(*yre)) != 0) {
        if (flags.f.range_error_ignore)
            *yre = inf < 0 ? NEG_HUGE_PHLOAT : POS_HUGE_PHLOAT;
        else
            return ERR_OUT_OF_RANGE;
    }
    if ((inf = p_isinf(*yim)) != 0) {
        if (flags.f.range_error_ignore)
            *yim = inf < 0 ? NEG_HUGE_PHLOAT : POS_HUGE_PHLOAT;
        else
            return ERR_OUT_OF_RANGE;
    }
    return ERR_NONE;
}

[achan001]

Hello, I see you picked this version.

tan(x+iy) = (sin(x)cos(x) + i cosh(y)sinh(y)) / (cos(x)^2 + sinh(y)^2)

For finite complex argument, tan(x+iy) should not return ERR_OUT_OF_RANGE*

>>> from math import *
>>> cos(pi/2) # size not much smaller than machine epsilon
6.123233995736766e-17

With correct angle reduction code, d = (cos(x)^2 + sinh(y)^2) > 0

• Test for d=0 is dead code.

For *yre, it is easier to explain with t = tan(x), S = sinh(y)

yre = t / (1 + SS (1 + tt))

With correct angle reduction code, |t| should never be huge, and |t| ≥ |*yre|

• Test for infinite *yre is dead code.

yim = (cosh(y)sinh(y)) / (cos(x)^2 + sinh(y)^2)

When |sinh(y)| is huge, it is numerically equal to cosh(y) For example, with Free42-Decimal:

sinh(40) = cosh(40) = 117692633418509992.7039499553745174

If denominator does not overflow, numerator will also not overflow. Since the code already handled overflowed d ...

• Test for infinite *yim is dead-code.

[achan001]

HP71B MATH-ROM also use same formula, from file cpr.a http://www.jeffcalc.hp41.eu/emu71/mathrom.html

* ************************
* doctan
* tan((x,y)) = ( sin(x)*cos(x)/d , sinh(y)*cosh(y)/d )
* d is sinh(y)^2+cos(x)^2
* other expression (used in HPBW):
* tan((x,y)) = ( sin(2*x)/d , sinh(2*y)/d )
* d is cos(2*x)+cosh(2*y)
* ************************

Expectedly, it does not check for zero d. Unexpected, it does not check for infinite d.

Since *yim numerator and denominator should overflow at the same time, it check for division exception.

m8285   GOSUBL dv15s  m7818    sinh(y)*cosh(y)/d = imag part
        GOSUBL orsb3  m7879
        ?XM=0                  exception?
        GOYES  m82A9             no
        A=0    W               yes, puts +/-1 in (A,B)

[achan001]

if ((inf = p_isinf(*yim)) != 0) {
    if (flags.f.range_error_ignore)
        *yim = inf < 0 ? NEG_HUGE_PHLOAT : POS_HUGE_PHLOAT;
    else
        return ERR_OUT_OF_RANGE;
}

This is not only dead-code, returns are wrong.

yim = (cosh(y)sinh(y)) / (cos(x)^2 + sinh(y)^2)

If yim not finite, it suggested both numerator and denominator overflowed. This only happens when |y| is big. Numerically, sinh(|y|) = cosh(y) = e^|y| / 2 It should have yim returned sign(y) = ±1

The only way to force |*yim| big is when denominator is tiny. With tiny y, sinh(y) ≈ y, cosh(y) ≈ 1

Let u = min(|cos(x)|), not much smaller than machine epsilon. With small denominator, we have:

|*yim| = |y| / (u^2 + y^2) = 1 / (u^2/|y| + |y|) ≤ 1 / (2u)

From previous post, we showed |*yre| ≤ max(tan(x))

tan(z) real and imag parts size not much bigger than reciprocal of machine epsilon

[achan001]

I guess it is OK to check d=0. Technically d≠0, but it might underflow to 0

This is highly unlikely, but possible. (I would love to see an example ...) From Accurate trigonometric functions for large arguments

X = 0x1.6ac5b262ca1ffp+849 = 5.3193726483265414167 ... * 10^255
X mod (2π) = 1.57079632679489661970003828406521420 ...
π / 2      = 1.57079632679489661923132169163975144 ...
diff       = 0.00000000000000000046871659242546276 ... ≈ 5*10^-19

This is amazing ! However, reduced angle is still not much smaller than machine epsilon.

>>> from gmpy2 import *
>>> x = mpfr('0x1.6ac5b262ca1ffp+849')
>>> y = cos(x)
>>> mpc(x, y)
mpc('5.3193726483265414e+255-4.6871659242546277e-19j')
>>> tan(_)
mpc('-1.066742692876852e+18-1.066742692876852e+18j')
>>> 1/(2*y)
mpfr('-1.066742692876852e+18')

For finite, non-zero d, infinity check for (yre, yim) may be skipped.

[thomasokken]

I guess it is OK to check d=0. Technically d≠0, but it might underflow to 0

OK... I'm also not sure if "correct angle reduction" is a safe assumption to make. I looked at the behavior of SIN on an HP calculator (I think it was the 42S, but I'm not sure, I can't find my notes at the moment), and for large arguments, it starts losing digits, I think it only used 3 extra digits of pi for the angle reduction. I don't think hardware FPUs or the Intel library handle arbitrarily large arguments correctly, either.

I understand that the check for infinite yre isn't really necessary, since the numerator is sin(xre)cos(xre), which should never be outside [-0.5, 0.5]. And the possibility of sinh(xim)*cosh(xim) getting large is already covered by the check for infinite d, and so the NaN caused by Inf/Inf should never happen.

So if I drop the checks for infinite yre and yim at the end, and leave the rest as it is, it should be OK?

[achan001]

So if I drop the checks for infinite yre and yim at the end, and leave the rest as it is, it should be OK?

YES

[Ajaja]

if (d == 0) {
        if (flags.f.range_error_ignore) {
            *yre = POS_HUGE_PHLOAT;
            *yim = POS_HUGE_PHLOAT;
            return ERR_NONE;
        } else
            return ERR_OUT_OF_RANGE;
    }

Shouldn't it be NEG_HUGE_PHLOAT/POS_HUGE_PHLOAT, depending on sinxre * cosxre and sinhxim * coshxim sign?

[thomasokken]

Strictly speaking, tan((n+1/2)*pi) is undefined, since the value of the limit depends on the direction you're approaching from. "Divide by 0" would be the appropriate error message. I went with infinity for consistency with TAN for real numbers on the HP-42S; Mathematica uses ComplexInfinity, i.e. undirected infinity.

[thomasokken]

I applied the changes to TAN and TANH, and uploaded new Plus42 builds to https://thomasokken.com/free42/download/test/