Python sympy 模块，Pow() 实例源码

def test_commutator():
    A = Operator('A')
    B = Operator('B')
    c = Commutator(A, B)
    c_tall = Commutator(A**2, B)
    assert str(c) == '[A,B]'
    assert pretty(c) == '[A,B]'
    assert upretty(c) == u('[A,B]')
    assert latex(c) == r'\left[A,B\right]'
    sT(c, "Commutator(Operator(Symbol('A')),Operator(Symbol('B')))")
    assert str(c_tall) == '[A**2,B]'
    ascii_str = \
"""\
[ 2  ]\n\
[A ,B]\
"""
    ucode_str = \
u("""\
? 2  ?\n\
?A ,B?\
""")
    assert pretty(c_tall) == ascii_str
    assert upretty(c_tall) == ucode_str
    assert latex(c_tall) == r'\left[\left(A\right)^{2},B\right]'
    sT(c_tall, "Commutator(Pow(Operator(Symbol('A')), Integer(2)),Operator(Symbol('B')))")

def _expand_powers(factors):
    """
    Helper function for normal_ordered_form and normal_order: Expand a
    power expression to a multiplication expression so that that the
    expression can be handled by the normal ordering functions.
    """

    new_factors = []
    for factor in factors.args:
        if isinstance(factor, Pow):
            for n in range(factor.args[1]):
                new_factors.append(factor.args[0])
        else:
            new_factors.append(factor)

    return new_factors

def _extractFunctionsFromSympy(self, expr, top=False):
        newVars = []
        newArgs = []
        isExpensive = False
        for arg in expr.args:
            (newVarsFromArg, newArg,dontcare) = self._extractFunctionsFromSympy(arg)
            newVars += newVarsFromArg
            newArgs.append(newArg)
        if newVars:
            expr = expr.func(*newArgs)
        if isinstance(expr.func, type(sympy.Function)) or (
                expr.func == sympy.Pow and (
                    not expr.exp.is_constant or
                    int(expr.exp) != expr.exp)):
            if top:
                isExpensive = True
            else:
                newSym = self.addSSA("_expensive_functions", expr)
                expr = newSym
                newVars.append(newSym)
        return (newVars, expr, isExpensive)

def _collect_factor_and_dimension(expr):

        if isinstance(expr, Quantity):
            return expr.scale_factor, expr.dimension
        elif isinstance(expr, Mul):
            factor = 1
            dimension = 1
            for arg in expr.args:
                arg_factor, arg_dim = Quantity._collect_factor_and_dimension(arg)
                factor *= arg_factor
                dimension *= arg_dim
            return factor, dimension
        elif isinstance(expr, Pow):
            factor, dim = Quantity._collect_factor_and_dimension(expr.base)
            return factor ** expr.exp, dim ** expr.exp
        elif isinstance(expr, Add):
            raise NotImplementedError
        else:
            return 1, 1

def _expand_powers(factors):
    """
    Helper function for normal_ordered_form and normal_order: Expand a
    power expression to a multiplication expression so that that the
    expression can be handled by the normal ordering functions.
    """

    new_factors = []
    for factor in factors.args:
        if (isinstance(factor, Pow)
                and isinstance(factor.args[1], Integer)
                and factor.args[1] > 0):
            for n in range(factor.args[1]):
                new_factors.append(factor.args[0])
        else:
            new_factors.append(factor)

    return new_factors

def _print_Pow(self, power):
        from sympy.simplify.simplify import fraction
        b, e = power.as_base_exp()
        if power.is_commutative:
            if e is S.NegativeOne:
                return prettyForm("1")/self._print(b)
            n, d = fraction(e)
            if n is S.One and d.is_Atom and not e.is_Integer:
                return self._print_nth_root(b, e)
            if e.is_Rational and e < 0:
                return prettyForm("1")/self._print(Pow(b, -e, evaluate=False))

        if b.is_Relational:
            return prettyForm(*self._print(b).parens()).__pow__(self._print(e))

        return self._print(b)**self._print(e)

def _real_to_rational(expr, tolerance=None):
    """
    Replace all reals in expr with rationals.

    >>> from sympy import nsimplify
    >>> from sympy.abc import x

    >>> nsimplify(.76 + .1*x**.5, rational=True)
    sqrt(x)/10 + 19/25

    """
    p = expr
    reps = {}
    reduce_num = None
    if tolerance is not None and tolerance < 1:
        reduce_num = ceiling(1/tolerance)
    for float in p.atoms(C.Float):
        key = float
        if reduce_num is not None:
            r = Rational(float).limit_denominator(reduce_num)
        elif (tolerance is not None and tolerance >= 1 and
                float.is_Integer is False):
            r = Rational(tolerance*round(float/tolerance)
                ).limit_denominator(int(tolerance))
        else:
            r = nsimplify(float, rational=False)
            # e.g. log(3).n() -> log(3) instead of a Rational
            if not r.is_Rational:
                if float < 0:
                    float = -float
                    d = Pow(10, int((mpmath.log(float)/mpmath.log(10))))
                    r = -Rational(str(float/d))*d
                elif float > 0:
                    d = Pow(10, int((mpmath.log(float)/mpmath.log(10))))
                    r = Rational(str(float/d))*d
                else:
                    r = Integer(0)
        reps[key] = r
    return p.subs(reps, simultaneous=True)

