Python simplify.powsimp函数代码示例

def mrv(e, x):
    "Returns a python set of  most rapidly varying (mrv) subexpressions of 'e'"
    e = powsimp(e, deep=True, combine='exp')
    assert isinstance(e, Basic)
    if not e.has(x):
        return set([])
    elif e == x:
        return set([x])
    elif e.is_Mul:
        a, b = e.as_two_terms()
        return mrv_max(mrv(a,x), mrv(b,x), x)
    elif e.is_Add:
        a, b = e.as_two_terms()
        return mrv_max(mrv(a,x), mrv(b,x), x)
    elif e.is_Pow:
        if e.exp.has(x):
            return mrv(exp(e.exp * log(e.base)), x)
        else:
            return mrv(e.base, x)
    elif e.func is log:
        return mrv(e.args[0], x)
    elif e.func is exp:
        if limitinf(e.args[0], x) in [oo,-oo]:
            return mrv_max(set([e]), mrv(e.args[0], x), x)
        else:
            return mrv(e.args[0], x)
    elif e.is_Function:
        if len(e.args) == 1:
            return mrv(e.args[0], x)
        #only functions of 1 argument currently implemented
        raise NotImplementedError("Functions with more arguments: '%s'" % e)
    raise NotImplementedError("Don't know how to calculate the mrv of '%s'" % e)

    def doit(self):
        prod = self._eval_product()

        if prod is not None:
            return powsimp(prod)
        else:
            return self

def mrv(e, x):
    """Returns a SubsSet of most rapidly varying (mrv) subexpressions of 'e',
       and e rewritten in terms of these"""
    e = powsimp(e, deep=True, combine='exp')
    assert isinstance(e, Basic)
    if not e.has(x):
        return SubsSet(), e
    elif e == x:
        s = SubsSet()
        return s, s[x]
    elif e.is_Mul or e.is_Add:
        i, d = e.as_independent(x)  # throw away x-independent terms
        if d.func != e.func:
            s, expr = mrv(d, x)
            return s, e.func(i, expr)
        a, b = d.as_two_terms()
        s1, e1 = mrv(a, x)
        s2, e2 = mrv(b, x)
        return mrv_max1(s1, s2, e.func(i, e1, e2), x)
    elif e.is_Pow:
        b, e = e.as_base_exp()
        if e.has(x):
            return mrv(exp(e * log(b)), x)
        else:
            s, expr = mrv(b, x)
            return s, expr**e
    elif e.func is log:
        s, expr = mrv(e.args[0], x)
        return s, log(expr)
    elif e.func is exp:
        # We know from the theory of this algorithm that exp(log(...)) may always
        # be simplified here, and doing so is vital for termination.
        if e.args[0].func is log:
            return mrv(e.args[0].args[0], x)
        if limitinf(e.args[0], x).is_unbounded:
            s1 = SubsSet()
            e1 = s1[e]
            s2, e2 = mrv(e.args[0], x)
            su = s1.union(s2)[0]
            su.rewrites[e1] = exp(e2)
            return mrv_max3(s1, e1, s2, exp(e2), su, e1, x)
        else:
            s, expr = mrv(e.args[0], x)
            return s, exp(expr)
    elif e.is_Function:
        l = [mrv(a, x) for a in e.args]
        l2 = [s for (s, _) in l if s != SubsSet()]
        if len(l2) != 1:
            # e.g. something like BesselJ(x, x)
            raise NotImplementedError("MRV set computation for functions in"
                                      " several variables not implemented.")
        s, ss = l2[0], SubsSet()
        args = [ss.do_subs(x[1]) for x in l]
        return s, e.func(*args)
    elif e.is_Derivative:
        raise NotImplementedError("MRV set computation for derviatives"
                                  " not implemented yet.")
        return mrv(e.args[0], x)
    raise NotImplementedError(
        "Don't know how to calculate the mrv of '%s'" % e)

    def doit(self, **hints):
        f = g = self.function
        for index, limit in enumerate(self.limits):
            i, a, b = limit
            dif = b - a
            if dif.is_Integer and dif < 0:
                return 1                
            g = self._eval_product(f, (i, a, b))
            if g is None:
                return Product(powsimp(f), *self.limits[index:])
            else:
                f = g

        if hints.get('deep', True):
            return f.doit(**hints)
        else:
            return powsimp(f)

