Use mpf_log1p in acos_asin() helper (implementing Hull et al algorithm)

Closes #787

mpmath/libmp/libelefun.py

@@ -723,6 +723,16 @@ def mpf_log(x, prec, rnd=round_fast):
         m -= n*ln2_fixed(wp)
     return from_man_exp(m, -wp, prec, rnd)
 
+def mpf_log1p(x, prec, rnd=round_fast):
+    """
+    Computes log(1+x) accurately.
+    """
+    wp = prec + 10
+    u = mpf_add(fone, x, wp*2)
+    return mpf_mul(mpf_log(u, wp),
+                   mpf_div(x, mpf_sub(u, fone, wp),
+                           wp), prec, rnd)
+
 def mpf_log_hypot(a, b, prec, rnd):
     """
     Computes log(sqrt(a^2+b^2)) accurately.

mpmath/libmp/libmpc.py

@@ -8,9 +8,9 @@ from .backend import MPZ
 from .libelefun import (mpf_acos, mpf_acosh, mpf_asin, mpf_atan, mpf_atan2,
                         mpf_cos, mpf_cos_pi, mpf_cos_sin, mpf_cos_sin_pi,
                         mpf_cosh, mpf_cosh_sinh, mpf_exp, mpf_fibonacci,
-                        mpf_log, mpf_log_hypot, mpf_nthroot, mpf_phi, mpf_pi,
-                        mpf_pow_int, mpf_sin, mpf_sin_pi, mpf_sinh, mpf_tan,
-                        mpf_tanh)
+                        mpf_log, mpf_log1p, mpf_log_hypot, mpf_nthroot,
+                        mpf_phi, mpf_pi, mpf_pow_int, mpf_sin, mpf_sin_pi,
+                        mpf_sinh, mpf_tan, mpf_tanh)
 from .libintmath import giant_steps, lshift, rshift
 from .libmpf import (ComplexResult, fhalf, finf, fnan, fninf, fnone, fone,
                      from_float, from_int, from_man_exp, ftwo, fzero, mpf_abs,
@@ -719,7 +719,7 @@ def acos_asin(z, prec, rnd, n):
             Am1 = mpf_shift(mpf_add(c1, c2, wp), -1)
         # im = log(1 + Am1 + sqrt(Am1*(alpha+1)))
         im = mpf_mul(Am1, mpf_add(alpha, fone, wp), wp)
-        im = mpf_log(mpf_add(fone, mpf_add(Am1, mpf_sqrt(im, wp), wp), wp), wp)
+        im = mpf_log1p(mpf_add(Am1, mpf_sqrt(im, wp), wp), wp)
     else:
         # im = log(alpha + sqrt(alpha*alpha - 1))
         im = mpf_sqrt(mpf_sub(mpf_mul(alpha, alpha, wp), fone, wp), wp)

mpmath/tests/test_functions.py

@@ -326,6 +326,12 @@ def test_asin():
     assert asin(mpc(+2, 0)).ae(mpc(+pi2, -log(2 + sqrt(3))))
     assert asin(mpc(0.5, 0)).ae(pi/6)
 
+    # issue 787
+    assert asin(mpc(0, 1e-22)).ae(1e-22j)
+    mp.prec = 700
+    assert asin(mpc(0, 1e-220)).ae(1e-220j)
+    mp.prec = 53
+
 def test_acos():
     pi4 = pi/4
     assert acos(mpc(+inf, +inf)) == mpc(+pi4, -inf)