def _get_min_qubits(a_gate):
    if isinstance(a_gate, Pow):
        return a_gate.base.min_qubits
    else:
        return a_gate.min_qubits

def _sympystr(self, printer, *args):
        from sympy.printing.str import sstr
        length = len(self.args)
        s = ''
        for i in range(length):
            if isinstance(self.args[i], (Add, Pow, Mul)):
                s = s + '('
            s = s + sstr(self.args[i])
            if isinstance(self.args[i], (Add, Pow, Mul)):
                s = s + ')'
            if i != length - 1:
                s = s + 'x'
        return s

def _print_Pow(self, power):
        from sympy.simplify.simplify import fraction
        b, e = power.as_base_exp()
        if power.is_commutative:
            if e is S.NegativeOne:
                return prettyForm("1")/self._print(b)
            n, d = fraction(e)
            if n is S.One and d.is_Atom and not e.is_Integer:
                return self._print_nth_root(b, e)
            if e.is_Rational and e < 0:
                return prettyForm("1")/self._print(C.Pow(b, -e, evaluate=False))

        return self._print(b)**self._print(e)

def _print_ProductSet(self, p):
        if len(p.sets) > 1 and not has_variety(p.sets):
            from sympy import Pow
            return self._print(Pow(p.sets[0], len(p.sets), evaluate=False))
        else:
            prod_char = u('\xd7')
            return self._print_seq(p.sets, None, None, ' %s ' % prod_char,
                parenthesize=lambda set: set.is_Union or set.is_Intersection)

def _eval_nseries(self, x, n, logx):
        # NOTE this is fairly inefficient
        from sympy import log, EulerGamma, Pow
        n += 1
        if self.args[0].subs(x, 0) != 0:
            return super(TrigonometricIntegral, self)._eval_nseries(x, n, logx)
        baseseries = self._trigfunc(x)._eval_nseries(x, n, logx)
        if self._trigfunc(0) != 0:
            baseseries -= 1
        baseseries = baseseries.replace(Pow, lambda t, n: t**n/n, simultaneous=False)
        if self._trigfunc(0) != 0:
            baseseries += EulerGamma + log(x)
        return baseseries.subs(x, self.args[0])._eval_nseries(x, n, logx)

def main():
    x = Pow(2, 50, evaluate=False)
    y = Pow(10, -50, evaluate=False)
    # A large, unevaluated expression
    m = Mul(x, y, evaluate=False)
    # Evaluating the expression
    e = S(2)**50/S(10)**50
    print("%s == %s" % (m, e))

def expr(self):
        """
        Process the atomic expression, by returning a sympy expression
        """
        if DEBUG:
            print("> Power ")
        expr = self.op.expr
        if self.exponent is None:
            return expr
        else:
            exponent = self.exponent.expr
            return Pow(expr, exponent)

def expr(self):
        """
        Process the term, by returning a sympy expression
        """
        if DEBUG:
            print("> Term ")

        ret = self.op[0].expr
        for operation, operand in zip(self.op[1::2], self.op[2::2]):
            if operation == '*':
                ret = Mul(ret, operand.expr)
            else:
                a   = Pow(operand.expr, -1)
                ret = Mul(ret, a)
        return ret

def powerProcess(self, x, y):
        if (not x.astUnit.isNull()) and y.sympy.is_constant():
            newUnit = x.astUnit ** float(y.sympy)
        else:
            newUnit = ASTUnit.null()
        return AST(sympy.Pow(x.sympy, y.sympy), newUnit)

    ############################################

def p_multiplicitiveExpr_impl(self, p):
        '''multiplicitiveExpr : multiplicitiveExpr '*' unaryExpr
                              | multiplicitiveExpr '/' unaryExpr
        '''
        lhs = p[1]
        rhs = p[3]
        if p[2] == '/':
            rhs = AST(sympy.Pow(rhs.sympy,sympy.Integer(-1)), rhs.astUnit ** -1)
        ret = AST(sympy.Mul(lhs.sympy,rhs.sympy), lhs.astUnit*rhs.astUnit)
        p[0] = ret

def test_Pow_is_algebraic():
    e = Symbol('e', algebraic=True)

    assert Pow(1, e, evaluate=False).is_algebraic
    assert Pow(0, e, evaluate=False).is_algebraic

    a = Symbol('a', algebraic=True)
    na = Symbol('na', algebraic=False)
    ia = Symbol('ia', algebraic=True, irrational=True)
    ib = Symbol('ib', algebraic=True, irrational=True)
    r = Symbol('r', rational=True)
    x = Symbol('x')
    assert (a**r).is_algebraic
    assert (a**x).is_algebraic is None
    assert (na**r).is_algebraic is None
    assert (ia**r).is_algebraic
    assert (ia**ib).is_algebraic is False

    assert (a**e).is_algebraic is None

    # Gelfond-Schneider constant:
    assert Pow(2, sqrt(2), evaluate=False).is_algebraic is False

    assert Pow(S.GoldenRatio, sqrt(3), evaluate=False).is_algebraic is False

    # issue 8649
    t = Symbol('t', real=True, transcendental=True)
    n = Symbol('n', integer=True)
    assert (t**n).is_algebraic is None
    assert (t**n).is_integer is None

