Repository: mpmath/mpmath (python)

Issue:
Poor result in log1p(complex)
A case where lop1p of a complex number gets a poor real part:

```python
>>> import mpmath
>>> mpmath.__version__
'1.3.0'
>>> from mpmath import mp, mpc, log, log1p, workprec
>>> c = 1.8370676479640493e-39 - 4.6885882517313053e-20j
>>> complex(log1p(c))
(1.0991429897136409e-39-4.6885882517313053e-20j)
>>> with workprec(500): complex(log1p(c))
(2.93621063767769e-39-4.6885882517313053e-20j)
```

With the higher precision, the real part of the result is entirely different. I don't understand the code:

```python
>>> mp.mag(c), mp.prec
(-63, 53)
```

Because -63 < -53, I expected the code to take the branch that uses a 2-term expansion, which works fine:

```python
>>> c - 0.5*c**2 # matches the high-precision result
(2.93621063767769e-39-4.6885882517313053e-20j)
```

But I must not have found the correct source code, since editing it myself has no visible effect.

The other branch in the code calls `log()` instead after adding 1 to c with double precision:

```python
>>> complex(log(mp.fadd(c, 1, prec=106)))
(1.0991429897136409e-39-4.6885882517313053e-20j) # reproduces the "bad" result
```

It's not enough extra precision to help. In this case,

```python
>>> mp.mag((2*c.real + c.real**2 + c.imag**2))
-127
>>> complex(log(mp.fadd(c, 1, prec=127+53))) # so this is enough
(2.93621063767769e-39-4.6885882517313053e-20j)
>>> complex(log(mp.fadd(c, 1, prec=127+53-2))) # but 2 bits less not quite enough
(2.9362106376776874e-39-4.6885882517313053e-20j)
```

To be clear, "the code" I'm talking about is this, from functions/functions.py:

```python
@defun_wrapped
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))
```

LATER: found the right source code. `ctx.mag(x)` and `-ctx.prec` in the above _both_ return -63 in this case, so that's why the "wrong" branch is taken. I have no idea how `ctx.prec` ended up at 63. As shown above, I never changed it from the default 53.

Pull request description:
Improve accuracy of log1p()
Patch, ported from Tim Peters version at:
https://github.com/mpmath/mpmath/issues/790#issuecomment-2084094631

And few tests from the issue thread

Closes #790

Patch:
diff --git a/.gitignore b/.gitignore
index 7698467..b434fda 100644
--- a/.gitignore
+++ b/.gitignore
@@ -62,3 +62,6 @@ htmlcov/
 
 # PyCharm files
 .idea/
+
+# Ignore files cached by Hypothesis (including examples directory so far)
+.hypothesis/
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 11db29a..0644870 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -6,6 +6,9 @@ build:
     os: ubuntu-22.04
     tools:
         python: "3"
+    jobs:
+        post_checkout:
+            - git fetch --unshallow
 python:
     install:
         - method: pip
diff --git a/README.rst b/README.rst
index ddde56d..3019c88 100644
--- a/README.rst
+++ b/README.rst
@@ -145,7 +145,7 @@ are human-readable, but can be compiled to prettier HTML using
 
 The most recent documentation is also available in HTML format:
 
-https://mpmath.org/doc/current/
+https://mpmath.readthedocs.io/
 
 3. Running tests
 ----------------
diff --git a/docs/basics.rst b/docs/basics.rst
index 680257b..3d40bd4 100644
--- a/docs/basics.rst
+++ b/docs/basics.rst
@@ -53,7 +53,7 @@ The ``mpf`` type is analogous to Python's built-in ``float``. It holds a real nu
     >>> mpf(mpf(2))
     mpf('2.0')
     >>> mpf("inf")
-    mpf('+inf')
+    mpf('inf')
 
 The ``mpc`` type represents a complex number in rectangular form as a pair of ``mpf`` instances. It can be constructed from a Python ``complex``, a real number, or a pair of real numbers:
 
diff --git a/docs/contexts.rst b/docs/contexts.rst
index 354e1a0..77ab898 100644
--- a/docs/contexts.rst
+++ b/docs/contexts.rst
@@ -193,12 +193,12 @@ Some transcendental functions are supported::
     >>> iv.exp(0)
     [1.0, 1.0]
     >>> iv.exp(['-inf','inf'])
-    [0.0, +inf]
+    [0.0, inf]
     >>>
     >>> iv.exp(['-inf',0])
     [0.0, 1.0]
     >>> iv.exp([0,'inf'])
