def test_simplify_ratio():
    # roots of x**3-3*x+5
    roots = ['(5/2 + 21**(1/2)/2)**(1/3)*(1/2 - I*3**(1/2)/2)'
             ' + 1/((1/2 - I*3**(1/2)/2)*(5/2 + 21**(1/2)/2)**(1/3))',
             '(5/2 + 21**(1/2)/2)**(1/3)*(1/2 + I*3**(1/2)/2)'
             ' + 1/((1/2 + I*3**(1/2)/2)*(5/2 + 21**(1/2)/2)**(1/3))',
             '-1/(5/2 + 21**(1/2)/2)**(1/3) - (5/2 + 21**(1/2)/2)**(1/3)']
    for r in roots:
        r = S(r)
        assert count_ops(simplify(r, ratio=1)) <= count_ops(r)
        # If ratio=oo, simplify() is always applied:
        assert simplify(r, ratio=oo) is not r

def expression_complexity(expr, complexity=None):
    '''
    Returns the complexity of an expression (either string or sympy)

    The complexity is defined as 1 for each arithmetic operation except divide which is 2,
    and all other operations are 20. This can be overridden using the complexity
    argument.

    Note: calling this on a statement rather than an expression is likely to lead to errors.

    Parameters
    ----------
    expr: `sympy.Expr` or str
        The expression.
    complexity: None or dict (optional)
        A dictionary mapping expression names to their complexity, to overwrite default behaviour.

    Returns
    -------
    complexity: int
        The complexity of the expression.
    '''
    subs = {'ADD':1, 'DIV':2, 'MUL':1, 'SUB':1}
    if complexity is not None:
        subs.update(complexity)
    ops = sympy.count_ops(expr, visual=True)
    for atom in ops.atoms():
        if hasattr(atom, 'name'):
            subs[atom.name] = 20 # unknown operations assumed to have a large cost
    return ops.evalf(subs=subs)

def subs_matrix_verbose(A_expr,subs,simultaneous=False):

    print "flashlight.sympy: subs_matrix_verbose(...) begin..."
    A_subs_expr = sympy.Matrix.zeros(A_expr.rows,A_expr.cols)
    for r in range(A_expr.rows):
        for c in range(A_expr.cols):
            print "    ",r,c,len(subs),sympy.count_ops(A_expr[r,c])
            A_subs_expr[r,c] = A_expr[r,c].subs(subs,simultaneous=simultaneous)
    print "flashlight.sympy: subs_matrix_verbose(...) end."
    return A_subs_expr

def test_issue_2827_trigsimp_methods():
    measure1 = lambda expr: len(str(expr))
    measure2 = lambda expr: -count_ops(expr)
                                       # Return the most complicated result
    expr = (x + 1)/(x + sin(x)**2 + cos(x)**2)
    ans = Matrix([1])
    M = Matrix([expr])
    assert trigsimp(M, method='fu', measure=measure1) == ans
    assert trigsimp(M, method='fu', measure=measure2) != ans
    # all methods should work with Basic expressions even if they
    # aren't Expr
    M = Matrix.eye(1)
    assert all(trigsimp(M, method=m) == M for m in
        'fu matching groebner old'.split())
    # watch for E in exptrigsimp, not only exp()
    eq = 1/sqrt(E) + E
    assert exptrigsimp(eq) == eq

def TestTraverseExpr(tmpExp):
    # Inspired from sympy-master/sympy/core/expr.py:

    assert isinstance(tmpExp, sympy.Expr) #fraction)  #Equality)
    assert tmpExp.is_Mul #fraction)  #Equality)
    print("Number of operators on the rhs exp: %s" % \
                sympy.count_ops(tmpExp, visual=True))
    args = tmpExp.as_ordered_factors() #order=order)
    print("args = %s" % args)
    print("args[1] = %s" % args[1])
    assert args[1].is_Pow
    args1 = args[1].args
    print args1
    print args1[0]
    assert args1[1] == -1
    #print("args1 (the denominator factors) = %s" % str(args1))
    #print type(args[1])
    assert args[1].is_Mul

    def render_node(self, node):
        expr = NodeRenderer(use_vectorisation_idx=False).render_node(node)

        if is_scalar_expression(expr, self.variables) and not has_non_float(expr,
                                                                            self.variables):
            if expr in self.optimisations:
                name = self.optimisations[expr]
            else:
                # Do not pull out very simple expressions (including constants
                # and numbers)
                sympy_expr = str_to_sympy(expr)
                if sympy.count_ops(sympy_expr, visual=False) < 2:
                    return expr
                self.n += 1
                name = '_lio_const_'+str(self.n)
                self.optimisations[expr] = name
            return name
        else:
            return NodeRenderer.render_node(self, node)

