def test_triangular():
    assert ask(Q.upper_triangular(X+Z.T+Identity(2)), Q.upper_triangular(X) &
            Q.lower_triangular(Z)) is True
    assert ask(Q.upper_triangular(X*Z.T), Q.upper_triangular(X) &
            Q.lower_triangular(Z)) is True
    assert ask(Q.lower_triangular(Identity(3))) is True
    assert ask(Q.lower_triangular(ZeroMatrix(3, 3))) is True

def test_assumptions():
    n = Symbol('n')
    A = MatrixSymbol('A', 1, n)
    P = PermutationMatrix(A)
    assert ask(Q.integer_elements(P))
    assert ask(Q.real_elements(P))
    assert ask(Q.complex_elements(P))

def test_non_trivial_implies():
    X = MatrixSymbol('X', 3, 3)
    Y = MatrixSymbol('Y', 3, 3)
    assert ask(Q.lower_triangular(X+Y), Q.lower_triangular(X) &
               Q.lower_triangular(Y))
    assert ask(Q.triangular(X), Q.lower_triangular(X))
    assert ask(Q.triangular(X+Y), Q.lower_triangular(X) &
               Q.lower_triangular(Y))

 def _eval_determinant(self):
     if self.blockshape == (2, 2):
         [[A, B],
          [C, D]] = self.blocks.tolist()
         if ask(Q.invertible(A)):
             return det(A)*det(D - C*A.I*B)
         elif ask(Q.invertible(D)):
             return det(D)*det(A - B*D.I*C)
     return Determinant(self)

def test_orthogonal():
    assert ask(Q.orthogonal(X), Q.orthogonal(X))
    assert ask(Q.orthogonal(X.T), Q.orthogonal(X)) is True
    assert ask(Q.orthogonal(X.I), Q.orthogonal(X)) is True
    assert ask(Q.orthogonal(Y)) is False
    assert ask(Q.orthogonal(X)) is None
    assert ask(Q.orthogonal(X*Z*X), Q.orthogonal(X) & Q.orthogonal(Z)) is True
    assert ask(Q.orthogonal(Identity(3))) is True
    assert ask(Q.orthogonal(ZeroMatrix(3, 3))) is False
    assert ask(Q.invertible(X), Q.orthogonal(X))
    assert not ask(Q.orthogonal(X+Z), Q.orthogonal(X) & Q.orthogonal(Z))

def test_invertible():
    assert ask(Q.invertible(X), Q.invertible(X))
    assert ask(Q.invertible(Y)) is False
    assert ask(Q.invertible(X*Y), Q.invertible(X)) is False
    assert ask(Q.invertible(X*Z), Q.invertible(X)) is None
    assert ask(Q.invertible(X*Z), Q.invertible(X) & Q.invertible(Z)) is True
    assert ask(Q.invertible(X.T)) is None
    assert ask(Q.invertible(X.T), Q.invertible(X)) is True
    assert ask(Q.invertible(X.I)) is True
    assert ask(Q.invertible(Identity(3))) is True
    assert ask(Q.invertible(ZeroMatrix(3, 3))) is False

def _test_orthogonal_unitary(predicate):
    assert ask(predicate(X), predicate(X))
    assert ask(predicate(X.T), predicate(X)) is True
    assert ask(predicate(X.I), predicate(X)) is True
    assert ask(predicate(Y)) is False
    assert ask(predicate(X)) is None
    assert ask(predicate(X*Z*X), predicate(X) & predicate(Z)) is True
    assert ask(predicate(Identity(3))) is True
    assert ask(predicate(ZeroMatrix(3, 3))) is False
    assert ask(Q.invertible(X), predicate(X))
    assert not ask(predicate(X + Z), predicate(X) & predicate(Z))

def test_fullrank():
    assert ask(Q.fullrank(X), Q.fullrank(X))
    assert ask(Q.fullrank(X**2), Q.fullrank(X))
    assert ask(Q.fullrank(X.T), Q.fullrank(X)) is True
    assert ask(Q.fullrank(X)) is None
    assert ask(Q.fullrank(Y)) is None
    assert ask(Q.fullrank(X*Z), Q.fullrank(X) & Q.fullrank(Z)) is True
    assert ask(Q.fullrank(Identity(3))) is True
    assert ask(Q.fullrank(ZeroMatrix(3, 3))) is False
    assert ask(Q.invertible(X), ~Q.fullrank(X)) == False

