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 factor() function. By default 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 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 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 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 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 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 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 Wester`s example was \(x^{100} - 1\), but was truncated for readability purpose. Note that this is not a big struggle for 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 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 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 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 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 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 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 sympy.polys is to expand all expressions given as arguments to polynomial manipulation functions and 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 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 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 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 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 together() and 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, https://www.math.unm.edu/~wester/cas/book/Wester.pdf