sympy/polys/algorithms.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.symbol import Symbol
   6 from sympy.core.basic import Basic, S
   7 from sympy.core.numbers import Integer
   8 from sympy.core.sympify import sympify
   9
  10 from sympy.core.numbers import igcd, ilcm
  11
  12 from polynomial import Poly, SymbolsError, MultivariatePolyError
  13
  14 from monomial import monomial_cmp, monomial_lcm, \
  15     monomial_gcd, monomial_mul, monomial_div
  16
  17 from sympy.utilities.iterables import all, any
  18
  19 def poly_div(f, g, *symbols):
  20     """Generalized polynomial division with remainder.
  21
  22        Given polynomial f and a set of polynomials g = (g_1, ..., g_n)
  23        compute a set of quotients q = (q_1, ..., q_n) and remainder r
  24        such that f = q_1*f_1 + ... + q_n*f_n + r, where r = 0 or r is
  25        a completely reduced polynomial with respect to g.
  26
  27        In particular g can be a tuple, list or a singleton. All g_i
  28        and f can be given as Poly class instances or as expressions.
  29
  30        For more information on the implemented algorithm refer to:
  31
  32        [1] D. Cox, J. Little, D. O'Shea, Ideals, Varieties and
  33            Algorithms, Springer, Second Edition, 1997, pp. 62
  34
  35        [2] I.A. Ajwa, Z. Liu, P.S. Wang, Groebner Bases Algorithm,
  36            http://citeseer.ist.psu.edu/ajwa95grbner.html, 1995
  37
  38     """
  39     if not isinstance(f, Poly):
  40         f = Poly(f, *symbols)
  41     elif symbols:
  42         raise SymbolsError("Redundant symbols were given")
  43
  44     f, g = f.unify_with(g)
  45
  46     symbols, flags = f.symbols, f.flags
  47
  48     r = Poly((), *symbols, **flags)
  49
  50     if isinstance(g, Basic):
  51         if g.is_constant:
  52             if g.is_zero:
  53                 raise ZeroDivisionError
  54             elif g.is_one:
  55                 return f, r
  56             else:
  57                 return f.div_term(g.LC), r
  58
  59         if g.is_monomial:
  60             LC, LM = g.lead_term
  61
  62             q_coeffs, q_monoms = [], []
  63             r_coeffs, r_monoms = [], []
  64
  65             for coeff, monom in f.iter_terms():
  66                 quotient = monomial_div(monom, LM)
  67
  68                 if quotient is not None:
  69                     coeff /= LC
  70
  71                     q_coeffs.append(Poly.cancel(coeff))
  72                     q_monoms.append(quotient)
  73                 else:
  74                     r_coeffs.append(coeff)
  75                     r_monoms.append(monom)
  76
  77             return (Poly((q_coeffs, q_monoms), *symbols, **flags),
  78                     Poly((r_coeffs, r_monoms), *symbols, **flags))
  79
  80         g, q = [g], [r]
  81     else:
  82         q = [r] * len(g)
  83
  84     while not f.is_zero:
  85         for i, h in enumerate(g):
  86             monom = monomial_div(f.LM, h.LM)
  87
  88             if monom is not None:
  89                 coeff = Poly.cancel(f.LC / h.LC)
  90
  91                 q[i] = q[i].add_term(coeff, monom)
  92                 f -= h.mul_term(coeff, monom)
  93
  94                 break
  95         else:
  96             r = r.add_term(*f.LT)
  97             f = f.kill_lead_term()
  98
  99     if len(q) != 1:
 100         return q, r
 101     else:
 102         return q[0], r
 103
 104 def poly_pdiv(f, g, *symbols):
 105     """Univariate polynomial pseudo-division with remainder.
 106
 107        Given univariate polynomials f and g over an integral domain D[x]
 108        applying classical division algorithm to LC(g)**(d + 1) * f and g
 109        where  d = max(-1, deg(f) - deg(g)),  compute polynomials q and r
 110        such that LC(g)**(d + 1)*f = g*q + r and r = 0 or deg(r) < deg(g).
 111        Polynomials q and r are called the pseudo-quotient of f by g and
 112        the pseudo-remainder of f modulo g respectively.
 113
 114        For more information on the implemented algorithm refer to:
 115
 116        [1] M. Bronstein, Symbolic Integration I: Transcendental
 117            Functions, Second Edition, Springer-Verlang, 2005
 118
 119     """
 120     if not isinstance(f, Poly):
 121         f = Poly(f, *symbols)
 122     elif symbols:
 123         raise SymbolsError("Redundant symbols were given")
 124
 125     f, g = f.unify_with(g)
 126
 127     if f.is_multivariate:
 128         raise MultivariatePolyError(f)
 129
 130     symbols, flags = f.symbols, f.flags
 131
 132     q, r = Poly((), *symbols, **flags), f
 133     coeff, N = g.LC, f.degree - g.degree + 1
 134
 135     while not r.is_zero:
 136         M = r.degree - g.degree
 137
 138         if M < 0:
 139             break
 140         else:
 141             T, N = (r.LC, (M,)), N - 1
 142
 143             q = q.mul_term(coeff).add_term(*T)
 144             r = r.mul_term(coeff)-g.mul_term(*T)
 145
 146     return (q.mul_term(coeff**N), r.mul_term(coeff**N))
 147
 148 def poly_groebner(f, *symbols, **flags):
 149     """Computes reduced Groebner basis for a set of polynomials.
 150
 151        Given a set of multivariate polynomials F, find another set G,
 152        such that Ideal F = Ideal G and G is a reduced Groebner basis.
 153
 154        The resulting basis is unique and has monic generators.
 155
 156        Groebner bases can be used to choose specific generators for a
 157        polynomial ideal. Because these bases are unique you can check
 158        for ideal equality by comparing the Groebner bases.  To see if
 159        one polynomial lies in an ideal, divide by the elements in the
 160        base and see if the remainder vanishes.
 161
 162        They can also be used to  solve systems of polynomial equations
 163        as,  by choosing lexicographic ordering,  you can eliminate one
 164        variable at a time, provided that the ideal is zero-dimensional
 165        (finite number of solutions).
 166
 167        >>> from sympy import *
 168        >>> x,y = symbols('xy')
 169
 170        >>> G = poly_groebner([x**2 + y**3, y**2-x], x, y, order='lex')
 171
 172        >>> [ g.as_basic() for g in G ]
 173        [x - y**2, y**3 + y**4]
 174
 175        For more information on the implemented algorithm refer to:
 176
 177        [1] N.K. Bose, B. Buchberger, J.P. Guiver, Multidimensional
 178            Systems Theory and Applications, Springer, 2003, pp. 98+
 179
 180        [2] A. Giovini, T. Mora, "One sugar cube, please" or Selection
 181            strategies in Buchberger algorithm, Proc. ISSAC '91, ACM
 182
 183        [3] I.A. Ajwa, Z. Liu, P.S. Wang, Groebner Bases Algorithm,
 184            http://citeseer.ist.psu.edu/ajwa95grbner.html, 1995
 185
 186        [4] D. Cox, J. Little, D. O'Shea, Ideals, Varieties and
 187            Algorithms, Springer, Second Edition, 1997, pp. 62
 188
 189     """
 190     if isinstance(f, (tuple, list, set)):
 191         f, g = f[0], list(f[1:])
 192
 193         if not isinstance(f, Poly):
 194             f = Poly(f, *symbols, **flags)
 195         elif symbols or flags:
 196             raise SymbolsError("Redundant symbols or flags were given")
 197
 198         f, g = f.unify_with(g)
 199
 200         symbols, flags = f.symbols, f.flags
 201     else:
 202         if not isinstance(f, Poly):
 203             f = Poly(f, *symbols, **flags)
 204         elif symbols or flags:
 205             raise SymbolsError("Redundant symbols or flags were given")
 206
 207         return [f.as_monic()]
 208
 209     compare = monomial_cmp(flags.get('order'))
 210
 211     f = [ h for h in [f] + g if h ]
 212
 213     if not f:
 214         return [Poly((), *symbols, **flags)]
 215
 216     R, P, G, B, F = set(), set(), set(), {}, {}
 217
 218     for i, h in enumerate(f):
 219         F[h] = i; R.add(i)
 220
 221     def normal(g, H):
 222         h = poly_div(g, [ f[i] for i in H ])[1]
 223
 224         if h.is_zero:
 225             return None
 226         else:
 227             if not F.has_key(h):
 228                 F[h] = len(f)
 229                 f.append(h)
 230
 231             return F[h], h.LM
 232
 233     def generate(R, P, G, B):
 234         while R:
 235             h = normal(f[R.pop()], G | P)
 236
 237             if h is not None:
 238                 k, LM = h
 239
 240                 G0 = set(g for g in G if monomial_div(f[g].LM, LM))
 241                 P0 = set(p for p in P if monomial_div(f[p].LM, LM))
 242
 243                 G, P, R = G - G0, P - P0 | set([k]), R | G0 | P0
 244
 245                 for i, j in set(B):
 246                     if i in G0 or j in G0:
 247                         del B[(i, j)]
 248
 249         G |= P
 250
 251         for i in G:
 252             for j in P:
 253                 if i == j:
 254                     continue
 255
 256                 if i < j:
 257                    k = (i, j)
 258                 else:
 259                    k = (j, i)
 260
 261                 if not B.has_key(k):
 262                     B[k] = monomial_lcm(f[i].LM, f[j].LM)
 263
 264         G = set([ normal(f[g], G - set([g]))[0] for g in G ])
 265
 266         return R, P, G, B
 267
 268     R, P, G, B = generate(R, P, G, B)
 269
 270     while B:
 271         k, M = B.items()[0]
 272
 273         for l, N in B.iteritems():
 274             if compare(M, N) == 1:
 275                 k, M = l, N
 276
 277         del B[k]
 278
 279         i, j = k[0], k[1]
 280         p, q = f[i], f[j]
 281
 282         p_LM, q_LM = p.LM, q.LM
 283
 284         if M == monomial_mul(p_LM, q_LM):
 285             continue
 286
 287         criterion = False
 288
 289         for g in G:
 290             if g == i or g == j:
 291                 continue
 292
 293             if not B.has_key((min(i, g), max(i, g))):
 294                 continue
 295
 296             if not B.has_key((min(j, g), max(j, g))):
 297                 continue
 298
 299             if not monomial_div(M, f[g].LM):
 300                 continue
 301
 302             criterion = True
 303             break
 304
 305         if criterion:
 306             continue
 307
 308         p = p.mul_term(1/p.LC, monomial_div(M, p_LM))
 309         q = q.mul_term(1/q.LC, monomial_div(M, q_LM))
 310
 311         h = normal(p - q, G)
 312
 313         if h is not None:
 314             k, LM = h
 315
 316             G0 = set(g for g in G if monomial_div(f[g].LM, LM))
 317
 318             R, P, G = G0, set([k]), G - G0
 319
 320             for i, j in set(B):
 321                 if i in G0 or j in G0:
 322                     del B[(i, j)]
 323
 324             R, P, G, B = generate(R, P, G, B)
 325
 326     G = [ f[g].as_monic() for g in G ]
 327
 328     G = sorted(G, compare, lambda p: p.LM)
 329
 330     return list(reversed(G))
 331
 332 def poly_lcm(f, g, *symbols):
 333     """Computes least common multiple of two polynomials.
 334
 335        Given two univariate polynomials,  the LCM is computed  via
 336        f*g = gcd(f, g)*lcm(f, g) formula. In multivariate case, we
 337        compute the unique generator of the intersection of the two
 338        ideals, generated by f and g.  This is done by computing  a
 339        Groebner basis, with respect to any lexicographic ordering,
 340        of t*f and (1 - t)*g, where t is an unrelated symbol and
 341        filtering out solution that does not contain t.
 342
 343        For more information on the implemented algorithm refer to:
 344
 345        [1] D. Cox, J. Little, D. O'Shea, Ideals, Varieties and
 346            Algorithms, Springer, Second Edition, 1997, pp. 187
 347
 348     """
 349     if not isinstance(f, Poly):
 350         f = Poly(f, *symbols)
 351     elif symbols:
 352         raise SymbolsError("Redundant symbols were given")
 353
 354     f, g = f.unify_with(g)
 355
 356     symbols, flags = f.symbols, f.flags
 357
 358     if f.is_monomial and g.is_monomial:
 359         monom = monomial_lcm(f.LM, g.LM)
 360
 361         fc, gc = f.LC, g.LC
 362
 363         if fc.is_Rational and gc.is_Rational:
 364             coeff = Integer(ilcm(fc.p, gc.p))
 365         else:
 366             coeff = S.One
 367
 368         return Poly((coeff, monom), *symbols, **flags)
 369
 370     fc, f = f.as_primitive()
 371     gc, g = g.as_primitive()
 372
 373     lcm = ilcm(int(fc), int(gc))
 374
 375     if f.is_multivariate:
 376         t = Symbol('t', dummy=True)
 377         lex = { 'order' : 'lex' }
 378
 379         f_monoms = [ (1,) + monom for monom in f.monoms ]
 380
 381         F = Poly((f.coeffs, f_monoms), t, *symbols, **lex)
 382
 383         g_monoms = [ (0,) + monom for monom in g.monoms ] + \
 384                    [ (1,) + monom for monom in g.monoms ]
 385
 386         g_coeffs = list(g.coeffs) + [ -coeff for coeff in g.coeffs ]
 387         G = Poly(dict(zip(g_monoms, g_coeffs)), t, *symbols, **lex)
 388
 389         def independent(h):
 390             return all(not monom[0] for monom in h.monoms)
 391
 392         H = [ h for h in poly_groebner((F, G)) if independent(h) ]
 393
 394         if lcm != 1:
 395             h_coeffs = [ coeff*lcm for coeff in H[0].coeffs ]
 396         else:
 397             h_coeffs = H[0].coeffs
 398
 399         h_monoms = [ monom[1:] for monom in H[0].monoms ]
 400
 401         return Poly(dict(zip(h_monoms, h_coeffs)), *symbols, **flags)
 402     else:
 403         h = poly_div(f * g, poly_gcd(f, g))[0]
 404
 405         if lcm != 1:
 406             return h.mul_term(lcm / h.LC)
 407         else:
 408             return h.as_monic()
 409
 410 def poly_gcd(f, g, *symbols):
 411     """Compute greatest common divisor of two polynomials.
 412
 413        Given two univariate polynomials, subresultants are used
 414        to compute the GCD.  In multivariate case Groebner basis
 415        approach is used together with f*g = gcd(f, g)*lcm(f, g)
 416        well known formula.
 417
 418        For more information on the implemented algorithm refer to:
 419
 420        [1] D. Cox, J. Little, D. O'Shea, Ideals, Varieties and
 421            Algorithms, Springer, Second Edition, 1997, pp. 187
 422
 423     """
 424     if not isinstance(f, Poly):
 425         f = Poly(f, *symbols)
 426     elif symbols:
 427         raise SymbolsError("Redundant symbols were given")
 428
 429     f, g = f.unify_with(g)
 430
 431     symbols, flags = f.symbols, f.flags
 432
 433     if f.is_zero and g.is_zero:
 434         return f
 435
 436     if f.is_constant:
 437         if f.is_zero:
 438             cont, g = g.as_primitive()
 439             return g.mul_term(cont / g.LC)
 440         if f.is_one:
 441             return f
 442
 443     if g.is_constant:
 444         if g.is_zero:
 445             cont, f = f.as_primitive()
 446             return f.mul_term(cont / f.LC)
 447         if g.is_one:
 448             return g
 449
 450     if f.is_monomial and g.is_monomial:
 451         monom = monomial_gcd(f.LM, g.LM)
 452
 453         fc, gc = f.LC, g.LC
 454
 455         if fc.is_Rational and gc.is_Rational:
 456             coeff = Integer(igcd(fc.p, gc.p))
 457         else:
 458             coeff = S.One
 459
 460         return Poly((coeff, monom), *symbols, **flags)
 461
 462     cf, f = f.as_primitive()
 463     cg, g = g.as_primitive()
 464
 465     gcd = igcd(int(cf), int(cg))
 466
 467     if f.is_multivariate:
 468         h = poly_div(f*g, poly_lcm(f, g))[0]
 469     else:
 470         h = poly_subresultants(f, g)[-1]
 471
 472     if gcd != 1:
 473         return h.mul_term(gcd / h.LC)
 474     else:
 475         return h.as_monic()
 476
 477 def poly_gcdex(f, g, *symbols):
 478     """Extended Euclidean algorithm.
 479
 480        Given univariate polynomials f and g over an Euclidean domain,
 481        computes polynomials s, t and h,  such that h = gcd(f, g) and
 482        s*f + t*g = h.
 483
 484        For more information on the implemented algorithm refer to:
 485
 486        [1] M. Bronstein, Symbolic Integration I: Transcendental
 487            Functions, Second Edition, Springer-Verlang, 2005
 488
 489     """
 490     s, h = poly_half_gcdex(f, g, *symbols)
 491     return s, poly_div(h - s*f, g)[0], h
 492
 493 def poly_half_gcdex(f, g, *symbols):
 494     """Half extended Euclidean algorithm.
 495
 496        Efficiently computes gcd(f, g)  and one of the coefficients
 497        in extended Euclidean algorithm. Formally, given univariate
 498        polynomials f and g over an Euclidean domain, computes s
 499        and h, such that h = gcd(f, g) and s*f = h (mod g).
 500
 501        For more information on the implemented algorithm refer to:
 502
 503        [1] M. Bronstein, Symbolic Integration I: Transcendental
 504            Functions, Second Edition, Springer-Verlang, 2005
 505
 506     """
 507     if not isinstance(f, Poly):
 508         f = Poly(f, *symbols)
 509     elif symbols:
 510         raise SymbolsError("Redundant symbols were given")
 511
 512     f, g = f.unify_with(g)
 513
 514     if f.is_multivariate:
 515         raise MultivariatePolyError(f)
 516
 517     symbols, flags = f.symbols, f.flags
 518
 519     a = Poly(S.One, *symbols, **flags)
 520     b = Poly((), *symbols, **flags)
 521
 522     while not g.is_zero:
 523         q, r = poly_div(f, g)
 524
 525         f, g = g, r
 526         c = a - q*b
 527         a, b = b, c
 528
 529     return a, f
 530
 531 def poly_resultant(f, g, *symbols):
 532     """Computes resultant of two univariate polynomials.
 533
 534        Resultants are a classical algebraic tool for determining if
 535        a  system of n polynomials in n-1 variables have common root
 536        without explicitly solving for the roots.
 537
 538        They are efficiently represented as  determinants of Bezout
 539        matrices whose entries are computed using O(n**2) additions
 540        and multiplications where n = max(deg(f), deg(g)).
 541
 542        >>> from sympy import *
 543        >>> x,y = symbols('xy')
 544
 545        Polynomials x**2-1 and (x-1)**2 have common root:
 546
 547        >>> poly_resultant(x**2-1, (x-1)**2, x)
 548        0
 549
 550        For more information on the implemented algorithm refer to:
 551
 552        [1] Eng-Wee Chionh, Fast Computation of the Bezout and Dixon
 553            Resultant Matrices, Journal of Symbolic Computation, ACM,
 554            Volume 33, Issue 1, January 2002, Pages 13-29
 555
 556     """
 557     if not isinstance(f, Poly):
 558         f = Poly(f, *symbols)
 559     elif symbols:
 560         raise SymbolsError("Redundant symbols were given")
 561
 562     f, g = f.unify_with(g)
 563
 564     if f.is_multivariate:
 565         raise MultivariatePolyError(f)
 566
 567     n, m = f.degree, g.degree
 568
 569     N = max(n, m)
 570
 571     if n < m:
 572         p = f.as_uv_dict()
 573         q = g.as_uv_dict()
 574     else:
 575         q = f.as_uv_dict()
 576         p = g.as_uv_dict()
 577
 578     import sympy.matrices
 579
 580     B = sympy.matrices.zeros(N)
 581
 582     for i in xrange(N):
 583         for j in xrange(i, N):
 584             if p.has_key(i) and q.has_key(j+1):
 585                 B[i, j] += p[i] * q[j+1]
 586
 587             if p.has_key(j+1) and q.has_key(i):
 588                 B[i, j] -= p[j+1] * q[i]
 589
 590     for i in xrange(1, N-1):
 591         for j in xrange(i, N-1):
 592             B[i, j] += B[i-1, j+1]
 593
 594     for i in xrange(N):
 595         for j in xrange(i+1, N):
 596             B[j, i] = B[i, j]
 597
 598     det = B.det()
 599
 600     if not det:
 601         return det
 602     else:
 603         if n >= m:
 604             det /= f.LC**(n-m)
 605         else:
 606             det /= g.LC**(m-n)
 607
 608         sign = (-1)**(n*(n-1)/2)
 609
 610         if det.is_Atom:
 611             return sign * det
 612         else:
 613             return sign * Poly.cancel(det)
 614
 615 def poly_subresultants(f, g, *symbols):
 616     """Computes subresultant PRS of two univariate polynomials.
 617
 618        Polynomial remainder sequence (PRS) is a fundamental tool in
 619        computer algebra as it gives as a sub-product the polynomial
 620        greatest common divisor (GCD), provided that the coefficient
 621        domain is an unique factorization domain.
 622
 623        There are several methods for computing PRS, eg.: Euclidean
 624        PRS, where the most famous algorithm is used, primitive PRS
 625        and, finally, subresultants which are implemented here.
 626
 627        The Euclidean approach is reasonably efficient but suffers
 628        severely from coefficient growth.  The primitive algorithm
 629        avoids this but requires a lot of coefficient computations.
 630
 631        Subresultants solve both problems and so it is efficient and
 632        have moderate coefficient growth. The current implementation
 633        uses pseudo-divisions  which is well suited for coefficients
 634        in integral domains or number fields.
 635
 636        Formally,  given univariate polynomials f and g over an UFD,
 637        then a sequence (R_0, R_1, ..., R_k, 0, ...) is a polynomial
 638        remainder sequence where R_0 = f, R_1 = g, R_k != 0 and R_k
 639        is similar to gcd(f, g).
 640
 641        For more information on the implemented algorithm refer to:
 642
 643        [1] M. Bronstein, Symbolic Integration I: Transcendental
 644            Functions, Second Edition, Springer-Verlang, 2005
 645
 646        [2] M. Keber, Division-Free computation of subresultants
 647            using Bezout matrices, Tech. Report MPI-I-2006-1-006,
 648            Saarbrucken, 2006
 649
 650     """
 651     if not isinstance(f, Poly):
 652         f = Poly(f, *symbols)
 653     elif symbols:
 654         raise SymbolsError("Redundant symbols were given")
 655
 656     f, g = f.unify_with(g)
 657
 658     if f.is_multivariate:
 659         raise MultivariatePolyError(f)
 660
 661     symbols, flags = f.symbols, f.flags
 662
 663     n, m = f.degree, g.degree
 664
 665     if n < m:
 666         f, g = g, f
 667         n, m = m, n
 668
 669     prs = [f, g]
 670
 671     d = n - m
 672
 673     b = (-1)**(d + 1)
 674
 675     h = poly_pdiv(f, g)[1]
 676     h = h.mul_term(b)
 677
 678     k = h.degree
 679
 680     c = S.NegativeOne
 681
 682     while not h.is_zero:
 683         prs.append(h)
 684
 685         coeff = g.LC
 686
 687         c = (-coeff)**d / c**(d-1)
 688
 689         b = -coeff * c**(m-k)
 690
 691         f, g, m, d = g, h, k, m-k
 692
 693         h = poly_pdiv(f, g)[1]
 694         h = h.div_term(b)
 695
 696         k = h.degree
 697
 698     return prs
 699
 700 def poly_sqf(f, *symbols):
 701     """Compute square-free decomposition of an univariate polynomial.
 702
 703        Given an univariate polynomial f over an unique factorization domain
 704        returns tuple (f_1, f_2, ..., f_n),  where all  A_i are co-prime and
 705        square-free polynomials and f = f_1 * f_2**2 * ... * f_n**n.
 706
 707        >>> from sympy import *
 708        >>> x,y = symbols('xy')
 709
 710        >>> p, q = poly_sqf(x*(x+1)**2, x)
 711
 712        >>> p.as_basic()
 713        x
 714        >>> q.as_basic()
 715        1 + x
 716
 717        For more information on the implemented algorithm refer to:
 718
 719        [1] M. Bronstein, Symbolic Integration I: Transcendental
 720            Functions, Second Edition, Springer-Verlang, 2005
 721
 722        [2] J. von zur Gathen, J. Gerhard, Modern Computer Algebra,
 723            Second Edition, Cambridge University Press, 2003
 724
 725     """
 726     if not isinstance(f, Poly):
 727         f = Poly(f, *symbols)
 728     elif symbols:
 729         raise SymbolsError("Redundant symbols were given")
 730
 731     if f.is_multivariate:
 732         raise MultivariatePolyError(f)
 733
 734     coeff, f = f.as_primitive()
 735
 736     sqf = []
 737
 738     h = f.diff()
 739
 740     g = poly_gcd(f, h)
 741
 742     p = poly_div(f, g)[0]
 743     q = poly_div(h, g)[0]
 744
 745     while True:
 746         h = q - p.diff()
 747
 748         if h.is_zero:
 749             break
 750
 751         g = poly_gcd(p, h)
 752
 753         sqf.append(g)
 754
 755         p = poly_div(p, g)[0]
 756         q = poly_div(h, g)[0]
 757
 758     sqf.append(p)
 759
 760     head, tail = sqf[0], sqf[1:]
 761     head = head.mul_term(coeff)
 762
 763     return [head] + tail
 764
 765 def poly_decompose(f, *symbols):
 766     """Computes functional decomposition of an univariate polynomial.
 767
 768        Besides factorization and square-free decomposition, functional
 769        decomposition is another important, but very different,  way of
 770        breaking down polynomials into simpler parts.
 771
 772        Formally given an univariate polynomial f with coefficients in a
 773        field of characteristic zero, returns tuple (f_1, f_2, ..., f_n)
 774        where f = f_1 o f_2 o ... f_n = f_1(f_2(... f_n)) and f_2, ...,
 775        f_n are monic and homogeneous polynomials of degree at least 2.
 776
 777        Unlike factorization, complete functional decompositions of
 778        polynomials are not unique, consider examples:
 779
 780         [1] f o g = f(x + b) o (g - b)
 781         [2] x**n o x**m = x**m o x**n
 782         [3] T_n o T_m = T_m o T_n
 783
 784        where T_n and T_m are Chebyshev polynomials.
 785
 786        >>> from sympy import *
 787        >>> x,y = symbols('xy')
 788
 789        >>> p, q = poly_decompose(x**4+2*x**2 + y, x)
 790
 791        >>> p.as_basic()
 792        y + 2*x + x**2
 793        >>> q.as_basic()
 794        x**2
 795
 796        For more information on the implemented algorithm refer to:
 797
 798        [1] D. Kozen, S. Landau, Polynomial decomposition algorithms,
 799            Journal of Symbolic Computation 7 (1989), pp. 445-456
 800
 801     """
 802     if not isinstance(f, Poly):
 803         f = Poly(f, *symbols)
 804     elif symbols:
 805         raise SymbolsError("Redundant symbols were given")
 806
 807     if f.is_multivariate:
 808         raise MultivariatePolyError(f)
 809
 810     symbols = f.symbols
 811     flags = f.flags
 812
 813     def right_factor(f, s):
 814         n, lc = f.degree, f.LC
 815
 816         f = f.as_uv_dict()
 817         q = { s : S.One }
 818
 819         r = n // s
 820
 821         for k in xrange(1, s):
 822             coeff = S.Zero
 823
 824             for j in xrange(0, k):
 825                 if not f.has_key(n+j-k):
 826                     continue
 827
 828                 if not q.has_key(s-j):
 829                     continue
 830
 831                 fc, qc = f[n+j-k], q[s-j]
 832
 833                 coeff += (k - r*j)*fc*qc
 834
 835             if coeff is not S.Zero:
 836                 q[s-k] = coeff / (k*r*lc)
 837
 838         return Poly(q, *symbols, **flags)
 839
 840     def left_factor(f, h):
 841         g, i = {}, 0
 842
 843         while not f.is_zero:
 844             q, r = poly_div(f, h)
 845
 846             if not r.is_constant:
 847                 return None
 848             else:
 849                 if r.LC is not S.Zero:
 850                     g[i] = r.LC
 851
 852                 f, i = q, i + 1
 853
 854         return Poly(g, *symbols, **flags)
 855
 856     def decompose(f):
 857         deg = f.degree
 858
 859         for s in xrange(2, deg):
 860             if deg % s != 0:
 861                 continue
 862
 863             h = right_factor(f, s)
 864
 865             if h is not None:
 866                 g = left_factor(f, h)
 867
 868                 if g is not None:
 869                     return (g, h)
 870
 871         return None
 872
 873     F = []
 874
 875     while True:
 876         result = decompose(f)
 877
 878         if result is not None:
 879             f, h = result
 880             F = [h] + F
 881         else:
 882             break
 883
 884     return [f] + F