def test_symmetric():
    assert ask(Q.symmetric(X), Q.symmetric(X))
    assert ask(Q.symmetric(X*Z), Q.symmetric(X)) is None
    assert ask(Q.symmetric(X*Z), Q.symmetric(X) & Q.symmetric(Z)) is True
    assert ask(Q.symmetric(X+Z), Q.symmetric(X) & Q.symmetric(Z)) is True
    assert ask(Q.symmetric(Y)) is False
    assert ask(Q.symmetric(Y*Y.T)) is True
    assert ask(Q.symmetric(Y.T*X*Y)) is None
    assert ask(Q.symmetric(Y.T*X*Y), Q.symmetric(X)) is True
    assert ask(Q.symmetric(X*X*X*X*X*X*X*X*X*X), Q.symmetric(X)) is True

def dtype_of(expr, *assumptions):
    if hasattr(expr, 'fortran_type'):
        return expr.fortran_type()

    with assuming(*assumptions):
        if ask(Q.integer(expr) | Q.integer_elements(expr)) or expr.is_integer:
            result = 'integer'
        elif ask(Q.real(expr) | Q.real_elements(expr)) or expr.is_real:
            result = 'real(kind=8)'
        elif ask(Q.complex(expr) | Q.complex_elements(expr)) or expr.is_complex:
            result = 'complex(kind=8)'
        else:
            raise TypeError('Could not infer type of %s'%str(expr))
    return result

def test_MatrixSlice():
    X = MatrixSymbol('X', 4, 4)
    B = MatrixSlice(X, (1, 3), (1, 3))
    C = MatrixSlice(X, (0, 3), (1, 3))
    assert ask(Q.symmetric(B), Q.symmetric(X))
    assert ask(Q.invertible(B), Q.invertible(X))
    assert ask(Q.diagonal(B), Q.diagonal(X))
    assert ask(Q.orthogonal(B), Q.orthogonal(X))
    assert ask(Q.upper_triangular(B), Q.upper_triangular(X))

    assert not ask(Q.symmetric(C), Q.symmetric(X))
    assert not ask(Q.invertible(C), Q.invertible(X))
    assert not ask(Q.diagonal(C), Q.diagonal(X))
    assert not ask(Q.orthogonal(C), Q.orthogonal(X))
    assert not ask(Q.upper_triangular(C), Q.upper_triangular(X))

 def doit(self, expand=False):
     if ask(Q.singular(self)):
         return S.Zero
     try:
         return self.arg._eval_determinant()
     except (AttributeError, NotImplementedError):
         return self

def isolate(alg, eps=None, fast=False):
    """Give a rational isolating interval for an algebraic number. """
    alg = sympify(alg)

    if alg.is_Rational:
        return (alg, alg)
    elif not ask(Q.real(alg)):
        raise NotImplementedError(
            "complex algebraic numbers are not supported")

    func = lambdify((), alg, modules="mpmath", printer=IntervalPrinter())

    poly = minpoly(alg, polys=True)
    intervals = poly.intervals(sqf=True)

    dps, done = mp.dps, False

    try:
        while not done:
            alg = func()

            for a, b in intervals:
                if a <= alg.a and alg.b <= b:
                    done = True
                    break
            else:
                mp.dps *= 2
    finally:
        mp.dps = dps

    if eps is not None:
        a, b = poly.refine_root(a, b, eps=eps, fast=fast)

    return (a, b)

def arctan_rule(integral):
    integrand, symbol = integral
    base, exp = integrand.as_base_exp()

    if sympy.simplify(exp + 1) == 0:
        a = sympy.Wild('a', exclude=[symbol])
        b = sympy.Wild('b', exclude=[symbol])
        match = base.match(a + b*symbol**2)
        if match:
            a, b = match[a], match[b]
            if ((isinstance(a, sympy.Number) and a < 0) or (isinstance(b, sympy.Number) and b < 0)):
                return
            if (sympy.ask(sympy.Q.negative(a) | sympy.Q.negative(b) | sympy.Q.is_true(a <= 0) | sympy.Q.is_true(b <= 0))):
                return
            #   /    dx       1  /   dx             1   /     dx                |                     |    1     1         /   du
            #  | --------- = -- | -------------- = --  | -------------------- = | sqrt(b/a)x = u      | =  -- ----------  | -------
            # /  a + bx^2    a /   1  + (b/a)x^2   a  /   1 + (sqrt(b/a)x)^2    | dx = du / sqrt(b/a) |    a   sqrt(b/a) /  1 + u^2
            if a == 1 and b == 1:
                return ArctanRule(integrand, symbol)
            if a == b:
                    constant = 1 / a
                    integrand_ = 1 / (1 + symbol**2)
                    substep = ArctanRule(integrand_, symbol)
                    return ConstantTimesRule(constant, integrand_, substep, integrand, symbol)
            u_var = new_symbol_(symbol)
            u_func = sympy.sqrt(sympy.sympify(b) / a) * symbol
            integrand_ = 1 / (1 + u_func**2)
            constant = 1 / sympy.sqrt(sympy.sympify(b) / a)
            substituted = 1 / (1 + u_var**2)
            substep = ArctanRule(substituted, u_var)
            substep = ConstantTimesRule(constant, substituted, substep, constant*substituted, u_var)
            substep = URule(u_var, u_func, constant, substep, constant*substituted, integrand_, symbol)
            return ConstantTimesRule(1/a, integrand_, substep, integrand, symbol)