def test_Pow_is_prime():
    from sympy import Pow
    x = Symbol('x', positive=True, integer=True)
    y = Symbol('y', positive=True, integer=True)
    assert (x**y).is_prime is None

    x = Symbol('x', positive=True)
    assert (x**y).is_prime is None

    assert Pow(6, S.One, evaluate=False).is_prime is False
    assert Pow(9, S.Half, evaluate=False).is_prime is True
    assert Pow(5, S.One, evaluate=False).is_prime is True

def test_dim_simplify_pow():
    assert dim_simplify(Pow(L, 2)) == L**2
    assert dim_simplify(L**2) == L**2

def dim_simplify(expr):
    """
    NOTE: this function could be deprecated in the future.

    Simplify expression by recursively evaluating the dimension arguments.

    This function proceeds to a very rough dimensional analysis. It tries to
    simplify expression with dimensions, and it deletes all what multiplies a
    dimension without being a dimension. This is necessary to avoid strange
    behavior when Add(L, L) be transformed into Mul(2, L).
    """

    if isinstance(expr, Dimension):
        return expr

    if isinstance(expr, Pow):
        return dim_simplify(expr.base)**dim_simplify(expr.exp)
    elif isinstance(expr, Function):
        return dim_simplify(expr.args[0])
    elif isinstance(expr, Add):
        if (all(isinstance(arg, Dimension) for arg in expr.args) or
            all(arg.is_dimensionless for arg in expr.args if isinstance(arg, Dimension))):
            return reduce(lambda x, y: x.add(y), expr.args)
        else:
            raise ValueError("Dimensions cannot be added: %s" % expr)
    elif isinstance(expr, Mul):
        return Dimension(Mul(*[dim_simplify(i).name for i in expr.args if isinstance(i, Dimension)]))

    raise ValueError("Cannot be simplifed: %s", expr)

def _get_min_qubits(a_gate):
    if isinstance(a_gate, Pow):
        return a_gate.base.min_qubits
    else:
        return a_gate.min_qubits

def test_anticommutator():
    A = Operator('A')
    B = Operator('B')
    ac = AntiCommutator(A, B)
    ac_tall = AntiCommutator(A**2, B)
    assert str(ac) == '{A,B}'
    assert pretty(ac) == '{A,B}'
    assert upretty(ac) == u'{A,B}'
    assert latex(ac) == r'\left\{A,B\right\}'
    sT(ac, "AntiCommutator(Operator(Symbol('A')),Operator(Symbol('B')))")
    assert str(ac_tall) == '{A**2,B}'
    ascii_str = \
"""\
/ 2  \\\n\
<A ,B>\n\
\\    /\
"""
    ucode_str = \
u("""\
? 2  ?\n\
?A ,B?\n\
?    ?\
""")
    assert pretty(ac_tall) == ascii_str
    assert upretty(ac_tall) == ucode_str
    assert latex(ac_tall) == r'\left\{A^{2},B\right\}'
    sT(ac_tall, "AntiCommutator(Pow(Operator(Symbol('A')), Integer(2)),Operator(Symbol('B')))")

def _sympystr(self, printer, *args):
        from sympy.printing.str import sstr
        length = len(self.args)
        s = ''
        for i in range(length):
            if isinstance(self.args[i], (Add, Pow, Mul)):
                s = s + '('
            s = s + sstr(self.args[i])
            if isinstance(self.args[i], (Add, Pow, Mul)):
                s = s + ')'
            if i != length - 1:
                s = s + 'x'
        return s

def _print_ProductSet(self, p):
        if len(p.sets) > 1 and not has_variety(p.sets):
            from sympy import Pow
            return self._print(Pow(p.sets[0], len(p.sets), evaluate=False))
        else:
            prod_char = u"\N{MULTIPLICATION SIGN}" if self._use_unicode else 'x'
            return self._print_seq(p.sets, None, None, ' %s ' % prod_char,
                                   parenthesize=lambda set: set.is_Union or
                                   set.is_Intersection or set.is_ProductSet)

def _eval_nseries(self, x, n, logx):
        # NOTE this is fairly inefficient
        from sympy import log, EulerGamma, Pow
        n += 1
        if self.args[0].subs(x, 0) != 0:
            return super(TrigonometricIntegral, self)._eval_nseries(x, n, logx)
        baseseries = self._trigfunc(x)._eval_nseries(x, n, logx)
        if self._trigfunc(0) != 0:
            baseseries -= 1
        baseseries = baseseries.replace(Pow, lambda t, n: t**n/n, simultaneous=False)
        if self._trigfunc(0) != 0:
            baseseries += EulerGamma + log(x)
        return baseseries.subs(x, self.args[0])._eval_nseries(x, n, logx)