-    [1.0, +inf]
+    [1.0, inf]
     >>> iv.exp([0,1])
     [1.0, 2.7182818284590455349]
     >>>
@@ -207,7 +207,7 @@ Some transcendental functions are supported::
     >>> iv.log([0,1])
     [-inf, 0.0]
     >>> iv.log([0,'inf'])
-    [-inf, +inf]
+    [-inf, inf]
     >>> iv.log(2)
     [0.69314718055994528623, 0.69314718055994539725]
     >>>
diff --git a/docs/setup.rst b/docs/setup.rst
index f489e6f..fd08ecf 100644
--- a/docs/setup.rst
+++ b/docs/setup.rst
@@ -17,7 +17,7 @@ of the Mpmath with pip::
 
 or some specific version with::
 
-    pip install mpmath==0.19
+    pip install mpmath==1.3.0
 
 You can install also extra dependencies, e.g. `gmpy
 <https://github.com/aleaxit/gmpy>`_ support::
@@ -84,7 +84,7 @@ Python interpreter and do the following::
 Using gmpy (optional)
 ---------------------
 
-If `gmpy <https://github.com/aleaxit/gmpy>`_ version 2.2.0rc1 or later is
+If `gmpy <https://github.com/aleaxit/gmpy>`_ version 2.2.0 or later is
 installed on your system, mpmath will automatically detect it and transparently
 use gmpy integers instead of Python integers.  This makes mpmath much faster,
 especially at high precision (approximately above 100 digits).
diff --git a/mpmath/ctx_mp_python.py b/mpmath/ctx_mp_python.py
index 2d28777..ba45e44 100644
--- a/mpmath/ctx_mp_python.py
+++ b/mpmath/ctx_mp_python.py
@@ -681,7 +681,7 @@ def convert(ctx, x, strings=True):
         prec, rounding = ctx._prec_rounding
         if isinstance(x, numbers.Rational):
             p, q = x.numerator, x.denominator
-            return ctx.make_mpf(from_rational(p, q, prec))
+            return ctx.make_mpf(from_rational(p, q, prec, rounding))
         if strings and isinstance(x, str):
             try:
                 _mpf_ = from_str(x, prec, rounding)
@@ -1095,7 +1095,7 @@ def mag(ctx, x):
             >>> mag(0.01), int(ceil(log(0.01,2)))
             (-6, -6)
             >>> mag(0), mag(inf), mag(-inf), mag(nan)
