Repository: sympy/sympy
Pull request: Add automatic evaluation for compositions of inverse hyperbolic functions
#### References to other Issues or PRs
Closes #17421. The more difficult issue with evaluation near branch cuts is not addressed, but there is another issue that covers that (#6137).

#### Brief description of what is fixed or changed
Add automatic evaluation of `atanh(tanh(x)), acosh(cosh(x)), asinh(sinh(x))` for suitable `x`.

#### Release Notes
<!-- BEGIN RELEASE NOTES -->
* functions
  * atanh(tanh(x)), acosh(cosh(x)) and asinh(sinh(x)) are automatically evaluated
<!-- END RELEASE NOTES -->


Problem statement:
Incorrect evaluation of acos(-I + acosh(cosh(cosh(1) + I)))
`-I + acosh(cosh(cosh(1) + I))` should be equivalent to `cosh(1)`, and this can be confirmed:

```
>>> N(-I + acosh(cosh(cosh(1) + I)))
1.54308063481524 + 0.e-20*I
>>> N(cosh(1))
1.54308063481524
```

Taking `acos`, the sign gets flipped:

```
>>> a = -I + acosh(cosh(cosh(1) + I))
>>> N(acos(a))
3.60377845017333e-22 - 1.0*I
>>> N(acos(cosh(1)))
1.0*I
```

Trying with explicit evaluation:

```
>>> a = -I + acosh(cosh(cosh(1) + I))
>>> N(acos(a))
3.60377845017333e-22 - 1.0*I
>>> N(a)
1.54308063481524 + 0.e-20*I
>>> acos(_)
1.15320910405546e-20 - 1.0*I
>>> acos(1.54308063481524 + 0.e-20*I)
0.999999999999997*I
>>> acos(1.54308063481524)
0.999999999999997*I
```

Hints from issue discussion:
Actually, I think this is more of a mpmath issue- will close unless anyone feels otherwise.
What makes you think it is an issue with mpmath rather than sympy?
Never mind. I guess the issue is that any positive imaginary part no matter how small is going to cause the sign to flip from `I` to `-I`.
perhaps related to #6137? See this [comment](https://github.com/sympy/sympy/issues/6137#issuecomment-37004430).
Maybe evalf should have some understanding of branch cuts at least so that it doesn't report excess precision near to them.
Or maybe just any way to simplify `acosh(cosh(x))`? As far as I know there's nothing in `trigsimp` (or anywhere) that deals with hyperbolic inverse compositions.

Something like this:

```
def acosh_simp(x):
    # simplify acosh(cosh(x))
    if x.is_real:
        if x.is_nonnegative:
            return x
        elif x.is_negative:
            return -x
    else:
        r, i = x.as_real_imag()
        if r.is_number and i.is_number:
            n1 = floor(i / pi)
            if n1 % 2 == 0:
                m = r + I * (i - pi*n1)
                if r.is_nonnegative:
                    return m
                else:
                    return -m
            else:
                m = r + I * (i - pi*n1 - pi)
                if r.is_nonpositive:
                    return -m
                else:
                    return m
```

Commit message: fix merge problem

Patch to review:
diff --git a/sympy/core/tests/test_evalf.py b/sympy/core/tests/test_evalf.py
index 52aa0a56eb35..63a034768552 100644
--- a/sympy/core/tests/test_evalf.py
+++ b/sympy/core/tests/test_evalf.py
@@ -1,7 +1,7 @@
 from sympy import (Abs, Add, atan, ceiling, cos, E, Eq, exp, factor,
     factorial, fibonacci, floor, Function, GoldenRatio, I, Integral,
     integrate, log, Mul, N, oo, pi, Pow, product, Product,
-    Rational, S, Sum, simplify, sin, sqrt, sstr, sympify, Symbol, Max, nfloat)
+    Rational, S, Sum, simplify, sin, sqrt, sstr, sympify, Symbol, Max, nfloat, cosh, acosh, acos)
 from sympy.core.numbers import comp
 from sympy.core.evalf import (complex_accuracy, PrecisionExhausted,
     scaled_zero, get_integer_part, as_mpmath, evalf)
@@ -567,3 +567,7 @@ def test_issue_13425():
     assert N('2**.5', 30) == N('sqrt(2)', 30)
     assert N('x - x', 30) == 0
     assert abs((N('pi*.1', 22)*10 - pi).n()) < 1e-22
