Combinatorics and Computer Algebra – Polynomial arithmetic

Combinatorics and Computer Algebra – Polynomial arithmetic#

Viviane Pons (based on the content of Florent Hivert)

references

  • J. von zur Gathen, J. Gerhard, Modern Computer Algebra, 3rd ed., Cambridge University Press, 2013.

  • A. Casamayou, N. Cohen, G. Connan, T. Dumont, L. Fousse, F. Maltey, M. Meulien, M. Mezzarobba, C. Pernet, N. M. Thiéry, P. Zimmermann Calcul mathématique avec Sage, 468 pages, 2013, ISBN: 978-14-811-9104-3. Link

  • A. Bostan, F. Chyzak, M. Giusti, R. Lebreton, G. Lecerf, B. Salvy and É. Schost, Algorithmes Efficaces en Calcul Formel, 2017, ISBN:979-10-699-0947-2. Link

Symbolic computation: data structure#

Question: What data structures can we use for mathematical calculations?

Examples of symbolic computations with SageMath

How do we compute:

\(\frac{\partial}{\partial a}\left(1 + F(a, 2ab) + b\right)= 2b\, \mathrm{D}_{1}\left(F\right)\left(a, 2 \, a b\right) + \mathrm{D}_{0}\left(F\right)\left(a, 2 \, a b\right)\)

(we use var function diff)

# Write your code here

# Write your code here

How do we compute:

\[\sum_{i=0}^{i=n} i^2 = \frac{n(n+1)(2n+1)}{6}\]
\[\sum_{n=1}^{\infty} 2^{-n} = 1\]
\[\sum_{i=1}^{\infty} \frac{1}{i^2} = \frac{\pi^2}{6}\]

We can use sum (which overrides the python sum function). Infinite is writen +oo

sum?

# Write your code here

# Write your code here

# Write your code here

How are these represented with SageMath?

var('a','b')
function('F')
e = 1 + F(a,2*a*b) + b
e

e.operands()

e.operator()

e2 = e.operands()[1]
e2

e2.operands()

e2.operator()

e3 = e2.operands()[1]
e3

e3.operands()

e3.operator()

The expression is seen as a syntax tree

Syntax trees are also used in compilers (look at syntax analysis and abstract syntax trees)

Advantages

  • very versatile data structure, easily extensible

  • programming with pattern matching

  • easy to implement in a functional language

  • very well suited to analytical computation (differentiation, integration…)

Disadvantages

  • very unstructured data

  • recursive data structure

  • a lot of dynamic allocation, bad data locality

  • little code reuse, few generic algorithms

We run into semantic problems

Example:

does the order of the operands matter?

.

.

.

.

.

Answer: it depends of the operator AND of the type of variables (examples: product of matrices). How can we deal with untyped variables?

Example: factorization of polynomials

x = var('x')
p = 54*x^4+36*x^3-102*x^2-72*x-12
p

factor(p)

Are we happy with this answer?

It is indeed a factorization of \(p\), but is it optimal?

.

.

.

.

.

It depends on context! Now Sage consideds \(p\) as a symobolic expression, it does not know if we want to factorize \(p\) as a polynomial with integer coefficient or with rational coefficient or something else.

We can change that, and tell Sage in which context we consider \(p\).

Let’s first look at \(p\) as polynomial with integer coefficient.

R = ZZ['x']; R

q = R(p); q

\(p\) and \(q\) look the same, but they are not the same as objects

p.parent()

q.parent()

Now the result of the factorization has no ambiguity

factor(q)

Let’s now change the context and work on rational numbers

R = QQ[x]
R

q = R(q)
q

factor(q)

In this new context, the factorization is again non ambiguous but different than the precedent one.

Note that Sage knows that \(R\) is an euclidian ring

R.category()

which means that the factorization is unique.

Let us now factorize using complex numbers using a 16 bits numerical approximation.

R = ComplexField(16)['x']; R

q = R(p); q

factor(q)

Another option is to enlarge the field of rational numbers by adding \(\sqrt{2}\)

R = QQ[sqrt(2)]['x']; R

q = R(p); q

factor(q)

And what if we want to consider coefficients taken modulo 5?

R = GF(5)['x']; R

q = R(p); q

factor(q)

Are all these factorization « equal »? How can we test if an expression is equal to another?

.

.

.

.

.

Theorem (Richardson-Matiyasevich)

