Improve accuracy of log1p() (#803)

Patch, ported from Tim Peters version at: https://github.com/mpmath/mpmath/issues/790#issuecomment-2084094631

Co-authored-by: Tim Peters <tim.peters@gmail.com>

mpmath/functions/functions.py

@@ -166,9 +166,80 @@ def expm1(ctx, x):
 def log1p(ctx, x):
     if not x:
         return ctx.zero
-    if ctx.mag(x) < -ctx.prec:
-        return x - 0.5*x**2
-    return ctx.log(ctx.fadd(1, x, prec=2*ctx.prec))
+    LOG1P_EXTRAPREC = 10  # ctx._wrap_specfun()
+    # Note that all cases could by handled by log(1+c) provided the
+    # add is done exactly. Our aim here is to be much faster than that,
+    # especially when |c| is small.
+    c = ctx.convert(x)
+    cmag = ctx.mag(c)
+    a, b = c.real, c.imag
+    wp = ctx.prec
+    if cmag >= -wp:
+        # |c| isn't very small. We call log(1+c) instead, but
+        # are careful about the precision used by the add. The
+        # real part of the result is log(|c+1|). That's
+        # determined by 1 + 2*a + a**2 + b**2, and the add has
+        # to preserve enough info so that no important bits of
+        # that sum are lost. It doesn't matter to this that 2*a,
+        # a**2, etc, are not computed explicitly here: we're
+        # deducing how many bits have to be present in the sum
+        # for log() to "reverse engineer" the value of 2*a +
+        # a**2 + b**2 to `prec` good bits,
+        if cmag < 4:
+            # |c| isn't very small, or large.
+            if ctx.mag(a) > ctx.mag(b):
+                # `a` already contributes the most to c's norm.
+                # After adding 1, it will utterly dominate it.
+                # We only need enough extra precision to avoid
+                # losing any of a's `prec` most significant bits
+                # when addiog, `b**2` is too small to matter.
+                wp *= 2
+            else:
+                # b**2 is the larger of the square terms. The
+                # smallest b can be is about 2**-prec, so the
+                # smallest b**2 can be is about 2**(-2*prec). So
+                # for a bit to matter compared to b**2, it has
+                # to be at least about 2**(-3*prec). Bits of 2*a
+                # (if any) >= 2**(-3*prec) will be preserved if
+                # we use 3*prec bits for the add.
+                wp *= 3
+        # Else (cmag >= 4), |c+1| >= |c| - 1 is so large that
+        # working precision is fine (although that takes some
+        # careful analysis for cmag=4, given that .mag() _may_
+        # return a rexult too large by 2), So leave wp alone.
+        arg = ctx.fadd(1.0, c, prec=wp)
+        result = ctx.log(arg)
+    else:
+        # Else c is "very small", and we use a series expansion,
+        # c - c**2/2. The real part of that is a+(b*b-a*a)/2,
+        # and the imag part b-a*b. Given that cmag < -prec, it
+        # can be shown that "a*b" is numerically insignifcant in
+        # the imag part, and _usually_ the "a*a/2" in the real
+        # part. What remains is cheap to compute. In the real
+        # part, though, if `a` is negative, the remaining
+        # a+b**2/2 can suffer massive cancellation - even total.
+        real = a + b*b*0.5 # usually the real part of the result
+        if (a < 0.0
+            and ctx.mag(real) <= ctx.mag(a) - LOG1P_EXTRAPREC):
+            # The guard bits were lost to cancellation. Rare. At
+            # the contrived
+            # -1.999999873062092e-40+1.999999936531045e-20j
+            # _all_ bits cancel out. Since a ~= -b*b/2 in this
+            # case, and |b| is at largest (worst case) about
+            # 2**-prec, |a| is about 2**(-2*prec), and the true
+            # result may be as small as a**2/2, which is about
+            # 2**(-4*prec), of which we want the leading prec
+            # bits. To get the leading prec bits starting at
+            # 2**(-4**prec) from addends starting at
+            # 2**-(2*prec), we need the subtraction to handle
+            # 3*prec bits (the first 2*prec of which may cancel
+            # to exactly 0).
+            a2 = a*a # only need at worst prec bits
+            b2 = ctx.fmul(b, b, prec=2*wp)
+            diff = ctx.fsub(b2, a2, prec=3*wp)
+            real = a + ctx.ldexp(diff, -1)
+        result = ctx.mpc(real, b)
+    return result
 
 @defun_wrapped
 def powm1(ctx, x, y):

mpmath/tests/test_functions.py

@@ -967,6 +967,13 @@ def test_log1p():
     assert log1p(inf) == inf
     assert log1p(1e-50).ae(1e-50)
     assert (log1p(1e-10)*1e10).ae(0.99999999995)
+    # issue 790
+    assert log1p(1.8370676479640493e-39-4.6885882517313053e-20j).ae(2.93621063767769e-39-4.6885882517313053e-20j, 0)
+    assert log1p(-2.0476815825463086e-80-2.0235857941734692e-40j).ae(-2.3184935597344513e-84-2.0235857941734692e-40j, 0)
+    assert log1p(-6.4922176418510124e-21+1.1394926627101214e-10j).ae(-1.4201199664289643e-37+1.1394926627101214e-10j, 0)
+    assert log1p(-1.430796815051627e-72+1.691624553529315e-36j).ae(4.6709293580298264e-91+1.691624553529315e-36j, 0)
+    assert log1p(-3.1061140011623543e-21+7.881768838480807e-11j).ae(4.3173401185459216e-38+7.881768838480807e-11j, 0)
+    assert log1p(-1.999999873062092e-40+1.999999936531045e-20j).ae(1.9999997461241924e-80+1.999999936531045e-20j)
 
 def test_powm1():
     assert powm1(2,3) == 7