.. _polys-wester: ============================== Examples from Wester's Article ============================== Introduction ============ In this tutorial we present examples from Wester's article concerning comparison and critique of mathematical abilities of several computer algebra systems (see [Wester1999]_). All the examples are related to polynomial and algebraic computations and SymPy specific remarks were added to all of them. Examples ======== All examples in this tutorial are computable, so one can just copy and paste them into a Python shell and do something useful with them. All computations were done using the following setup:: >>> from sympy import * >>> init_printing(use_unicode=True) >>> var('x,y,z,s,c') (x, y, z, s, c) Simple univariate polynomial factorization ------------------------------------------ To obtain a factorization of a polynomial use :func:~.factor function. By default :func:~.factor returns the result in unevaluated form, so the content of the input polynomial is left unexpanded, as in the following example:: >>> factor(6*x - 10) 2⋅(3⋅x - 5) To achieve the same effect in a more systematic way use :func:~.primitive function, which returns the content and the primitive part of the input polynomial:: >>> primitive(6*x - 10) (2, 3⋅x - 5) .. note:: The content and the primitive part can be computed only over a ring. To simplify coefficients of a polynomial over a field use :func:~.monic. Univariate GCD, resultant and factorization ------------------------------------------- Consider univariate polynomials f, g and h over integers:: >>> f = 64*x**34 - 21*x**47 - 126*x**8 - 46*x**5 - 16*x**60 - 81 >>> g = 72*x**60 - 25*x**25 - 19*x**23 - 22*x**39 - 83*x**52 + 54*x**10 + 81 >>> h = 34*x**19 - 25*x**16 + 70*x**7 + 20*x**3 - 91*x - 86 We can compute the greatest common divisor (GCD) of two polynomials using :func:~.gcd function:: >>> gcd(f, g) 1 We see that f and g have no common factors. However, f*h and g*h have an obvious factor h:: >>> gcd(expand(f*h), expand(g*h)) - h 0 The same can be verified using the resultant of univariate polynomials:: >>> resultant(expand(f*h), expand(g*h)) 0 Factorization of large univariate polynomials (of degree 120 in this case) over integers is also possible:: >>> factor(expand(f*g)) ⎛ 60 47 34 8 5 ⎞ ⎛ 60 52 39 25 23 10 ⎞ -⎝16⋅x + 21⋅x - 64⋅x + 126⋅x + 46⋅x + 81⎠⋅⎝72⋅x - 83⋅x - 22⋅x - 25⋅x - 19⋅x + 54⋅x + 81⎠ Multivariate GCD and factorization ---------------------------------- What can be done in univariate case, can be also done for multivariate polynomials. Consider the following polynomials f, g and h in \mathbb{Z}[x,y,z]:: >>> f = 24*x*y**19*z**8 - 47*x**17*y**5*z**8 + 6*x**15*y**9*z**2 - 3*x**22 + 5 >>> g = 34*x**5*y**8*z**13 + 20*x**7*y**7*z**7 + 12*x**9*y**16*z**4 + 80*y**14*z >>> h = 11*x**12*y**7*z**13 - 23*x**2*y**8*z**10 + 47*x**17*y**5*z**8 As previously, we can verify that f and g have no common factors:: >>> gcd(f, g) 1 However, f*h and g*h have an obvious factor h:: >>> gcd(expand(f*h), expand(g*h)) - h 0 Multivariate factorization of large polynomials is also possible:: >>> factor(expand(f*g)) 7 ⎛ 9 9 3 7 6 5 12 7⎞ ⎛ 22 17 5 8 15 9 2 19 8 ⎞ -2⋅y ⋅z⋅⎝6⋅x ⋅y ⋅z + 10⋅x ⋅z + 17⋅x ⋅y⋅z + 40⋅y ⎠⋅⎝3⋅x + 47⋅x ⋅y ⋅z - 6⋅x ⋅y ⋅z - 24⋅x⋅y ⋅z - 5⎠ Support for symbols in exponents -------------------------------- Polynomial manipulation functions provided by :mod:sympy.polys are mostly used with integer exponents. However, it's perfectly valid to compute with symbolic exponents, e.g.:: >>> n = var('n') >>> gcd(x**n - x**(2*n), x**n) n x Results may depend on powers being expanded (which will depend on assumptions of the base): >>> gcd(x**(n + 4), x**(n + 1) + 3*x**n) 1 >>> x = var('x', positive=True) >>> gcd(x**(n + 4), x**(n + 1) + 3*x**n) n x Testing if polynomials have common zeros ---------------------------------------- To test if two polynomials have a root in common we can use :func:~.resultant function. The theory says that the resultant of two polynomials vanishes if there is a common zero of those polynomials. For example:: >>> x = var('x') >>> resultant(3*x**4 + 3*x**3 + x**2 - x - 2, x**3 - 3*x**2 + x + 5) 0 We can visualize this fact by factoring the polynomials:: >>> factor(3*x**4 + 3*x**3 + x**2 - x - 2) ⎛ 3 ⎞ (x + 1)⋅⎝3⋅x + x - 2⎠ >>> factor(x**3 - 3*x**2 + x + 5) ⎛ 2 ⎞ (x + 1)⋅⎝x - 4⋅x + 5⎠ In both cases we obtained the factor x + 1 which tells us that the common root is x = -1. Normalizing simple rational functions ------------------------------------- To remove common factors from the numerator and the denominator of a rational function the elegant way, use :func:~.cancel function. For example:: >>> cancel((x**2 - 4)/(x**2 + 4*x + 4)) x - 2 ───── x + 2 Expanding expressions and factoring back ---------------------------------------- One can work easily we expressions in both expanded and factored forms. Consider a polynomial f in expanded form. We differentiate it and factor the result back:: >>> f = expand((x + 1)**20) >>> g = diff(f, x) >>> factor(g) 19 20⋅(x + 1) The same can be achieved in factored form:: >>> diff((x + 1)**20, x) 19 20⋅(x + 1) Factoring in terms of cyclotomic polynomials -------------------------------------------- SymPy can very efficiently decompose polynomials of the form x^n \pm 1 in terms of cyclotomic polynomials:: >>> factor(x**15 - 1) ⎛ 2 ⎞ ⎛ 4 3 2 ⎞ ⎛ 8 7 5 4 3 ⎞ (x - 1)⋅⎝x + x + 1⎠⋅⎝x + x + x + x + 1⎠⋅⎝x - x + x - x + x - x + 1⎠ The original Westers example was x^{100} - 1, but was truncated for readability purpose. Note that this is not a big struggle for :func:~.factor to decompose polynomials of degree 1000 or greater. Univariate factoring over Gaussian numbers ------------------------------------------ Consider a univariate polynomial f with integer coefficients:: >>> f = 4*x**4 + 8*x**3 + 77*x**2 + 18*x + 153 We want to obtain a factorization of f over Gaussian numbers. To do this we use :func:~.factor as previously, but this time we set gaussian keyword to True:: >>> factor(f, gaussian=True) ⎛ 3⋅ⅈ⎞ ⎛ 3⋅ⅈ⎞ 4⋅⎜x - ───⎟⋅⎜x + ───⎟⋅(x + 1 - 4⋅ⅈ)⋅(x + 1 + 4⋅ⅈ) ⎝ 2 ⎠ ⎝ 2 ⎠ As the result we got a splitting factorization of f with monic factors (this is a general rule when computing in a field with SymPy). The gaussian keyword is useful for improving code readability, however the same result can be computed using more general syntax:: >>> factor(f, extension=I) ⎛ 3⋅ⅈ⎞ ⎛ 3⋅ⅈ⎞ 4⋅⎜x - ───⎟⋅⎜x + ───⎟⋅(x + 1 - 4⋅ⅈ)⋅(x + 1 + 4⋅ⅈ) ⎝ 2 ⎠ ⎝ 2 ⎠ Computing with automatic field extensions ----------------------------------------- Consider two univariate polynomials f and g:: >>> f = x**3 + (sqrt(2) - 2)*x**2 - (2*sqrt(2) + 3)*x - 3*sqrt(2) >>> g = x**2 - 2 We would like to reduce degrees of the numerator and the denominator of a rational function f/g. To do this we employ :func:~.cancel function:: >>> cancel(f/g) 3 2 2 x - 2⋅x + √2⋅x - 3⋅x - 2⋅√2⋅x - 3⋅√2 ─────────────────────────────────────── 2 x - 2 Unfortunately nothing interesting happened. This is because by default SymPy treats \sqrt{2} as a generator, obtaining a bivariate polynomial for the numerator. To make :func:~.cancel recognize algebraic properties of \sqrt{2}, one needs to use extension keyword:: >>> cancel(f/g, extension=True) 2 x - 2⋅x - 3 ──────────── x - √2 Setting extension=True tells :func:~.cancel to find minimal algebraic number domain for the coefficients of f/g. The automatically inferred domain is \mathbb{Q}(\sqrt{2}). If one doesn't want to rely on automatic inference, the same result can be obtained by setting the extension keyword with an explicit algebraic number:: >>> cancel(f/g, extension=sqrt(2)) 2 x - 2⋅x - 3 ──────────── x - √2 Univariate factoring over various domains ----------------------------------------- Consider a univariate polynomial f with integer coefficients:: >>> f = x**4 - 3*x**2 + 1 With :mod:sympy.polys we can obtain factorizations of f over different domains, which includes: * rationals:: >>> factor(f) ⎛ 2 ⎞ ⎛ 2 ⎞ ⎝x - x - 1⎠⋅⎝x + x - 1⎠ * finite fields:: >>> factor(f, modulus=5) 2 2 (x - 2) ⋅(x + 2) * algebraic numbers:: >>> alg = AlgebraicNumber((sqrt(5) - 1)/2, alias='alpha') >>> factor(f, extension=alg) (x - α)⋅(x + α)⋅(x - 1 - α)⋅(x + α + 1) Factoring polynomials into linear factors ----------------------------------------- Currently SymPy can factor polynomials into irreducibles over various domains, which can result in a splitting factorization (into linear factors). However, there is currently no systematic way to infer a splitting field (algebraic number field) automatically. In future the following syntax will be implemented:: >>> factor(x**3 + x**2 - 7, split=True) Traceback (most recent call last): ... NotImplementedError: 'split' option is not implemented yet Note this is different from extension=True, because the later only tells how expression parsing should be done, not what should be the domain of computation. One can simulate the split keyword for several classes of polynomials using :func:~.solve function. Advanced factoring over finite fields ------------------------------------- Consider a univariate polynomial f with integer coefficients:: >>> f = x**11 + x + 1 We can factor f over a large finite field F_{65537}:: >>> factor(f, modulus=65537) ⎛ 2 ⎞ ⎛ 9 8 6 5 3 2 ⎞ ⎝x + x + 1⎠⋅⎝x - x + x - x + x - x + 1⎠ and expand the resulting factorization back:: >>> expand(_) 11 x + x + 1 obtaining polynomial f. This was done using symmetric polynomial representation over finite fields The same thing can be done using non-symmetric representation:: >>> factor(f, modulus=65537, symmetric=False) ⎛ 2 ⎞ ⎛ 9 8 6 5 3 2 ⎞ ⎝x + x + 1⎠⋅⎝x + 65536⋅x + x + 65536⋅x + x + 65536⋅x + 1⎠ As with symmetric representation we can expand the factorization to get the input polynomial back. This time, however, we need to truncate coefficients of the expanded polynomial modulo 65537:: >>> trunc(expand(_), 65537) 11 x + x + 1 Working with expressions as polynomials --------------------------------------- Consider a multivariate polynomial f in \mathbb{Z}[x,y,z]:: >>> f = expand((x - 2*y**2 + 3*z**3)**20) We want to compute factorization of f. To do this we use factor as usually, however we note that the polynomial in consideration is already in expanded form, so we can tell the factorization routine to skip expanding f:: >>> factor(f, expand=False) 20 ⎛ 2 3⎞ ⎝x - 2⋅y + 3⋅z ⎠ The default in :mod:sympy.polys is to expand all expressions given as arguments to polynomial manipulation functions and :class:~.Poly class. If we know that expanding is unnecessary, then by setting expand=False we can save quite a lot of time for complicated inputs. This can be really important when computing with expressions like:: >>> g = expand((sin(x) - 2*cos(y)**2 + 3*tan(z)**3)**20) >>> factor(g, expand=False) 20 ⎛ 2 3 ⎞ ⎝-sin(x) + 2⋅cos (y) - 3⋅tan (z)⎠ Computing reduced Gröbner bases ------------------------------- To compute a reduced Gröbner basis for a set of polynomials use the :func:~sympy.polys.polytools.groebner function. The function accepts various monomial orderings, e.g.: lex, grlex and grevlex, or a user defined one, via order keyword. The lex ordering is the most interesting because it has elimination property, which means that if the system of polynomial equations to :func:~sympy.polys.polytools.groebner is zero-dimensional (has finite number of solutions) the last element of the basis is a univariate polynomial. Consider the following example:: >>> f = expand((1 - c**2)**5 * (1 - s**2)**5 * (c**2 + s**2)**10) >>> groebner([f, c**2 + s**2 - 1]) ⎛⎡ 2 2 20 18 16 14 12 10⎤ ⎞ GroebnerBasis⎝⎣c + s - 1, c - 5⋅c + 10⋅c - 10⋅c + 5⋅c - c ⎦, s, c, domain=ℤ, order=lex⎠ The result is an ordinary Python list, so we can easily apply a function to all its elements, for example we can factor those elements:: >>> list(map(factor, _)) ⎡ 2 2 10 5 5⎤ ⎣c + s - 1, c ⋅(c - 1) ⋅(c + 1) ⎦ From the above we can easily find all solutions of the system of polynomial equations. Or we can use :func:~.solve to achieve this in a more systematic way:: >>> solve([f, s**2 + c**2 - 1], c, s) [(-1, 0), (0, -1), (0, 1), (1, 0)] Multivariate factoring over algebraic numbers --------------------------------------------- Computing with multivariate polynomials over various domains is as simple as in univariate case. For example consider the following factorization over \mathbb{Q}(\sqrt{-3}):: >>> factor(x**3 + y**3, extension=sqrt(-3)) ⎛ ⎛ 1 √3⋅ⅈ⎞⎞ ⎛ ⎛ 1 √3⋅ⅈ⎞⎞ (x + y)⋅⎜x + y⋅⎜- ─ - ────⎟⎟⋅⎜x + y⋅⎜- ─ + ────⎟⎟ ⎝ ⎝ 2 2 ⎠⎠ ⎝ ⎝ 2 2 ⎠⎠ .. note:: Currently multivariate polynomials over finite fields aren't supported. Partial fraction decomposition ------------------------------ Consider a univariate rational function f with integer coefficients:: >>> f = (x**2 + 2*x + 3)/(x**3 + 4*x**2 + 5*x + 2) To decompose f into partial fractions use :func:~.apart function:: >>> apart(f) 3 2 2 ───── - ───── + ──────── x + 2 x + 1 2 (x + 1) To return from partial fractions to the rational function use a composition of :func:~.together and :func:~.cancel:: >>> cancel(together(_)) 2 x + 2⋅x + 3 ─────────────────── 3 2 x + 4⋅x + 5⋅x + 2 Literature ========== .. [Wester1999] Michael J. Wester, A Critique of the Mathematical Abilities of CA Systems, 1999, `_