def rewrite(e, Omega, x, wsym):
    """e(x) ... the function
    Omega ... the mrv set
    wsym ... the symbol which is going to be used for w

    Returns the rewritten e in terms of w and log(w). See test_rewrite1()
    for examples and correct results.
    """
    assert isinstance(Omega, SubsSet)
    assert len(Omega) != 0
    #all items in Omega must be exponentials
    for t in Omega.keys():
        assert t.func is exp
    rewrites = Omega.rewrites
    Omega = Omega.items()

    nodes = build_expression_tree(Omega, rewrites)
    Omega.sort(key=lambda x: nodes[x[1]].ht(), reverse=True)

    g, _ = Omega[-1] #g is going to be the "w" - the simplest one in the mrv set
    sig = sign(g.args[0], x)
    if sig == 1:
        wsym = 1/wsym #if g goes to oo, substitute 1/w
    elif sig != -1:
        raise NotImplementedError('Result depends on the sign of %s' % sig)
    #O2 is a list, which results by rewriting each item in Omega using "w"
    O2 = []
    for f, var in Omega:
        c = limitinf(f.args[0]/g.args[0], x)
        arg = f.args[0]
        if var in rewrites:
            assert rewrites[var].func is exp
            arg = rewrites[var].args[0]
        O2.append((var, exp((arg - c*g.args[0]).expand())*wsym**c))

    #Remember that Omega contains subexpressions of "e". So now we find
    #them in "e" and substitute them for our rewriting, stored in O2

    # the following powsimp is necessary to automatically combine exponentials,
    # so that the .subs() below succeeds:
    # TODO this should not be necessary
    f = powsimp(e, deep=True, combine='exp')
    for a, b in O2:
        f = f.subs(a, b)

    for _, var in Omega:
        assert not f.has(var)

    #finally compute the logarithm of w (logw).
    logw = g.args[0]
    if sig == 1:
        logw = -logw     #log(w)->log(1/w)=-log(w)

    return f, logw

    def doit(self, **hints):
        f = self.function

        for index, limit in enumerate(self.limits):
            i, a, b = limit
            dif = b - a
            if dif.is_Integer and dif < 0:
                a, b = b + 1, a - 1
                f = 1 / f
            if isinstance(i, Idx):
                i = i.label

            g = self._eval_product(f, (i, a, b))
            if g in (None, S.NaN):
                return self.func(powsimp(f), *self.limits[index:])
            else:
                f = g

        if hints.get('deep', True):
            return f.doit(**hints)
        else:
            return powsimp(f)

    def doit(self, **hints):
        term = self.term
        lower = self.lower
        upper = self.upper
        if hints.get('deep', True):
            term = term.doit(**hints)
            lower = lower.doit(**hints)
            upper = upper.doit(**hints)

        prod = self._eval_product(lower, upper, term)

        if prod is not None:
            return powsimp(prod)
        else:
            return self

