# 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, wrap_line=False)
>>> 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#

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