sympy/polys/polynomial.py

   1
   2 from sympy.core.add import Add
   3 from sympy.core.mul import Mul
   4 from sympy.core.power import Pow
   5 from sympy.core.sympify import sympify
   6 from sympy.core.basic import Basic, S, C
   7 from sympy.core.symbol import Symbol, Wild
   8 from sympy.core.numbers import Integer, igcd, ilcm
   9
  10 from sympy.utilities import all, any
  11
  12 from sympy.polys.monomial import monomial_cmp, monomial_mul, \
  13     monomial_div, monomial_max, monomial_min, monomial_as_basic
  14
  15 import sympy.polys
  16
  17 import math
  18
  19 class SymbolsError(Exception):
  20     pass
  21
  22 class PolynomialError(Exception):
  23     pass
  24
  25 class CoefficientError(PolynomialError):
  26     pass
  27
  28 class UnivariatePolyError(PolynomialError):
  29
  30     def __init__(self, poly):
  31         self.message = "%s is an univariate polynomial" % poly
  32
  33     def __str__(self):
  34         return self.message
  35
  36 class MultivariatePolyError(PolynomialError):
  37
  38     def __init__(self, poly):
  39         self.message = "%s is a multivariate polynomial" % poly
  40
  41     def __str__(self):
  42         return self.message
  43
  44 ##
  45 ## TODO:
  46 ##
  47 ##  [!] Copy + Write tests for everything.
  48 ##  [2] Analyze dense MV multiplication.
  49 ##  [~] Create MV poly with int coeffs.
  50 ##  [4] Think on GF polys (for factoring).
  51 ##  [5] Improve coefficients analysis.
  52 ##  [6] Implement FactoredPoly.
  53 ##  [7] Concept of Monomial.
  54 ##
  55
  56 class Poly(Basic):
  57     """Represents polynomials with symbolic coefficients.
  58
  59        Polynomials are internally represented as two lists containing
  60        coefficients and monomials (tuples of exponents) respectively.
  61        Stored are only terms with non-zero coefficients, so generally
  62        all polynomials are considered sparse. However algorithms will
  63        detect dense polynomials and use the appropriate method to
  64        solve the given problem in the most efficient way.
  65
  66        The most common way to initialize a polynomial instance is to
  67        provide a valid expression together witch a set of ordered
  68        symbols and, additionally, monomial ordering:
  69
  70             Poly(expression, x_1, x_2, ..., x_n, order='grlex')
  71
  72        If the given expression is not a polynomial with respect to
  73        the given variables, an exception is raised. Alternatively
  74        you can use Basic.as_poly to avoid exception handling.
  75
  76        By default ordering of monomials can be omitted. In this case
  77        graded lexicographic order will be used. Anyway remember that
  78        'order' is a keyword argument.
  79
  80        Currently there are supported four standard orderings:
  81
  82            [1] lex       -> lexicographic order
  83            [2] grlex     -> graded lex order
  84            [3] grevlex   -> reversed grlex order
  85            [4] 1-el      -> first elimination order
  86
  87        Polynomial can be also constructed explicitly by specifying
  88        a collection of coefficients and monomials. This can be done
  89        in at least three different ways:
  90
  91            [1] [(c_1, M_1), (c_2, M_2), ..., (c_1, M_1)]
  92
  93            [2] (c_1, c2, ..., c_n), (M_1, M_2, ..., M_n)
  94
  95            [3] { M_1 : c_1, M_2 : c_2, ..., M_n : c_n }
  96
  97        Although all three representation look similar, they are
  98        designed for different tasks and have specific properties:
  99
 100            [1] All coefficients and monomials  are validated before
 101                polynomial instance is created. Monomials are sorted
 102                with respect to the given order.
 103
 104            [2] No validity checking, no sorting, just a raw copy.
 105
 106            [3] Also no validity checking however monomials are
 107                sorted with respect to the given order.
 108
 109        For interactive usage choose [1] as it's safe and relatively
 110        fast. Use [2] or [3] internally for time critical algorithms,
 111        when you know that coefficients and monomials will be valid.
 112
 113        Implemented methods:
 114
 115            U --> handles uniformly int | long and Basic coefficients
 116                  (other methods need to be overridden in IntegerPoly)
 117
 118            P --> a property (with appropriate decorator)
 119
 120            - --> otherwise
 121
 122        [1] Representation conversion:
 123
 124           [1.1] [--] as_basic   --> converts to a valid sympy expression
 125           [1.2] [U-] as_dict    --> returns { M_1 : c_1, ..., M_1 : c_1 }
 126           [1.3] [U-] as_uv_dict --> returns dict with integer keys
 127
 128        [2] Polynomial internals:
 129
 130           [2.1] [UP] coeffs  --> list of coefficients
 131           [2.2] [UP] monoms  --> list of monomials
 132           [2.3] [UP] symbols --> an ordered list of symbols
 133           [2.4] [UP] order   --> the ordering of monomials
 134           [2.5] [UP] stamp   --> frozen set of polynomial's variables
 135           [2.6] [UP] domain  --> coefficient domain, see [8] for details
 136           [2.7] [UP] args    --> required for printing and hashing
 137           [2.8] [UP] flags   --> dict of keyword arguments to Poly
 138
 139        [3] General characteristics:
 140
 141           [3.1] [UP] norm         --> computes oo-norm of coefficients
 142           [3.2] [UP] length       --> counts the total number of terms
 143           [3.3] [UP] degree       --> returns degree of the leading monomial
 144           [3.4] [UP] total_degree --> returns true total degree of a polynomial
 145
 146        [4] Items generators:
 147
 148           [4.1] [U-] iter_coeffs     --> iterate over coefficients
 149           [4.2] [U-] iter_monoms     --> iterate over monomials
 150           [4.3] [U-] iter_terms      --> iterate over terms
 151
 152           [4.4] [U-] iter_all_coeffs --> iterate over all coefficients
 153           [4.5] [U-] iter_all_monoms --> iterate over all monomials
 154           [4.6] [U-] iter_all_terms  --> iterate over all terms
 155
 156        [5] Leading and trailing items:
 157
 158           [5.1] [UP] lead_coeff or LC --> returns the leading coefficient
 159           [5.2] [UP] lead_monom or LM --> returns the leading monomial
 160           [5.3] [UP] lead_term  or LT --> returns the leading term
 161
 162           [5.4] [UP] tail_coeff or TC --> returns the trailing coefficient
 163           [5.5] [UP] tail_monom or TM --> returns the trailing monomial
 164           [5.6] [UP] tail_term  or TT --> returns the trailing term
 165
 166        [6] Arithmetic operators:
 167
 168           [6.1]  [U-] __abs__     --> for all i, abs(c_i)
 169           [6.2]  [U-] __neg__     --> polynomial negation
 170           [6.3]  [U-] __add__     --> polynomial addition
 171           [6.4]  [U-] __sub__     --> polynomial subtraction
 172           [6.5]  [--] __mul__     --> polynomial multiplication
 173           [6.6]  [U-] __pow__     --> polynomial exponentiation
 174           [6.7]  [U-] __div__     --> polynomial quotient
 175           [6.8]  [U-] __mod__     --> polynomial remainder
 176           [6.9]  [U-] __divmod__  --> both 'div' and 'mod'
 177           [6.10] [U-] __truediv__ --> the same as 'div'
 178
 179        [7] Polynomial properties:
 180
 181           [7.1]  [UP] is_zero          --> only one term c_0 == 0
 182           [7.2]  [UP] is_one           --> only one term c_0 == 1
 183           [7.3]  [-P] is_number        --> only numeric constant term
 184           [7.4]  [UP] is_constant      --> only arbitrary constant term
 185           [7.5]  [UP] is_monomial      --> number of terms == 1
 186           [7.6]  [UP] is_univariate    --> number of variables == 1
 187           [7.7]  [UP] is_multivariate  --> number of variables > 1
 188           [7.8]  [UP] is_homogeneous   --> has constant term
 189           [7.9]  [UP] is_inhomogeneous --> no constant term
 190           [7.10] [UP] is_sparse        --> filled with less than 90% of terms
 191           [7.11] [UP] is_dense         --> filled with more than 90% of terms
 192           [7.12] [UP] is_monic         --> returns True if leading coeff is one
 193           [7.13] [UP] is_primitive     --> returns True if content is one
 194           [7.14] [UP] is_squarefree    --> no factors of multiplicity >= 2
 195           [7.15] [UP] is_linear        --> all monomials of degree <= 1
 196
 197        [8] Coefficients properties:
 198
 199           [8.1] [UP] has_integer_coeffs  --> for all i, c_i in Z
 200           [8.2] [UP] has_rational_coeffs --> for all i, c_i in Q
 201           [8.3] [UP] has_radical_coeffs  --> for all i, c_i in R
 202           [8.4] [UP] has_complex_coeffs  --> for all i, c_i in C
 203           [8.5] [UP] has_symbolic_coeffs --> otherwise
 204
 205        [9] Special functionality:
 206
 207           [9.1]  [--] as_monic      --> divides all coefficients by LC
 208           [9.2]  [--] as_integer    --> returns polynomial with integer coeffs
 209
 210           [9.3]  [-P] content       --> returns GCD of polynomial coefficients
 211           [9.4]  [--] as_primitive  --> returns content and primitive part
 212
 213           [9.5]  [U-] as_squarefree --> returns square-free part of a polynomial
 214
 215           [9.6]  [U-] as_reduced    --> extract GCD of polynomial monomials
 216
 217           [9.7]  [--] map_coeffs    --> applies a function to all coefficients
 218           [9.8]  [U-] coeff         --> returns coefficient of the given monomial
 219
 220           [9.9]  [U-] unify_with    --> returns polys with a common set of symbols
 221
 222           [9.10] [U-] diff          --> efficient polynomial differentiation
 223           [9.11] [U-] integrate     --> efficient polynomial integration
 224
 225        [10] Operations on terms:
 226
 227           [10.1] [U-] add_term --> add a term to a polynomial
 228           [10.2] [U-] sub_term --> subtract a term from a polynomial
 229
 230           [10.3] [--] mul_term --> multiply a polynomial with a term
 231           [10.4] [--] div_term --> divide a polynomial by a term
 232
 233           [10.5] [U-] remove_lead_term --> remove leading term (up to order)
 234           [10.6] [U-] remove_last_term --> remove last term (up to order)
 235
 236        [11] Substitution and evaluation:
 237
 238           [11.1] [U-] __call__    --> evaluates poly a the given point
 239           [11.2] [U-] evaluate    --> evaluates poly for specific vars
 240           [11.3] [--] _eval_subs  --> efficiently substitute variables
 241
 242        [12] Comparison functions:
 243
 244           [12.1] [U-] __eq__      --> P1 == P, up to order of monomials
 245           [12.2] [U-] __ne__      --> P1 != P, up to order of monomials
 246           [12.3] [U-] __nonzero__ --> check if zero polynomial
 247
 248        [13] String representation functions:
 249
 250           [13.1] [U-] repr      --> returns represantation of 'self'
 251           [13.3] [U-] str       --> returns pretty representation
 252
 253        [14] Other (derived from Basic) functionality:
 254
 255           [14.1] [U-] _eval_is_polynomial --> checks if poly is a poly
 256
 257        For general information on polynomials and algorithms, refer to:
 258
 259        [1] D. E. Knuth, The Art of Computer Programming: Seminumerical
 260            Algorithms, v.2, Addison-Wesley Professional, 1997
 261
 262        [2] J. von zur Gathen, J. Gerhard, Modern Computer Algebra,
 263            Second Edition, Cambridge University Press, 2003
 264
 265        [3] K. Geddes, S. R. Czapor, G. Labahn, Algorithms for
 266            Computer Algebra, Springer, 1992
 267
 268        [4] D. Bini, V. Y.Pan, Polynomial and Matrix Computations:
 269            Volume 1: Fundamental Algorithms, Birkhauser Boston, 1994
 270
 271        [5] R. T. Moenck, Practical Fast Polynomial Multiplication,
 272            Proceedings of the third ACM symposium on Symbolic and
 273            algebraic computation, ACM Press, Yorktown Heights,
 274            New York, USA, 1976, pp. 136-148
 275
 276        [6] M. Bronstein, Symbolic Integration I: Transcendental
 277            Functions, Second Edition, Springer-Verlang, 2005
 278
 279        See also documentation of particular methods.
 280
 281     """
 282
 283     def __new__(cls, poly, *symbols, **flags):
 284         N, terms = len(symbols), {}
 285         coeffs, monoms = (), ()
 286
 287         order = flags.pop('order', 'grlex')
 288
 289         if N == 0:
 290             if isinstance(poly, Poly):
 291                 if order != poly.order:
 292                     symbols = poly.symbols
 293                     N = len(poly.symbols)
 294                     poly = poly.as_dict()
 295                 else:
 296                     return poly
 297             else:
 298                 raise SymbolsError("No symbols were given")
 299
 300         stamp = frozenset(symbols)
 301
 302         if len(stamp) != N:
 303             raise SymbolsError("Duplicate symbols: %s" % (symbols,))
 304
 305         symbols = tuple(map(sympify, symbols))
 306
 307         if any(not s.is_Symbol for s in symbols):
 308             raise SymbolsError("Invalid symbols: %s" % (symbols,))
 309
 310         if any(not s.is_commutative for s in symbols):
 311             raise SymbolsError("Non-commutative symbols: %s" % (symbols,))
 312
 313         # { M1: c1, M2: c2, ... }
 314         if type(poly) is dict:
 315             if not poly:
 316                 coeffs = (S.Zero,)
 317                 monoms = ((0,) * N,)
 318             else:
 319                 terms = poly
 320
 321         # ((c1, c2, ...), (M1, M2, ...))
 322         elif type(poly) is tuple:
 323             if not poly:
 324                 coeffs = (S.Zero,)
 325                 monoms = ((0,) * N,)
 326             else:
 327                 if len(poly) == 2:
 328                     coeffs, monoms = poly[0], poly[1]
 329
 330                     if isinstance(coeffs, (tuple, list)):
 331                         if not (coeffs or monoms):
 332                             coeffs = (S.Zero,)
 333                             monoms = ((0,) * N,)
 334                     else:
 335                         coeffs = (coeffs,)
 336                         monoms = (monoms,)
 337                 else:
 338                     raise PolynomialError("Invalid arguments," \
 339                         " should get '(coeffs, monoms)' tuple")
 340
 341         # [ (c1,M1), (c2,M2), ... ]
 342         elif type(poly) is list:
 343             if not poly:
 344                 coeffs = (S.Zero,)
 345                 monoms = ((0,) * N,)
 346             else:
 347                 for coeff, monom in poly:
 348                     coeff = sympify(coeff)
 349
 350                     if coeff.has_any_symbols(*symbols):
 351                         raise CoefficientError("%s coefficient is dependent" \
 352                             " of polynomial's symbols %s" % (coeff, symbols))
 353
 354                     if type(monom) is int:
 355                         monom = (monom,)
 356
 357                     if len(monom) != N:
 358                         raise PolynomialError("Dimension of %s monomial does" \
 359                             " not match the number of symbols (%d)" % (monom, N))
 360
 361                     if not coeff:
 362                         continue
 363
 364                     if terms.has_key(monom):
 365                         coeff += terms[monom]
 366
 367                         if not coeff:
 368                             del terms[monom]
 369                             continue
 370
 371                     terms[monom] = coeff
 372
 373                 if not terms:
 374                     coeffs = (S.Zero,)
 375                     monoms = ((0,) * N,)
 376         elif isinstance(poly, Poly):
 377             if stamp <= poly.stamp:
 378                 if symbols == poly.symbols:
 379                     if N > 1 and order != poly.order:
 380                         terms = poly.as_dict()
 381                     else:
 382                         return poly
 383                 else:
 384                     terms = Poly._permute(poly, *symbols)
 385             else:
 386                 if any(coeff.has_any_symbols(*symbols) for coeff in poly.coeffs):
 387                     terms = Poly._decompose(poly.as_basic(), *symbols)
 388                 else:
 389                     K = len(poly.symbols)
 390
 391                     if symbols[:K] == poly.symbols:
 392                         coeffs, T = poly.coeffs, (0,) * (N - K)
 393                         monoms = tuple(M + T for M in poly.monoms)
 394
 395                         if order != poly.order:
 396                             terms = dict(zip(monoms, coeffs))
 397                     else:
 398                         terms = Poly._permute(poly, *symbols)
 399         else:
 400             terms = Poly._decompose(poly, *symbols)
 401
 402         if terms:
 403             if N == 1 and type(terms.keys()[0]) is not tuple:
 404                 keys = tuple(reversed(sorted(terms.keys())))
 405
 406                 coeffs = [ terms[i] for i in keys ]
 407                 monoms = [ (i,) for i in keys ]
 408             else:
 409                 f = monomial_cmp(order)
 410
 411                 monoms = terms.keys()
 412                 monoms.sort(f, reverse=True)
 413
 414                 coeffs = [ terms[monom] for monom in monoms ]
 415
 416         args = (tuple(coeffs), tuple(monoms),
 417                 symbols, order, stamp, None)
 418
 419         return Basic.__new__(cls, *args, **flags)
 420
 421     @staticmethod
 422     def _permute(poly, *symbols):
 423         terms, N = {}, len(symbols)
 424
 425         if symbols == poly.symbols[:N]:
 426             symbols = poly.symbols[N:]
 427
 428             for coeff, monom in poly.iter_terms():
 429                 coeff *= monomial_as_basic(monom[N:], *symbols)
 430
 431                 monom = monom[:N]
 432
 433                 if terms.has_key(monom):
 434                     coeff += terms[monom]
 435
 436                     if not coeff:
 437                         del terms[monom]
 438                         continue
 439
 440                 terms[monom] = coeff
 441         else:
 442             new_symbols, indices = list(poly.symbols), []
 443
 444             for s in new_symbols[:]:
 445                 for i, t in enumerate(symbols):
 446                     if s is t:
 447                         new_symbols.remove(s)
 448                         indices.append(i)
 449                         break
 450                 else:
 451                     indices.append(None)
 452
 453             for coeff, monom in poly.iter_terms():
 454                 exponents, j = [0] * N, 0
 455
 456                 for exp, i in zip(monom, indices):
 457                     if i is None:
 458                         coeff *= new_symbols[j]**exp
 459                         j += 1
 460                     else:
 461                         exponents[i] = exp
 462
 463                 monom = tuple(exponents)
 464
 465                 if terms.has_key(monom):
 466                     coeff += terms[monom]
 467
 468                     if not coeff:
 469                         del terms[monom]
 470                         continue
 471
 472                 terms[monom] = coeff
 473
 474         return terms
 475
 476     @staticmethod
 477     def _decompose(poly, *symbols):
 478         """Converts basic expression to dictionary representation.
 479
 480            This procedure will expand given expression and scan it,
 481            extracting coefficients and monomials according to a set
 482            of symbols, if possible. At this point there is no need
 483            for specifying ordering information. If conversion is
 484            not possible an exception is raised.
 485
 486            This method should be left for internal use. To convert
 487            expression to a polynomial use Polynomial constructor
 488            explicitly or Basic.as_poly method.
 489
 490            >>> from sympy import *
 491            >>> x, y = symbols('xy')
 492
 493            >>> Poly._decompose(y*x**2 + 3, x)
 494            {(0,): 3, (2,): y}
 495
 496            >>> Poly._decompose(y*x**2 + 3, x, y)
 497            {(0, 0): 3, (2, 1): 1}
 498
 499         """
 500         result, indices, N = {}, {}, len(symbols)
 501
 502         for i, sym in enumerate(symbols):
 503             indices[sym] = i
 504
 505         poly = sympify(poly).expand()
 506
 507         if poly.is_Add:
 508             terms = poly.args
 509         else:
 510             if poly.is_Number:
 511                 return { (0,) * N : poly }
 512             else:
 513                 terms = [ poly ]
 514
 515         for term in terms:
 516             if not term.has_any_symbols(*symbols):
 517                 coeff, monom = term, (0,) * N
 518             else:
 519                 if term.is_Mul:
 520                     factors = term.args
 521                 else:
 522                     factors = [ term ]
 523
 524                 coeff, monom = S.One, [0] * N
 525
 526                 for factor in factors:
 527                     if factor.has_any_symbols(*symbols):
 528                         if factor.is_Pow:
 529                             if factor.exp.is_Integer:
 530                                 b, e = factor.base, factor.exp.p
 531
 532                                 if b.is_Symbol and e > 0:
 533                                     monom[indices[b]] += e
 534                                     continue
 535                         elif factor.is_Symbol:
 536                             monom[indices[factor]] += 1
 537                             continue
 538
 539                         raise PolynomialError("Can't decompose %s" % factor)
 540                     else:
 541                         coeff *= factor
 542
 543                 monom = tuple(monom)
 544
 545             if result.has_key(monom):
 546                 coeff += result[monom]
 547
 548                 if not coeff:
 549                     del result[monom]
 550                     continue
 551
 552             result[monom] = coeff
 553
 554         if not result:
 555             return { (0,) * N : S.One }
 556         else:
 557             return result
 558
 559     @staticmethod
 560     def cancel(f, *symbols):
 561         """Cancel common factors in a fractional expression.
 562
 563            Given a quotient of polynomials, cancel common factors from
 564            the numerator and the denominator.  The result is formed in
 565            an expanded form, even if there was no cancellation.
 566
 567            To improve the  speed of calculations a list of symbols can
 568            be specified. This can be particularly useful in univariate
 569            case with additional parameters, to avoid Groebner basis
 570            computations.
 571
 572            The input fraction can be given as a single expression or
 573            as a tuple consisting of the numerator and denominator.
 574
 575            >>> from sympy import *
 576            >>> x,y = symbols('xy')
 577
 578            >>> Poly.cancel((x**2-y**2)/(x-y), x, y)
 579            x + y
 580
 581            >>> Poly.cancel((x**2-y**2, x-y), x, y)
 582            x + y
 583
 584         """
 585         if type(f) is tuple:
 586             numer, denom = f
 587         else:
 588             if f.is_Atom:
 589                 return f
 590
 591             if not symbols:
 592                 symbols = f.atoms(Symbol)
 593
 594                 if not symbols:
 595                     symbols = None
 596                 else:
 597                     symbols = sorted(symbols)
 598
 599             numer, denom = f.as_numer_denom()
 600
 601         numer = numer.expand()
 602
 603         if numer is S.Zero:
 604             return numer
 605
 606         denom = denom.expand()
 607
 608         if denom.is_Number:
 609             return numer / denom
 610
 611         if symbols is None:
 612             if denom.is_Add:
 613                 a, b, c = map(Wild, 'abc')
 614
 615                 r = denom.match(a + b*c**S.Half)
 616
 617                 if r is not None:
 618                     if r[b] is S.Zero:
 619                         denom = match[a]
 620                     else:
 621                         a, b, c = r[a], r[b], r[c]
 622
 623                         numer *= a-b*c**S.Half
 624                         numer = numer.expand()
 625
 626                         denom = a**2 - c*b**2
 627
 628             return numer / denom
 629
 630         try:
 631             p = Poly(numer, *symbols)
 632             q = Poly(denom, *symbols)
 633         except PolynomialError:
 634             return numer / denom
 635
 636         g = sympy.polys.algorithms.poly_gcd(p, q)
 637
 638         p = sympy.polys.algorithms.poly_div(p, g)[0]
 639         q = sympy.polys.algorithms.poly_div(q, g)[0]
 640
 641         return p.as_basic() / q.as_basic()
 642
 643     def as_basic(self):
 644         """Converts polynomial to a valid sympy expression.
 645
 646            Takes list representation of a polynomial and constructs
 647            expression using symbols used during creation of a given
 648            polynomial. If you wish to have different symbols in the
 649            resulting expressions, use subs() first.
 650
 651            Note that the result is always returned in expanded form.
 652
 653            >>> from sympy import *
 654            >>> x,y = symbols('xy')
 655
 656            >>> p = Poly(x**2+3, x)
 657            >>> p.as_basic()
 658            3 + x**2
 659
 660         """
 661         multinomial = dict(zip(self.monoms, self.coeffs))
 662         return multinomial_as_basic(multinomial, *self.symbols)
 663
 664     def as_dict(self):
 665         """Convert list representation to dictionary representation.
 666
 667            Due to hashing of immutable expressions in Python, we can
 668            not store a dictionary directly in a polynomial instance.
 669            Also it would be inconvenient in case of ordering of
 670            monomials (unnecessary sorting).
 671
 672            This method creates dictionary efficiently, so it can be
 673            used in several algorithms which perform best when dicts
 674            are being used, eg. addition, multiplication etc.
 675
 676            >>> from sympy import *
 677            >>> x,y = symbols('xy')
 678
 679            >>> p = Poly(x**2+3, x)
 680            >>> p.as_dict()
 681            {(0,): 3, (2,): 1}
 682
 683         """
 684         return dict(zip(self.monoms, self.coeffs))
 685
 686     def as_uv_dict(self):
 687         """Return dictionary representation with integer keys.
 688
 689            In case of univariate polynomials it is inefficient to
 690            to handle exponents as singleton tuples, as an example
 691            see multiplication algorithm.  This method will return
 692            dictionary representation with all tuples converted to
 693            explicit integers.
 694
 695            >>> from sympy import *
 696            >>> x,y = symbols('xy')
 697
 698            >>> p = Poly(x**2+3, x)
 699            >>> p.as_uv_dict()
 700            {0: 3, 2: 1}
 701
 702         """
 703         if self.is_univariate:
 704             return dict(zip([ M[0] for M in self.monoms ], self.coeffs))
 705         else:
 706             raise MultivariatePolyError(self)
 707
 708     @property
 709     def coeffs(self):
 710         return self._args[0]
 711
 712     @property
 713     def monoms(self):
 714         return self._args[1]
 715
 716     @property
 717     def symbols(self):
 718         return self._args[2]
 719
 720     @property
 721     def order(self):
 722         return self._args[3]
 723
 724     @property
 725     def stamp(self):
 726         return self._args[4]
 727
 728     @property
 729     def domain(self):
 730         return self._args[5]
 731
 732     @property
 733     def args(self):
 734         return self._args[0:4]
 735
 736     @property
 737     def flags(self):
 738         return { 'order' : self.order }
 739
 740     @property
 741     def length(self):
 742         return len(self.monoms)
 743
 744     @property
 745     def degree(self):
 746         """Returns degree of the leading term. """
 747         return sum(self.monoms[0])
 748
 749     @property
 750     def total_degree(self):
 751         """Returns true total degree of a polynomial. """
 752         return max([ sum(monom) for monom in self.monoms ])
 753
 754     @property
 755     def norm(self):
 756         """Computes oo-norm of polynomial's coefficients. """
 757         return max([ abs(coeff) for coeff in self.coeffs ])
 758
 759     def iter_basic_args(self):
 760         for coeff in self.coeffs:
 761             yield coeff
 762
 763         for i, symbol in enumerate(self.symbols):
 764             if any(monom[i] != 0 for monom in self.monoms):
 765                 yield symbol
 766
 767     def iter_coeffs(self):
 768         for coeff in self.coeffs:
 769             yield coeff
 770
 771     def iter_monoms(self):
 772         for monom in self.monoms:
 773             yield monom
 774
 775     def iter_terms(self):
 776         for coeff, monom in zip(self.coeffs, self.monoms):
 777             yield coeff, monom
 778
 779     def iter_all_coeffs(self):
 780         if self.is_univariate:
 781             terms = self.as_uv_dict()
 782
 783             for i in xrange(self.degree, -1, -1):
 784                 if terms.has_key(i):
 785                     yield terms[i]
 786                 else:
 787                     yield S.Zero
 788         else:
 789             raise MultivariatePolyError(self)
 790
 791     def iter_all_monoms(self):
 792         if self.is_univariate:
 793             for i in xrange(self.degree, -1, -1):
 794                 yield (i,)
 795         else:
 796             raise MultivariatePolyError(self)
 797
 798     def iter_all_terms(self):
 799         if self.is_univariate:
 800             terms = self.as_uv_dict()
 801
 802             for i in xrange(self.degree, -1, -1):
 803                 if terms.has_key(i):
 804                     yield terms[i], (i,)
 805                 else:
 806                     yield S.Zero, (i,)
 807         else:
 808             raise MultivariatePolyError(self)
 809
 810     @property
 811     def lead_coeff(self):
 812         return self.coeffs[0]
 813
 814     @property
 815     def lead_monom(self):
 816         return self.monoms[0]
 817
 818     @property
 819     def lead_term(self):
 820         return (self.coeffs[0], self.monoms[0])
 821
 822     LC = lead_coeff
 823     LM = lead_monom
 824     LT = lead_term
 825
 826     @property
 827     def tail_coeff(self):
 828         return self.coeffs[-1]
 829
 830     @property
 831     def tail_monom(self):
 832         return self.monoms[-1]
 833
 834     @property
 835     def tail_term(self):
 836         return (self.coeffs[-1], self.monoms[-1])
 837
 838     TC = tail_coeff
 839     TM = tail_monom
 840     TT = tail_term
 841
 842     def __abs__(self):
 843         return self.__class__(([ abs(coeff) for coeff in self.coeffs ],
 844             self.monoms), *self.symbols, **self.flags)
 845
 846     def __neg__(self):
 847         return self.__class__(([ -coeff for coeff in self.coeffs ],
 848             self.monoms), *self.symbols, **self.flags)
 849
 850     def __add__(self, other):
 851         try:
 852             self, poly = self.unify_with(other)
 853         except PolynomialError:
 854             return self.as_basic() + other.as_basic()
 855
 856         if poly.is_zero:
 857             return self
 858
 859         if self.is_zero:
 860             return poly
 861
 862         if poly.is_monomial:
 863             return self.add_term(*poly.LT)
 864
 865         if self.is_monomial:
 866             return poly.add_term(*self.LT)
 867
 868         terms = self.as_dict()
 869
 870         for coeff, monom in poly.iter_terms():
 871             if terms.has_key(monom):
 872                 coeff += terms[monom]
 873
 874                 coeff = Poly.cancel(coeff)
 875
 876                 if not coeff:
 877                     del terms[monom]
 878                     continue
 879
 880             terms[monom] = coeff
 881
 882         return self.__class__(terms, *self.symbols, **self.flags)
 883
 884     def __sub__(self, other):
 885         try:
 886             self, poly = self.unify_with(other)
 887         except PolynomialError:
 888             return self.as_basic() - other.as_basic()
 889
 890         if poly.is_zero:
 891             return self
 892
 893         if self.is_zero:
 894             return -poly
 895
 896         if poly.is_monomial:
 897             return self.sub_term(*poly.LT)
 898
 899         if self.is_monomial:
 900             return (-poly).add_term(*self.LT)
 901
 902         terms = self.as_dict()
 903
 904         for coeff, monom in poly.iter_terms():
 905             if terms.has_key(monom):
 906                 coeff -= terms[monom]
 907
 908                 coeff = Poly.cancel(coeff)
 909
 910                 if not coeff:
 911                     del terms[monom]
 912                     continue
 913
 914             terms[monom] = -coeff
 915
 916         return self.__class__(terms, *self.symbols, **self.flags)
 917
 918     def __mul__(self, other):
 919         """Efficiently multiply sparse and dense polynomials.
 920
 921            Given polynomials P and Q, if both of them are dense, then
 922            in univariate case perform dense multiplication, otherwise
 923            apply Kronecker 'trick', mapping multivariate polynomials
 924            to univariate ones in last variable and then multiply.
 925
 926            If any of the polynomials is sparse then in both univariate
 927            and multivariate cases use naive multiplication algorithm.
 928
 929            For more information on implemented algorithms refer to:
 930
 931            [1] R. T. Moenck, Practical Fast Polynomial Multiplication,
 932                Proceedings of the third ACM symposium on Symbolic and
 933                algebraic computation, ACM Press, Yorktown Heights,
 934                New York, USA, 1976, pp. 136-148
 935         """
 936         try:
 937             self, poly = self.unify_with(other)
 938         except PolynomialError:
 939             return self.as_basic() * other.as_basic()
 940
 941         if self.is_constant:
 942             if self.is_zero:
 943                 return self
 944             elif self.is_one:
 945                 return poly
 946             else:
 947                 return poly.mul_term(self.LC)
 948
 949         if poly.is_constant:
 950             if poly.is_zero:
 951                 return poly
 952             elif poly.is_one:
 953                 return self
 954             else:
 955                 return self.mul_term(poly.LC)
 956
 957         if self.is_monomial:
 958             return poly.mul_term(*self.LT)
 959
 960         if poly.is_monomial:
 961             return self.mul_term(*poly.LT)
 962
 963         if self.is_dense and poly.is_dense:
 964             if self.is_multivariate:
 965                 a = monomial_max(*self.monoms)
 966                 b = monomial_max(*poly.monoms)
 967
 968                 degs, p, q = [1], {}, {}
 969
 970                 for x, y in reversed(zip(a, b)[1:]):
 971                     degs.insert(0, degs[0]*(x+y+1))
 972
 973                 N = sum([ x*y for x, y in zip(self.LM, degs) ])
 974                 M = sum([ x*y for x, y in zip(poly.LM, degs) ])
 975
 976                 for coeff, monom in self.iter_terms():
 977                     p[sum([ x*y for x, y in zip(monom, degs) ])] = coeff
 978
 979                 for coeff, monom in poly.iter_terms():
 980                     q[sum([ x*y for x, y in zip(monom, degs) ])] = coeff
 981             else:
 982                 p, N = self.as_uv_dict(), self.degree
 983                 q, M = poly.as_uv_dict(), poly.degree
 984
 985             coeffs, monoms = [], []
 986
 987             for k in xrange(N+M, -1, -1):
 988                 coeff = 0
 989
 990                 for i in xrange(k+1):
 991                     if p.has_key(i) and q.has_key(k-i):
 992                         coeff += Poly.cancel(p[i] * q[k-i])
 993
 994                 if coeff:
 995                     coeffs.append(coeff)
 996                     monoms.append((k,))
 997
 998             if self.is_multivariate:
 999                 terms, degs = {}, degs[:-1]