def rewrite(e, Omega, x, wsym):
    """e(x) ... the function
    Omega ... the mrv set
    wsym ... the symbol which is going to be used for w

    Returns the rewritten e in terms of w and log(w). See test_rewrite1()
    for examples and correct results.
    """
    assert isinstance(Omega, set)
    assert len(Omega) != 0
    # all items in Omega must be exponentials
    for t in Omega:
        assert t.func is exp

    def cmpfunc(a, b):
        return -cmp(len(mrv(a, x)), len(mrv(b, x)))

    # sort Omega (mrv set) from the most complicated to the simplest ones
    # the complexity of "a" from Omega: the length of the mrv set of "a"
    Omega = list(Omega)
    Omega.sort(cmp=cmpfunc)
    g = Omega[-1]  # g is going to be the "w" - the simplest one in the mrv set
    sig = sign(g.args[0], x) == 1
    if sig:
        wsym = 1 / wsym  # if g goes to oo, substitute 1/w
    # O2 is a list, which results by rewriting each item in Omega using "w"
    O2 = []
    for f in Omega:
        c = mrv_leadterm(f.args[0] / g.args[0], x)
        # the c is a constant, because both f and g are from Omega:
        assert c[1] == 0
        O2.append(exp((f.args[0] - c[0] * g.args[0]).expand()) * wsym ** c[0])

    # Remember that Omega contains subexpressions of "e". So now we find
    # them in "e" and substitute them for our rewriting, stored in O2

    # the following powsimp is necessary to automatically combine exponentials,
    # so that the .subs() below succeeds:
    f = powsimp(e, deep=True, combine="exp")
    for a, b in zip(Omega, O2):
        f = f.subs(a, b)

    # finally compute the logarithm of w (logw).
    logw = g.args[0]
    if sig:
        logw = -logw  # log(w)->log(1/w)=-log(w)

    return f, logw

    def doit(self, **hints):
        term = self.term
        if term == 0:
            return S.Zero
        elif term == 1:
            return S.One
        lower = self.lower
        upper = self.upper
        if hints.get('deep', True):
            term = term.doit(**hints)
            lower = lower.doit(**hints)
            upper = upper.doit(**hints)
        dif = upper - lower
        if dif.is_Number and dif < 0:
            upper, lower = lower, upper

        prod = self._eval_product(lower, upper, term)

        if prod is not None:
            return powsimp(prod)
        else:
            return self

def mrv(e, x):
    "Returns a python set of  most rapidly varying (mrv) subexpressions of 'e'"
    e = powsimp(e, deep=True, combine="exp")
    assert isinstance(e, Basic)
    if not e.has(x):
        return set([])
    elif e == x:
        return set([x])
    elif e.is_Mul or e.is_Add:
        while 1:
            i, d = e.as_independent(x)  # throw away x-independent terms
            if d.func != e.func and (d.is_Add or d.is_Mul):
                e = d
                continue
            break
        if d.func != e.func:
            return mrv(d, x)
        a, b = d.as_two_terms()
        return mrv_max(mrv(a, x), mrv(b, x), x)
    elif e.is_Pow:
        b, e = e.as_base_exp()
        if e.has(x):
            return mrv(exp(e * log(b)), x)
        else:
            return mrv(b, x)
    elif e.func is log:
        return mrv(e.args[0], x)
    elif e.func is exp:
        if limitinf(e.args[0], x).is_unbounded:
            return mrv_max(set([e]), mrv(e.args[0], x), x)
        else:
            return mrv(e.args[0], x)
    elif e.is_Function:
        return reduce(lambda a, b: mrv_max(a, b, x), [mrv(a, x) for a in e.args])
    elif e.is_Derivative:
        return mrv(e.args[0], x)
    raise NotImplementedError("Don't know how to calculate the mrv of '%s'" % e)

def _quintic_simplify(expr):
    expr = powsimp(expr)
    expr = cancel(expr)
    return together(expr)