def _separatevars(expr, force):
    if len(expr.free_symbols) == 1:
        return expr
    # don't destroy a Mul since much of the work may already be done
    if expr.is_Mul:
        args = list(expr.args)
        changed = False
        for i, a in enumerate(args):
            args[i] = separatevars(a, force)
            changed = changed or args[i] != a
        if changed:
            expr = expr.func(*args)
        return expr

    # get a Pow ready for expansion
    if expr.is_Pow:
        expr = Pow(separatevars(expr.base, force=force), expr.exp)

    # First try other expansion methods
    expr = expr.expand(mul=False, multinomial=False, force=force)

    _expr, reps = posify(expr) if force else (expr, {})
    expr = factor(_expr).subs(reps)

    if not expr.is_Add:
        return expr

    # Find any common coefficients to pull out
    args = list(expr.args)
    commonc = args[0].args_cnc(cset=True, warn=False)[0]
    for i in args[1:]:
        commonc &= i.args_cnc(cset=True, warn=False)[0]
    commonc = Mul(*commonc)
    commonc = commonc.as_coeff_Mul()[1]  # ignore constants
    commonc_set = commonc.args_cnc(cset=True, warn=False)[0]

    # remove them
    for i, a in enumerate(args):
        c, nc = a.args_cnc(cset=True, warn=False)
        c = c - commonc_set
        args[i] = Mul(*c)*Mul(*nc)
    nonsepar = Add(*args)

    if len(nonsepar.free_symbols) > 1:
        _expr = nonsepar
        _expr, reps = posify(_expr) if force else (_expr, {})
        _expr = (factor(_expr)).subs(reps)

        if not _expr.is_Add:
            nonsepar = _expr

    return commonc*nonsepar

def __new__(cls, *args):
        """ Construct a Trace object.

        Parameters
        ==========
        args = sympy expression
        indices = tuple/list if indices, optional

        """

        # expect no indices,int or a tuple/list/Tuple
        if (len(args) == 2):
            if not isinstance(args[1], (list, Tuple, tuple)):
                indices = Tuple(args[1])
            else:
                indices = Tuple(*args[1])

            expr = args[0]
        elif (len(args) == 1):
            indices = Tuple()
            expr = args[0]
        else:
            raise ValueError("Arguments to Tr should be of form "
                             "(expr[, [indices]])")

        if isinstance(expr, Matrix):
            return expr.trace()
        elif hasattr(expr, 'trace') and callable(expr.trace):
            #for any objects that have trace() defined e.g numpy
            return expr.trace()
        elif isinstance(expr, Add):
            return Add(*[Tr(arg, indices) for arg in expr.args])
        elif isinstance(expr, Mul):
            c_part, nc_part = expr.args_cnc()
            if len(nc_part) == 0:
                return Mul(*c_part)
            else:
                obj = Expr.__new__(cls, Mul(*nc_part), indices )
                #this check is needed to prevent cached instances
                #being returned even if len(c_part)==0
                return Mul(*c_part)*obj if len(c_part) > 0 else obj
        elif isinstance(expr, Pow):
            if (_is_scalar(expr.args[0]) and
                    _is_scalar(expr.args[1])):
                return expr
            else:
                return Expr.__new__(cls, expr, indices)
        else:
            if (_is_scalar(expr)):
                return expr

            return Expr.__new__(cls, expr, indices)

def gate_sort(circuit):
    """Sorts the gates while keeping track of commutation relations

    This function uses a bubble sort to rearrange the order of gate
    application. Keeps track of Quantum computations special commutation
    relations (e.g. things that apply to the same Qubit do not commute with
    each other)

    circuit is the Mul of gates that are to be sorted.
    """
    # Make sure we have an Add or Mul.
    if isinstance(circuit, Add):
        return sum(gate_sort(t) for t in circuit.args)
    if isinstance(circuit, Pow):
        return gate_sort(circuit.base)**circuit.exp
    elif isinstance(circuit, Gate):
        return circuit
    if not isinstance(circuit, Mul):
        return circuit

    changes = True
    while changes:
        changes = False
        circ_array = circuit.args
        for i in xrange(len(circ_array) - 1):
            # Go through each element and switch ones that are in wrong order
            if isinstance(circ_array[i], (Gate, Pow)) and \
                    isinstance(circ_array[i + 1], (Gate, Pow)):
                # If we have a Pow object, look at only the base
                first_base, first_exp = circ_array[i].as_base_exp()
                second_base, second_exp = circ_array[i + 1].as_base_exp()

                # Use sympy's hash based sorting. This is not mathematical
                # sorting, but is rather based on comparing hashes of objects.
                # See Basic.compare for details.
                if first_base.compare(second_base) > 0:
                    if Commutator(first_base, second_base).doit() == 0:
                        new_args = (circuit.args[:i] + (circuit.args[i + 1],) +
                                   (circuit.args[i],) + circuit.args[i + 2:])
                        circuit = Mul(*new_args)
                        circ_array = circuit.args
                        changes = True
                        break
                    if AntiCommutator(first_base, second_base).doit() == 0:
                        new_args = (circuit.args[:i] + (circuit.args[i + 1],) +
                                   (circuit.args[i],) + circuit.args[i + 2:])
                        sign = Integer(-1)**(first_exp*second_exp)
                        circuit = sign*Mul(*new_args)
                        circ_array = circuit.args
                        changes = True
                        break
    return circuit


#-----------------------------------------------------------------------------
# Utility functions
#-----------------------------------------------------------------------------