+
+
+def test_issue_17421():
+    assert N(acos(-I + acosh(cosh(cosh(1) + I)))) == 1.0*I
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index 5a5ee285715b..4cd84f026f7f 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -491,11 +491,11 @@ def _eval_rewrite_as_Pow(self, arg, **kwargs):
                 return Pow(logs[0].args[0], arg.coeff(logs[0]))
 
 
-def _match_real_imag(expr):
+def match_real_imag(expr):
     """
     Try to match expr with a + b*I for real a and b.
 
-    ``_match_real_imag`` returns a tuple containing the real and imaginary
+    ``match_real_imag`` returns a tuple containing the real and imaginary
     parts of expr or (None, None) if direct matching is not possible. Contrary
     to ``re()``, ``im()``, ``as_real_imag()``, this helper won't force things
     by returning expressions themselves containing ``re()`` or ``im()`` and it
@@ -607,7 +607,7 @@ def eval(cls, arg, base=None):
         if isinstance(arg, exp) and arg.args[0].is_extended_real:
             return arg.args[0]
         elif isinstance(arg, exp) and arg.args[0].is_number:
-            r_, i_ = _match_real_imag(arg.args[0])
+            r_, i_ = match_real_imag(arg.args[0])
             if i_ and i_.is_comparable:
                 i_ %= 2*S.Pi
                 if i_ > S.Pi:
diff --git a/sympy/functions/elementary/hyperbolic.py b/sympy/functions/elementary/hyperbolic.py
index dde790eb5562..fa315d480312 100644
--- a/sympy/functions/elementary/hyperbolic.py
+++ b/sympy/functions/elementary/hyperbolic.py
@@ -1,11 +1,12 @@
 from __future__ import print_function, division
 
-from sympy.core import S, sympify, cacheit
+from sympy.core import S, sympify, cacheit, pi, I
 from sympy.core.add import Add
 from sympy.core.function import Function, ArgumentIndexError, _coeff_isneg
 from sympy.functions.combinatorial.factorials import factorial, RisingFactorial
-from sympy.functions.elementary.exponential import exp, log
+from sympy.functions.elementary.exponential import exp, log, match_real_imag
 from sympy.functions.elementary.miscellaneous import sqrt
+from sympy.functions.elementary.integers import floor
 
 
 def _rewrite_hyperbolics_as_exp(expr):
@@ -926,6 +927,21 @@ def eval(cls, arg):
                 if _coeff_isneg(arg):
                     return -cls(-arg)
 
+        if isinstance(arg, sinh):
+            x = arg.args[0]
+            if x.is_number:
+                if x.is_real:
+                    return x
+                r, i = match_real_imag(x)
+                if r is not None and i is not None:
+                    f = floor((i + pi/2)/pi)
+                    if f.is_number and f.is_integer:
+                        m = r + I*i - I*pi*f
+                        if f % 2 == 0:
+                            return m
+                        else:
+                            return -m
+
     @staticmethod
     @cacheit
     def taylor_term(n, x, *previous_terms):
@@ -1033,6 +1049,29 @@ def eval(cls, arg):
         if arg == -S.ImaginaryUnit*S.Infinity:
             return S.Infinity - S.ImaginaryUnit*S.Pi/2
 
+        if isinstance(arg, cosh):
+            x = arg.args[0]
+            if x.is_number:
+                if x.is_real:
+                    from sympy.functions.elementary.complexes import Abs
+                    return Abs(x)
+                r, i = match_real_imag(x)
+                if r is not None and i is not None:
+                    f = floor(i/pi)
+                    if f.is_number and f.is_integer:
+                        m = r + I*i - I*pi*f
+                        if f % 2 == 0:
+                            if r.is_nonnegative:
+                                return m
+                            elif r.is_negative:
+                                return -m
+                        else:
+                            m -= I*pi
+                            if r.is_nonpositive:
+                                return -m
+                            elif r.is_positive:
+                                return m
+
     @staticmethod
     @cacheit
     def taylor_term(n, x, *previous_terms):
@@ -1121,6 +1160,20 @@ def eval(cls, arg):
                 if _coeff_isneg(arg):
                     return -cls(-arg)
 
+        if isinstance(arg, tanh):
+            x = arg.args[0]
+            if x.is_number:
+                if x.is_real:
+                    return x
+                r, i = match_real_imag(x)
+                if r is not None and i is not None:
+                    f = floor(2*i/pi)
+                    if f.is_number and f.is_integer:
+                        if f % 2 == 0:
+                            return r + I*i - I*f*pi/2
+                        else:
+                            return r + I*i - I*(f + 1)*pi/2
+
     @staticmethod
     @cacheit
     def taylor_term(n, x, *previous_terms):
diff --git a/sympy/functions/elementary/tests/test_exponential.py b/sympy/functions/elementary/tests/test_exponential.py
index c7794d845e31..9a18a0988a2f 100644
--- a/sympy/functions/elementary/tests/test_exponential.py
+++ b/sympy/functions/elementary/tests/test_exponential.py
@@ -3,7 +3,7 @@
     LambertW, sqrt, Rational, expand_log, S, sign, conjugate, refine,
     sin, cos, sinh, cosh, tanh, exp_polar, re, Function, simplify,
     AccumBounds, MatrixSymbol, Pow, gcd, Sum, Product)
-from sympy.functions.elementary.exponential import _match_real_imag
+from sympy.functions.elementary.exponential import match_real_imag
 from sympy.abc import x, y, z
 from sympy.core.expr import unchanged
 from sympy.core.function import ArgumentIndexError
@@ -221,16 +221,16 @@ def test_log_values():
 def test_match_real_imag():
     x, y = symbols('x,y', real=True)
     i = Symbol('i', imaginary=True)
-    assert _match_real_imag(S.One) == (1, 0)
-    assert _match_real_imag(I) == (0, 1)
-    assert _match_real_imag(3 - 5*I) == (3, -5)
-    assert _match_real_imag(-sqrt(3) + S.Half*I) == (-sqrt(3), S.Half)
-    assert _match_real_imag(x + y*I) == (x, y)
-    assert _match_real_imag(x*I + y*I) == (0, x + y)
-    assert _match_real_imag((x + y)*I) == (0, x + y)
-    assert _match_real_imag(-S(2)/3*i*I) == (None, None)
-    assert _match_real_imag(1 - 2*i) == (None, None)
-    assert _match_real_imag(sqrt(2)*(3 - 5*I)) == (None, None)
+    assert match_real_imag(S.One) == (1, 0)
+    assert match_real_imag(I) == (0, 1)
+    assert match_real_imag(3 - 5*I) == (3, -5)
+    assert match_real_imag(-sqrt(3) + S.Half*I) == (-sqrt(3), S.Half)
+    assert match_real_imag(x + y*I) == (x, y)
+    assert match_real_imag(x*I + y*I) == (0, x + y)
+    assert match_real_imag((x + y)*I) == (0, x + y)
+    assert match_real_imag(-S(2)/3*i*I) == (None, None)
+    assert match_real_imag(1 - 2*i) == (None, None)
+    assert match_real_imag(sqrt(2)*(3 - 5*I)) == (None, None)
 
 
 def test_log_exact():
diff --git a/sympy/functions/elementary/tests/test_hyperbolic.py b/sympy/functions/elementary/tests/test_hyperbolic.py
index dd21df540d65..237533224ea2 100644
--- a/sympy/functions/elementary/tests/test_hyperbolic.py
+++ b/sympy/functions/elementary/tests/test_hyperbolic.py
@@ -1,7 +1,7 @@
 from sympy import (symbols, Symbol, sinh, nan, oo, zoo, pi, asinh, acosh, log,
     sqrt, coth, I, cot, E, tanh, tan, cosh, cos, S, sin, Rational, atanh, acoth,
     Integer, O, exp, sech, sec, csch, asech, acsch, acos, asin, expand_mul,
-    AccumBounds, im, re)
+    AccumBounds, im, re, Abs)
 
 from sympy.core.expr import unchanged
 from sympy.core.function import ArgumentIndexError
@@ -508,6 +508,20 @@ def test_asinh():
     # Symmetry
     assert asinh(-S.Half) == -asinh(S.Half)
 
+    # inverse composition
+    assert unchanged(asinh, sinh(Symbol('v1')))
+
+    assert asinh(sinh(0, evaluate=False)) == 0
+    assert asinh(sinh(-3, evaluate=False)) == -3
+    assert asinh(sinh(2, evaluate=False)) == 2
+    assert asinh(sinh(I, evaluate=False)) == I
+    assert asinh(sinh(-I, evaluate=False)) == -I
+    assert asinh(sinh(5*I, evaluate=False)) == -2*I*pi + 5*I
+    assert asinh(sinh(15 + 11*I)) == 15 - 4*I*pi + 11*I
+    assert asinh(sinh(-73 + 97*I)) == 73 - 97*I + 31*I*pi
+    assert asinh(sinh(-7 - 23*I)) == 7 - 7*I*pi + 23*I
+    assert asinh(sinh(13 - 3*I)) == -13 - I*pi + 3*I
+
 
 def test_asinh_rewrite():
     x = Symbol('x')
@@ -570,6 +584,22 @@ def test_acosh():
     assert str(acosh(5*I).n(6)) == '2.31244 + 1.5708*I'
     assert str(acosh(-5*I).n(6)) == '2.31244 - 1.5708*I'
 
+    # inverse composition
+    assert unchanged(acosh, Symbol('v1'))
+
+    assert acosh(cosh(-3, evaluate=False)) == 3
+    assert acosh(cosh(3, evaluate=False)) == 3
+    assert acosh(cosh(0, evaluate=False)) == 0
+    assert acosh(cosh(I, evaluate=False)) == I
+    assert acosh(cosh(-I, evaluate=False)) == I
+    assert acosh(cosh(7*I, evaluate=False)) == -2*I*pi + 7*I
+    assert acosh(cosh(1 + I)) == 1 + I
+    assert acosh(cosh(3 - 3*I)) == 3 - 3*I
+    assert acosh(cosh(-3 + 2*I)) == 3 - 2*I
+    assert acosh(cosh(-5 - 17*I)) == 5 - 6*I*pi + 17*I
+    assert acosh(cosh(-21 + 11*I)) == 21 - 11*I + 4*I*pi
+    assert acosh(cosh(cosh(1) + I)) == cosh(1) + I
+
 
 def test_acosh_rewrite():
     x = Symbol('x')
@@ -780,6 +810,24 @@ def test_atanh():
     # Symmetry
     assert atanh(-S.Half) == -atanh(S.Half)
 
+    # inverse composition
+    assert unchanged(atanh, tanh(Symbol('v1')))
+
+    assert atanh(tanh(-5, evaluate=False)) == -5
+    assert atanh(tanh(0, evaluate=False)) == 0
+    assert atanh(tanh(7, evaluate=False)) == 7
+    assert atanh(tanh(I, evaluate=False)) == I
+    assert atanh(tanh(-I, evaluate=False)) == -I
+    assert atanh(tanh(-11*I, evaluate=False)) == -11*I + 4*I*pi
+    assert atanh(tanh(3 + I)) == 3 + I
+    assert atanh(tanh(4 + 5*I)) == 4 - 2*I*pi + 5*I
+    assert atanh(tanh(pi/2)) == pi/2
+    assert atanh(tanh(pi)) == pi
+    assert atanh(tanh(-3 + 7*I)) == -3 - 2*I*pi + 7*I
+    assert atanh(tanh(9 - 2*I/3)) == 9 - 2*I/3
+    assert atanh(tanh(-32 - 123*I)) == -32 - 123*I + 39*I*pi
+
+
 def test_atanh_rewrite():
     x = Symbol('x')
     assert atanh(x).rewrite(log) == (log(1 + x) - log(1 - x)) / 2


Code review comments:
diff_hunk: @@ -926,6 +927,21 @@ def eval(cls, arg):
                 if _coeff_isneg(arg):
                     return -cls(-arg)
 
+        if isinstance(arg, sinh):
+            x = arg.args[0]
+            if x.is_number:
path: sympy/functions/elementary/hyperbolic.py
text: @user1:
You can reduce the indentation here by combining the conditions:
```python
if isinstance(arg, sinh) and arg.args[0].is_number:
    ...
```
Or:
```python
if isinstance(arg, sinh):
    x = arg.args[0]
    if not x.is_number:
        return
    ...
```

diff_hunk: @@ -926,6 +927,21 @@ def eval(cls, arg):
                 if _coeff_isneg(arg):
                     return -cls(-arg)
 
+        if isinstance(arg, sinh):
+            x = arg.args[0]
+            if x.is_number:
+                if x.is_real:
+                    return x
+                r, i = match_real_imag(x)
+                if r is not None and i is not None:
+                    f = floor((i + pi/2)/pi)
+                    if f.is_number and f.is_integer:
+                        m = r + I*i - I*pi*f
+                        if f % 2 == 0:
+                            return m
path: sympy/functions/elementary/hyperbolic.py
text: @user1:
We already know that f.is_number because f comes from x and we checked that x.is_number. Likewise if f is floor(...) then f.is_integer *should* always be True (although it might sometimes return None). It is still possible though that f.is_number and f.is_integer but we can't calculate f%2. I would say that the check would be better as:
```python
even = f.is_even
if even is True:
    return m
elif even is False:
    return -m
# Note the 3rd case even is None where we don't evaluate
```