1000
1001                 for i, d in enumerate(monoms):
1002                     monom, d = [], d[0]
1003
1004                     for j in degs:
1005                         k, d = divmod(d, j)
1006                         monom.append(k)
1007
1008                     terms[tuple(monom)+(d,)] = coeffs[i]
1009             else:
1010                 terms = coeffs, monoms
1011         else:
1012             terms = {}
1013
1014             for coeff_p, M in self.iter_terms():
1015                 for coeff_q, N in poly.iter_terms():
1016                     monom = monomial_mul(M, N)
1017
1018                     coeff = coeff_p * coeff_q
1019                     coeff = Poly.cancel(coeff)
1020
1021                     if terms.has_key(monom):
1022                         coeff += terms[monom]
1023
1024                         if not coeff:
1025                             del terms[monom]
1026                             continue
1027
1028                     terms[monom] = coeff
1029
1030         return self.__class__(terms, *self.symbols, **self.flags)
1031
1032     def __pow__(self, other):
1033         """Polynomial exponentiation using binary method.
1034
1035            Given a polynomial P and integer exponent N, raise P to power
1036            N in floor(lg(N)) + 1(N) multiplications,  where lg(N) stands
1037            for binary logarithm and 1(N) is the number of ones in binary
1038            representation of N.
1039
1040            For more information on the implemented algorithm refer to:
1041
1042            [1] D. E. Knuth, The Art of Computer Programming: Seminumerical
1043                Algorithms, v.2, Addison-Wesley Professional, 1997, pp. 461
1044
1045         """
1046         other = sympify(other)
1047
1048         if isinstance(other, Poly):
1049             if other.is_constant:
1050                 other = sympify(other.lead_coeff)
1051             else:
1052                 return self.as_basic()**other.as_basic()
1053
1054         if other.is_Integer:
1055             if other is S.One:
1056                 return self
1057
1058             P = self.__class__(S.One, *self.symbols, **self.flags)
1059
1060             if other is S.Zero:
1061                 return P
1062
1063             N, Q = abs(int(other)), self
1064
1065             while True:
1066                 N, M = N/2, N
1067
1068                 if M & 1:
1069                     P *= Q
1070
1071                     if N == 0:
1072                         break
1073
1074                 Q *= Q
1075
1076             if other.is_positive:
1077                 return P
1078             else:
1079                 return Pow(P, S.NegativeOne)
1080         else:
1081             return self.as_basic()**other
1082
1083     def __div__(self, other):
1084         try:
1085             return sympy.polys.algorithms.poly_div(self, other)[0]
1086         except PolynomialError:
1087             return self.as_basic() / other.as_basic()
1088
1089     def __mod__(self, other):
1090         try:
1091             return sympy.polys.algorithms.poly_div(self, other)[1]
1092         except PolynomialError:
1093             return S.Zero
1094
1095     def __divmod__(self, other):
1096         try:
1097             return sympy.polys.algorithms.poly_div(self, other)
1098         except PolynomialError:
1099             return (self.as_basic() / other.as_basic(), S.Zero)
1100
1101     __truediv__ = __div__
1102
1103     @property
1104     def is_zero(self):
1105         return self.coeffs in ((S.Zero,), (0,))
1106
1107     @property
1108     def is_one(self):
1109         return self.coeffs in ((S.One,), (1,)) and \
1110             all(e == 0 for e in self.monoms[0])
1111
1112     @property
1113     def is_number(self):
1114         return self.is_constant and self.coeffs[0].is_number
1115
1116     @property
1117     def is_constant(self):
1118         return len(self.monoms) == 1 and \
1119             all(e == 0 for e in self.monoms[0])
1120
1121     @property
1122     def is_monomial(self):
1123         return len(self.monoms) == 1
1124
1125     @property
1126     def is_univariate(self):
1127         return len(self.symbols) == 1
1128
1129     @property
1130     def is_multivariate(self):
1131         return len(self.symbols) != 1
1132
1133     @property
1134     def is_homogeneous(self):
1135         return any(e != 0 for e in self.monoms[-1])
1136
1137     @property
1138     def is_inhomogeneous(self):
1139         return all(e == 0 for e in self.monoms[-1])
1140
1141     @property
1142     def is_sparse(self):
1143         return not self.is_dense
1144
1145     @property
1146     def is_dense(self):
1147         """Returns True if 'self' is dense.
1148
1149            Let 'self' be a polynomial in M variables of a total degree N
1150            and having L terms (with non-zero coefficients). Let K be the
1151            number of monomials in M variables of degree at most N.  Then
1152            'self' is considered dense if it's at least 90% filled.
1153
1154            The total number of monomials is given by (M + N)! / (M! N!),
1155            where the factorials are estimated using improved Stirling's
1156            approximation:
1157
1158                        n! = sqrt((2 n + 1/3) pi) n**n exp(-n)
1159
1160            For more information on this subject refer to:
1161
1162              http://mathworld.wolfram.com/StirlingsApproximation.html
1163
1164         """
1165         def enum(M, N):
1166             U = float(M+N)
1167             V = (6*M+1)*(6*N+1)
1168             S = 2.4*math.sqrt(U/V)
1169             A = math.pow(U/M, M)
1170             B = math.pow(U/N, N)
1171             return S * A * B
1172
1173         V, N = len(self.symbols), self.total_degree
1174
1175         return (self.length / enum(V, N)) > 0.9
1176
1177     @property
1178     def is_primitive(self):
1179         return self.content in (S.One, 1)
1180
1181     @property
1182     def is_monic(self):
1183         return self.lead_coeff in (S.One, 1)
1184
1185     @property
1186     def is_squarefree(self):
1187         """Returns True if 'self' has no factors of multiplicity >= 2.
1188
1189            >>> from sympy import *
1190            >>> x,y = symbols('xy')
1191
1192            >>> Poly(x-1, x).is_squarefree
1193            True
1194
1195            >>> Poly((x-1)**2, x).is_squarefree
1196            False
1197
1198         """
1199         return self.as_squarefree() is self
1200
1201     @property
1202     def is_linear(self):
1203         return all(sum(monom) <= 1 for monom in self.iter_monoms())
1204
1205     @property
1206     def has_integer_coeffs(self):
1207         return self.domain == 'Z'
1208
1209     @property
1210     def has_rational_coeffs(self):
1211         return self.domain == 'Q'
1212
1213     @property
1214     def has_radical_coeffs(self):
1215         return self.domain == 'R'
1216
1217     @property
1218     def has_complex_coeffs(self):
1219         return self.domain == 'C'
1220
1221     @property
1222     def has_symbolic_coeffs(self):
1223         return self.domain == 'S'
1224
1225     def as_monic(self):
1226         """Returns a monic polynomial.
1227
1228            >>> from sympy import *
1229            >>> x,y = symbols('xy')
1230
1231            >>> Poly(x**2 + 4, x).as_monic()
1232            Poly(x**2 + 4, x)
1233
1234            >>> Poly(2*x**2 + 4, x).as_monic()
1235            Poly(x**2 + 2, x)
1236
1237            >>> Poly(y*x**2 + y**2 + y, x).as_monic()
1238            Poly(x**2 + 1 + y, x)
1239
1240         """
1241         LC = self.lead_coeff
1242
1243         if not self or LC is S.One:
1244             return self
1245
1246         coeffs = [ coeff / LC for coeff in self.coeffs ]
1247
1248         for i, coeff in enumerate(coeffs):
1249             coeffs[i] = Poly.cancel(coeff)
1250
1251         return self.__class__((coeffs, self.monoms),
1252             *self.symbols, **self.flags)
1253
1254     def as_integer(self):
1255         """Makes all coefficients integers if possible.
1256
1257            Given a polynomial with rational coefficients,  returns a tuple
1258            consisting of a common denominator of those coefficients and an
1259            integer polynomial. Otherwise an exception is being raised.
1260
1261            >>> from sympy import *
1262            >>> x,y = symbols("xy")
1263
1264            >>> Poly(3*x**2 + x, x).as_integer()
1265            (1, Poly(3*x**2 + x, x))
1266
1267            >>> Poly(3*x**2 + x/2, x).as_integer()
1268            (2, Poly(6*x**2 + x, x))
1269
1270         """
1271         denom = 1
1272
1273         for coeff in self.iter_coeffs():
1274             if coeff.is_Rational:
1275                 if not coeff.is_Integer:
1276                     denom = ilcm(denom, coeff.q)
1277             else:
1278                 raise CoefficientError("%s is not a rational number" % coeff)
1279
1280         denom = sympify(denom)
1281
1282         if denom is S.One:
1283             return denom, self
1284         else:
1285             coeffs = [ coeff * denom for coeff in self.iter_coeffs() ]
1286
1287             return denom, self.__class__((coeffs, self.monoms),
1288                 *self.symbols, **self.flags)
1289
1290     @property
1291     def content(self):
1292         """Returns integer GCD of all the coefficients.
1293
1294            >>> from sympy import *
1295            >>> x,y,z = symbols('xyz')
1296
1297            >>> p = Poly(2*x + 5*x*y, x, y)
1298            >>> p.content
1299            1
1300
1301            >>> p = Poly(2*x + 4*x*y, x, y)
1302            >>> p.content
1303            2
1304
1305            >>> p = Poly(2*x + z*x*y, x, y)
1306            >>> p.content
1307            1
1308
1309         """
1310         if self.is_zero:
1311             return S.One
1312         else:
1313             content = 0
1314
1315             for coeff in self.coeffs:
1316                 if coeff.is_Rational:
1317                     content = igcd(content, coeff.p)
1318                 else:
1319                     return S.One
1320
1321             return Integer(content)
1322
1323     def as_primitive(self):
1324         """Returns content and primitive part of a polynomial.
1325
1326            >>> from sympy import *
1327            >>> x,y = symbols('xy')
1328
1329            >>> p = Poly(4*x**2 + 2*x, x)
1330            >>> p.as_primitive()
1331            (2, Poly(2*x**2 + x, x))
1332
1333         """
1334         content = self.content
1335
1336         if content is S.Zero or content is S.One:
1337             return content, self
1338         else:
1339             coeffs = [ coeff / content for coeff in self.coeffs ]
1340
1341             return content, self.__class__((coeffs,
1342                 self.monoms), *self.symbols, **self.flags)
1343
1344     def as_squarefree(self):
1345         """Returns square-free part of a polynomial.
1346
1347            >>> from sympy import *
1348            >>> x,y = symbols('xy')
1349
1350            >>> Poly((x-1)**2, x).as_squarefree()
1351            Poly(x - 1, x)
1352
1353         """
1354         f, A = self, sympy.polys.algorithms
1355
1356         if f.is_univariate:
1357             F = f.as_primitive()[1]
1358
1359             h = f.diff()
1360
1361             g = A.poly_gcd(F, h)
1362             r = A.poly_div(f, g)
1363
1364             return r[0]
1365         else:
1366             raise MultivariatePolyError(f)
1367
1368     def as_reduced(self):
1369         """Remove GCD of monomials from 'self'.
1370
1371            >>> from sympy import *
1372            >>> x,y = symbols('xy')
1373
1374            >>> Poly(x**3 + x, x).as_reduced()
1375            ((1,), Poly(x**2 + 1, x))
1376
1377            >>> Poly(x**3*y+x**2*y**2, x, y).as_reduced()
1378            ((2, 1), Poly(x + y, x, y))
1379
1380         """
1381         if self.is_inhomogeneous:
1382             return (0,)*len(self.symbols), self
1383         else:
1384             if self.is_univariate:
1385                 gcd = self.monoms[-1]
1386
1387                 terms = self.coeffs, [ (monom - gcd[0],)
1388                     for (monom,) in self.monoms ]
1389             else:
1390                 gcd = monomial_min(*self.monoms)
1391
1392                 if all(not n for n in gcd):
1393                     return gcd, self
1394
1395                 terms = {}
1396
1397                 for coeff, monom in self.iter_terms():
1398                     terms[monomial_div(monom, gcd)] = coeff
1399
1400             return gcd, Poly(terms, *self.symbols, **self.flags)
1401
1402     def map_coeffs(self, f, *args, **kwargs):
1403         """Apply a function to all the coefficients.
1404
1405            >>> from sympy import *
1406            >>> x,y,u,v = symbols('xyuv')
1407
1408            >>> p = Poly(x**2 + 2*x*y, x, y)
1409            >>> q = p.map_coeffs(lambda c: 2*c)
1410
1411            >>> q.as_basic()
1412            4*x*y + 2*x**2
1413
1414            >>> p = Poly(u*x**2 + v*x*y, x, y)
1415            >>> q = p.map_coeffs(expand, complex=True)
1416
1417            >>> q.as_basic()
1418            x**2*(I*im(u) + re(u)) + x*y*(I*im(v) + re(v))
1419
1420         """
1421         terms = []
1422
1423         for coeff, monom in self.iter_terms():
1424             coeff = f(coeff, *args, **kwargs)
1425
1426             if coeff.has_any_symbols(*self.symbols):
1427                 raise CoefficientError("%s coefficient is dependent" \
1428                     " of polynomial's symbols %s" % (coeff, self.symbols))
1429             elif coeff is not S.Zero:
1430                 terms.append((coeff, monom))
1431
1432         return self.__class__(terms, *self.symbols, **self.flags)
1433
1434     def coeff(self, *monom):
1435         """Returns coefficient of a particular monomial.
1436
1437            >>> from sympy import *
1438            >>> x,y = symbols('xy')
1439
1440            >>> p = Poly(3*x**2*y + 4*x*y**2 + 1, x, y)
1441
1442            >>> p.coeff(2, 1)
1443            3
1444            >>> p.coeff(1, 2)
1445            4
1446            >>> p.coeff(1, 1)
1447            0
1448
1449            If len(monom) == 0 then returns coeff of the constant term:
1450
1451            >>> p.coeff(0, 0)
1452            1
1453            >>> p.coeff()
1454            1
1455
1456         """
1457         if not monom:
1458             if all(e == 0 for e in self.monoms[-1]):
1459                 return self.coeffs[-1]
1460         else:
1461             for i in xrange(len(self.monoms)):
1462                 if self.monoms[i] == monom:
1463                     return self.coeffs[i]
1464
1465         return S.Zero
1466
1467     def unify_with(self, other):
1468         """Unify 'self' with a polynomial or a set of polynomials.
1469
1470            This method will return polynomials of the same type,  dominated
1471            by 'self', with a common set of symbols (which is a union of all
1472            symbols from all polynomials) and with common monomial ordering.
1473
1474            You can pass a polynomial  or an expression to this method, or
1475            a list or a tuple of polynomials or expressions. If any of the
1476            inputs would be an expression then it will be converted to a
1477            polynomial.
1478
1479         """
1480         symbols, stamp = self.symbols, self.stamp
1481         flags, cls = self.flags, self.__class__
1482
1483         if isinstance(other, (tuple, list, set)):
1484             for poly in other:
1485                 if isinstance(poly, Poly):
1486                     stamp |= poly.stamp
1487                 #elif atoms:
1488                 #    stamp |= poly.atoms(Symbol)
1489
1490             if not (stamp <= self.stamp):
1491                 symbols = tuple(sorted(stamp))
1492                 self = cls(self, *symbols, **flags)
1493
1494             other = other.__class__( cls(poly,
1495                 *symbols, **flags) for poly in other )
1496         else:
1497             other = sympify(other)
1498
1499             if isinstance(other, Poly):
1500                 stamp |= other.stamp
1501             #elif atoms:
1502             #    stamp |= other.atoms(Symbol)
1503
1504             if not (stamp <= self.stamp):
1505                 symbols = tuple(sorted(stamp))
1506                 self = cls(self, *symbols, **flags)
1507
1508             other = cls(other, *symbols, **flags)
1509
1510         return self, other
1511
1512     def diff(self, *symbols):
1513         """Efficiently differentiate polynomials.
1514
1515            Differentiate a polynomial with respect to a set of symbols.  If
1516            a symbol is polynomial's variable, then monomial differentiation
1517            is performed and coefficients updated with integer factors. In
1518            other case each of the coefficients is differentiated.
1519
1520            Additionally, for each of the symbols you can specify a single
1521            positive integer which will indicate the number of times to
1522            perform differentiation using this symbol.
1523
1524            >>> from sympy import *
1525            >>> x,y,z = symbols('xyz')
1526
1527            >>> p = Poly(x**2*y + z**2, x, y)
1528
1529            >>> p.diff(x)
1530            Poly(2*x*y, x, y)
1531
1532            >>> p.diff(x, 2)
1533            Poly(2*y, x, y)
1534
1535            >>> p.diff(x, 2, y)
1536            Poly(2, x, y)
1537
1538            >>> p.diff(z)
1539            Poly(2*z, x, y)
1540
1541         """
1542         if self.is_zero:
1543             return self
1544
1545         new_symbols = []
1546
1547         for s in symbols:
1548             s = sympify(s)
1549
1550             if s.is_Symbol:
1551                 new_symbols.append((s, 1))
1552             elif s.is_Integer:
1553                 sym, _ = new_symbols.pop()
1554                 new_symbols.append((sym, int(s)))
1555             else:
1556                 raise TypeError
1557
1558         if not new_symbols:
1559             if self.is_univariate:
1560                 new_symbols = [(self.symbols[0], 1)]
1561             else:
1562                 return self
1563
1564         indices, symbols = {}, self.stamp
1565
1566         for i, s in enumerate(self.symbols):
1567             indices[s] = i
1568
1569         poly = self.as_dict()
1570
1571         for s, k in new_symbols:
1572             new_poly = {}
1573
1574             if s in symbols:
1575                 i = indices[s]
1576
1577                 for M, coeff in poly.iteritems():
1578                     n = M[i] - k
1579
1580                     if n >= 0:
1581                         monom = M[:i]+(n,)+M[i+1:]
1582
1583                         for j in xrange(n, M[i]):
1584                             coeff *= j+1
1585
1586                         if new_poly.has_key(monom):
1587                             coeff += new_poly[monom]
1588
1589                             if not coeff:
1590                                 del new_poly[monom]
1591                                 continue
1592
1593                         new_poly[monom] = coeff
1594             elif not isinstance(self, IntegerPoly):
1595                 for monom, coeff in poly.iteritems():
1596                     if coeff.has_any_symbols(s):
1597                         coeff = coeff.diff(*([s]*k))
1598
1599                         if coeff:
1600                             new_poly[monom] = coeff
1601
1602             poly = new_poly
1603
1604             if not poly:
1605                 break
1606
1607         return self.__class__(poly, *self.symbols, **self.flags)
1608
1609     def integrate(self, *symbols):
1610         """Efficiently integrate polynomials.
1611
1612            Integrate a polynomial with respect a set of symbols. If a
1613            symbol is polynomial's variable, then monomial integration
1614            is performed and coefficients updated with integer factors.
1615            In other case each of the coefficients is integrated.
1616
1617            Additionally, for each of the symbols you can specify a
1618            single positive integer which will indicate the number
1619            of times to perform integration using this symbol.
1620
1621            >>> from sympy import *
1622            >>> x,y,z = symbols('xyz')
1623
1624            >>> p = Poly(x**2*y + z**2, x, y)
1625
1626            >>> p.integrate(x)
1627            Poly(1/3*x**3*y + z**2*x, x, y)
1628
1629            >>> p.integrate(x, 2)
1630            Poly(1/12*x**4*y + 1/2*z**2*x**2, x, y)
1631
1632            >>> p.integrate(x, 2, y)
1633            Poly(1/24*x**4*y**2 + 1/2*z**2*x**2*y, x, y)
1634
1635            >>> p.integrate(z)
1636            Poly(z*x**2*y + 1/3*z**3, x, y)
1637
1638         """
1639         if self.is_zero:
1640             return self
1641
1642         new_symbols = []
1643
1644         for s in symbols:
1645             s = sympify(s)
1646
1647             if s.is_Symbol:
1648                 new_symbols.append((s, 1))
1649             elif s.is_Integer:
1650                 sym, _ = new_symbols.pop()
1651                 new_symbols.append((sym, int(s)))
1652             else:
1653                 raise TypeError
1654
1655         indices, symbols = {}, self.stamp
1656
1657         for i, s in enumerate(self.symbols):
1658             indices[s] = i
1659
1660         poly = self.as_dict()
1661
1662         for s, k in new_symbols:
1663             new_poly = {}
1664
1665             if s in symbols:
1666                 i = indices[s]
1667
1668                 for M, coeff in poly.iteritems():
1669                     n = M[i] + k
1670
1671                     monom = M[:i]+(n,)+M[i+1:]
1672
1673                     for j in xrange(M[i], n):
1674                         coeff /= j+1
1675
1676                     if new_poly.has_key(monom):
1677                         coeff += new_poly[monom]
1678
1679                         if not coeff:
1680                             del new_poly[monom]
1681                             continue
1682
1683                     new_poly[monom] = coeff
1684             else:
1685                 for M, coeff in poly.iteritems():
1686                     new_poly[M] = C.Integral(coeff, *([s]*k)).doit()
1687
1688             if not new_poly:
1689                 break
1690             else:
1691                 poly = new_poly
1692
1693         return self.__class__(poly, *self.symbols, **self.flags)
1694
1695     def add_term(self, coeff, monom):
1696         """Efficiently add a single term to 'self'.
1697
1698            The list of monomials is sorted at initialization time, this
1699            motivates usage of binary search algorithm to find an index
1700            of an existing monomial or a suitable place for a new one.
1701            This gives O(lg(n)) complexity, where 'n' is the initial
1702            number of terms, superior to naive approach.
1703
1704            For more information on the implemented algorithm refer to:
1705
1706            [1] D. E. Knuth, The Art of Computer Programming: Sorting
1707                and Searching, v.1, Addison-Wesley Professional, 1998
1708
1709         """
1710         if self.is_zero:
1711             coeffs = (coeff,)
1712             monoms = (monom,)
1713         else:
1714             coeffs = list(self.coeffs)
1715             monoms = list(self.monoms)
1716
1717             compare = monomial_cmp(self.order)
1718
1719             if compare(monom, monoms[0]) > 0:
1720                 coeffs.insert(0, coeff)
1721                 monoms.insert(0, monom)
1722             elif compare(monom, monoms[-1]) < 0:
1723                 coeffs.append(coeff)
1724                 monoms.append(monom)
1725             else:
1726                 lo, hi = 0, len(monoms)-1
1727
1728                 while lo <= hi:
1729                     i = (lo + hi) // 2
1730
1731                     k = compare(monom, monoms[i])
1732
1733                     if not k:
1734                         coeff += coeffs[i]
1735
1736                         if coeff:
1737                             coeffs[i] = coeff
1738                         else:
1739                             del coeffs[i]
1740                             del monoms[i]
1741
1742                         break
1743                     else:
1744                         if k > 0:
1745                             hi = i - 1
1746                         else:
1747                             lo = i + 1
1748                 else:
1749                     coeffs.insert(i, coeff)
1750                     monoms.insert(i, monom)
1751
1752         return self.__class__((coeffs, monoms),
1753             *self.symbols, **self.flags)
1754
1755     def sub_term(self, coeff, monom):
1756         """Efficiently subtract a single term from 'self'.
1757
1758            The list of monomials is sorted at initialization time, this
1759            motivates usage of binary search algorithm to find an index
1760            of an existing monomial or a suitable place for a new one.
1761            This gives O(lg(n)) complexity, where 'n' is the initial
1762            number of terms, superior to naive approach.
1763
1764            For more information on the implemented algorithm refer to:
1765
1766            [1] D. E. Knuth, The Art of Computer Programming: Sorting
1767                and Searching, v.2, Addison-Wesley Professional, 1998
1768
1769         """
1770         if self.is_zero:
1771             coeffs = (-coeff,)
1772             monoms = (monom,)
1773         else:
1774             coeffs = list(self.coeffs)
1775             monoms = list(self.monoms)
1776
1777             compare = monomial_cmp(self.order)
1778
1779             if compare(monom, monoms[0]) > 0:
1780                 coeffs.insert(0, -coeff)
1781                 monoms.insert(0, monom)
1782             elif compare(monom, monoms[-1]) < 0:
1783                 coeffs.append(-coeff)
1784                 monoms.append(monom)
1785             else:
1786                 lo, hi = 0, len(monoms)-1
1787
1788                 while lo <= hi:
1789                     i = (lo + hi) // 2
1790
1791                     k = compare(monom, monoms[i])
1792
1793                     if not k:
1794                         coeff -= coeffs[i]
1795
1796                         if coeff:
1797                             coeffs[i] = -coeff
1798                         else:
1799                             del coeffs[i]
1800                             del monoms[i]
1801
1802                         break
1803                     else:
1804                         if k > 0:
1805                             hi = i - 1
1806                         else:
1807                             lo = i + 1
1808                 else:
1809                     coeffs.insert(i, -coeff)
1810                     monoms.insert(i, monom)
1811
1812         return self.__class__((coeffs, monoms),
1813             *self.symbols, **self.flags)
1814
1815     def mul_term(self, coeff, monom=None):
1816         """Efficiently multiply 'self' with a single term. """
1817         coeff = sympify(coeff)
1818
1819         if coeff is S.Zero:
1820             terms = ()
1821         else:
1822             if coeff is S.One:
1823                 if monom is None:
1824                     return self
1825                 else:
1826                     coeffs = self.coeffs
1827             else:
1828                 coeffs = [ p_coeff * coeff for p_coeff in self.coeffs ]
1829
1830                 for i, coeff in enumerate(coeffs):
1831                     coeffs[i] = Poly.cancel(coeff)
1832
1833             if monom is None:
1834                 monoms = self.monoms
1835             else:
1836                 monoms = [ monomial_mul(p_monom, monom)
1837                     for p_monom in self.monoms ]
1838
1839             terms = coeffs, monoms
1840
1841         return self.__class__(terms, *self.symbols, **self.flags)
1842
1843     def div_term(self, coeff, monom=None):
1844         """Efficiently divide 'self' by a single term. """
1845         coeff = sympify(coeff)
1846
1847         if coeff is S.Zero:
1848             raise ZeroDivisionError
1849
1850         if coeff is S.One:
1851             if monom is None:
1852                 return self
1853             else:
1854                 coeffs = self.coeffs
1855         else:
1856             coeffs = [ p_coeff / coeff for p_coeff in self.coeffs ]
1857
1858             for i, coeff in enumerate(coeffs):
1859                 coeffs[i] = Poly.cancel(coeff)
1860
1861         if monom is None:
1862             monoms = self.monoms
1863         else:
1864             monoms = [ monomial_div(p_monom, monom)
1865                 for p_monom in self.monoms ]
1866
1867             if any(monom is None for monom in monoms):
1868                 raise PolynomialError("%s monomial must divide" \
1869                     " exactly the given polynomial" % (monom,))
1870
1871         return self.__class__((coeffs, monoms),
1872             *self.symbols, **self.flags)
1873
1874     def kill_lead_term(self):
1875         """Removes leading term from 'self'. """
1876         terms = self.coeffs[1:], self.monoms[1:]
1877         return self.__class__(terms, *self.symbols, **self.flags)
1878
1879     def kill_last_term(self):
1880         """Removes last term from 'self'. """
1881         terms = self.coeffs[:-1], self.monoms[:-1]
1882         return self.__class__(terms, *self.symbols, **self.flags)
1883
1884     def __call__(self, *point):
1885         """Efficiently evaluate polynomial at a given point.
1886
1887            Evaluation is always done using Horner scheme.  In multivariate
1888            case a greedy algorithm is used to obtain a sequence of partial
1889            evaluations which minimizes the total number of multiplications
1890            required to perform this evaluation. This strategy is efficient
1891            for most of multivariate polynomials.
1892
1893            Note that evaluation is done for all variables, which means the
1894            dimension of the given point must match the number of symbols.
1895
1896            If you wish to efficiently evaluate polynomial for a subset of
1897            symbols use 'evaluate' method instead. Alternatively one can
1898            use Basic.subs() for this purpose.
1899
1900            >>> from sympy import *
1901            >>> x,y = symbols('xy')
1902
1903            >>> p = Poly(x**2 + 2*x + 1, x)
1904
1905            >>> p(2)
1906            9
1907            >>> p(y)
1908            1 + y*(2 + y)
1909
1910            For more information on the implemented algorithm refer to:
1911
1912            [1] M. Ceberio, V. Kreinovich, Greedy Algorithms for Optimizing
1913                Multivariate Horner Schemes, ACM SIGSAM Bulletin, Volume 38,
1914                Issue 1, 2004, pp. 8-15
1915
1916         """
1917         N = len(self.symbols)
1918
1919         if len(point) != N:
1920             raise ValueError("Dimension of %s does not match" \
1921                 " the number of symbols (%d)" % (point, N))
1922
1923         if self.is_univariate:
1924             terms = self.as_uv_dict()
1925
1926             point, result = point[0], 0
1927
1928             for k in xrange(self.degree, -1, -1):
1929                 result *= point
1930
1931                 if terms.has_key(k):
1932                     result += terms[k]
1933
1934             return result
1935         else:
1936             def evaluate(terms):
1937                 count = [0] * N
1938
1939                 for monom in terms.iterkeys():
1940                     for i, M in enumerate(monom):
1941                         if M != 0:
1942                             count[i] += 1
1943
1944                 K = max(count)
1945
1946                 if K <= 1:
1947                     result = 0
1948
1949                     for monom, coeff in terms.iteritems():
1950                         for base, exp in zip(point, monom):
1951                             if exp != 0:
1952                                 if exp == 1:
1953                                     coeff *= base
1954                                 else:
1955                                     coeff *= base**exp
1956
1957                         result += coeff
1958
1959                     return result
1960                 else:
1961                     k, indeps, depend = count.index(K), {}, {}
1962
1963                     n = min([ M[k] for M in terms.iterkeys() if M[k] ])
1964
1965                     for M, coeff in terms.iteritems():
1966                         if M[k] != 0:
1967                             depend[M[:k]+(M[k]-n,)+M[k+1:]] = coeff
1968                         else:
1969                             indeps[M] = coeff
1970
1971                     result = point[k]**n * evaluate(depend)
1972
1973                     if indeps:
1974                         return result + evaluate(indeps)
1975                     else:
1976                         return result
1977
1978             return evaluate(self.as_dict())
1979
1980     def evaluate(self, pattern):
1981         """Evaluates polynomial for a given set of symbols. """
1982         symbols = list(self.symbols)
1983
1984         if type(pattern) is dict:
1985             pattern = pattern.items()
1986         elif type(pattern) is tuple:
1987             pattern = [pattern]
1988
1989         poly = self.as_dict()
1990
1991         for s, value in pattern:
1992             if s not in self.stamp:
1993                 raise PolynomialError("%s symbol does" \
1994                     " not belong to %s" % (s, self.symbols))
1995             else:
1996                 i = symbols.index(s)
1997
1998             # 'value' might be int | long
1999             if isinstance(value, Symbol):
2000                 if value in self.stamp:
2001                     raise PolynomialError("%s symbol must" \
2002                         " not belong to %s" % (s, self.symbols))
2003
2004             terms = {}
2005
2006             for M, coeff in poly.iteritems():
2007                 monom = M[:i] + M[i+1:]
2008                 coeff *= value ** M[i]
2009
2010                 coeff = Poly.cancel(coeff)
2011
2012                 if terms.has_key(monom):
2013                     coeff += terms[monom]
2014
2015                     if not coeff:
2016                         del terms[monom]
2017                         continue
2018                 elif not coeff:
2019                     continue
2020
2021                 terms[monom] = coeff
2022
2023             del symbols[i]
2024             poly = terms
2025
2026         if not symbols:
2027             if not poly:
2028                 return S.Zero
2029             else:
2030                 return poly.popitem()[1]
2031         else:
2032             return self.__class__(poly, *symbols, **self.flags)
2033
2034     def _eval_subs(self, old, new):
2035         symbols = list(self.symbols)
2036
2037         if old in self.stamp:
2038             terms, i = {}, symbols.index(old)
2039
2040             if new.is_Symbol:
2041                 if new in self.stamp:
2042                     j = symbols.index(new)
2043
2044                     for coeff, M in self.iter_terms():
2045                         n, m = M[i], M[j]
2046
2047                         if i < j:
2048                             monom = M[:i]+M[i+1:j]+(m+n,)+M[j+1:]
2049                         else:
2050                             monom = M[:j]+(m+n,)+M[j+1:i]+M[i+1:]
2051
2052                         if terms.has_key(monom):
2053                             coeff += terms[monom]
2054
2055                             if not coeff:
2056                                 del terms[monom]
2057                                 continue
2058
2059                         terms[monom] = coeff
2060
2061                     del symbols[i]
2062
2063                     return self.__class__(terms, *symbols, **self.flags)
2064                 else:
2065                     for coeff in self.coeffs:
2066                         if coeff.has_any_symbols(new):
2067                             break
2068                     else:
2069                         symbols[i], terms = new, (self.coeffs, self.monoms)
2070                         return self.__class__(terms, *symbols, **self.flags)
2071             elif new.is_number:
2072                 if len(symbols) == 1:
2073                     return self(new)
2074                 else:
2075                     return self.evaluate((old, new))
2076         elif not new.has_any_symbols(*symbols):
2077             coeffs = [ sympify(coeff).subs(old, new) for coeff in self.coeffs ]
2078             return self.__class__((coeffs, self.monoms), *symbols, **self.flags)
2079
2080         result = self.as_basic().subs(old, new)
2081
2082         try:
2083             return self.__class__(result, *symbols, **self.flags)
2084         except PolynomialError:
2085             return result
2086
2087     def __eq__(self, other):
2088         """Compare polynomials up to order of symbols and monomials. """
2089         try:
2090             poly = self.__class__(other, *self.symbols, **self.flags)
2091         except PolynomialError:
2092             return False
2093
2094         if self.length != poly.length:
2095             return False
2096
2097         if hash(self.monoms) != hash(poly.monoms):
2098             return False
2099
2100         if hash(self.coeffs) != hash(poly.coeffs):
2101             for a, b in zip(self.coeffs, poly.coeffs):
2102                 if a != b:
2103                     return False
2104
2105         return True
2106
2107     def __ne__(self, other):
2108         """Compare polynomials up to order of symbols and monomials. """
2109         try:
2110             poly = self.__class__(other, *self.symbols, **self.flags)
2111         except PolynomialError:
2112             return True
2113
2114         if self.length != poly.length:
2115             return True
2116
2117         if hash(self.monoms) != hash(poly.monoms):
2118             return True
2119
2120         if hash(self.coeffs) != hash(poly.coeffs):
2121             for a, b in zip(self.coeffs, poly.coeffs):
2122                 if a != b:
2123                     return True
2124
2125         return False
2126
2127     def __nonzero__(self):
2128         return self.coeffs not in ((S.Zero,), (0,))
2129
2130     def _eval_is_polynomial(self, symbols):
2131         try:
2132             self.__class__(self, *symbols, **self.flags)
2133         except PolynomialError:
2134             return False
2135
2136         return True
2137
2138 class IntegerPoly(Poly):
2139     pass
2140
2141 def multinomial_as_basic(multinomial, *symbols):
2142     """
2143     Converts the multinomial to Add/Mul/Pow instances.
2144
2145     multinomial is a dict of {powers: coefficient} pairs, powers is a tuple of
2146     python integers, coefficient is a python integer.
2147     """
2148     l = []
2149     for powers, k in multinomial.iteritems():
2150         term = [k]
2151         for i,e in enumerate(powers):
2152             term.append(Pow(symbols[i], powers[i]))
2153         l.append(Mul(*term))
2154     result = Add(*l)
2155     return result