Repository: sympy/sympy

Issue:
Inconsistent block_collapse of BlockDiagMatrix's
``` Python
>>> from sympy import ImmutableMatrix, BlockDiagMatrix, block_collapse
>>> from sympy.abc import x, y, z
>>> bm1 = BlockDiagMatrix(ImmutableMatrix([x]), ImmutableMatrix([y]))
>>> bm2 = BlockDiagMatrix(ImmutableMatrix([y]), ImmutableMatrix([z]))
>>> 
>>> type(block_collapse(bm1))
sympy.matrices.expressions.blockmatrix.BlockDiagMatrix
>>> type(block_collapse(bm1 + bm2))
sympy.matrices.expressions.blockmatrix.BlockDiagMatrix
>>> 
>>> type(bm1.T)
sympy.matrices.expressions.blockmatrix.BlockMatrix
>>> type(block_collapse(bm1 - bm2))
sympy.matrices.expressions.blockmatrix.BlockMatrix
```


Discussion:

Patch:
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index 32fd694088..4ac1397a3b 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -4,8 +4,7 @@
 from sympy.core.add import Add
 from sympy.core.cache import cacheit
 from sympy.core.compatibility import range
-from sympy.core.function import (Function, ArgumentIndexError, _coeff_isneg,
-        expand_mul)
+from sympy.core.function import Function, ArgumentIndexError, _coeff_isneg
 from sympy.core.logic import fuzzy_not
 from sympy.core.mul import Mul
 from sympy.core.numbers import Integer
@@ -13,7 +12,6 @@
 from sympy.core.singleton import S
 from sympy.core.symbol import Wild, Dummy
 from sympy.functions.combinatorial.factorials import factorial
-from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.ntheory import multiplicity, perfect_power
 
 # NOTE IMPORTANT
@@ -419,7 +417,7 @@ def _eval_is_extended_positive(self):
     def _eval_nseries(self, x, n, logx):
         # NOTE Please see the comment at the beginning of this file, labelled
         #      IMPORTANT.
-        from sympy import limit, oo, Order, powsimp, Wild, expand_complex
+        from sympy import limit, oo, Order, powsimp
         arg = self.args[0]
         arg_series = arg._eval_nseries(x, n=n, logx=logx)
         if arg_series.is_Order:
@@ -434,12 +432,7 @@ def _eval_nseries(self, x, n, logx):
         r = exp(arg0)*exp_series.subs(t, arg_series - arg0)
         r += Order(o.expr.subs(t, (arg_series - arg0)), x)
         r = r.expand()
-        r = powsimp(r, deep=True, combine='exp')
-        # powsimp may introduce unexpanded (-1)**Rational; see PR #17201
-        simplerat = lambda x: x.is_Rational and x.q in [3, 4, 6]
-        w = Wild('w', properties=[simplerat])
-        r = r.replace((-1)**w, expand_complex((-1)**w))
-        return r
+        return powsimp(r, deep=True, combine='exp')
 
     def _taylor(self, x, n):
         from sympy import Order
@@ -494,7 +487,6 @@ def _eval_rewrite_as_Pow(self, arg, **kwargs):
 class log(Function):
     r"""
     The natural logarithm function `\ln(x)` or `\log(x)`.
-
     Logarithms are taken with the natural base, `e`. To get
     a logarithm of a different base ``b``, use ``log(x, b)``,
     which is essentially short-hand for ``log(x)/log(b)``.
@@ -502,19 +494,16 @@ class log(Function):
     Examples
     ========
 
-    >>> from sympy import log, sqrt, S, I
+    >>> from sympy import log, S
     >>> log(8, 2)
     3
     >>> log(S(8)/3, 2)
     -log(3)/log(2) + 3
-    >>> log(-1 + I*sqrt(3))
-    log(2) + 2*I*pi/3
 
     See Also
     ========
 
     exp
-
     """
 
     def fdiff(self, argindex=1):
