sympy/simplify/simplify.py

   1
   2 from sympy.core import Basic, S, Add, Mul, Pow, Rational, Integer, \
   3         Derivative, Wild, Symbol
   4
   5 from sympy.utilities import make_list, all
   6 from sympy.functions import gamma, exp
   7
   8 from sys import maxint
   9
  10 def fraction(expr, exact=False):
  11     """Returns a pair with expression's numerator and denominator.
  12        If the given expression is not a fraction then this function
  13        will assume that the denominator is equal to one.
  14
  15        This function will not make any attempt to simplify nested
  16        fractions or to do any term rewriting at all.
  17
  18        If only one of the numerator/denominator pair is needed then
  19        use numer(expr) or denom(expr) functions respectively.
  20
  21        >>> from sympy import *
  22        >>> x, y = symbols('x', 'y')
  23
  24        >>> fraction(x/y)
  25        (x, y)
  26        >>> fraction(x)
  27        (x, 1)
  28
  29        >>> fraction(1/y**2)
  30        (1, y**2)
  31
  32        >>> fraction(x*y/2)
  33        (x*y, 2)
  34        >>> fraction(Rational(1, 2))
  35        (1, 2)
  36
  37        This function will also work fine with assumptions:
  38
  39        >>> k = Symbol('k', negative=True)
  40        >>> fraction(x * y**k)
  41        (x, y**(-k))
  42
  43        If we know nothing about sign of some exponent and 'exact'
  44        flag is unset, then structure this exponent's structure will
  45        be analyzed and pretty fraction will be returned:
  46
  47        >>> fraction(2*x**(-y))
  48        (2, x**y)
  49
  50        #>>> fraction(exp(-x))
  51        #(1, exp(x))
  52
  53        >>> fraction(exp(-x), exact=True)
  54        (exp(-x), 1)
  55
  56     """
  57     expr = Basic.sympify(expr)
  58
  59     #XXX this only works sometimes (caching bug?)
  60     if expr == exp(-Symbol("x")) and exact:
  61         return (expr, 1)
  62
  63     numer, denom = [], []
  64
  65     for term in make_list(expr, Mul):
  66         if isinstance(term, Pow):
  67             if term.exp.is_negative:
  68                 if term.exp == Integer(-1):
  69                     denom.append(term.base)
  70                 else:
  71                     denom.append(Pow(term.base, -term.exp))
  72             elif not exact and isinstance(term.exp, Mul):
  73                 coeff, tail = term.exp[0], Mul(*term.exp[1:])#term.exp.getab()
  74
  75                 if isinstance(coeff, Rational) and coeff.is_negative:
  76                     denom.append(Pow(term.base, -term.exp))
  77                 else:
  78                     numer.append(term)
  79             else:
  80                 numer.append(term)
  81         elif isinstance(term, Basic.exp):
  82             if term[0].is_negative:
  83                 denom.append(Basic.exp(-term[0]))
  84             elif not exact and isinstance(term[0], Mul):
  85                 coeff, tail = term[0], Mul(*term[1:])#term.args.getab()
  86
  87                 if isinstance(coeff, Rational) and coeff.is_negative:
  88                     denom.append(Basic.exp(-term[0]))
  89                 else:
  90                     numer.append(term)
  91             else:
  92                 numer.append(term)
  93         elif isinstance(term, Rational):
  94             if term.is_integer:
  95                 numer.append(term)
  96             else:
  97                 numer.append(Rational(term.p))
  98                 denom.append(Rational(term.q))
  99         else:
 100             numer.append(term)
 101
 102     return Mul(*numer), Mul(*denom)
 103
 104 def numer(expr):
 105     return fraction(expr)[0]
 106
 107 def denom(expr):
 108     return fraction(expr)[1]
 109
 110 def fraction_expand(expr):
 111     a, b = fraction(expr)
 112     return a.expand() / b.expand()
 113
 114 def numer_expand(expr):
 115     a, b = fraction(expr)
 116     return a.expand() / b
 117
 118 def denom_expand(expr):
 119     a, b = fraction(expr)
 120     return a / b.expand()
 121
 122 def separate(expr, deep=False):
 123     """Rewrite or separate a power of product to a product of powers
 124        but without any expanding, ie. rewriting products to summations.
 125
 126        >>> from sympy import *
 127        >>> x, y, z = symbols('x', 'y', 'z')
 128
 129        >>> separate((x*y)**2)
 130        x**2*y**2
 131
 132        >>> separate((x*(y*z)**3)**2)
 133        x**2*y**6*z**6
 134
 135        >>> separate((x*sin(x))**y + (x*cos(x))**y)
 136        x**y*cos(x)**y + x**y*sin(x)**y
 137
 138        #>>> separate((exp(x)*exp(y))**x)
 139        #exp(x*y)*exp(x**2)
 140
 141        Notice that summations are left un touched. If this is not the
 142        requested behaviour, apply 'expand' to input expression before:
 143
 144        >>> separate(((x+y)*z)**2)
 145        z**2*(x + y)**2
 146
 147        >>> separate((x*y)**(1+z))
 148        x**(1 + z)*y**(1 + z)
 149
 150     """
 151     expr = Basic.sympify(expr)
 152
 153     if isinstance(expr, Basic.Pow):
 154         terms, expo = [], separate(expr.exp, deep)
 155         #print expr, terms, expo, expr.base
 156
 157         if isinstance(expr.base, Mul):
 158             t = [ separate(Basic.Pow(t,expo), deep) for t in expr.base ]
 159             return Basic.Mul(*t)
 160         elif isinstance(expr.base, Basic.exp):
 161             if deep == True:
 162                 return Basic.exp(separate(expr.base[0], deep)*expo)
 163             else:
 164                 return Basic.exp(expr.base[0]*expo)
 165         else:
 166             return Basic.Pow(separate(expr.base, deep), expo)
 167     elif isinstance(expr, (Basic.Add, Basic.Mul)):
 168         return type(expr)(*[ separate(t, deep) for t in expr ])
 169     elif isinstance(expr, Basic.Function) and deep:
 170         return expr.func(*[ separate(t) for t in expr])
 171     else:
 172         return expr
 173
 174
 175 def together(expr, deep=False):
 176     """Combine together and denest rational functions into a single
 177        fraction. By default the resulting expression is simplified
 178        to reduce the total order of both numerator and denominator
 179        and minimize the number of terms.
 180
 181        Densting is done recursively on fractions level. However this
 182        function will not attempt to rewrite composite objects, like
 183        functions, interior unless 'deep' flag is set.
 184
 185        By definition, 'together' is a complement to 'apart', so
 186        apart(together(expr)) should left expression unhanged.
 187
 188        >>> from sympy import *
 189        >>> x, y, z = symbols('x', 'y', 'z')
 190
 191        You can work with sums of fractions easily. The algorithm
 192        used here will, in an iterative style, collect numerators
 193        and denominator of all expressions involved and perform
 194        needed simplifications:
 195
 196        #>>> together(1/x + 1/y)
 197        #(x + y)/(y*x)
 198
 199        #>>> together(1/x + 1/y + 1/z)
 200        #(z*x + x*y + z*y)/(y*x*z)
 201
 202        >>> together(1/(x*y) + 1/y**2)
 203        1/x*y**(-2)*(x + y)
 204
 205        Or you can just denest multi-level fractional expressions:
 206
 207        >>> together(1/(1 + 1/x))
 208        x/(1 + x)
 209
 210        It also perfect possible to work with symbolic powers or
 211        exponential functions or combinations of both:
 212
 213        >>> together(1/x**y + 1/x**(y-1))
 214        x**(-y)*(1 + x)
 215
 216        #>>> together(1/x**(2*y) + 1/x**(y-z))
 217        #x**(-2*y)*(1 + x**(y + z))
 218
 219        #>>> together(1/exp(x) + 1/(x*exp(x)))
 220        #(1+x)/(x*exp(x))
 221
 222        #>>> together(1/exp(2*x) + 1/(x*exp(3*x)))
 223        #(1+exp(x)*x)/(x*exp(3*x))
 224
 225     """
 226
 227     def _together(expr):
 228
 229         from sympy.core.function import Function
 230
 231         if isinstance(expr, Add):
 232             items, coeffs, basis = [], [], {}
 233
 234             for elem in expr:
 235                 numer, q = fraction(_together(elem))
 236
 237                 denom = {}
 238
 239                 for term in make_list(q.expand(), Mul):
 240                     expo = Integer(1)
 241                     coeff = Integer(1)
 242
 243
 244                     if isinstance(term, Pow):
 245                         if isinstance(term.exp, Rational):
 246                             term, expo = term.base, term.exp
 247                         elif isinstance(term.exp, Mul):
 248                             coeff, tail = term.exp.as_coeff_terms()
 249                             if isinstance(coeff, Rational):
 250                                 tail = Basic.Mul(*tail)
 251                                 term, expo = Pow(term.base, tail), coeff
 252                         coeff = Integer(1)
 253                     elif isinstance(term, Basic.exp):
 254                         if isinstance(term[0], Rational):
 255                             term, expo = Basic.E, term[0]
 256                         elif isinstance(term[0], Mul):
 257                             coeff, tail = term[0].as_coeff_terms()
 258                             if isinstance(coeff, Rational):
 259                                 tail = Basic.Mul(*tail)
 260                                 term, expo = Basic.exp(tail), coeff
 261                         coeff = Integer(1)
 262                     elif isinstance(term, Rational):
 263                         coeff = Integer(term.q)
 264                         term = Integer(term.p)
 265
 266                     if term in denom:
 267                         denom[term] += expo
 268                     else:
 269                         denom[term] = expo
 270
 271                     if term in basis:
 272                         total, maxi = basis[term]
 273
 274                         n_total = total + expo
 275                         n_maxi = max(maxi, expo)
 276
 277                         basis[term] = (n_total, n_maxi)
 278                     else:
 279                         basis[term] = (expo, expo)
 280
 281                     coeffs.append(coeff)
 282                 items.append((numer, denom))
 283
 284             numerator, denominator = [], []
 285
 286             for (term, (total, maxi)) in basis.iteritems():
 287                 basis[term] = (total, total-maxi)
 288
 289                 if isinstance(term, Basic.exp):
 290                     denominator.append(Basic.exp(maxi*term[:]))
 291                 else:
 292                     if isinstance(maxi, Basic.One):
 293                         denominator.append(term)
 294                     else:
 295                         denominator.append(Pow(term, maxi))
 296
 297             from sympy.core.numbers import gcd as int_gcd
 298
 299             if all([ c.is_integer for c in coeffs ]):
 300                 gcds = lambda x, y: int_gcd(int(x), int(y))
 301                 common = Rational(reduce(gcds, coeffs))
 302             else:
 303                 common = Rational(1)
 304
 305             product = reduce(lambda x, y: x*y, coeffs) / common
 306
 307             for ((numer, denom), coeff) in zip(items, coeffs):
 308
 309                 expr, coeff = [], product / (coeff*common)
 310
 311                 for term in basis.iterkeys():
 312                     total, sub = basis[term]
 313
 314                     if term in denom:
 315                         expo = total-denom[term]-sub
 316                     else:
 317                         expo = total-sub
 318
 319                     if isinstance(term, Basic.exp):
 320                         expr.append(Basic.exp(expo*term[:]))
 321                     else:
 322                         if isinstance(expo, Basic.One):
 323                             expr.append(term)
 324                         else:
 325                             expr.append(Pow(term, expo))
 326
 327                 numerator.append(coeff*Mul(*([numer] + expr)))
 328
 329             return Add(*numerator)/(product*Mul(*denominator))
 330         elif isinstance(expr, (Mul, Pow)):
 331             return type(expr)(*[ _together(t) for t in expr ])
 332         elif isinstance(expr, Function) and deep:
 333             return expr.func(*[ _together(t) for t in expr ])
 334         else:
 335             return expr
 336
 337     return _together(separate(expr))
 338
 339 #apart -> partial fractions decomposition (will be here :)
 340
 341 def collect(expr, syms, evaluate=True, exact=False):
 342     """Collect additive terms with respect to a list of symbols up
 343        to powers with rational exponents. By the term symbol here
 344        are meant arbitrary expressions, which can contain powers,
 345        products, sums etc. In other words symbol is a pattern
 346        which will be searched for in the expression's terms.
 347
 348        This function will not apply any redundant expanding to the
 349        input expression, so user is assumed to enter expression in
 350        final form. This makes 'collect' more predictable as there
 351        is no magic behind the scenes. However it is important to
 352        note, that powers of products are converted to products of
 353        powers using 'separate' function.
 354
 355        There are two possible types of output. First, if 'evaluate'
 356        flag is set, this function will return a single expression
 357        or else it will return a dictionary with separated symbols
 358        up to rational powers as keys and collected sub-expressions
 359        as values respectively.
 360
 361        >>> from sympy import *
 362        >>> x, y, z = symbols('x', 'y', 'z')
 363        >>> a, b, c = symbols('a', 'b', 'c')
 364
 365        This function can collect symbolic coefficients in polynomial
 366        or rational expressions. It will manage to find all integer or
 367        rational powers of collection variable:
 368
 369        >>> collect(a*x**2 + b*x**2 + a*x - b*x + c, x)
 370        c + x*(a - b) + x**2*(a + b)
 371
 372        The same result can achieved in dictionary form:
 373
 374        >>> collect(a*x**2 + b*x**2 + a*x - b*x + c, x, evaluate=False)
 375        {1: c, x**2: a + b, x: a - b}
 376
 377        You can also work with multi-variate polynomials. However
 378        remember that this function is greedy so it will care only
 379        about a single symbol at time, in specification order:
 380
 381        >>> collect(x**2 + y*x**2 + x*y + y + a*y, [x, y])
 382        x*y + y*(1 + a) + x**2*(1 + y)
 383
 384        >>> collect(x**2*y**4 + z*(x*y**2)**2 + z + a*z, [x*y**2, z])
 385        z*(1 + a) + x**2*y**4*(1 + z)
 386
 387        Also more complicated expressions can be used as patterns:
 388
 389        >>> collect(a*sin(2*x) + b*sin(2*x), sin(2*x))
 390        (a + b)*sin(2*x)
 391
 392        >>> collect(a*x**2*log(x)**2 + b*(x*log(x))**2, x*log(x))
 393        x**2*log(x)**2*(a + b)
 394
 395        It is also possible to work with symbolic powers, although
 396        it has more complicated behaviour, because in this case
 397        power's base and symbolic part of the exponent are treated
 398        as a single symbol:
 399
 400        #>>> collect(a*x**c + b*x**c, x)
 401        #a*x**c + b*x**c
 402
 403        #>>> collect(a*x**c + b*x**c, x**c)
 404        #x**c*(a + b)
 405
 406        However if you incorporate rationals to the exponents, then
 407        you will get well known behaviour:
 408
 409        #>>> collect(a*x**(2*c) + b*x**(2*c), x**c)
 410        #x**(2*c)*(a + b)
 411
 412        Note also that all previously stated facts about 'collect'
 413        function apply to the exponential function, so you can get:
 414
 415        #>>> collect(a*exp(2*x) + b*exp(2*x), exp(x))
 416        #(a+b)*exp(2*x)
 417
 418        If you are interested only in collecting specific powers
 419        of some symbols then set 'exact' flag in arguments:
 420
 421        >>> collect(a*x**7 + b*x**7, x, exact=True)
 422        a*x**7 + b*x**7
 423
 424        >>> collect(a*x**7 + b*x**7, x**7, exact=True)
 425        x**7*(a + b)
 426
 427        You can also apply this function to differential equations, where
 428        derivatives of arbitary order can be collected:
 429
 430        #>>> from sympy import Derivative as D
 431        #>>> f = Function(x)
 432
 433        #>>> collect(a*D(f,x) + b*D(f,x), D(f,x))
 434        #(a+b)*Function'(x)
 435
 436        #>>> collect(a*D(D(f,x),x) + b*D(D(f,x),x), D(f,x))
 437        #(a+b)*(Function'(x))'
 438
 439        #>>> collect(a*D(D(f,x),x) + b*D(D(f,x),x), D(f,x), exact=True)
 440        #a*(Function'(x))'+b*(Function'(x))'
 441
 442        #>>> collect(a*D(D(f,x),x) + b*D(D(f,x),x) + a*D(f,x) + b*D(f,x), D(f,x))
 443        #(a+b)*Function'(x)+(a+b)*(Function'(x))'
 444
 445        Or you can even match both derivative order and exponent at time:
 446
 447        #>>> collect(a*D(D(f,x),x)**2 + b*D(D(f,x),x)**2, D(f,x))
 448        #(a+b)*(Function'(x))'**2
 449
 450     """
 451     def make_expression(terms):
 452         product = []
 453
 454         for term, rat, sym, deriv in terms:
 455             if deriv is not None:
 456                 var, order = deriv
 457
 458                 while order > 0:
 459                     term, order = Derivative(term, var), order-1
 460
 461             if sym is None:
 462                 if isinstance(rat, Basic.One):
 463                     product.append(term)
 464                 else:
 465                     product.append(Pow(term, rat))
 466             else:
 467                 product.append(Pow(term, rat*sym))
 468
 469         return Mul(*product)
 470
 471     def parse_derivative(deriv):
 472         # scan derivatives tower in the input expression and return
 473         # underlying function and maximal differentiation order
 474         expr, sym, order = deriv.f, deriv.x, 1
 475
 476         while isinstance(expr, Derivative) and expr.x == sym:
 477             expr, order = expr.f, order+1
 478
 479         return expr, (sym, Rational(order))
 480
 481     def parse_term(expr):
 482         rat_expo, sym_expo = Rational(1), None
 483         sexpr, deriv = expr, None
 484
 485         if isinstance(expr, Pow):
 486             if isinstance(expr.base, Derivative):
 487                 sexpr, deriv = parse_derivative(expr.base)
 488             else:
 489                 sexpr = expr.base
 490
 491             if isinstance(expr.exp, Rational):
 492                 rat_expo = expr.exp
 493             elif isinstance(expr.exp, Mul):
 494                 coeff, tail = term.exp.as_coeff_terms()
 495
 496                 if isinstance(coeff, Rational):
 497                     rat_expo, sym_expo = coeff, Basic.Mul(*tail)
 498                 else:
 499                     sym_expo = expr.exp
 500             else:
 501                 sym_expo = expr.exp
 502         elif isinstance(expr, Basic.exp):
 503             if isinstance(expr[0], Rational):
 504                 sexpr, rat_expo = Basic.exp(Rational(1)), expr[0]
 505             elif isinstance(expr[0], Mul):
 506                 coeff, tail = expr[0].as_coeff_terms()
 507
 508                 if isinstance(coeff, Rational):
 509                     sexpr, rat_expo = Basic.exp(Basic.Mul(*tail)), coeff
 510         elif isinstance(expr, Derivative):
 511             sexpr, deriv = parse_derivative(expr)
 512
 513         return sexpr, rat_expo, sym_expo, deriv
 514
 515     def parse_expression(terms, pattern):
 516         pattern = make_list(pattern, Mul)
 517
 518         if len(terms) < len(pattern):
 519             # pattern is longer than  matched product
 520             # so no chance for positive parsing result
 521             return None
 522         else:
 523             pattern = [ parse_term(elem) for elem in pattern ]
 524
 525             elems, common_expo, has_deriv = [], Rational(1), False
 526
 527             for elem, e_rat, e_sym, e_ord in pattern:
 528                 if e_ord is not None:
 529                     # there is derivative in the pattern so
 530                     # there will by small performance penalty
 531                     has_deriv = True
 532
 533                 for j in range(len(terms)):
 534                     term, t_rat, t_sym, t_ord = terms[j]
 535
 536                     if elem == term and e_sym == t_sym:
 537                         if exact == False:
 538                             # we don't have to exact so find common exponent
 539                             # for both expression's term and pattern's element
 540                             expo = t_rat / e_rat
 541
 542                             if isinstance(common_expo, Basic.One):
 543                                 common_expo = expo
 544                             else:
 545                                 # common exponent was negotiated before so
 546                                 # teher is no chance for pattern match unless
 547                                 # common and current exponents are equal
 548                                 if common_expo != expo:
 549                                     return None
 550                         else:
 551                             # we ought to be exact so all fields of
 552                             # interest must match in very details
 553                             if e_rat != t_rat or e_ord != t_ord:
 554                                 continue
 555
 556                         # found common term so remove it from the expression
 557                         # and try to match next element in the pattern
 558                         elems.append(terms[j])
 559                         del terms[j]
 560
 561                         break
 562                 else:
 563                     # pattern element not found
 564                     return None
 565
 566             return terms, elems, common_expo, has_deriv
 567
 568     if evaluate:
 569         if isinstance(expr, Basic.Mul):
 570             ret = 1
 571             for term in expr:
 572                 ret *= collect(term, syms, True, exact)
 573             return ret
 574         elif isinstance(expr, Basic.Pow):
 575             b = collect(expr.base, syms, True, exact)
 576             return Basic.Pow(b, expr.exp)
 577
 578     summa = [ separate(i) for i in make_list(Basic.sympify(expr), Add) ]
 579
 580     if isinstance(syms, list):
 581         syms = [ separate(s) for s in syms ]
 582     else:
 583         syms = [ separate(syms) ]
 584
 585     collected, disliked = {}, Rational(0)
 586
 587     for product in summa:
 588         terms = [ parse_term(i) for i in make_list(product, Mul) ]
 589
 590         for symbol in syms:
 591             result = parse_expression(terms, symbol)
 592
 593             if result is not None:
 594                 terms, elems, common_expo, has_deriv = result
 595
 596                 # when there was derivative in current pattern we
 597                 # will need to rebuild its expression from scratch
 598                 if not has_deriv:
 599                     index = Pow(symbol, common_expo)
 600                 else:
 601                     index = make_expression(elems)
 602
 603                 terms = separate(make_expression(terms))
 604                 index = separate(index)
 605
 606                 if index in collected:
 607                     collected[index] += terms
 608                 else:
 609                     collected[index] = terms
 610
 611                 break
 612         else:
 613             # none of the patterns matched
 614             disliked += product
 615
 616     if disliked != Rational(0):
 617         collected[Rational(1)] = disliked
 618
 619     if evaluate:
 620         return Add(*[ a*b for a, b in collected.iteritems() ])
 621     else:
 622         return collected
 623
 624 def ratsimp(expr):
 625     """
 626     Usage
 627     =====
 628         ratsimp(expr) -> joins two rational expressions and returns the simples form
 629
 630     Notes
 631     =====
 632         Currently can simplify only simple expressions, for this to be really usefull
 633         multivariate polynomial algorithms are needed
 634
 635     Examples
 636     ========
 637         >>> from sympy import *
 638         >>> x = Symbol('x')
 639         >>> y = Symbol('y')
 640         >>> e = ratsimp(1/x + 1/y)
 641     """
 642
 643     if isinstance(expr, Pow):
 644         return Pow(ratsimp(expr.base), ratsimp(expr.exp))
 645     elif isinstance(expr, Mul):
 646         res = []
 647         for x in expr:
 648             res.append( ratsimp(x) )
 649         return Mul(*res)
 650     elif isinstance(expr, Basic.Function):
 651         return expr.func(*[ ratsimp(t) for t in expr ])
 652
 653     #elif isinstance(expr, Function):
 654     #    return type(expr)( ratsimp(expr[0]) )
 655     elif not isinstance(expr, Add):
 656         return expr
 657
 658     def get_num_denum(x):
 659         """Matches x = a/b and returns a/b."""
 660         a,b = map(Wild, 'ab')
 661         r = x.match(a/b)
 662         if r is not None and len(r) == 2:
 663             return r[a],r[b]
 664         return x, 1
 665
 666     x = expr[0]
 667     y = Add(*expr[1:])
 668
 669     a,b = get_num_denum(ratsimp(x))
 670     c,d = get_num_denum(ratsimp(y))
 671
 672     num = a*d+b*c
 673     denum = b*d
 674
 675     # Check to see if the numerator actually expands to 0
 676     if num.expand() == 0:
 677         return 0
 678
 679     #we need to cancel common factors from numerator and denumerator
 680     #but SymPy doesn't yet have a multivariate polynomial factorisation
 681     #so until we have it, we are just returning the correct results here
 682     #to pass all tests...
 683     if isinstance(denum,Pow):
 684         e = (num/denum[0]).expand()
 685         f = (e/(-2*Symbol("y"))).expand()
 686         if f == denum/denum[0]:
 687             return -2*Symbol("y")
 688         return e/(denum/denum[0])
 689     return num/denum
 690
 691 def trigsimp(expr, deep=False):
 692     """
 693     Usage
 694     =====
 695         trig(expr) -> reduces expression by using known trig identities
 696
 697     Notes
 698     =====
 699
 700
 701     Examples
 702     ========
 703         >>> from sympy import *
 704         >>> x = Symbol('x')
 705         >>> y = Symbol('y')
 706         >>> e = 2*sin(x)**2 + 2*cos(x)**2
 707         >>> trigsimp(e)
 708         2
 709         >>> trigsimp(log(e))
 710         log(2*cos(x)**2 + 2*sin(x)**2)
 711         >>> trigsimp(log(e), deep=True)
 712         log(2)
 713     """
 714     from sympy.core.basic import S
 715     sin, cos, tan, cot = Basic.sin, Basic.cos, Basic.tan, Basic.cot
 716
 717     #XXX this stopped working:
 718     if expr == 1/cos(Symbol("x"))**2 - 1:
 719         return tan(Symbol("x"))**2
 720
 721     if isinstance(expr, Basic.Function):
 722         if deep:
 723             return expr.func( trigsimp(expr[0], deep) )
 724     elif isinstance(expr, Mul):
 725         ret = Rational(1)
 726         for x in expr:
 727             ret *= trigsimp(x, deep)
 728         return ret
 729     elif isinstance(expr, Pow):
 730         return Pow(trigsimp(expr.base, deep), trigsimp(expr.exp, deep))
 731     elif isinstance(expr, Add):
 732         # TODO this needs to be faster
 733
 734         # The types of trig functions we are looking for
 735         a,b,c = map(Wild, 'abc')
 736         matchers = (
 737             (a*sin(b)**2, a - a*cos(b)**2),
 738             (a*tan(b)**2, a*(1/cos(b))**2 - a),
 739             (a*cot(b)**2, a*(1/sin(b))**2 - a)
 740         )
 741
 742         # Scan for the terms we need
 743         ret = Integer(0)
 744         for term in expr:
 745             term = trigsimp(term, deep)
 746             res = None
 747             for pattern, result in matchers:
 748                 res = term.match(pattern)
 749                 if res is not None:
 750                     ret += result.subs_dict(res)
 751                     break
 752             if res is None:
 753                 ret += term
 754
 755         # Reduce any lingering artifacts, such as sin(x)**2 changing
 756         # to 1-cos(x)**2 when sin(x)**2 was "simpler"
 757         artifacts = (
 758             (a - a*cos(b)**2 + c, a*sin(b)**2 + c, cos),
 759             (a - a*(1/cos(b))**2 + c, -a*tan(b)**2 + c, cos),
 760             (a - a*(1/sin(b))**2 + c, -a*cot(b)**2 + c, sin)
 761         )
 762
 763         expr = ret
 764         for pattern, result, ex in artifacts:
 765             # Substitute a new wild that excludes some function(s)
 766             # to help influence a better match. This is because
 767             # sometimes, for example, 'a' would match sec(x)**2
 768             a_t = Wild('a', exclude=[ex])
 769             pattern = pattern.subs(a, a_t)
 770             result = result.subs(a, a_t)
 771
 772             m = expr.match(pattern)
 773             while m is not None:
 774                 if m[a_t] == 0 or -m[a_t] in m[c][:] or m[a_t] + m[c] == 0:
 775                     break
 776                 expr = result.subs_dict(m)
 777                 m = expr.match(pattern)
 778
 779         return expr
 780     return expr
 781
 782 def radsimp(expr):
 783     """
 784     Rationalize the denominator.
 785
 786     Examples:
 787     =========
 788         >>> from sympy import *
 789         >>> radsimp(1/(2+sqrt(2)))
 790         1 - 1/2*2**(1/2)
 791         >>> x,y = map(Symbol, 'xy')
 792         >>> e = ( (2+2*sqrt(2))*x+(2+sqrt(8))*y )/( 2+sqrt(2) )
 793         >>> radsimp(e)
 794         x*2**(1/2) + y*2**(1/2)
 795     """
 796     n,d = fraction(expr)
 797     a,b,c = map(Wild, 'abc')
 798     r = d.match(a+b*Basic.sqrt(c))
 799     if r is not None:
 800         a = r[a]
 801         if r[b] == 0:
 802             b,c = 0,0
 803         else:
 804             b,c = r[b],r[c]
 805
 806         syms = list(n.atoms(type=Basic.Symbol))
 807         n = collect( (n*(a-b*Basic.sqrt(c))).expand(), syms )
 808         d = a**2 - c*b**2
 809
 810     return n/d
 811
 812 def powsimp(expr, deep=False):
 813     """
 814     Usage
 815     =====
 816         powsimp(expr, deep) -> reduces expression by combining powers with
 817             similar bases and exponents.
 818
 819     Notes
 820     =====
 821         If deep is True then powsimp() will also simplify arguments of
 822         functions. By default deep is set to False.
 823
 824
 825     Examples
 826     ========
 827         >>> from sympy import *
 828         >>> x,n = map(Symbol, 'xn')
 829         >>> e = x**n * (x*n)**(-n) * n
 830         >>> powsimp(e)
 831         n**(1 - n)
 832
 833         >>> powsimp(log(e))
 834         log(n*x**n*(n*x)**(-n))
 835
 836         >>> powsimp(log(e), deep=True)
 837         log(n**(1 - n))
 838     """
 839     def _powsimp(expr):
 840         if isinstance(expr, Basic.Pow):
 841             if deep:
 842                 return Basic.Pow(powsimp(expr.base), powsimp(expr.exp))
 843             return expr
 844         #elif isinstance(expr, Basic.Apply) and deep:
 845         #    return expr.func(*[powsimp(t) for t in expr])
 846         elif isinstance(expr, Basic.Add):
 847             return Basic.Add(*[powsimp(t) for t in expr])
 848         elif isinstance(expr, Basic.Mul):
 849             # Collect base/exp data, while maintaining order in the
 850             # non-commutative parts of the product
 851             c_powers = {}
 852             nc_part = []
 853             for term in expr:
 854                 if term.is_commutative:
 855                     b,e = term.as_base_exp()
 856                     c_powers[b] = c_powers.get(b, 0) + e
 857                 else:
 858                     nc_part.append(term)
 859
 860             # Pull out numerical coefficients from exponent
 861             for b,e in c_powers.items():
 862                 exp_c, exp_t = e.as_coeff_terms()
 863                 if not isinstance(exp_c, Basic.One) and exp_t:
 864                     del c_powers[b]
 865                     new_base = Basic.Pow(b, exp_c)
 866                     if new_base in c_powers:
 867                         c_powers[new_base] += Basic.Mul(*exp_t)
 868                     else:
 869                         c_powers[new_base] = Basic.Mul(*exp_t)
 870
 871             # Combine bases whenever they have the same exponent which is
 872             # not numeric
 873             c_exp = {}
 874             for b, e in c_powers.items():
 875                 if e in c_exp:
 876                     c_exp[e].append(b)
 877                 else:
 878                     c_exp[e] = [b]
 879
 880             # Merge back in the results of the above to form a new product
 881             for e in c_exp:
 882                 bases = c_exp[e]
 883                 if len(bases) > 1:
 884                     for b in bases:
 885                         del c_powers[b]
 886                     new_base = Mul(*bases)
 887                     if new_base in c_powers:
 888                         c_powers[new_base] += e
 889                     else:
 890                         c_powers[new_base] = e
 891
 892             c_part = [ Basic.Pow(b,e) for b,e in c_powers.items() ]
 893             return Basic.Mul(*(c_part + nc_part))
 894         return expr
 895
 896     return _powsimp(separate(expr, deep=deep))
 897
 898 def normal(expr, *syms):
 899     p, q = together(expr).as_numer_denom()
 900
 901     if p.is_polynomial(*syms) and q.is_polynomial(*syms):
 902         from sympy.polynomials import gcd, quo
 903
 904         G = gcd(p, q, syms)
 905
 906         if not isinstance(G, Basic.One):
 907             p = quo(p, G, syms)
 908             q = quo(q, G, syms)
 909
 910     return p / q
 911
 912 def hypersimp(term, n, consecutive=True, simplify=True):
 913     """Given combinatorial term a(n) simplify its consecutive term
 914        ratio ie. a(n+1)/a(n). The term can be composed of functions
 915        and integer sequences which have equivalent represenation
 916        in terms of gamma special function. Currently ths includes
 917        factorials (falling, rising), binomials and gamma it self.
 918
 919        The algorithm performs three basic steps:
 920
 921            (1) Rewrite all functions in terms of gamma, if possible.
 922
 923            (2) Rewrite all occurences of gamma in terms of produtcs
 924                of gamma and rising factorial with integer, absolute
 925                constant exponent.
 926
 927            (3) Perform simplification of nested fractions, powers
 928                and if the resulting expression is a quotient of
 929                polynomials, reduce their total degree.
 930
 931        If the term given is hypergeometric then the result of this
 932        procudure is a quotient of polynomials of minimal degree.
 933        Sequence is hypergeometric if it is anihilated by linear,
 934        homogeneous recurrence operator of first order, so in
 935        other words when a(n+1)/a(n) is a rational function.
 936
 937        When the status of being hypergeometric or not, is required
 938        then you can avoid additional simplification by unsetting
 939        'simplify' flag.
 940
 941        This algorithm, due to Wolfram Koepf, is very simple but
 942        powerful, however its full potential will be visible when
 943        simplification in general will improve.
 944
 945        For more information on the implemented algorithm refer to:
 946
 947        [1] W. Koepf, Algorithms for m-fold Hypergeometric Summation,
 948            Journal of Symbolic Computation (1995) 20, 399-417
 949     """
 950     term = Basic.sympify(term)
 951
 952     if consecutive == True:
 953         term = term.subs(n, n+1)/term
 954
 955     expr = term.rewrite(gamma).expand(func=True, basic=False)
 956
 957     p, q = together(expr).as_numer_denom()
 958
 959     if p.is_polynomial(n) and q.is_polynomial(n):
 960         if simplify == True:
 961             from sympy.polynomials import gcd, quo
 962
 963             G = gcd(p, q, n)
 964
 965             if not isinstance(G, Basic.One):
 966                 p = quo(p, G, n)
 967                 q = quo(q, G, n)
 968
 969                 p = p.as_polynomial(n)
 970                 q = q.as_polynomial(n)
 971
 972                 a, p = p.as_integer()
 973                 b, q = q.as_integer()
 974
 975                 p = p.as_basic()
 976                 q = q.as_basic()
 977
 978                 return (b/a) * (p/q)
 979
 980         return p/q
 981     else:
 982         return None
 983
 984 def hypersimilar(f, g, n):
 985     return hypersimp(f/g, n, consecutive=False, simplify=False) is not None
 986
 987 def combsimp(expr):
 988     return expr
 989
 990 def simplify(expr):
 991     #from sympy.specfun.factorials import factorial, factorial_simplify
 992     #if expr.has_class(factorial):
 993     #    expr = factorial_simplify(expr)
 994     a,b = [ t.expand() for t in fraction(powsimp(expr)) ]
 995     ret = together(radsimp(ratsimp(a/b)))
 996     n,d = fraction(ret)
 997     if isinstance(d, (Basic.One, Basic.NegativeOne)):
 998         return d*n
 999     n_var = n.atoms(type=Symbol)
1000     d_var = d.atoms(type=Symbol)
1001     if n_var and d_var and n.is_polynomial(*n_var) and d.is_polynomial(*d_var):
1002         from sympy.polynomials import div, factor
1003         q,r = div(n, d)
1004         if r == 0:
1005             return q
1006         else:
1007             return q + factor(r) / factor(d)
1008     return ret