In the class of expressions formed from a variable \(X\) and the constant \(1\) using the ring operations \(+\), \(-\), \(\times\), composition and the \(\sin(\cdot)\) and absolute value \(|\cdot|\) functions, the test for equivalence to \(0\) is undecidable.

Solution To control the computational semantics and equality problems, we limit ourselves to a specific domain of expressions.

Especially, we are often looking for a normal form.

In this lecture, we study the specific domain of: polynomials.

  • simple, but already very general

  • similar to arbitrary-precision integers.

Multi-precision integers#

Reminder Machine integers have limited precision.

C++ examples taken from Florent’s Hivert machine:

type

bits

min

max

int

32

-2 147 483 648

2 147 483 647

unsigned int

32

0

4 294 967 295

long

64

-9 223 372 036 854 775 808

9 223 372 036 854 775 807

unsigned long

64

0

18 446 744 073 709 551 615

In python: arbitrary long integers

factorial(100)

a = factorial(100)
len(a.bits())

Question How can we represent numbers of arbitrary sizes?

Solution

We choose a base \(B\) that fits in a machine word, and we express the number in this base: $\( n = (a_k,\dots a_2a_1a_0)_B = \sum_{i=0}^{k} a_i B^i\,. \)$

Example: GMP (for GNU Multi Precision)

\(B = 2^{64}\) on the x86_64 architecture.

A number written in base \(B\) is a polynomial: $\( n = P(B) \text{ where } P= \sum_{i=0}^{k} a_i X^i\,. \)$

Roughly speaking: integer arithmetic = polynomial arithmetic + carry

Polynomials#

How to represent polynomials?#

.

.

.

.

.

.

Basically 2 solutions:

  • dense representation: we store all the coefficients up to the degree in an array

  • sparse representation: we only store non-zero coefficients with their degrees (association table, ordered list)

TODO: implement a basic class of dense polynomials in python

What are the basic operations? What complexity?#

(on dense polynomials)

.

.

.

.

.

  • addition

  • evaluation

  • product

Complexity of addition: \(O(d)\) for a polynomial of degree \(d\)

Complexity of evaluation? Of product?

Evaluation of a dense Polynomial#

\(P(X) = \sum_{i=0}^d a_i X^i\)

How many products and additions to compute P(a)?

Example on degree \(3\):

\(P(X) = a_3 X^3 + a_2 X^2 + a_1 X + a_0\)

Try to find to most efficient way to compute \(P(a)\).

.

.

.

.

.

Solution Horner’s method

\(P(x)=(\cdots ((a_dx+a_{d-1})x+a_{n-2})x+\cdots +a_1)x+a_0\)

Example in degree \(3\):

\(P(X) = ((a_3 X + a_2)X + a_1)X + a_0\)

Result: \(d\) products and \(d\) additions, complexity \(O(d)\)

This is optimal if we need to calculate the value at \(1\) point, but not if we want to calculate values at more points.

  • integers \(0\dots n\): Newton potential operator \(\Delta\) (see TD)

  • For roots of unity: fast Fourier transform

  • For a large number of points: Euclidean division + divide and conquer. (see next class)

TODO: implement the Horner’s method to evaluate polynomials on your python class

Product: what complexity?#

Naive multiplication of \(A(X)\) of degree \(d\) and \(B(X)\) of degree \(d'\).

\[A(X) = \sum_{i=0}^{d} a_i X^i\]
\[B(X) = \sum_{i=0}^{d'} b_i X^i\]
\[AB(X) = \sum_{k=0}^{d+d'} \left( \sum_{i+j = k} a_i b_j \right) X^k\]

Complexity: \(O(d^2)\)

We can do better:

  • Karatsuba \(O(n^{\log_2(3)})=O(n^{1.58})\) (see next class)

  • Toom-Cook \(O(n^{\log_3(5)})=O(n^{1.46})\)

  • Schönhage-Strassen (Fourier transform) \(O(n\log(n))\)

The Schönhage-Strassen algorithm uses the following principle:

  1. Evaluate \(A(a)\) and \(B(a)\) for many values of \(a\)

  2. Compute \(AB(a)\) for all these values

  3. interpolate to find \(P=AB\) from the values

In the case of the Schönhage-Strassen, the values \(a\) are the roots of unity.

Interpolation of polynomials#

With \(d+1\) values \(P(a)\), we can recover the coefficient of the polynomial \(P\) (see TD).