def test_operator():
    a = Operator('A')
    b = Operator('B', Symbol('t'), S(1)/2)
    inv = a.inv()
    f = Function('f')
    x = symbols('x')
    d = DifferentialOperator(Derivative(f(x), x), f(x))
    op = OuterProduct(Ket(), Bra())
    assert str(a) == 'A'
    assert pretty(a) == 'A'
    assert upretty(a) == u('A')
    assert latex(a) == 'A'
    sT(a, "Operator(Symbol('A'))")
    assert str(inv) == 'A**(-1)'
    ascii_str = \
"""\
 -1\n\
A  \
"""
    ucode_str = \
u("""\
 -1\n\
A  \
""")
    assert pretty(inv) == ascii_str
    assert upretty(inv) == ucode_str
    assert latex(inv) == r'\left(A\right)^{-1}'
    sT(inv, "Pow(Operator(Symbol('A')), Integer(-1))")
    assert str(d) == 'DifferentialOperator(Derivative(f(x), x),f(x))'
    ascii_str = \
"""\
                    /d            \\\n\
DifferentialOperator|--(f(x)),f(x)|\n\
                    \dx           /\
"""
    ucode_str = \
u("""\
                    ?d            ?\n\
DifferentialOperator???(f(x)),f(x)?\n\
                    ?dx           ?\
""")
    assert pretty(d) == ascii_str
    assert upretty(d) == ucode_str
    assert latex(d) == \
        r'DifferentialOperator\left(\frac{d}{d x} f{\left (x \right )},f{\left (x \right )}\right)'
    sT(d, "DifferentialOperator(Derivative(Function('f')(Symbol('x')), Symbol('x')),Function('f')(Symbol('x')))")
    assert str(b) == 'Operator(B,t,1/2)'
    assert pretty(b) == 'Operator(B,t,1/2)'
    assert upretty(b) == u('Operator(B,t,1/2)')
    assert latex(b) == r'Operator\left(B,t,\frac{1}{2}\right)'
    sT(b, "Operator(Symbol('B'),Symbol('t'),Rational(1, 2))")
    assert str(op) == '|psi><psi|'
    assert pretty(op) == '|psi><psi|'
    assert upretty(op) == u('??????')
    assert latex(op) == r'{\left|\psi\right\rangle }{\left\langle \psi\right|}'
    sT(op, "OuterProduct(Ket(Symbol('psi')),Bra(Symbol('psi')))")

def tensor_product_simp(e, **hints):
    """Try to simplify and combine TensorProducts.

    In general this will try to pull expressions inside of ``TensorProducts``.
    It currently only works for relatively simple cases where the products have
    only scalars, raw ``TensorProducts``, not ``Add``, ``Pow``, ``Commutators``
    of ``TensorProducts``. It is best to see what it does by showing examples.

    Examples
    ========

    >>> from sympy.physics.quantum import tensor_product_simp
    >>> from sympy.physics.quantum import TensorProduct
    >>> from sympy import Symbol
    >>> A = Symbol('A',commutative=False)
    >>> B = Symbol('B',commutative=False)
    >>> C = Symbol('C',commutative=False)
    >>> D = Symbol('D',commutative=False)

    First see what happens to products of tensor products:

    >>> e = TensorProduct(A,B)*TensorProduct(C,D)
    >>> e
    AxB*CxD
    >>> tensor_product_simp(e)
    (A*C)x(B*D)

    This is the core logic of this function, and it works inside, powers, sums,
    commutators and anticommutators as well:

    >>> tensor_product_simp(e**2)
    (A*C)x(B*D)**2

    """
    if isinstance(e, Add):
        return Add(*[tensor_product_simp(arg) for arg in e.args])
    elif isinstance(e, Pow):
        return tensor_product_simp(e.base) ** e.exp
    elif isinstance(e, Mul):
        return tensor_product_simp_Mul(e)
    elif isinstance(e, Commutator):
        return Commutator(*[tensor_product_simp(arg) for arg in e.args])
    elif isinstance(e, AntiCommutator):
        return AntiCommutator(*[tensor_product_simp(arg) for arg in e.args])
    else:
        return e

def find_substitutions(integrand, symbol, u_var):
    results = []

    def test_subterm(u, u_diff):
        substituted = integrand / u_diff
        if symbol not in substituted.free_symbols:
            # replaced everything already
            return False

        substituted = substituted.subs(u, u_var).cancel()
        if symbol not in substituted.free_symbols:
            return substituted.as_independent(u_var, as_Add=False)

        return False

    def possible_subterms(term):
        if isinstance(term, (TrigonometricFunction,
                             sympy.asin, sympy.acos, sympy.atan,
                             sympy.exp, sympy.log, sympy.Heaviside)):
            return [term.args[0]]
        elif isinstance(term, sympy.Mul):
            r = []
            for u in term.args:
                r.append(u)
                r.extend(possible_subterms(u))
            return r
        elif isinstance(term, sympy.Pow):
            if term.args[1].is_constant(symbol):
                return [term.args[0]]
            elif term.args[0].is_constant(symbol):
                return [term.args[1]]
        elif isinstance(term, sympy.Add):
            r = []
            for arg in term.args:
                r.append(arg)
                r.extend(possible_subterms(arg))
            return r
        return []

    for u in possible_subterms(integrand):
        if u == symbol:
            continue
        u_diff = manual_diff(u, symbol)
        new_integrand = test_subterm(u, u_diff)
        if new_integrand is not False:
            constant, new_integrand = new_integrand
            substitution = (u, constant, new_integrand)
            if substitution not in results:
                results.append(substitution)

    return results