def expression_complexity(expr, complexity=None):
    '''
    Returns the complexity of an expression (either string or sympy)

    The complexity is defined as 1 for each arithmetic operation except divide which is 2,
    and all other operations are 20. This can be overridden using the complexity
    argument.

    Note: calling this on a statement rather than an expression is likely to lead to errors.

    Parameters
    ----------
    expr: `sympy.Expr` or str
        The expression.
    complexity: None or dict (optional)
        A dictionary mapping expression names to their complexity, to overwrite default behaviour.

    Returns
    -------
    complexity: int
        The complexity of the expression.
    '''
    if isinstance(expr, str):
        # we do this because sympy.count_ops doesn't handle inequalities (TODO: handle sympy as well str)
        for op in ['<=', '>=', '==', '<', '>']:
            expr = expr.replace(op, '+')
        # work around bug with rand() and randn() (TODO: improve this)
        expr = expr.replace('rand()', 'rand(0)')
        expr = expr.replace('randn()', 'randn(0)')
    subs = {'ADD':1, 'DIV':2, 'MUL':1, 'SUB':1}
    if complexity is not None:
        subs.update(complexity)
    ops = sympy.count_ops(expr, visual=True)
    for atom in ops.atoms():
        if hasattr(atom, 'name'):
            subs[atom.name] = 20 # unknown operations assumed to have a large cost
    return ops.evalf(subs=subs)

def Read(fileName):
    fin = open(fileName)

    while (True):
        str = fin.readline()
        if str:
            #pass
            str = str.rstrip(rstripStr)
        else:
            break

        if str.startswith("#"):
            continue

        for index, v in enumerate(defVars):
            cmpStr = v + " = "
            #print("cmpStr = %s" % cmpStr)
            #" - v"
            if str.startswith(cmpStr):
                #TestTraverseExpr(tmpExp)

                tmpExp = Parse(str[len(cmpStr) : ])
                lhs = sympy.Symbol(v)
                print (tmpExp - lhs)
                lhs, rhs = PreprocessEq(lhs, tmpExp)

                #assert False

                #eqList[index] = Parse(str[len(cmpStr) : ] + " - " + v)
                eqList[index] = rhs - lhs

                #print("eqList[index] = %s" % str(eqList[index]))
                print("eqList[%d] = %s" % (index, eqList[index]))
                #print(separatevars(expr=eqList[index], symbols=[xdst, xsrc]))
                print("Number of operators: %s" % \
                        sympy.count_ops(eqList[index], visual=True))
                GenCode(eqList[index])

 def count(val):
     return count_ops(val, visual=True)

 def count_ops(self, visual=None):
     """wrapper for count_ops that returns the operation count."""
     from sympy import count_ops
     return count_ops(self, visual)
     return sum(a.count_ops(visual) for a in self.args)