@@ -537,7 +526,6 @@ def eval(cls, arg, base=None):
         from sympy import unpolarify
         from sympy.calculus import AccumBounds
         from sympy.sets.setexpr import SetExpr
-        from sympy.functions.elementary.complexes import Abs, re, im
 
         arg = sympify(arg)
 
@@ -590,9 +578,8 @@ def eval(cls, arg, base=None):
             return arg._eval_func(cls)
 
         if arg.is_number:
-            I = S.ImaginaryUnit
             if arg.is_negative:
-                return S.Pi * I + cls(-arg)
+                return S.Pi * S.ImaginaryUnit + cls(-arg)
             elif arg is S.ComplexInfinity:
                 return S.ComplexInfinity
             elif arg is S.Exp1:
@@ -613,63 +600,6 @@ def eval(cls, arg, base=None):
                     else:
                         return -S.Pi * S.ImaginaryUnit * S.Half + cls(-coeff)
 
-        if arg.is_number and arg.is_algebraic:
-            I = S.ImaginaryUnit
-            # Match arg = coeff*(r_ + i_*I) with coeff>0, r_ and i_ real.
-            coeff, arg_ = arg.as_independent(I, as_Add=False)
-            if coeff.is_negative:
-                coeff *= -1
-                arg_ *= -1
-            arg_ = expand_mul(arg_, deep=False)
-            r_, i_ = arg_.as_independent(I, as_Add=True)
-            i_ = i_.as_coefficient(I)
-            if coeff.is_real and i_ and i_.is_real and r_.is_real:
-                if r_.is_zero:
-                    if i_.is_positive:
-                        return S.Pi * I * S.Half + cls(coeff * i_)
-                    elif i_.is_negative:
-                        return -S.Pi * I * S.Half + cls(coeff * -i_)
-                    elif i_.is_zero:
-                        return zoo
-                else:
-                    from sympy.simplify import ratsimp
-                    # Check for arguments involving rational multiples of pi
-                    t = (i_/r_).cancel()
-                    atan_table = {
-                        # first quadrant only
-                        sqrt(3): S.Pi/3,
-                        1: S.Pi/4,
-                        sqrt(5 - 2*sqrt(5)): S.Pi/5,
-                        sqrt(2)*sqrt(5 - sqrt(5))/(1 + sqrt(5)): S.Pi/5,
-                        sqrt(5 + 2*sqrt(5)): 2*S.Pi/5,
-                        sqrt(2)*sqrt(sqrt(5) + 5)/(-1 + sqrt(5)): 2*S.Pi/5,
-                        sqrt(3)/3: S.Pi/6,
-                        sqrt(2) - 1: S.Pi/8,
-                        sqrt(2 - sqrt(2))/sqrt(sqrt(2) + 2): S.Pi/8,
-                        sqrt(2) + 1: 3*S.Pi/8,
-                        sqrt(sqrt(2) + 2)/sqrt(2 - sqrt(2)): 3*S.Pi/8,
-                        sqrt(1 - 2*sqrt(5)/5): S.Pi/10,
-                        (-sqrt(2) + sqrt(10))/(2*sqrt(sqrt(5) + 5)): S.Pi/10,
-                        sqrt(1 + 2*sqrt(5)/5): 3*S.Pi/10,
-                        (sqrt(2) + sqrt(10))/(2*sqrt(5 - sqrt(5))): 3*S.Pi/10,
-                        2 - sqrt(3): S.Pi/12,
-                        (-1 + sqrt(3))/(1 + sqrt(3)): S.Pi/12,
-                        2 + sqrt(3): 5*S.Pi/12,
-                        (1 + sqrt(3))/(-1 + sqrt(3)): 5*S.Pi/12
-                    }
-                    if t in atan_table:
-                        modulus = ratsimp(coeff * Abs(arg_))
-                        if r_.is_positive:
-                            return cls(modulus) + I * atan_table[t]
-                        else:
-                            return cls(modulus) + I * (atan_table[t] - S.Pi)
-                    elif -t in atan_table:
-                        modulus = ratsimp(coeff * Abs(arg_))
-                        if r_.is_positive:
-                            return cls(modulus) + I * (-atan_table[-t])
-                        else:
-                            return cls(modulus) + I * (S.Pi - atan_table[-t])
-
     def as_base_exp(self):
         """
         Returns this function in the form (base, exponent).
diff --git a/sympy/matrices/expressions/blockmatrix.py b/sympy/matrices/expressions/blockmatrix.py
index 4989278284..6f648c3955 100644
--- a/sympy/matrices/expressions/blockmatrix.py
+++ b/sympy/matrices/expressions/blockmatrix.py
@@ -327,6 +327,9 @@ def colblocksizes(self):
     def _eval_inverse(self, expand='ignored'):
         return BlockDiagMatrix(*[mat.inverse() for mat in self.args])
 
+    def _eval_transpose(self):
+        return BlockDiagMatrix(*[mat.transpose() for mat in self.args])
+
     def _blockmul(self, other):
         if (isinstance(other, BlockDiagMatrix) and
                 self.colblocksizes == other.rowblocksizes):
@@ -419,11 +422,21 @@ def bc_block_plus_ident(expr):
 def bc_dist(expr):
     """ Turn  a*[X, Y] into [a*X, a*Y] """
     factor, mat = expr.as_coeff_mmul()
-    if factor != 1 and isinstance(unpack(mat), BlockMatrix):
-        B = unpack(mat).blocks
-        return BlockMatrix([[factor * B[i, j] for j in range(B.cols)]
-                                              for i in range(B.rows)])
-    return expr
+    if factor == 1:
+        return expr
+
+    unpacked = unpack(mat)
+
+    if isinstance(unpacked, BlockDiagMatrix):
+        B = unpacked.diag
+        new_B = [factor * mat for mat in B]
+        return BlockDiagMatrix(*new_B)
+    elif isinstance(unpacked, BlockMatrix):
+        B = unpacked.blocks
+        new_B = [
+            [factor * B[i, j] for j in range(B.cols)] for i in range(B.rows)]
+        return BlockMatrix(new_B)
+    return unpacked
 
 
 def bc_matmul(expr):
@@ -452,10 +465,14 @@ def bc_matmul(expr):
     return MatMul(factor, *matrices).doit()
 
 def bc_transpose(expr):
-    return BlockMatrix(block_collapse(expr.arg).blocks.applyfunc(transpose).T)
+    collapse = block_collapse(expr.arg)
+    return collapse._eval_transpose()
 
 
 def bc_inverse(expr):
+    if isinstance(expr.arg, BlockDiagMatrix):
+        return expr._eval_inverse()
+
     expr2 = blockinverse_1x1(expr)
     if expr != expr2:
         return expr2
diff --git a/sympy/printing/mathematica.py b/sympy/printing/mathematica.py
index ff33f738d0..0ffb0f1377 100644
--- a/sympy/printing/mathematica.py
+++ b/sympy/printing/mathematica.py
@@ -205,9 +205,7 @@ def _print_GoldenRatio(self, expr):
         return 'GoldenRatio'
 
     def _print_TribonacciConstant(self, expr):
-        expanded = expr.expand(func=True)
-        PREC = precedence(expr)
-        return self.parenthesize(expanded, PREC)
+        return self.doprint(expr._eval_expand_func())
 
     def _print_EulerGamma(self, expr):
         return 'EulerGamma'


Test patch:
diff --git a/sympy/functions/elementary/tests/test_exponential.py b/sympy/functions/elementary/tests/test_exponential.py
index b1f094223d..55ee1a4b5e 100644
--- a/sympy/functions/elementary/tests/test_exponential.py
+++ b/sympy/functions/elementary/tests/test_exponential.py
@@ -1,8 +1,8 @@
 from sympy import (
     symbols, log, ln, Float, nan, oo, zoo, I, pi, E, exp, Symbol,
     LambertW, sqrt, Rational, expand_log, S, sign, conjugate, refine,
-    sin, cos, sinh, cosh, tanh, exp_polar, re, Function, simplify,
-    AccumBounds, MatrixSymbol, Pow, gcd)
+    sin, cos, sinh, cosh, tanh, exp_polar, re, simplify,
+    AccumBounds, MatrixSymbol, Pow)
 from sympy.abc import x, y, z
 from sympy.core.expr import unchanged
 from sympy.core.function import ArgumentIndexError
@@ -217,48 +217,6 @@ def test_log_values():
     assert log(2*3**2).func is log
 
 
-def test_log_exact():
-    # check for pi/2, pi/3, pi/4, pi/6, pi/8, pi/12; pi/5, pi/10:
-    for n in range(-23, 24):
-        if gcd(n, 24) != 1:
-            assert log(exp(n*I*pi/24).rewrite(sqrt)) == n*I*pi/24
-        for n in range(-9, 10):
-            assert log(exp(n*I*pi/10).rewrite(sqrt)) == n*I*pi/10
-
-    assert log(S.Half - I*sqrt(3)/2) == -I*pi/3
-    assert log(-S.Half + I*sqrt(3)/2) == I*2*pi/3
-    assert log(-sqrt(2)/2 - I*sqrt(2)/2) == -I*3*pi/4
-    assert log(-sqrt(3)/2 - I*S.Half) == -I*5*pi/6
-
-    assert log(-S(1)/4 + sqrt(5)/4 - I*sqrt(sqrt(5)/8 + S(5)/8)) == -I*2*pi/5
-    assert log(sqrt(S(5)/8 - sqrt(5)/8) + I*(S(1)/4 + sqrt(5)/4)) == I*3*pi/10
-    assert log(-sqrt(sqrt(2)/4 + S(1)/2) + I*sqrt(S(1)/2 - sqrt(2)/4)) == I*7*pi/8
-    assert log(-sqrt(6)/4 - sqrt(2)/4 + I*(-sqrt(6)/4 + sqrt(2)/4)) == -I*11*pi/12
-
-    assert log(-1 + I*sqrt(3)) == log(2) + I*2*pi/3
-    assert log(5 + 5*I) == log(5*sqrt(2)) + I*pi/4
-    assert log(sqrt(-12)) == log(2*sqrt(3)) + I*pi/2
-    assert log(-sqrt(6) + sqrt(2) - I*sqrt(6) - I*sqrt(2)) == log(4) - I*7*pi/12
-    assert log(-sqrt(6-3*sqrt(2)) - I*sqrt(6+3*sqrt(2))) == log(2*sqrt(3)) - 5*I*pi/8
-    assert log(1 + I*sqrt(2-sqrt(2))/sqrt(2+sqrt(2))) == log(2/sqrt(sqrt(2) + 2)) + I*pi/8
-    assert log(cos(7*pi/12) + I*sin(7*pi/12)) == 7*I*pi/12
-    assert log(cos(6*pi/5) + I*sin(6*pi/5)) == -4*I*pi/5
-
-    assert log(5*(1 + I)/sqrt(2)) == log(5) + I*pi/4
-    assert log(sqrt(2)*(-sqrt(3) + 1 - sqrt(3)*I - I)) == log(4) - I*7*pi/12
-    assert log(-sqrt(2)*(1 - I*sqrt(3))) == log(2*sqrt(2)) + 2*I*pi/3
-    assert log(sqrt(3)*I*(-sqrt(6 - 3*sqrt(2)) - I*sqrt(3*sqrt(2) + 6))) == log(6) - I*pi/8
-
-    zero = (1 + sqrt(2))**2 - 3 - 2*sqrt(2)
-    assert log(zero - I*sqrt(3)) == log(sqrt(3)) - I*pi/2
-    assert unchanged(log, zero + I*zero) or log(zero + zero*I) == zoo
-
-    # bail quickly if no obvious simplification is possible:
-    assert unchanged(log, (sqrt(2)-1/sqrt(sqrt(3)+I))**1000)
-    # beware of non-real coefficients
-    assert unchanged(log, sqrt(2-sqrt(5))*(1 + I))
-
-
 def test_log_base():
     assert log(1, 2) == 0
     assert log(2, 2) == 1
diff --git a/sympy/matrices/expressions/tests/test_blockmatrix.py b/sympy/matrices/expressions/tests/test_blockmatrix.py
index bc162b2824..5de357ce2c 100644
--- a/sympy/matrices/expressions/tests/test_blockmatrix.py
+++ b/sympy/matrices/expressions/tests/test_blockmatrix.py
@@ -1,6 +1,7 @@
-from sympy.matrices.expressions.blockmatrix import (block_collapse, bc_matmul,
-        bc_block_plus_ident, BlockDiagMatrix, BlockMatrix, bc_dist, bc_matadd,
-        bc_transpose, blockcut, reblock_2x2, deblock)
+from sympy.matrices.expressions.blockmatrix import (
+    block_collapse, bc_matmul, bc_block_plus_ident, BlockDiagMatrix,
+    BlockMatrix, bc_dist, bc_matadd, bc_transpose, bc_inverse,
+    blockcut, reblock_2x2, deblock)
 from sympy.matrices.expressions import (MatrixSymbol, Identity,
         Inverse, trace, Transpose, det)
 from sympy.matrices import Matrix, ImmutableMatrix, ones
@@ -203,3 +204,14 @@ def test_deblock():
                     for i in range(4)])
 
     assert deblock(reblock_2x2(B)) == B
+
+def test_block_collapse_type():
+    bm1 = BlockDiagMatrix(ImmutableMatrix([1]), ImmutableMatrix([2]))
+    bm2 = BlockDiagMatrix(ImmutableMatrix([3]), ImmutableMatrix([4]))
+
+    assert bm1.T.__class__ == BlockDiagMatrix
+    assert block_collapse(bm1 - bm2).__class__ == BlockDiagMatrix
+    assert block_collapse(Inverse(bm1)).__class__ == BlockDiagMatrix
+    assert block_collapse(Transpose(bm1)).__class__ == BlockDiagMatrix
+    assert bc_transpose(Transpose(bm1)).__class__ == BlockDiagMatrix
+    assert bc_inverse(Inverse(bm1)).__class__ == BlockDiagMatrix
diff --git a/sympy/printing/tests/test_mathematica.py b/sympy/printing/tests/test_mathematica.py
index 0e88788885..a2aa5cc16a 100644
--- a/sympy/printing/tests/test_mathematica.py
+++ b/sympy/printing/tests/test_mathematica.py
@@ -127,11 +127,8 @@ def test_constants():
     assert mcode(pi) == "Pi"
     assert mcode(S.GoldenRatio) == "GoldenRatio"
     assert mcode(S.TribonacciConstant) == \
-        "(1/3 + (1/3)*(19 - 3*33^(1/2))^(1/3) + " \
-        "(1/3)*(3*33^(1/2) + 19)^(1/3))"
-    assert mcode(2*S.TribonacciConstant) == \
-        "2*(1/3 + (1/3)*(19 - 3*33^(1/2))^(1/3) + " \
-        "(1/3)*(3*33^(1/2) + 19)^(1/3))"
+        "1/3 + (1/3)*(19 - 3*33^(1/2))^(1/3) + " \
+        "(1/3)*(3*33^(1/2) + 19)^(1/3)"
     assert mcode(S.EulerGamma) == "EulerGamma"
     assert mcode(S.Catalan) == "Catalan"
 
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index 367e6d1654..7ffc5f335a 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -1345,10 +1345,16 @@ def test_failing_assumptions():
 
 def test_issue_6056():
     assert solve(tanh(x + 3)*tanh(x - 3) - 1) == []
-    assert solve(tanh(x - 1)*tanh(x + 1) + 1) == \
-            [-3*I*pi/4, -I*pi/4, I*pi/4, 3*I*pi/4]
-    assert solve((tanh(x + 3)*tanh(x - 3) + 1)**2) == \
-            [-3*I*pi/4, -I*pi/4, I*pi/4, 3*I*pi/4]
+    assert set([simplify(w) for w in solve(tanh(x - 1)*tanh(x + 1) + 1)]) == set([
+        -log(2)/2 + log(1 - I),
+        -log(2)/2 + log(-1 - I),
+        -log(2)/2 + log(1 + I),
+        -log(2)/2 + log(-1 + I),])
+    assert set([simplify(w) for w in solve((tanh(x + 3)*tanh(x - 3) + 1)**2)]) == set([
+        -log(2)/2 + log(1 - I),
+        -log(2)/2 + log(-1 - I),
+        -log(2)/2 + log(1 + I),
+        -log(2)/2 + log(-1 + I),])
 
 
 def test_issue_5673():
@@ -1856,16 +1862,19 @@ def test_issue_8828():
 @slow
 def test_issue_2840_8155():
     assert solve(sin(3*x) + sin(6*x)) == [
-        0, -5*pi/3, -4*pi/3, -pi, -2*pi/3, -pi/3, pi/3, 2*pi/3, pi, 4*pi/3,
-        14*pi/9, 5*pi/3, 16*pi/9, 2*pi, -2*I*log(-(-1)**(S(1)/9)),
-        -2*I*log(-(-1)**(S(2)/9)), -2*I*log(-sin(pi/18) - I*cos(pi/18)),
-        -2*I*log(-sin(pi/18) + I*cos(pi/18)),
-        -2*I*log(sin(pi/18) - I*cos(pi/18)),
-        -2*I*log(sin(pi/18) + I*cos(pi/18)),
-        -2*I*log(exp(-2*I*pi/9)), -2*I*log(exp(-I*pi/9)),
-        -2*I*log(exp(I*pi/9)), -2*I*log(exp(2*I*pi/9))]
+        0, -pi, pi, 14*pi/9, 16*pi/9, 2*pi, 2*I*(log(2) - log(-1 - sqrt(3)*I)),
+        2*I*(log(2) - log(-1 + sqrt(3)*I)), 2*I*(log(2) - log(1 - sqrt(3)*I)),
+        2*I*(log(2) - log(1 + sqrt(3)*I)), 2*I*(log(2) - log(-sqrt(3) - I)),
+        2*I*(log(2) - log(-sqrt(3) + I)), 2*I*(log(2) - log(sqrt(3) - I)),
+        2*I*(log(2) - log(sqrt(3) + I)), -2*I*log(-(-1)**(S(1)/9)), -2*I*log(
+        -(-1)**(S(2)/9)), -2*I*log(-sin(pi/18) - I*cos(pi/18)), -2*I*log(-sin(
+        pi/18) + I*cos(pi/18)), -2*I*log(sin(pi/18) - I*cos(pi/18)), -2*I*log(
+        sin(pi/18) + I*cos(pi/18)), -2*I*log(exp(-2*I*pi/9)), -2*I*log(exp(
+        -I*pi/9)), -2*I*log(exp(I*pi/9)), -2*I*log(exp(2*I*pi/9))]
     assert solve(2*sin(x) - 2*sin(2*x)) == [
-        0, -5*pi/3, -pi, -pi/3, pi/3, pi, 5*pi/3]
+        0, -pi, pi, 2*I*(log(2) - log(-sqrt(3) - I)), 2*I*(log(2) -
+        log(-sqrt(3) + I)), 2*I*(log(2) - log(sqrt(3) - I)), 2*I*(log(2) -
+        log(sqrt(3) + I))]
 
 
 def test_issue_9567():