def processAssignment(self, var, operand, rhs):
        subsystem = self.currentSubsystem()
        scope = self.currentScope()
        #error if this is a diffvar
        if var in subsystem.diffvars:
            raise XXXSyntaxError("Can't assign to differential variable '%s'." % var)
        if var in subsystem.inputs:
            raise XXXSyntaxError("Can't assign to shared variable '%s'." % var)
        if var == self.timeVar:
            raise XXXSyntaxError("Can't assign to integration variable '%s'." % var)
        if var in subsystem.accums:
            if operand != '+=' and operand != '-=':
                raise XXXSyntaxError("Can only += or -= accumulation variable '%s'."%var)
        if var not in subsystem.outputs and not scope.hasSymbol(var):
            (exists, dontcare) = self.searchForJunction(var)
            if exists:
                raise XXXSyntaxError("Can't assign to shared variable '%s'." % var)

        if operand != '=':
            if not scope.hasSymbol(var):
                if var in subsystem.accums:
                    lhs = textToAST("0",ASTUnit(scope.getUnit(var),False))
                else:
                    raise XXXSyntaxError("'%s' used before assignment."%var)
            else:
                lhs = self.readAccessVar(var)

            if operand == '+=' or operand == '-=':
                if operand == '-=':
                    rhs = AST(sympy.Mul(sympy.Integer(-1),rhs.sympy), rhs.astUnit)
                rhs = AST(sympy.Add(lhs.sympy,rhs.sympy), self.checkExactUnits(lhs.astUnit,rhs.astUnit))
            elif operand == '*=' or operand == '/=':
                if operand == "/=":
                    rhs = AST(sympy.Pow(rhs.sympy,sympy.Integer(-1)), rhs.astUnit ** -1)
                rhs = AST(sympy.Mul(lhs.sympy,rhs.sympy), lhs.astUnit*rhs.astUnit)
            elif operand == '^=':
                rhs = self.powerProcess(lhs, rhs)
            else:
                assert(0)

        #Time to do final unit checks
        if scope.hasUnit(var):
            #print var
            #print scope.getUnit(var)
            #print rhs
            rhs = self.checkExplicitCast(scope.getUnit(var), rhs)

        #ok, we're ready to do the assignment!
        if var in subsystem.frozen:
            XXXSyntaxError('"%s" cannot be assigned once it has been read.' % var)
        symbol = Symbol(var)
        subsystem.ssa[symbol] = rhs
        scope.addInstruction(symbol)
        scope.setSymbol(var, symbol)

        #parameter checking code.
        #FIXME
        nameDepend = set([symbol.name for symbol in rhs.dependencies()])
        newlyReadParams = (set(subsystem.attributeMap.keys()) & nameDepend) - subsystem.frozen
        newlyReadParams -= set([var])
        subsystem.frozen |= newlyReadParams

def __eq__(self, other):
        """Return a boolean indicating whether a == b on the basis of
        their symbolic trees.

        This is the same as a.compare(b) == 0 but faster.

        Notes
        =====

        If a class that overrides __eq__() needs to retain the
        implementation of __hash__() from a parent class, the
        interpreter must be told this explicitly by setting __hash__ =
        <ParentClass>.__hash__. Otherwise the inheritance of __hash__()
        will be blocked, just as if __hash__ had been explicitly set to
        None.

        References
        ==========

        from http://docs.python.org/dev/reference/datamodel.html#object.__hash__
        """
        from sympy import Pow
        if self is other:
            return True

        from .function import AppliedUndef, UndefinedFunction as UndefFunc

        if isinstance(self, UndefFunc) and isinstance(other, UndefFunc):
            if self.class_key() == other.class_key():
                return True
            else:
                return False
        if type(self) is not type(other):
            # issue 6100 a**1.0 == a like a**2.0 == a**2
            if isinstance(self, Pow) and self.exp == 1:
                return self.base == other
            if isinstance(other, Pow) and other.exp == 1:
                return self == other.base
            try:
                other = _sympify(other)
            except SympifyError:
                return False    # sympy != other

            if isinstance(self, AppliedUndef) and isinstance(other,
                                                             AppliedUndef):
                if self.class_key() != other.class_key():
                    return False
            elif type(self) is not type(other):
                return False

        return self._hashable_content() == other._hashable_content()