def substitution_rule(integral):
    integrand, symbol = integral

    u_var = sympy.Dummy("u")
    substitutions = find_substitutions(integrand, symbol, u_var)
    if substitutions:
        ways = []
        for u_func, c, substituted in substitutions:
            subrule = integral_steps(substituted, u_var)
            if contains_dont_know(subrule):
                continue

            if sympy.simplify(c - 1) != 0:
                _, denom = c.as_numer_denom()
                subrule = ConstantTimesRule(c, substituted, subrule, substituted, symbol)

                if denom.free_symbols:
                    piecewise = []
                    could_be_zero = []

                    if isinstance(denom, sympy.Mul):
                        could_be_zero = denom.args
                    else:
                        could_be_zero.append(denom)

                    for expr in could_be_zero:
                        if not sympy.ask(~sympy.Q.zero(expr)):
                            substep = integral_steps(integrand.subs(expr, 0), symbol)

                            if substep:
                                piecewise.append((
                                    substep,
                                    sympy.Eq(expr, 0)
                                ))
                    piecewise.append((subrule, True))
                    subrule = PiecewiseRule(piecewise, substituted, symbol)

            ways.append(URule(u_var, u_func, c,
                              subrule,
                              integrand, symbol))

        if len(ways) > 1:
            return AlternativeRule(ways, integrand, symbol)
        elif ways:
            return ways[0]

    elif integrand.has(sympy.exp):
        u_func = sympy.exp(symbol)
        c = 1
        substituted = integrand / u_func.diff(symbol)
        substituted = substituted.subs(u_func, u_var)

        if symbol not in substituted.free_symbols:
            return URule(u_var, u_func, c,
                         integral_steps(substituted, u_var),
                         integrand, symbol)

def power_rule(integral):
    integrand, symbol = integral
    base, exp = integrand.as_base_exp()

    if symbol not in exp.free_symbols and isinstance(base, sympy.Symbol):
        if sympy.simplify(exp + 1) == 0:
            return LogRule(base, integrand, symbol)
        return PowerRule(base, exp, integrand, symbol)
    elif symbol not in base.free_symbols and isinstance(exp, sympy.Symbol):
        rule = ExpRule(base, exp, integrand, symbol)

        if sympy.ask(~sympy.Q.zero(sympy.log(base))):
            return rule
        elif sympy.ask(sympy.Q.zero(sympy.log(base))):
            return ConstantRule(1, 1, symbol)

        return PiecewiseRule([
            (ConstantRule(1, 1, symbol), sympy.Eq(sympy.log(base), 0)),
            (rule, True)
        ], integrand, symbol)

def test_DiagonalMatrix():
    x = MatrixSymbol('x', n, m)
    D = DiagonalMatrix(x)
    assert D.diagonal_length is None
    assert D.shape == (n, m)

    x = MatrixSymbol('x', n, n)
    D = DiagonalMatrix(x)
    assert D.diagonal_length == n
    assert D.shape == (n, n)
    assert D[1, 2] == 0
    assert D[1, 1] == x[1, 1]
    i = Symbol('i')
    j = Symbol('j')
    x = MatrixSymbol('x', 3, 3)
    ij = DiagonalMatrix(x)[i, j]
    assert ij != 0
    assert ij.subs({i:0, j:0}) == x[0, 0]
    assert ij.subs({i:0, j:1}) == 0
    assert ij.subs({i:1, j:1}) == x[1, 1]
    assert ask(Q.diagonal(D))  # affirm that D is diagonal

    x = MatrixSymbol('x', n, 3)
    D = DiagonalMatrix(x)
    assert D.diagonal_length == 3
    assert D.shape == (n, 3)
    assert D[2, m] == KroneckerDelta(2, m)*x[2, m]
    assert D[3, m] == 0
    raises(IndexError, lambda: D[m, 3])

    x = MatrixSymbol('x', 3, n)
    D = DiagonalMatrix(x)
    assert D.diagonal_length == 3
    assert D.shape == (3, n)
    assert D[m, 2] == KroneckerDelta(m, 2)*x[m, 2]
    assert D[m, 3] == 0
    raises(IndexError, lambda: D[3, m])

    x = MatrixSymbol('x', n, m)
    D = DiagonalMatrix(x)
    assert D.diagonal_length is None
    assert D.shape == (n, m)
    assert D[m, 4] != 0

    x = MatrixSymbol('x', 3, 4)
    assert [DiagonalMatrix(x)[i] for i in range(12)] == [
        x[0, 0], 0, 0, 0, 0, x[1, 1], 0, 0, 0, 0, x[2, 2], 0]

    # shape is retained, issue 12427
    assert (
        DiagonalMatrix(MatrixSymbol('x', 3, 4))*
        DiagonalMatrix(MatrixSymbol('x', 4, 2))).shape == (3, 2)