-            (-inf, +inf, +inf, nan)
+            (-inf, inf, inf, nan)
 
         """
         if hasattr(x, "_mpf_"):
diff --git a/mpmath/function_docs.py b/mpmath/function_docs.py
index db483ee..5963b33 100644
--- a/mpmath/function_docs.py
+++ b/mpmath/function_docs.py
@@ -420,7 +420,7 @@
     >>> sqrt(1+1j)
     (1.09868411346781 + 0.455089860562227j)
     >>> sqrt(inf)
-    +inf
+    inf
 
 Square root evaluation is fast at huge precision::
 
@@ -489,7 +489,7 @@
     >>> exp(-1)
     0.3678794411714423215955238
     >>> exp(inf)
-    +inf
+    inf
     >>> exp(-inf)
     0.0
 
@@ -565,7 +565,7 @@
     >>> cosh(1)
     1.543080634815243778477906
     >>> cosh(-inf), cosh(+inf)
-    (+inf, +inf)
+    (inf, inf)
 
 The hyperbolic cosine is an even, convex function with
 a global minimum at `x = 0`, having a Maclaurin series
@@ -595,7 +595,7 @@
     >>> sinh(1)
     1.175201193643801456882382
     >>> sinh(-inf), sinh(+inf)
-    (-inf, +inf)
+    (-inf, inf)
 
 The hyperbolic sine is an odd function, with a Maclaurin
 series that starts::
@@ -718,7 +718,7 @@
     >>> iv.tan([0,1])
     [0.0, 1.55740772465490223050697482944]
     >>> iv.tan([0,2])  # Interval includes a singularity
-    [-inf, +inf]
+    [-inf, inf]
 """
 
 sec = r"""
@@ -746,7 +746,7 @@
     >>> iv.sec([0,1])
     [1.0, 1.85081571768092561791175326276]
     >>> iv.sec([0,2])  # Interval includes a singularity
-    [-inf, +inf]
+    [-inf, inf]
 """
 
 csc = r"""
@@ -771,9 +771,9 @@
 
     >>> iv.dps = 25; iv.pretty = True
     >>> iv.csc([0,1])  # Interval includes a singularity
-    [1.18839510577812121626159943988, +inf]
+    [1.18839510577812121626159943988, inf]
     >>> iv.csc([0,2])
-    [1.0, +inf]
+    [1.0, inf]
 """
 
 cot = r"""
@@ -799,9 +799,9 @@
 
     >>> iv.dps = 25; iv.pretty = True
     >>> iv.cot([0,1])  # Interval includes a singularity
-    [0.642092615934330703006419974862, +inf]
+    [0.642092615934330703006419974862, inf]
     >>> iv.cot([1,2])
-    [-inf, +inf]
+    [-inf, inf]
 """
 
 acos = r"""
@@ -1627,7 +1627,7 @@
     >>> log(0)
     -inf
     >>> log(inf)
-    +inf
+    inf
 
 The natural logarithm is the antiderivative of `1/x`::
 
@@ -1957,7 +1957,7 @@
     4 24.0
     5 120.0
     >>> fac(inf)
-    +inf
+    inf
     >>> fac(0.5), sqrt(pi)/2
     (0.886226925452758, 0.886226925452758)
 
@@ -2013,7 +2013,7 @@
     4 6.0
     5 24.0
     >>> gamma(inf)
-    +inf
+    inf
     >>> gamma(0)
     Traceback (most recent call last):
       ...
@@ -2100,7 +2100,7 @@
 while higher orders tend to zero::
 
     >>> psi(0,inf), psi(1,inf), psi(2,inf)
-    (+inf, 0.0, 0.0)
+    (inf, 0.0, 0.0)
 
 Evaluation for a complex argument::
 
@@ -2181,7 +2181,7 @@
 The infinite harmonic series `1 + 1/2 + 1/3 + \ldots` diverges::
 
     >>> harmonic(inf)
-    +inf
+    inf
 
 :func:`~mpmath.harmonic` is evaluated using the digamma function rather
 than by summing the harmonic series term by term. It can therefore
@@ -2462,7 +2462,7 @@
     >>> beta(-3, 3)
     -0.333333333333333
     >>> beta(-2, 3)
-    +inf
+    inf
     >>> beta(inf, 1)
     0.0
     >>> beta(inf, 0)
@@ -2868,7 +2868,7 @@
     >>> hyper([1,2,3,4], [5,6,7], 1)   # 4F3 at finite-valued branch point
     1.141783505526870731311423
     >>> hyper([4,5,6,7], [1,2,3], 1)   # 4F3 at pole
-    +inf
+    inf
     >>> hyper([1,2,3,4,5], [6,7,8,9], 10)    # 5F4
     (1.543998916527972259717257 - 0.5876309929580408028816365j)
     >>> hyper([1,2,3,4,5,6], [7,8,9,10,11], 1j)   # 6F5
@@ -2879,7 +2879,7 @@
     >>> hyper(['1/3',1,'3/2',2], ['1/5','11/6','41/8'], 1)
     2.219433352235586121250027
     >>> hyper(['1/3',1,'3/2',2], ['1/5','11/6','5/4'], 1)
-    +inf
+    inf
     >>> eps1 = extradps(6)(lambda: 1 - mpf('1e-6'))()
     >>> hyper(['1/3',1,'3/2',2], ['1/5','11/6','5/4'], eps1)
     2923978034.412973409330956
@@ -3268,7 +3268,7 @@
     >>> hyp2f1(-2.5, 3, 4, 1)
     0.06926406926406926406926407
     >>> hyp2f1(2, 3, 4, 1)
-    +inf
+    inf
 
 Evaluation for huge arguments::
 
@@ -3395,7 +3395,7 @@
 The `U` function may be singular at `z = 0`::
 
     >>> hyperu(1.5, 2, 0)
-    +inf
+    inf
     >>> hyperu(1.5, -2, 0)
     0.1719434921288400112603671
 
@@ -3588,7 +3588,7 @@
     >>> gammainc(-3, -4, -5)
     (-0.2214577048967798566234192 + 0.0j)
     >>> lower_gamma(-3, 5)
-    +inf
+    inf
 
 If `z` is an integer, the recurrence reduces the incomplete gamma
 function to `P(a) \exp(-a) + Q(b) \exp(-b)` where `P` and
@@ -3768,7 +3768,7 @@
     >>> erfi(-1)
     -1.65042575879754
     >>> erfi(inf)
-    +inf
+    inf
     >>> erfi(-inf)
     -inf
 
@@ -3820,7 +3820,7 @@
     >>> erfinv(0)
     0.0
     >>> erfinv(1)
-    +inf
+    inf
     >>> erfinv(-1)
     -inf
 
@@ -4023,7 +4023,7 @@
     >>> ei(1)
     1.89511781635594
     >>> ei(inf)
-    +inf
+    inf
     >>> ei(-inf)
     0.0
 
@@ -4130,7 +4130,7 @@
     >>> findroot(li, 2)
     1.45136923488338105028396848589
     >>> li(inf)
-    +inf
+    inf
     >>> li(2, offset=True)
     0.0
     >>> li(1, offset=True)
@@ -4335,7 +4335,7 @@
     >>> chi(1)
     0.8378669409802082408946786
     >>> chi(inf)
-    +inf
+    inf
     >>> findroot(chi, 0.5)
     0.5238225713898644064509583
     >>> chi(2+3j)
@@ -4367,7 +4367,7 @@
     >>> shi(-1)
     -1.057250875375728514571842
     >>> shi(inf)
-    +inf
+    inf
     >>> shi(2+3j)
     (-0.1931890762719198291678095 + 2.645432555362369624818525j)
 
@@ -4650,11 +4650,11 @@
     0.3333333332991690159427932
     0.3333333333333333333333333
     >>> airyai(+inf,-2); airyai(+inf,-3)
-    +inf
-    +inf
+    inf
+    inf
     >>> airyai(-1000000,-2); airyai(-inf,-2)
     666666.4078472650651209742
-    +inf
+    inf
     >>> airyai(-1000000,-3); airyai(-inf,-3)
     -333333074513.7520264995733
     -inf
@@ -4714,7 +4714,7 @@
     >>> airybi(-1)
     0.10399738949694461188869
     >>> airybi(inf); airybi(-inf)
-    +inf
+    inf
     0.0
 
 Evaluation is supported for large magnitudes of the argument::
@@ -4832,10 +4832,10 @@
     0.0
     >>> airybi(10,-1); airybi(+inf,-1)
     147809803.1074067161675853
-    +inf
+    inf
     >>> airybi(+inf,-2); airybi(+inf,-3)
-    +inf
-    +inf
+    inf
+    inf
     >>> airybi(-1000000,-2); airybi(-inf,-2)
     0.4482883750599908479851085
     0.4482883573538263579148237
@@ -4846,7 +4846,7 @@
     -inf
     >>> airybi(-100000,-4); airybi(-inf,-4)
     2241411040.437759489540248
-    +inf
+    inf
 
 """
 