def __new__(cls, *args):
        """ Construct a Trace object.

        Parameters
        ==========
        args = sympy expression
        indices = tuple/list if indices, optional

        """

        # expect no indices,int or a tuple/list/Tuple
        if (len(args) == 2):
            if not isinstance(args[1], (list, Tuple, tuple)):
                indices = Tuple(args[1])
            else:
                indices = Tuple(*args[1])

            expr = args[0]
        elif (len(args) == 1):
            indices = Tuple()
            expr = args[0]
        else:
            raise ValueError("Arguments to Tr should be of form "
                             "(expr[, [indices]])")

        if isinstance(expr, Matrix):
            return expr.trace()
        elif hasattr(expr, 'trace') and callable(expr.trace):
            #for any objects that have trace() defined e.g numpy
            return expr.trace()
        elif isinstance(expr, Add):
            return Add(*[Tr(arg, indices) for arg in expr.args])
        elif isinstance(expr, Mul):
            c_part, nc_part = expr.args_cnc()
            if len(nc_part) == 0:
                return Mul(*c_part)
            else:
                obj = Expr.__new__(cls, Mul(*nc_part), indices )
                #this check is needed to prevent cached instances
                #being returned even if len(c_part)==0
                return Mul(*c_part)*obj if len(c_part) > 0 else obj
        elif isinstance(expr, Pow):
            if (_is_scalar(expr.args[0]) and
                    _is_scalar(expr.args[1])):
                return expr
            else:
                return Expr.__new__(cls, expr, indices)
        else:
            if (_is_scalar(expr)):
                return expr

            return Expr.__new__(cls, expr, indices)

def gate_sort(circuit):
    """Sorts the gates while keeping track of commutation relations

    This function uses a bubble sort to rearrange the order of gate
    application. Keeps track of Quantum computations special commutation
    relations (e.g. things that apply to the same Qubit do not commute with
    each other)

    circuit is the Mul of gates that are to be sorted.
    """
    # Make sure we have an Add or Mul.
    if isinstance(circuit, Add):
        return sum(gate_sort(t) for t in circuit.args)
    if isinstance(circuit, Pow):
        return gate_sort(circuit.base)**circuit.exp
    elif isinstance(circuit, Gate):
        return circuit
    if not isinstance(circuit, Mul):
        return circuit

    changes = True
    while changes:
        changes = False
        circ_array = circuit.args
        for i in range(len(circ_array) - 1):
            # Go through each element and switch ones that are in wrong order
            if isinstance(circ_array[i], (Gate, Pow)) and \
                    isinstance(circ_array[i + 1], (Gate, Pow)):
                # If we have a Pow object, look at only the base
                first_base, first_exp = circ_array[i].as_base_exp()
                second_base, second_exp = circ_array[i + 1].as_base_exp()

                # Use sympy's hash based sorting. This is not mathematical
                # sorting, but is rather based on comparing hashes of objects.
                # See Basic.compare for details.
                if first_base.compare(second_base) > 0:
                    if Commutator(first_base, second_base).doit() == 0:
                        new_args = (circuit.args[:i] + (circuit.args[i + 1],) +
                                   (circuit.args[i],) + circuit.args[i + 2:])
                        circuit = Mul(*new_args)
                        circ_array = circuit.args
                        changes = True
                        break
                    if AntiCommutator(first_base, second_base).doit() == 0:
                        new_args = (circuit.args[:i] + (circuit.args[i + 1],) +
                                   (circuit.args[i],) + circuit.args[i + 2:])
                        sign = Integer(-1)**(first_exp*second_exp)
                        circuit = sign*Mul(*new_args)
                        circ_array = circuit.args
                        changes = True
                        break
    return circuit


#-----------------------------------------------------------------------------
# Utility functions
#-----------------------------------------------------------------------------

def test_operator():
    a = Operator('A')
    b = Operator('B', Symbol('t'), S(1)/2)
    inv = a.inv()
    f = Function('f')
    x = symbols('x')
    d = DifferentialOperator(Derivative(f(x), x), f(x))
    op = OuterProduct(Ket(), Bra())
    assert str(a) == 'A'
    assert pretty(a) == 'A'
    assert upretty(a) == u'A'
    assert latex(a) == 'A'
    sT(a, "Operator(Symbol('A'))")
    assert str(inv) == 'A**(-1)'
    ascii_str = \
"""\
 -1\n\
A  \
"""
    ucode_str = \
u("""\
 -1\n\
A  \
""")
    assert pretty(inv) == ascii_str
    assert upretty(inv) == ucode_str
    assert latex(inv) == r'A^{-1}'
    sT(inv, "Pow(Operator(Symbol('A')), Integer(-1))")
    assert str(d) == 'DifferentialOperator(Derivative(f(x), x),f(x))'
    ascii_str = \
"""\
                    /d            \\\n\
DifferentialOperator|--(f(x)),f(x)|\n\
                    \dx           /\
"""
    ucode_str = \
u("""\
                    ?d            ?\n\
DifferentialOperator???(f(x)),f(x)?\n\
                    ?dx           ?\
""")
    assert pretty(d) == ascii_str
    assert upretty(d) == ucode_str
    assert latex(d) == \
        r'DifferentialOperator\left(\frac{d}{d x} f{\left (x \right )},f{\left (x \right )}\right)'
    sT(d, "DifferentialOperator(Derivative(Function('f')(Symbol('x')), Symbol('x')),Function('f')(Symbol('x')))")
    assert str(b) == 'Operator(B,t,1/2)'
    assert pretty(b) == 'Operator(B,t,1/2)'
    assert upretty(b) == u'Operator(B,t,1/2)'
    assert latex(b) == r'Operator\left(B,t,\frac{1}{2}\right)'
    sT(b, "Operator(Symbol('B'),Symbol('t'),Rational(1, 2))")
    assert str(op) == '|psi><psi|'
    assert pretty(op) == '|psi><psi|'
    assert upretty(op) == u'??????'
    assert latex(op) == r'{\left|\psi\right\rangle }{\left\langle \psi\right|}'
    sT(op, "OuterProduct(Ket(Symbol('psi')),Bra(Symbol('psi')))")

