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.

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, wrap_line=False, no_global=True)
>>> var('x,y,z,s,c,n')
(x, y, z, s, c, n)
```

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()`.

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⎠
```

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⎠
```

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.:

```
>>> gcd(2*x**(n + 4) - x**(n + 2), 4*x**(n + 1) + 3*x**n)
n
x
```

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:

```
>>> 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\).

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
```

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)
```

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.

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 ⎠
```

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`. Do 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
```

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)

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.

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
```

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)⎠
```

To compute a reduced Gröbner basis for a set of polynomials use
`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)]
```

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.

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
```

[Wester1999] | Michael J. Wester, A Critique of the Mathematical Abilities of CA Systems, 1999, http://www.math.unm.edu/~wester/cas/book/Wester.pdf |