def _mellin_transform(f, x, s_, integrator=_default_integrator, simplify=True):
    """ Backend function to compute mellin transforms. """
    from sympy import re, Max, Min, count_ops
    # We use a fresh dummy, because assumptions on s might drop conditions on
    # convergence of the integral.
    s = _dummy('s', 'mellin-transform', f)
    F = integrator(x**(s-1) * f, x)

    if not F.has(Integral):
        return _simplify(F.subs(s, s_), simplify), (-oo, oo), True

    if not F.is_Piecewise:
        raise IntegralTransformError('Mellin', f, 'could not compute integral')

    F, cond = F.args[0]
    if F.has(Integral):
        raise IntegralTransformError('Mellin', f, 'integral in unexpected form')

    def process_conds(cond):
        """
        Turn ``cond`` into a strip (a, b), and auxiliary conditions.
        """
        a = -oo
        b = oo
        aux = True
        conds = conjuncts(to_cnf(cond))
        t = Dummy('t', real=True)
        for c in conds:
            a_ = oo
            b_ = -oo
            aux_ = []
            for d in disjuncts(c):
                d_ = d.replace(re, lambda x: x.as_real_imag()[0]).subs(re(s), t)
                if not d.is_Relational or (d.rel_op != '<' and d.rel_op != '<=') \
                   or d_.has(s) or not d_.has(t):
                    aux_ += [d]
                    continue
                soln = _solve_inequality(d_, t)
                if not soln.is_Relational or \
                   (soln.rel_op != '<' and soln.rel_op != '<='):
                    aux_ += [d]
                    continue
                if soln.lhs == t:
                    b_ = Max(soln.rhs, b_)
                else:
                    a_ = Min(soln.lhs, a_)
            if a_ != oo and a_ != b:
                a = Max(a_, a)
            elif b_ != -oo and b_ != a:
                b = Min(b_, b)
            else:
                aux = And(aux, Or(*aux_))
        return a, b, aux

    conds = [process_conds(c) for c in disjuncts(cond)]
    conds = filter(lambda x: x[2] is not False, conds)
    conds.sort(key=lambda x: (x[0]-x[1], count_ops(x[2])))

    if not conds:
        raise IntegralTransformError('Mellin', f, 'no convergence found')

    a, b, aux = conds[0]
    return _simplify(F.subs(s, s_), simplify), (a, b), aux

from sympy import simplify, cos, sin, trigsimp, cancel
from sympy import sqrt, count_ops, oo, symbols, log
from sympy.abc import x, y

expr = (2*x + 3*x**2)/(4*x*sin(y)**2 + 2*x*cos(y)**2)
expr
simplify(expr)

trigsimp(expr)
cancel(_)

root = 4/(sqrt(2)+3)
simplify(root, ratio=1) == root
count_ops(simplify(root, ratio=oo)) > count_ops(root)
x, y = symbols('x y', positive=True)
expr2 = log(x) + log(y) + log(x)*log(1/y)

expr3 = simplify(expr2)
expr3
count_ops(expr2)
count_ops(expr3)
print(count_ops(expr2, visual=True))
print(count_ops(expr3, visual=True))

	def intersect(self, other, **flags):
		if not isinstance(other, Circle):
			raise BaryException()
		if sympy.count_ops(self.eqn()) > sympy.count_ops(other.eqn()):
			raise BaryException()
		return [Point(*sol) for sol in hsolve([self.eqn(), ((self._eqn - other._eqn) / (x + y + z)).simplify()], x, y, z, **flags)]

 def count_ops(self, visual=None):
     """wrapper for count_ops that returns the operation count."""
     from sympy import count_ops
     return count_ops(self, visual)