@@ -5013,7 +5013,7 @@
     >>> ellipk(-inf)
     0.0
     >>> ellipk(1)
-    +inf
+    inf
     >>> ellipk(-1)
     1.31102877714605990523242
     >>> ellipk(2)
@@ -6158,9 +6158,9 @@
 The point `x = 0` is a singularity (logarithmic if `n = 0`)::
 
     >>> besselk(0,0)
-    +inf
+    inf
     >>> besselk(1,0)
-    +inf
+    inf
     >>> for n in range(-4, 5):
     ...     print(besselk(n, '1e-1000'))
     ...
@@ -6326,17 +6326,17 @@
     >>> lambertw(e)
     1.0
     >>> lambertw(inf)
-    +inf
+    inf
     >>> lambertw(0, k=-1)
     -inf
     >>> lambertw(0, k=3)
     -inf
     >>> lambertw(inf, k=2)
-    (+inf + 12.56637061435917295385057j)
+    (inf + 12.56637061435917295385057j)
     >>> lambertw(inf, k=3)
-    (+inf + 18.84955592153875943077586j)
+    (inf + 18.84955592153875943077586j)
     >>> lambertw(-inf, k=3)
-    (+inf + 21.9911485751285526692385j)
+    (inf + 21.9911485751285526692385j)
 
 The `k = 0` and `k = -1` branches join at `z = -1/e` where
 `W(z) = -1` for both branches. Since `-1/e` can only be represented
@@ -6401,7 +6401,7 @@
     >>> barnesg(8)
     24883200.0
     >>> barnesg(inf)
-    +inf
+    inf
     >>> barnesg(0), barnesg(-1), barnesg(-2)
     (0.0, 0.0, 0.0)
 
@@ -6730,7 +6730,7 @@
     1.200973602347074224816022
     1.200973602347074224816022
     >>> loggamma(inf)
-    +inf
+    inf
 
 Huge arguments are permitted::
 
@@ -6802,7 +6802,7 @@
     >>> siegeltheta(0)
     0.0
     >>> siegeltheta(inf)
-    +inf
+    inf
     >>> siegeltheta(-inf)
     -inf
     >>> siegeltheta(1)