#.........这里部分代码省略.........

       None is returned if checksol() could not conclude.

       flags:
           'numerical=True (default)'
               do a fast numerical check if f has only one symbol.
           'minimal=True (default is False)'
               a very fast, minimal testing.
           'warning=True (default is False)'
               print a warning if checksol() could not conclude.
           'simplified=True (default)'
               solution should be simplified before substituting into function
               and function should be simplified after making substitution.
           'force=True (default is False)'
               make positive all symbols without assumptions regarding sign.
    """

    if sol is not None:
        sol = {symbol: sol}
    elif isinstance(symbol, dict):
        sol = symbol
    else:
        msg = 'Expecting sym, val or {sym: val}, None but got %s, %s'
        raise ValueError(msg % (symbol, sol))

    if hasattr(f, '__iter__') and hasattr(f, '__len__'):
        if not f:
            raise ValueError('no functions to check')
        rv = set()
        for fi in f:
            check = checksol(fi, sol, **flags)
            if check is False:
                return False
            rv.add(check)
        if None in rv: # rv might contain True and/or None
            return None
        assert len(rv) == 1 # True
        return True

    if isinstance(f, Poly):
        f = f.as_expr()
    elif isinstance(f, Equality):
        f = f.lhs - f.rhs

    if not f:
        return True

    if not f.has(*sol.keys()):
        return False

    attempt = -1
    numerical = flags.get('numerical', True)
    while 1:
        attempt += 1
        if attempt == 0:
            val = f.subs(sol)
        elif attempt == 1:
            if not val.atoms(Symbol) and numerical:
                # val is a constant, so a fast numerical test may suffice
                if val not in [S.Infinity, S.NegativeInfinity]:
                    # issue 2088 shows that +/-oo chops to 0
                    val = val.evalf(36).n(30, chop=True)
        elif attempt == 2:
            if flags.get('minimal', False):
                return
            # the flag 'simplified=False' is used in solve to avoid
            # simplifying the solution. So if it is set to False there
            # the simplification will not be attempted here, either. But
            # if the simplification is done here then the flag should be
            # set to False so it isn't done again there.
            # FIXME: this can't work, since `flags` is not passed to
            # `checksol()` as a dict, but as keywords.
            # So, any modification to `flags` here will be lost when returning
            # from `checksol()`.
            if flags.get('simplified', True):
                for k in sol:
                    sol[k] = simplify(sympify(sol[k]))
                flags['simplified'] = False
                val = simplify(f.subs(sol))
            if flags.get('force', False):
                val = posify(val)[0]
        elif attempt == 3:
            val = powsimp(val)
        elif attempt == 4:
            val = cancel(val)
        elif attempt == 5:
            val = val.expand()
        elif attempt == 6:
            val = together(val)
        elif attempt == 7:
            val = powsimp(val)
        else:
            break
        if val.is_zero:
            return True
        elif attempt > 0 and numerical and val.is_nonzero:
            return False

    if flags.get('warning', False):
        print("\n\tWarning: could not verify solution %s." % sol)

def rewrite(e, Omega, x, wsym):
    """e(x) ... the function
    Omega ... the mrv set
    wsym ... the symbol which is going to be used for w

    Returns the rewritten e in terms of w and log(w). See test_rewrite1()
    for examples and correct results.
    """
    from sympy import ilcm

    if not isinstance(Omega, SubsSet):
        raise TypeError("Omega should be an instance of SubsSet")
    if len(Omega) == 0:
        raise ValueError("Length can not be 0")
    # all items in Omega must be exponentials
    for t in Omega.keys():
        if not t.func is exp:
            raise ValueError("Value should be exp")
    rewrites = Omega.rewrites
    Omega = list(Omega.items())

    nodes = build_expression_tree(Omega, rewrites)
    Omega.sort(key=lambda x: nodes[x[1]].ht(), reverse=True)

    # make sure we know the sign of each exp() term; after the loop,
    # g is going to be the "w" - the simplest one in the mrv set
    for g, _ in Omega:
        sig = sign(g.args[0], x)
        if sig != 1 and sig != -1:
            raise NotImplementedError("Result depends on the sign of %s" % sig)
    if sig == 1:
        wsym = 1 / wsym  # if g goes to oo, substitute 1/w
    # O2 is a list, which results by rewriting each item in Omega using "w"
    O2 = []
    denominators = []
    for f, var in Omega:
        c = limitinf(f.args[0] / g.args[0], x)
        if c.is_Rational:
            denominators.append(c.q)
        arg = f.args[0]
        if var in rewrites:
            if not rewrites[var].func is exp:
                raise ValueError("Value should be exp")
            arg = rewrites[var].args[0]
        O2.append((var, exp((arg - c * g.args[0]).expand()) * wsym ** c))

    # Remember that Omega contains subexpressions of "e". So now we find
    # them in "e" and substitute them for our rewriting, stored in O2

    # the following powsimp is necessary to automatically combine exponentials,
    # so that the .subs() below succeeds:
    # TODO this should not be necessary
    f = powsimp(e, deep=True, combine="exp")
    for a, b in O2:
        f = f.subs(a, b)

    for _, var in Omega:
        assert not f.has(var)

    # finally compute the logarithm of w (logw).
    logw = g.args[0]
    if sig == 1:
        logw = -logw  # log(w)->log(1/w)=-log(w)

    # Some parts of sympy have difficulty computing series expansions with
    # non-integral exponents. The following heuristic improves the situation:
    exponent = reduce(ilcm, denominators, 1)
    f = f.subs(wsym, wsym ** exponent)
    logw /= exponent

    return f, logw