#.........这里部分代码省略.........
    tpf = (xFj - xPj)/(xTj - xPj)

    xP = [(((xF[i]/ppf)*(beta[i]**(NT+1) - 1))/(beta[i]**(NT+1) - beta[i]**(-NP))) \
                                                                            for i in r]
    xT = [(((xF[i]/tpf)*(1 - beta[i]**(-NP)))/(beta[i]**(NT+1) - beta[i]**(-NP))) \
                                                                            for i in r]
    rfeed = xFj / xF[k]
    rprod = xPj / xP[k]
    rtail = xTj / xT[k]

    # setup constraint equations
    numer = [ppf*xP[i]*log(rprod) + tpf*xT[i]*log(rtail) - xF[i]*log(rfeed) for i in r]
    denom = [log(beta[j]) * ((beta[i] - 1.0)/(beta[i] + 1.0)) for i in r]
    LoverF = sum([n/d for n, d in zip(numer, denom)])
    SWUoverF = -1.0 * sum(numer)
    SWUoverP = SWUoverF / ppf

    prod_constraint = (xPj/xFj)*ppf - (beta[j]**(NT+1) - 1)/\
                      (beta[j]**(NT+1) - beta[j]**(-NP))
    tail_constraint = (xTj/xFj)*(sum(xT)) - (1 - beta[j]**(-NP))/\
                      (beta[j]**(NT+1) - beta[j]**(-NP))
    #xp_constraint = 1.0 - sum(xP)
    #xf_constraint = 1.0 - sum(xF)
    #xt_constraint = 1.0 - sum(xT)

    # This is NT(NP,...) and is correct!
    #nt_closed = solve(prod_constraint, NT)[0] 

    # However, this is NT(NP,...) rewritten (by hand) to minimize the number of NP 
    # and M* instances in the expression.  Luckily this is only depends on the key 
    # component and remains general no matter the number of components.
    nt_closed = (-MW[0]*log(alpha) + Mstar*log(alpha) + log(xTj) + log((-1.0 + xPj/\
        xF[0])/(xPj - xTj)) - log(alpha**(NP*(MW[0] - Mstar))*(xF[0]*xPj - xPj*xTj)/\
        (-xF[0]*xPj + xF[0]*xTj) + 1))/((MW[0] - Mstar)*log(alpha))

    # new expression for normalized flow rate
    # NOTE: not needed, solved below
    #loverf = LoverF.xreplace({NT: nt_closed})

    # Define the constraint equation with which to solve NP. This is chosen such to 
    # minimize the number of ops in the derivatives (and thus np_closed).  Other, 
    # more verbose possibilities are commented out.
    #np_constraint = (xP[j]/sum(xP) - xPj).xreplace({NT: nt_closed})
    #np_constraint = (xP[j]- sum(xP)*xPj).xreplace({NT: nt_closed})
    #np_constraint = (xT[j]/sum(xT) - xTj).xreplace({NT: nt_closed})
    np_constraint = (xT[j] - sum(xT)*xTj).xreplace({NT: nt_closed})

    # get closed form approximation of NP via symbolic derivatives
    stat = _aggstatus(stat, "  order-{0} NP approximation".format(nporder), aggstat)
    d0NP = np_constraint.xreplace({NP: NP0})
    d1NP = diff(np_constraint, NP, 1).xreplace({NP: NP0})
    if 1 == nporder:
        np_closed = NP0 - d1NP / d0NP
    elif 2 == nporder:
        d2NP = diff(np_constraint, NP, 2).xreplace({NP: NP0})/2.0
        # taylor series polynomial coefficients, grouped by order
        # f(x) = ax**2 + bx + c
        a = d2NP
        b = d1NP - 2*NP0*d2NP
        c = d0NP - NP0*d1NP + NP0*NP0*d2NP
        # quadratic eq. (minus only)
        #np_closed = (-b - sqrt(b**2 - 4*a*c)) / (2*a)
        # However, we need to break up this expr as follows to prevent 
        # a floating point arithmetic bug if b**2 - 4*a*c is very close
        # to zero but happens to be negative.  LAME!!!
        np_2a = 2*a
        np_sqrt_base = b**2 - 4*a*c
        np_closed = (-NP_b - sqrt(NP_sqrt_base)) / (NP_2a)
    else:
        raise ValueError("nporder must be 1 or 2")

    # generate cse for writing out
    msg = "  minimizing ops by eliminating common sub-expressions"
    stat = _aggstatus(stat, msg, aggstat)
    exprstages = [Eq(NP_b, b), Eq(NP_2a, np_2a), 
                  # fix for floating point sqrt() error
                  Eq(NP_sqrt_base, np_sqrt_base), Eq(NP_sqrt_base, Abs(NP_sqrt_base)), 
                  Eq(NP1, np_closed), Eq(NT1, nt_closed).xreplace({NP: NP1})]
    cse_stages = cse(exprstages, numbered_symbols('n'))
    exprothers = [Eq(LpF, LoverF), Eq(PpF, ppf), Eq(TpF, tpf), 
                  Eq(SWUpF, SWUoverF), Eq(SWUpP, SWUoverP)] + \
                 [Eq(*z) for z in zip(xPi, xP)] + [Eq(*z) for z in zip(xTi, xT)]
    exprothers = [e.xreplace({NP: NP1, NT: NT1}) for e in exprothers]
    cse_others = cse(exprothers, numbered_symbols('g'))
    exprops = count_ops(exprstages + exprothers)
    cse_ops = count_ops(cse_stages + cse_others)
    msg = "    reduced {0} ops to {1}".format(exprops, cse_ops)
    stat = _aggstatus(stat, msg, aggstat)

    # create function body
    ccode, repnames = cse_to_c(*cse_stages, indent=6, debug=debug)
    ccode_others, repnames_others = cse_to_c(*cse_others, indent=6, debug=debug)
    ccode += ccode_others
    repnames |= repnames_others

    msg = "  completed in {0:.3G} s".format(time.time() - start_time)
    stat = _aggstatus(stat, msg, aggstat)
    if aggstat:
        print(stat)
    return ccode, repnames, stat