@@ -7174,7 +7174,7 @@
 Euler constants::
 
     >>> primezeta(1)
-    +inf
+    inf
     >>> extradps(25)(lambda x: primezeta(1+x)+log(x))(+eps)
     -0.31571845205389007685
     >>> mertens-euler
@@ -8104,7 +8104,7 @@
 Evaluation at zero::
 
     >>> whitm(1,-1,0); whitm(1,-0.5,0); whitm(1,0,0)
-    +inf
+    inf
     nan
     0.0
 
@@ -8158,11 +8158,11 @@
     >>> for m in [-1, -0.5, 0, 0.5, 1]:
     ...     whitw(1, m, 0)
     ...
-    +inf
+    inf
     nan
     0.0
     nan
-    +inf
+    inf
 
 We can verify that :func:`~mpmath.whitw` numerically satisfies the
 differential equation for arbitrarily chosen values::
@@ -9154,7 +9154,7 @@
     1.202056903159594285399738
     1.202056903159594285399738
     >>> dirichlet(1, [1])
-    +inf
+    inf
 
 The alternating zeta function::
 
@@ -9753,7 +9753,7 @@
     >>> eulerpoly(5, -inf)
     -inf
     >>> eulerpoly(6, -inf)
-    +inf
+    inf
 
 Computing Euler numbers::
 
@@ -9997,7 +9997,7 @@
     0.2988589049025509052765491
     0.2988589049025509052765491
     >>> scorerhi(+inf); scorerhi(-inf)
-    +inf
+    inf
     0.0
     >>> scorerhi(1)
     0.9722051551424333218376886
diff --git a/mpmath/functions/elliptic.py b/mpmath/functions/elliptic.py
index 6e7e906..6c0c0ed 100644
--- a/mpmath/functions/elliptic.py
+++ b/mpmath/functions/elliptic.py
@@ -241,7 +241,7 @@ def kfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
         >>> kfrom(q=1)
         1
         >>> kfrom(q=-1)
