--- jupytext: text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.16.0 kernelspec: display_name: SageMath 10.2 language: sage name: sagemath --- # Combinatorics and Computer Algebra -- Polynomial arithmetic # Polynomial multiplication and division Reminder: basic operations * addition $O(d)$ * evaluation $O(d)$ * naive multiplication $O(d^2)$ ## Naive multiplication To multiply two integers or polynomials, we proceed as we did in primary school, multiplying the first number by the digits of the second number and shifting, then adding everything together: 36 * 25 ---- 36*5 : 180 36*2 : +72 ==== 900 3x + 6 * 2x + 5 --------- (3x+6)*5 : 15x + 30 (3x+6)*2x : 6x^2 + 12x =============== 6x^2 + 27x + 30 **TODO: implement naive multiplication on polynomials** **Complexity** To multiply two numbers with $d_1$ and $d_2$ digits, we need to multiply digits $d_1d_2$ times. To multiply two polynomials of degree $d_1$ and $d_2$, we need to multiply coefficients $(d_1+1)(d_2+1)$ times. (reminder: a polynomial of degree $d$ has $d+1$ coefficients) **Can we do better than quadratic?** Yes! * Karatsuba $O(n^{\log_2(3)})=O(n^{1.58})$ * Toom-Cook $O(n^{\log_3(5)})=O(n^{1.46})$ * Schönhage-Strassen (Fourier transform) $O(n\log(n))$ ## Karatsuba algorithm ### The algorithm Idea: $$(a X^k + b)(c X^k + d) = ac X^{2k} + (ad + bc) X^k + bd$$ Requires to compute 4 products: $ac$, $ad$, $bc$, and $bd$ We can *rewrite* it as $$(a X^k + b)(c X^k + d) = ac X^{2k} + (ac + bd - (a - b)(c - d)) X^k + bd$$ **only 3 products** How can we use this principle to get a better muliplication algorithm? . . . . . . **Divide and conquer** Computation of $PQ$ 1. We divide $P$ and $Q$ into two $P = A X^k+B$ and $Q = C X^k +D$ 2. We recursively compute $AC$, $BD$ and $(A - B)(C - D)$ Works on both polynomials and integers **Example on integers** Product $1237 \times 2587$ $$1237 = 12 \times 10^2 + 37$$ $$2587 = 25 \times 10^2 + 87$$ | $A_0$ | $B_0$ | $C_0$ | $D_0$ | |--------|--------|--------|--------| | $12$ | $37$ | $25$ | $87$ | Product of $A_0 C_0 = 12 \times 25$ $$12 = 1 \times 10 + 2$$ $$25 = 2 \times 10 + 5$$ | $A_1$ | $B_1$ | $C_1$ | $D_1$ | |--------|--------|--------|--------| | $1$ | $2$ | $2$ | $5$ | $$A_1 C_1 = 1 \times 2 = 2 $$ $$B_1 D_1 = 2 \times 5 = 10 $$ $$(A_1 - B_1)(C_1 - D_1) = (-1) \times (-3) = 3$$ We find $$12 \times 25 = A_1 C_1 10^2 + (A_1 C_1 + B_1 D_1 - (A_1 - B_1)(C_1 - D_1)) 10 + B_1 D_1 = 2 \times 10^2 + 9 \times 10 + 10 = 300$$ We use a similar process to compute $$B_0 D_0 = 37 \times 87 = 3219$$ $$(A_0 - B_0)(C_0 - D_0) = (-25) \times (-62) = 1550$$ We obtain the final result | | | | | | | | | |----------------------------------------|---|---|---|---|---|---|---| | $A_0 C_0 \times 10^4$ | 3 | 0 | 0 | 0 | 0 | 0 | 0 | | $+ A_0 C_0 \times 10^2$ | | | 3 | 0 | 0 | 0 | 0 | | $+ B_0 D_0 \times 10^2$ | | 3 | 2 | 1 | 9 | 0 | 0 | | $- (A_0 - B_0)(C_0 - D_0) \times 10^2$ | | 1 | 5 | 5 | 0 | 0 | 0 | | $+ B_0 D_0$ | | | | 3 | 2 | 1 | 9 | | $=$ | 3 | 2 | 0 | 0 | 1 | 1 | 9 | ```{code-cell} ipython3 1237 * 2587 ``` **Example on polynomials:** $$ \begin{aligned} P(X) &= 3 X^3 - X^2 +1 = (3X - 1) \times X^2 + 1 \\ Q(X) &= 2 X^2 + X + 5 \\ \\ A(X) &= 3X+1 \\ B(X) &= 1 \\ C(X) &= 2 \\ D(X) &= X + 5 \\ \\ AC(X) &= 2(3X +1) = 6X -2 \\ BD(X) &= X + 5 \\ (A - B)(C - D)(X) &= (3X - 2)(-X -3) \\ &= -3X^2 + (-3 + 6 -5 \times 2) X + 6 \\ &= -3 X^2 -7X +6 \\ \\ PQ(X) &= (6X -2) X^4 + (6X - 2 + X + 5 + 3X^2 + 7X -6)X^2 + X + 5 \\ &= 6X^5 + X^4 + 14 X^3 - 3 X^2 + X + 5 \end{aligned} $$ ```{code-cell} ipython3 R = ZZ['x'] x = R.gen() P = 3*x^3 - x^2 +1 Q = 2*x^2 + x + 5 P*Q ``` +++ {"user_expressions": []} **TODO: implement Karatsuba algorithm for polynomial multiplication and compare efficiency with naive multiplication** +++ ### Complexity What is the complexity of the Karatsuba algorithm? . . . . . For simplicity, suppose that $P$ and $Q$ have degree $d = 2^k$. To make **one** product of degree $d$, we need to recursively perform **three** products of degree $d/2$. That gives | $2^k$ | $2^{k-1}$ | $2^{k-2}$ | $\dots$ | $1$ | |-------|-----------|-----------|---------|-------| | $1$ | $3$ | $3^2$ | $\dots$ | $3^k$ | as $k = \log_2(d)$, we get $$3^{\log_2(d)} = \exp\left(\log_2(d)\ln(3)\right) = \exp\left(\frac{\ln(d)}{\ln(2)}\ln(3)\right) = d^{\log_2(3)} \approx d^{1.58}$$ Note: this is an example of the [Master Theorem](https://en.wikipedia.org/wiki/Master_theorem_(analysis_of_algorithms)) to compute complexity of divide and conqueer algorithms. **Question: is $O(n^{\log_2(3)})$ much better than $O(n^2)$ ?** Illustration with an image ![Graphical illustration of the number of operation needed: the classical case is 4 operations drawn in a square. 1 iteration of the algo removes one square, 2 iterations removes a smaller square in each of the remaining squares, etc: after many iterations, we see a sierpinski triangle and most of the squares are empty](kara.png) +++ ## Polynomial Euclidian division Let us look at another useful operation on polynomials. We start with a definition ### Definition: integral and Euclidian ring **Definition** Let $(A,0,+,1,\times)$ be a unitary commutative ring. We say that $A$ is *integral* if it has no divisor of $0$ for any $a, b\in A$. In other words $ab = 0$ implies $a=0$ or $b=0$ An integral ring is said to be *Euclidean* if we can find an application $v:A/\{0\}\mapsto \mathbb{N}$ verifying that * for any $a,b\in A$ where $b\neq 0$, there exists $q,r\in A$ such that $a=bq+r$ and $r=0$ or $v(r) < v(b)$. * for any $a,b\in A/\{0\}$, $v(b)\leq v(ab)$. Note 1: it is convenient to suppose that $v(0)=-1$ or $v(0)=-\infty$. Note 2: the pair $q,r$ is not necessarily unique. **Example of Euclidian ring** * The integers $\mathbb{Z}$ with $v(n)=|n|$. * If $K$ is a field, the ring $K[X]$ of univariate polynomials with coefficients in $K$ with $v(P) = d^o(P)$. * The Gaussian integers (complex numbers $a+ib$ where $a,b\in\mathbb{Z}$) with $v(a+ib)=a^2+b^2$. In the first two cases, the pair $q, r$ is unique. ### Example of Euclidian division on polynomials $A(X) = X^5 - X^4 - X^3 +3 X^2 - 2X$ $B(X) = X^2 - X + 1$ | $X^5$ | $-X^4$ | $-X^3$ | $+3X^2$ | $-2X$ | $+0$ | $\vert$ | $X^2$ | $-X$ | $+1$ | |--------|--------|---------|---------|-------|------|---------|-------|-------|------| | $-X^5$ | $+X^4$ | $-X^3$ | | | | $\vert$ | $X^3$ | | | | | | $-2X^3$ | $+3X^2$ | $-2X$ | $+0$ | $\vert$ | | | | | | | $+2X^3$ | $-2X^2$ | $+2X$ | | $\vert$ | | $-2X$ | | | | | | $X^2$ | $+0$ | $+0$ | $\vert$ | | | | | | | | $-X^2$ | $+X$ | $-1$ | $\vert$ | | | $+1$ | | | | | | $X$ | $-1$ | $\vert$ | | | | $Q(X) = X^3 -2X +1$ $R(X) = X - 1$ $A(X) = (X^2 - X + 1)(X^3 -2X +1) + X - 1$ ```{code-cell} ipython3 R = ZZ['X'] X = R.gen() A = X^5 - X^4 - X^3 +3 *X^2 - 2*X B = X^2 - X + 1 Q = X^3 -2*X +1 R = X - 1 ``` ```{code-cell} ipython3 B*Q + R ``` ```{code-cell} ipython3 A /B ``` ```{code-cell} ipython3 A % B ``` ```{code-cell} ipython3 A // B ``` **Applications**: GCD, decomposition (factorization) ```{code-cell} ipython3 factor(A) ``` ```{code-cell} ipython3 factor(Q) ``` ```{code-cell} ipython3 gcd(A,Q) ``` **Question: what does it mean to be invertible? What are invertible polynomials?** **Can we still "invert" a polynomial? Where does the inverse live?** **TODO: implement euclidian division** +++ ### Complexity How many product of coefficients are needed to compute the Euclidien division? Can we do better? (See TD) ```{code-cell} ipython3 ```