sympy/simplify/simplify.py

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