-        (0.0 + +infj)
+        (0.0 + infj)
     """
     if k is not None:
         return ctx.convert(k)
@@ -662,10 +662,10 @@ def elliprf(ctx, x, y, z):
         >>> elliprf(2,2,2)**2
         0.5
         >>> elliprf(1,0,0); elliprf(0,0,1); elliprf(0,1,0); elliprf(0,0,0)
-        +inf
-        +inf
-        +inf
-        +inf
+        inf
+        inf
+        inf
+        inf
 
     Representing complete elliptic integrals in terms of `R_F`::
 
@@ -783,7 +783,7 @@ def elliprc(ctx, x, y, pv=True):
         3.141592653589793238462643
         3.141592653589793238462643
         >>> elliprc(1,0)
-        +inf
+        inf
         >>> elliprc(5,5)**2
         0.2
         >>> elliprc(1,inf); elliprc(inf,1); elliprc(inf,inf)
@@ -856,9 +856,9 @@ def elliprj(ctx, x, y, z, p, integration=1):
         1.380226776765915172432054
         1.380226776765915172432054
         >>> elliprj(1,3,2,0); elliprj(0,1,1,0); elliprj(0,0,0,0)
-        +inf
-        +inf
-        +inf
+        inf
+        inf
+        inf
         >>> elliprj(1,inf,1,0); elliprj(1,1,1,inf)
         0.0
         0.0
@@ -1050,8 +1050,8 @@ def ellipf(ctx, phi, m):
         1.415737208425956198892166
         1.415737208425956198892166
         >>> ellipf(pi/2+eps, 1); ellipf(-pi/2-eps, 1)
-        +inf
-        +inf
+        inf
+        inf
         >>> ellipf(1.5, 1)
         3.340677542798311003320813
 
@@ -1166,9 +1166,9 @@ def ellipe(ctx, *args):
         >>> ellipe(2)
         (0.5990701173677961037199612 + 0.5990701173677961037199612j)
         >>> ellipe(inf)
-        (0.0 + +infj)
+        (0.0 + infj)
         >>> ellipe(-inf)
-        +inf
+        inf
 
     Verifying the defining integral and hypergeometric
     representation::
@@ -1325,9 +1325,9 @@ def ellippi(ctx, *args):
         >>> ellippi(2,inf)
         0.0
         >>> abs(ellippi(1,5))
-        +inf
+        inf
         >>> abs(ellippi(0.25,1))
-        +inf
+        inf
 
     Evaluation in terms of simpler functions::
 
@@ -1373,13 +1373,13 @@ def ellippi(ctx, *args):
     Some degenerate cases::
 
         >>> ellippi(1,1)
-        +inf
+        inf
         >>> ellippi(1,0)
-        +inf
+        inf
         >>> ellippi(1,2,0)
-        +inf
+        inf
         >>> ellippi(1,2,1)
-        +inf
+        inf
         >>> ellippi(1,0,1)
         0.0
 
diff --git a/mpmath/functions/functions.py b/mpmath/functions/functions.py
index 7bc8bca..593ac30 100644
--- a/mpmath/functions/functions.py
+++ b/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):
diff --git a/mpmath/functions/zeta.py b/mpmath/functions/zeta.py
index 4960ff9..9aab919 100644
--- a/mpmath/functions/zeta.py
+++ b/mpmath/functions/zeta.py
@@ -953,7 +953,7 @@ def secondzeta(ctx, s, a = 0.015, **kwargs):
         >>> secondzeta(-8)
         -0.67236328125
         >>> secondzeta(-7)
-        +inf
+        inf
 
     **Implementation notes**
 
diff --git a/mpmath/libmp/libmpf.py b/mpmath/libmp/libmpf.py
index 7cabdb2..ddf81f4 100644
--- a/mpmath/libmp/libmpf.py
+++ b/mpmath/libmp/libmpf.py
@@ -1103,28 +1103,39 @@ def to_digits_exp(s, dps, base=10):
     exponent += len(digits) - fixdps - 1
     return sign, digits, exponent
 
-def round_digits(digits, dps, base):
+def round_digits(digits, dps, base, rounding=round_nearest):
     '''
     Returns the rounded digits, and the number of places the decimal point was
     shifted.
+
+    Supports three kinds of rounding: up, down, or nearest.
     '''
 
     assert len(digits) > dps
+    assert rounding in (round_nearest, round_up, round_down)
+
     exponent = 0
 
-    rnd_digs = stddigits[(base//2 + base % 2):base]
+    if rounding == round_down:
+        return digits[:dps], 0
+    elif rounding == round_nearest:
+        rnd_digs = stddigits[(base//2 + base % 2):base]
+    else:
+        rnd_digs = stddigits[:base]
+
+    tie_up = False
 
-    round_up = False
     # The first digit after dps is a 5.
     if digits[dps] == rnd_digs[0]:
         for i in range(dps+1, len(digits)):
             if digits[i] != '0':
-                round_up = True
+                tie_up = True
                 break
         if digits[dps-1] in stddigits[1:base:2]:
-            round_up = True
+            tie_up = True
+
+    if not tie_up:
 
-    if not round_up:
         digits = digits[:dps] + stddigits[int(digits[dps], base) - 1]
 
     # Rounding up kills some instances of "...99999"
@@ -1195,7 +1206,7 @@ def to_str(s, dps, strip_zeros=True, min_fixed=None, max_fixed=None,
             if show_zero_exponent:
                 t += sep + '+0'
             return prefix + t
-        if s == finf: return '+inf'
+        if s == finf: return 'inf'
         if s == fninf: return '-inf'
         if s == fnan: return 'nan'
         raise ValueError
@@ -1375,9 +1386,16 @@ def from_str(x, prec=0, rnd=round_fast, base=0):
     (?P<width>0|[1-9][0-9]*)?
     (?P<thousands_separators>[,_])?
     (?:\.(?P<precision>0|[1-9][0-9]*))?
+    (?P<rounding>[UDZN])?
     (?P<type>[eEfFgG])
 """, re.DOTALL | re.VERBOSE).fullmatch
 
+_GMPY_ROUND_CHAR_DICT = {
+        'U': round_ceiling,
+        'D': round_floor,
+        'Z': round_down,
+        'N': round_nearest
+        }
 
 def calc_padding(nchars, width, align):
     '''
@@ -1415,6 +1433,7 @@ def read_format_spec(format_spec):
         'thousands_separators': '',
         'width': -1,
         'precision': 6,
+        'rounding': round_nearest,
         'type': 'f'
         }
 
@@ -1429,8 +1448,12 @@ def read_format_spec(format_spec):
             or format_dict['thousands_separators']
         format_dict['width'] = int(match['width'] or format_dict['width'])
         format_dict['precision'] = int(match['precision'] or format_dict['precision'])
+        rounding_char = match['rounding']
         format_dict['type'] = match['type'] or format_dict['type']
 
+        if rounding_char is not None:
+            format_dict['rounding'] = _GMPY_ROUND_CHAR_DICT[rounding_char]
+
         if match['zeropad'] and match['fill_char']:
             raise ValueError('Cannot specify both 0-padding and a fill '
                              'character')
@@ -1453,7 +1476,8 @@ def format_fixed(s,
                  sep_range=3,
                  base=10,
                  alternate=False,
-                 no_neg_0=False):
+                 no_neg_0=False,
+                 rounding=round_nearest):
     '''
     Format a number into fixed point.
     Returns the sign character, and the string that represents the number in