def tensor_product_simp(e, **hints):
    """Try to simplify and combine TensorProducts.

    In general this will try to pull expressions inside of ``TensorProducts``.
    It currently only works for relatively simple cases where the products have
    only scalars, raw ``TensorProducts``, not ``Add``, ``Pow``, ``Commutators``
    of ``TensorProducts``. It is best to see what it does by showing examples.

    Examples
    ========

    >>> from sympy.physics.quantum import tensor_product_simp
    >>> from sympy.physics.quantum import TensorProduct
    >>> from sympy import Symbol
    >>> A = Symbol('A',commutative=False)
    >>> B = Symbol('B',commutative=False)
    >>> C = Symbol('C',commutative=False)
    >>> D = Symbol('D',commutative=False)

    First see what happens to products of tensor products:

    >>> e = TensorProduct(A,B)*TensorProduct(C,D)
    >>> e
    AxB*CxD
    >>> tensor_product_simp(e)
    (A*C)x(B*D)

    This is the core logic of this function, and it works inside, powers, sums,
    commutators and anticommutators as well:

    >>> tensor_product_simp(e**2)
    (A*C)x(B*D)**2

    """
    if isinstance(e, Add):
        return Add(*[tensor_product_simp(arg) for arg in e.args])
    elif isinstance(e, Pow):
        return tensor_product_simp(e.base) ** e.exp
    elif isinstance(e, Mul):
        return tensor_product_simp_Mul(e)
    elif isinstance(e, Commutator):
        return Commutator(*[tensor_product_simp(arg) for arg in e.args])
    elif isinstance(e, AntiCommutator):
        return AntiCommutator(*[tensor_product_simp(arg) for arg in e.args])
    else:
        return e

def qsimplify_pauli(e):
    """
    Simplify an expression that includes products of pauli operators.

    Parameters
    ==========

    e : expression
        An expression that contains products of Pauli operators that is
        to be simplified.

    Examples
    ========

    >>> from sympy.physics.quantum.pauli import SigmaX, SigmaY
    >>> from sympy.physics.quantum.pauli import qsimplify_pauli
    >>> sx, sy = SigmaX(), SigmaY()
    >>> sx * sy
    SigmaX()*SigmaY()
    >>> qsimplify_pauli(sx * sy)
    I*SigmaZ()
    """
    if isinstance(e, Operator):
        return e

    if isinstance(e, (Add, Pow, exp)):
        t = type(e)
        return t(*(qsimplify_pauli(arg) for arg in e.args))

    if isinstance(e, Mul):

        c, nc = e.args_cnc()

        nc_s = []
        while nc:
            curr = nc.pop(0)

            while (len(nc) and
                   isinstance(curr, SigmaOpBase) and
                   isinstance(nc[0], SigmaOpBase) and
                   curr.name == nc[0].name):

                x = nc.pop(0)
                y = _qsimplify_pauli_product(curr, x)
                c1, nc1 = y.args_cnc()
                curr = Mul(*nc1)
                c = c + c1

            nc_s.append(curr)

        return Mul(*c) * Mul(*nc_s)

    return e

def find_substitutions(integrand, symbol, u_var):
    results = []

    def test_subterm(u, u_diff):
        substituted = integrand / u_diff
        if symbol not in substituted.free_symbols:
            # replaced everything already
            return False

        substituted = substituted.subs(u, u_var).cancel()
        if symbol not in substituted.free_symbols:
            return substituted.as_independent(u_var, as_Add=False)

        return False

    def possible_subterms(term):
        if isinstance(term, (TrigonometricFunction,
                             sympy.asin, sympy.acos, sympy.atan,
                             sympy.exp, sympy.log, sympy.Heaviside)):
            return [term.args[0]]
        elif isinstance(term, sympy.Mul):
            r = []
            for u in term.args:
                r.append(u)
                r.extend(possible_subterms(u))
            return r
        elif isinstance(term, sympy.Pow):
            if term.args[1].is_constant(symbol):
                return [term.args[0]]
            elif term.args[0].is_constant(symbol):
                return [term.args[1]]
        elif isinstance(term, sympy.Add):
            r = []
            for arg in term.args:
                r.append(arg)
                r.extend(possible_subterms(arg))
            return r
        return []

    for u in possible_subterms(integrand):
        if u == symbol:
            continue
        u_diff = manual_diff(u, symbol)
        new_integrand = test_subterm(u, u_diff)
        if new_integrand is not False:
            constant, new_integrand = new_integrand
            if new_integrand == integrand.subs(symbol, u_var):
                continue
            substitution = (u, constant, new_integrand)
            if substitution not in results:
                results.append(substitution)

    return results