def canonise_log(equation):
    expanded_log = sympy.expand_log(equation, force=True)

    terms = expanded_log.as_ordered_terms()

    a, b = sympy.Wild('a'), sympy.Wild('b')

    total_interior = 1
    for term in terms:

        if sympy.ask(sympy.Q.complex(term)):
            term_interior *= -1
        else:
            term_interior = term.match(sympy.log(a) / b)[a]

            if term.could_extract_minus_sign():
                total_interior /= term_interior
            else:
                total_interior *= term_interior

    if isinstance(total_interior, sympy.Add):
        invert = False
    elif isinstance(total_interior, sympy.Mul):

        match = total_interior.together().match(x / b)  # for some reason, (x/3).match(a/b) gives {a: 1/3, b: 1/x} so we have to use a workaround
        if match is not None:
            invert = False
        else:
            match = total_interior.together().match(a / b)

            degree_numerator = 0 if isinstance(match[a], sympy.Rational) else match[a].as_poly().degree()
            degree_denominator = 0 if isinstance(match[b], sympy.Rational) else match[b].as_poly().degree()

            if degree_numerator < degree_denominator:
                invert = True
            else:
                invert = False
    elif isinstance(total_interior, sympy.Pow):
        index = total_interior.as_base_exp()[1]
        if index < 0:
            invert = True
        else:
            invert = False

    else:  # for debugging - wtf kind of c-c-c-class is it???
        print(total_interior, type(total_interior))

    if invert:
        return -sympy.log((1 / total_interior).together(), evaluate=False) / terms[0].as_coeff_Mul()[0].q
    else:
        return sympy.log(total_interior.together(), evaluate=False) / terms[0].as_coeff_Mul()[0].q

def test_diagonal():
    assert ask(Q.diagonal(X + Z.T + Identity(2)), Q.diagonal(X) &
               Q.diagonal(Z)) is True
    assert ask(Q.diagonal(ZeroMatrix(3, 3)))
    assert ask(Q.lower_triangular(X) & Q.upper_triangular(X), Q.diagonal(X))
    assert ask(Q.diagonal(X), Q.lower_triangular(X) & Q.upper_triangular(X))
    assert ask(Q.symmetric(X), Q.diagonal(X))
    assert ask(Q.triangular(X), Q.diagonal(X))

def arctan_rule(integral):
    integrand, symbol = integral
    base, exp = integrand.as_base_exp()

    if sympy.simplify(exp + 1) == 0:
        a = sympy.Wild('a', exclude=[symbol])
        b = sympy.Wild('b', exclude=[symbol])
        match = base.match(a + b*symbol**2)
        if match:
            a, b = match[a], match[b]

            if ((isinstance(a, sympy.Number) and a < 0) or
                (isinstance(b, sympy.Number) and b < 0)):
                return
            if (sympy.ask(sympy.Q.negative(a) | sympy.Q.negative(b) |
                          sympy.Q.is_true(a <= 0) | sympy.Q.is_true(b <= 0))):
                return

            if a != 1 or b != 1:
                b_condition = b >= 0
                u_var = sympy.Dummy("u")
                rewritten = (sympy.Integer(1) / a) * (base / a) ** (-1)
                u_func = sympy.sqrt(sympy.sympify(b) / a) * symbol
                constant = 1 / sympy.sqrt(sympy.sympify(b) / a)
                substituted = rewritten.subs(u_func, u_var)

                if a == b:
                    substep = ArctanRule(integrand, symbol)
                else:
                    subrule = ArctanRule(substituted, u_var)
                    if constant != 1:
                        b_condition = b > 0
                        subrule = ConstantTimesRule(
                            constant, substituted, subrule,
                            substituted, symbol)

                    substep = URule(u_var, u_func, constant,
                                    subrule,
                                    integrand, symbol)

                if a != 1:
                    other = (base / a) ** (-1)
                    substep = ConstantTimesRule(
                        sympy.Integer(1) / a, other, substep,
                        integrand, symbol)

                return PiecewiseRule([
                    (substep, sympy.And(a > 0, b_condition))
                ], integrand, symbol)

            return ArctanRule(integrand, symbol)