@@ -1490,7 +1514,7 @@ def format_fixed(s,
         if no_neg_0:
             sign = '' if sign_spec == '-' else sign_spec
     else:
-        digits, exp_add = round_digits(digits, dps, base)
+        digits, exp_add = round_digits(digits, dps, base, rounding)
         exponent += exp_add
 
         # Here we prepend the corresponding 0s to the digits string, according
@@ -1543,7 +1567,8 @@ def format_scientific(s,
                       base=10,
                       capitalize=False,
                       alternate=False,
-                      no_neg_0=False):
+                      rounding=round_nearest):
+
 
     sep = 'E' if capitalize else 'e'
 
@@ -1555,7 +1580,7 @@ def format_scientific(s,
     if sign != '-' and sign_spec != '-':
         sign = sign_spec
 
-    digits, exp_add = round_digits(digits, dps, base)
+    digits, exp_add = round_digits(digits, dps, base, rounding)
     exponent += exp_add
 
     if strip_zeros:
@@ -1593,6 +1618,13 @@ def format_mpf(num, format_spec):
     strip_last_zero = False
     strip_zeros = False
 
+    if format_dict['rounding'] == round_ceiling:
+        rounding = round_down if num[0] else round_up
+    elif format_dict['rounding'] == round_floor:
+        rounding = round_up if num[0] else round_down
+    else:
+        rounding = format_dict['rounding']
+
     if fmt_type == 'g':
         if not format_dict['alternate']:
             strip_last_zero = True
@@ -1620,7 +1652,8 @@ def format_mpf(num, format_spec):
                 sep_range=3,
                 base=10,
                 alternate=format_dict['alternate'],
-                no_neg_0=format_dict['no_neg_0']
+                no_neg_0=format_dict['no_neg_0'],
+                rounding=rounding
                 )
     else:  # The format type is scientific
         sign, digits = format_scientific(
@@ -1631,7 +1664,7 @@ def format_mpf(num, format_spec):
                 base=10,
                 capitalize=capitalize,
                 alternate=format_dict['alternate'],
-                no_neg_0=format_dict['no_neg_0']
+                rounding=rounding
                 )
 
     nchars = len(digits) + len(sign)


Test patch:
diff --git a/mpmath/tests/test_convert.py b/mpmath/tests/test_convert.py
index 678af08..73c84ca 100644
--- a/mpmath/tests/test_convert.py
+++ b/mpmath/tests/test_convert.py
@@ -41,7 +41,7 @@ def test_basic_string():
     assert str(mpf("-2163048125l")) == '-2163048125.0'
     assert str(mpf("-2163048125L/1088391168")) == '-1.98738118113799'
     assert str(mpf("2163048125/1088391168l")) == '1.98738118113799'
-    assert str(mpf('inf')) == '+inf'
+    assert str(mpf('inf')) == 'inf'
 
     # issue 613
     assert str(mpf('2_5_0_0.0')) == '2500.0'
diff --git a/mpmath/tests/test_format.py b/mpmath/tests/test_format.py
index ecec9b7..4379296 100644
--- a/mpmath/tests/test_format.py
+++ b/mpmath/tests/test_format.py
@@ -559,6 +559,91 @@ def test_mpf_fmt():
         assert f"{mp.mpf('0.24'):=+020.2e}" == '+000000000002.40e-01'
         assert f"{mp.mpf('0.24'):=+020.2g}" == '+0000000000000000.24'
 
