mp.asin is inaccurate for small complex inputs

[pearu]

trafficstars

As in the title.

Reproducer:

>>> mpmath.mp.asin(mpmath.mpc(0, 1e-22))
mpc(real='0.0', imag='0.0')
>>> mpmath.fp.asin(mpmath.mpc(0, 1e-22))  # expected result
1e-22j

A workaround is to increase precision:

>>> mpmath.mp.prec=80
>>> mpmath.mp.asin(mpmath.mpc(0, 1e-22))
mpc(real='0.0', imag='9.9999996826552253889673883e-23')

Points z=x+yj in the following (approximate) region

|x| < |y| < 1e-19

of the complex plane suffer from the defect reported here.

[skirpichev]

Perhaps, we could just increase working precision:

diff --git a/mpmath/libmp/libmpc.py b/mpmath/libmp/libmpc.py
index 2e22b22..b0893d4 100644
--- a/mpmath/libmp/libmpc.py
+++ b/mpmath/libmp/libmpc.py
@@ -609,7 +609,7 @@ def acos_asin(z, prec, rnd, n):
     and by abs(a) <= 1, in order to improve the numerical accuracy.
     """
     a, b = z
-    wp = prec + 10
+    wp = prec + 25
     # infinite real component
     if a in (finf, fninf):
         if b in (finf, fninf):

[pearu]

Yes, but notice that the increase by 25 is likely insufficient or it must depend on the input. For instance,

>>> mpmath.mp.prec=700
>>> mpmath.mp.asin(mpmath.mpc(0, 1e-220))
mpc(real='0.0', imag='0.0')
>>> mpmath.mp.prec=800
>>> mpmath.mp.asin(mpmath.mpc(0, 1e-220))
mpc(real='0.0', imag='9.9999999999999999239355998681967174227375122814278010784257107926662288959827321849441348805654645781222578141680149580443831699278821999049728767619913188664117273184328232545396962216463347269815350517239219609176864224305536231460588661161e-221')

So, I think the algorithm in mpmath itself should be fixed because the same Hull et al algorithm is used elsewhere with double FP without having the defect reported here.

[pearu]

Looking at the code, I see that mpmath implements the Hull et al algorithm partially on the so-called "safe" region (see p 307) but it does not handle the non-safe regions (see p 315). While some parts of the non-safe regions are ok to dismiss because of MP math, then handling the non-safe region 2 in the paper corresponds to resolving this issue. In short, asin(z) in region 2 should read:

# assume z.imag >= 0, z.real => 0
if z.imag < eps * abs(z.real - 1):
    if z.real < 1:
        abs_real = arcsin(z.real)
        abs_imag = z.imag / sqrt((1 - z.real) * (z.real + 1))
    else:
        abs_real = pi / 2
        abs_imag = log1p(z.real - 1 + sqrt((z.real - 1) * (z.real + 1)))

where eps = 2 ** (-prec). HTH