+        # Tests for different kinds of rounding
+        num = mp.mpf('-1.23456789999901234567')
+        assert f"{num:=.2Uf}" == "-1.23"
+        assert f"{num:=.2Df}" == "-1.24"
+        assert f"{num:=.2Zf}" == "-1.23"
+        assert f"{num:=.2Nf}" == "-1.23"
+
+        assert f"{num:=.3Uf}" == "-1.234"
+        assert f"{num:=.3Df}" == "-1.235"
+        assert f"{num:=.3Zf}" == "-1.234"
+        assert f"{num:=.3Nf}" == "-1.235"
+
+        assert f"{num:=.10Uf}" == "-1.2345678999"
+        assert f"{num:=.10Df}" == "-1.2345679000"
+        assert f"{num:=.10Zf}" == "-1.2345678999"
+        assert f"{num:=.10Nf}" == "-1.2345679000"
+
+        num = mp.mpf('1.23456789999901234567')
+        assert f"{num:=.2Uf}" == "1.24"
+        assert f"{num:=.2Df}" == "1.23"
+        assert f"{num:=.2Zf}" == "1.23"
+        assert f"{num:=.2Nf}" == "1.23"
+
+        assert f"{num:=.3Uf}" == "1.235"
+        assert f"{num:=.3Df}" == "1.234"
+        assert f"{num:=.3Zf}" == "1.234"
+        assert f"{num:=.3Nf}" == "1.235"
+
+        assert f"{num:=.10Uf}" == "1.2345679000"
+        assert f"{num:=.10Df}" == "1.2345678999"
+        assert f"{num:=.10Zf}" == "1.2345678999"
+        assert f"{num:=.10Nf}" == "1.2345679000"
+
+        num = mp.mpf('-123.456789999901234567')
+        assert f"{num:=.2Ue}" == "-1.23e+02"
+        assert f"{num:=.2De}" == "-1.24e+02"
+        assert f"{num:=.2Ze}" == "-1.23e+02"
+        assert f"{num:=.2Ne}" == "-1.23e+02"
+
+        assert f"{num:=.3Ue}" == "-1.234e+02"
+        assert f"{num:=.3De}" == "-1.235e+02"
+        assert f"{num:=.3Ze}" == "-1.234e+02"
+        assert f"{num:=.3Ne}" == "-1.235e+02"
+
+        assert f"{num:=.10Ue}" == "-1.2345678999e+02"
+        assert f"{num:=.10De}" == "-1.2345679000e+02"
+        assert f"{num:=.10Ze}" == "-1.2345678999e+02"
+        assert f"{num:=.10Ne}" == "-1.2345679000e+02"
+
+        num = mp.mpf('123456.789999901234567')
+        assert f"{num:=.2Ue}" == "1.24e+05"
+        assert f"{num:=.2De}" == "1.23e+05"
+        assert f"{num:=.2Ze}" == "1.23e+05"
+        assert f"{num:=.2Ne}" == "1.23e+05"
+
+        assert f"{num:=.3Ue}" == "1.235e+05"
+        assert f"{num:=.3De}" == "1.234e+05"
+        assert f"{num:=.3Ze}" == "1.234e+05"
+        assert f"{num:=.3Ne}" == "1.235e+05"
+
+        assert f"{num:=.10Ue}" == "1.2345679000e+05"
+        assert f"{num:=.10De}" == "1.2345678999e+05"
+        assert f"{num:=.10Ze}" == "1.2345678999e+05"
+        assert f"{num:=.10Ne}" == "1.2345679000e+05"
+
+        assert f"{mp.mpf('123.456'):.2Ug}" == "1.3e+02"
+        assert f"{mp.mpf('123.456'):.2Dg}" == "1.2e+02"
+        assert f"{mp.mpf('123.456'):.2Zg}" == "1.2e+02"
+        assert f"{mp.mpf('123.456'):.2Ng}" == "1.2e+02"
+
+        assert f"{mp.mpf('-123.456'):.2Ug}" == "-1.2e+02"
+        assert f"{mp.mpf('-123.456'):.2Dg}" == "-1.3e+02"
+        assert f"{mp.mpf('-123.456'):.2Zg}" == "-1.2e+02"
+        assert f"{mp.mpf('-123.456'):.2Ng}" == "-1.2e+02"
+
+        assert f"{mp.mpf('123.456'):.5Ug}" == "123.46"
+        assert f"{mp.mpf('123.456'):.5Dg}" == "123.45"
+        assert f"{mp.mpf('123.456'):.5Zg}" == "123.45"
+        assert f"{mp.mpf('123.456'):.5Ng}" == "123.46"
+
+        assert f"{mp.mpf('-123.456'):.5Ug}" == "-123.45"
+        assert f"{mp.mpf('-123.456'):.5Dg}" == "-123.46"
+        assert f"{mp.mpf('-123.456'):.5Zg}" == "-123.45"
+        assert f"{mp.mpf('-123.456'):.5Ng}" == "-123.46"
+
 
 def test_errors():
     with pytest.raises(ValueError):
diff --git a/mpmath/tests/test_functions.py b/mpmath/tests/test_functions.py
index aca69f0..01914f7 100644
--- a/mpmath/tests/test_functions.py
+++ b/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