.. _polys-internals: =============================================== Internals of the Polynomial Manipulation Module =============================================== All polynomial manipulations are relative to a *ground domain*. For example, when factoring a polynomial like `x^{10} - 1`, one has to decide what ring the coefficients are supposed to belong to, or less trivially, what coefficients are allowed to appear in the factorization. This choice of coefficients is called a ground domain. Typical choices include the integers `\mathbb{Z}`, the rational numbers `\mathbb{Q}` or various related rings and fields. But it is perfectly legitimate (although in this case uninteresting) to factorize over polynomial rings such as `k[Y]`, where `k` is some fixed field. Thus the polynomial manipulation algorithms (both complicated ones like factoring, and simpler ones like addition or multiplication) have to rely on other code to manipulate the coefficients. In the polynomial manipulation module, such code is encapsulated in so-called :mod:`~diofant.domains`. A domain is basically a factory object: it takes various representations of data, and converts them into objects with unified interface. Every object created by a domain has to implement the arithmetic operations `+`, `-` and `\times`. Other operations are accessed through the domain, e.g. as in ``ZZ.quo(ZZ(4), ZZ(2))``. Manipulation of sparse, distributed polynomials =============================================== Dense representations quickly require infeasible amounts of storage and computation time if the number of variables increases. For this reason, there is code to manipulate polynomials in a *sparse* representation. .. currentmodule:: diofant.polys.rings Sparse polynomials are represented as dictionaries. .. autofunction:: ring .. autoclass:: PolyElement :members: :special-members: .. currentmodule:: diofant.polys.univar .. autoclass:: UnivarPolyElement :members: Polynomial factorization algorithms =================================== Many variants of Euclid's algorithm: .. module:: diofant.polys.euclidtools Classical remainder sequence ---------------------------- Let `K` be a field, and consider the ring `K[X]` of polynomials in a single indeterminate `X` with coefficients in `K`. Given two elements `f` and `g` of `K[X]` with `g\neq 0` there are unique polynomials `q` and `r` such that `f = qg + r` and `\deg(r) < \deg(g)` or `r = 0`. They are denoted by `\mathrm{quo}(f,g)` (*quotient*) and `\mathrm{rem}(f,g)` (*remainder*), so we have the *division identity* .. math:: f = \mathrm{quo}(f,g)g + \mathrm{rem}(f,g). It follows that every ideal `I` of `K[X]` is a principal ideal, generated by any element `\neq 0` of minimum degree (assuming `I` non-zero). In fact, if `g` is such a polynomial and `f` is any element of `I`, `\mathrm{rem}(f,g)` belongs to `I` as a linear combination of `f` and `g`, hence must be zero; therefore `f` is a multiple of `g`. Using this result it is possible to find a `greatest common divisor `_ (gcd) of any polynomials `f,g,\ldots` in `K[X]`. If `I` is the ideal formed by all linear combinations of the given polynomials with coefficients in `K[X]`, and `d` is its generator, then every common divisor of the polynomials also divides `d`. On the other hand, the given polynomials are multiples of the generator `d`; hence `d` is a gcd of the polynomials, denoted `\mathrm{gcd}(f,g,\ldots)`. An algorithm for the gcd of two polynomials `f` and `g` in `K[X]` can now be obtained as follows. By the division identity, `r = \mathrm{rem}(f,g)` is in the ideal generated by `f` and `g`, as well as `f` is in the ideal generated by `g` and `r`. Hence the ideals generated by the pairs `(f,g)` and `(g,r)` are the same. Set `f_0 = f`, `f_1 = g`, and define recursively `f_i = \mathrm{rem}(f_{i-2},f_{i-1})` for `i\ge 2`. The recursion ends after a finite number of steps with `f_{k+1}=0`, since the degrees of the polynomials are strictly decreasing. By the above remark, all the pairs `(f_{i-1},f_i)` generate the same ideal. In particular, the ideal generated by `f` and `g` is generated by `f_k` alone as `f_{k+1} = 0`. Hence `d = f_k` is a gcd of `f` and `g`. The sequence of polynomials `f_0`, `f_1,\ldots, f_k` is called the *Euclidean polynomial remainder sequence* determined by `(f,g)` because of the analogy with the classical `Euclidean algorithm `_ for the gcd of natural numbers. The algorithm may be extended to obtain an expression for `d` in terms of `f` and `g` by using the full division identities to write recursively each `f_i` as a linear combination of `f` and `g`. This leads to an equation .. math:: d = uf + vg\qquad (u,v \in K[X]) analogous to `Bézout's identity `_ in the case of integers. Simplified remainder sequences ------------------------------ Assume, as is usual, that the coefficient field `K` is the field of fractions of an integral domain `A`. In this case the coefficients (numerators and denominators) of the polynomials in the Euclidean remainder sequence tend to grow very fast. If `A` is a unique factorization domain, the coefficients may be reduced by cancelling common factors of numerators and denominators. Further reduction is possible noting that a gcd of polynomials in `K[X]` is not unique: it may be multiplied by any (non-zero) constant factor. Any polynomial `f` in `K[X]` can be simplified by extracting the denominators and common factors of the numerators of its coefficients. This yields the representation `f = cF` where `c\in K` is the *content* of `f` and `F` is a *primitive* polynomial, i.e., a polynomial in `A[X]` with coprime coefficients. It is possible to start the algorithm by replacing the given polynomials `f` and `g` with their primitive parts. This will only modify `\mathrm{rem}(f,g)` by a constant factor. Replacing it with its primitive part and continuing recursively we obtain all the primitive parts of the polynomials in the Euclidean remainder sequence, including the primitive `\mathrm{gcd}(f,g)`. This sequence is the *primitive polynomial remainder sequence*. It is an example of *general polynomial remainder sequences* where the computed remainders are modified by constant multipliers (or divisors) in order to simplify the results. Subresultant sequence --------------------- The coefficients of the primitive polynomial sequence do not grow exceedingly, but the computation of the primitive parts requires extra processing effort. Besides, the method only works with fraction fields of unique factorization domains, excluding, for example, the general number fields. Collins :cite:`Collins1967subr` realized that the so-called *subresultant polynomials* of a pair of polynomials also form a generalized remainder sequence. The coefficients of these polynomials are expressible as determinants in the coefficients of the given polynomials. Hence (the logarithm of) their size only grows linearly. In addition, if the coefficients of the given polynomials are in the subdomain `A`, so are those of the subresultant polynomials. This means that the subresultant sequence is comparable to the primitive remainder sequence without relying on unique factorization in `A`. To see how subresultants are associated with remainder sequences recall that all polynomials `h` in the sequence are linear combinations of the given polynomials `f` and `g` .. math:: h = uf+vg with polynomials `u` and `v` in `K[X]`. Moreover, as is seen from the extended Euclidean algorithm, the degrees of `u` and `v` are relatively low, with limited growth from step to step. Let `n = \deg(f)`, and `m = \deg(g)`, and assume `n\ge m`. If `\deg(h) = j < m`, the coefficients of the powers `X^k` (`k > j`) in the products `uf` and `vg` cancel each other. In particular, the products must have the same degree, say, `l`. Then `\deg(u) = l - n` and `\deg(v) = l - m` with a total of `2l -n - m + 2` coefficients to be determined. On the other hand, the equality `h = uf + vg` implies that `l - j` linear combinations of the coefficients are zero, those associated with the powers `X^i` (`j < i \leq l`), and one has a given non-zero value, namely the leading coefficient of `h`. To satisfy these `l - j + 1` linear equations the total number of coefficients to be determined cannot be lower than `l - j + 1`, in general. This leads to the inequality `l \ge n + m - j - 1`. Taking `l = n + m - j - 1`, we obtain `\deg(u) = m - j - 1` and `\deg(v) = n - j - 1`. In the case `j = 0` the matrix of the resulting system of linear equations is the `Sylvester matrix `_ `S(f,g)` associated to `f` and `g`, an `(n+m)\times (n+m)` matrix with coefficients of `f` and `g` as entries. Its determinant is the `resultant `_ `\mathrm{res}(f,g)` of the pair `(f,g)`. It is non-zero if and only if `f` and `g` are relatively prime. For any `j` in the interval from `0` to `m` the matrix of the linear system is an `(n+m-2j)\times (n+m-2j)` submatrix of the Sylvester matrix. Its determinant `s_j(f,g)` is called the `j` th *scalar subresultant* of `f` and `g`. If `s_j(f,g)` is not zero, the associated equation `h = uf + vg` has a unique solution where `\deg(h) = j` and the leading coefficient of `h` has any given value; the one with leading coefficient `s_j(f,g)` is the `j` th *subresultant polynomial* or, briefly, *subresultant* of the pair `(f,g)`, and denoted `S_j(f,g)`. This choice guarantees that the remainining coefficients are also certain subdeterminants of the Sylvester matrix. In particular, if `f` and `g` are in `A[X]`, so is `S_j(f,g)` as well. This construction of subresultants applies to any `j` between `0` and `m` regardless of the value of `s_j(f,g)`; if it is zero, then `\deg(S_j(f,g)) < j`. The properties of subresultants are as follows. Let `n_0 = \deg(f)`, `n_1 = \deg(g)`, `n_2, \ldots, n_k` be the decreasing sequence of degrees of polynomials in a remainder sequence. Let `0 \le j \le n_1`; then - `s_j(f,g)\ne 0` if and only if `j = n_i` for some `i`. - `S_j(f,g)\ne 0` if and only if `j = n_i` or `j = n_i - 1` for some `i`. Normally, `n_{i-1} - n_i = 1` for `1 < i \le k`. If `n_{i-1} - n_i > 1` for some `i` (the *abnormal* case), then `S_{n_{i-1}-1}(f,g)` and `S_{n_i}(f,g)` are constant multiples of each other. Hence either one could be included in the polynomial remainder sequence. The former is given by smaller determinants, so it is expected to have smaller coefficients. Collins defined the *subresultant remainder sequence* by setting .. math:: f_i = S_{n_{i-1}-1}(f,g) \qquad (2\le i \le k). In the normal case, these are the same as the `S_{n_i}(f,g)`. He also derived expressions for the constants `\gamma_i` in the remainder formulas .. math:: \gamma_i f_i = \mathrm{rem}(f_{i-2},f_{i-1}) in terms of the leading coefficients of `f_1,\ldots,f_{i-1}`, working in the field `K`. Brown and Traub :cite:`Brown1971subres` later developed a recursive procedure for computing the coefficients `\gamma_i`. Their algorithm deals with elements of the domain `A` exclusively (assuming `f,g\in A[X]`). However, in the abnormal case there was a problem, a division in `A` which could only be conjectured to be exact. This was subsequently justified by Brown :cite:`Brown1978prs` who showed that the result of the division is, in fact, a scalar subresultant. More specifically, the constant appearing in the computation of `f_i` is `s_{n_{i-2}}(f,g)` (Theorem 3). The implication of this discovery is that the scalar subresultants are computed as by-products of the algorithm, all but `s_{n_k}(f,g)` which is not needed after finding `f_{k+1} = 0`. Completing the last step we obtain all non-zero scalar subresultants, including the last one which is the resultant if this does not vanish. Polynomial factorization in characteristic zero: .. automodule:: diofant.polys.factortools :members: Gröbner basis algorithms ======================== Gröbner bases can be used to answer many problems in computational commutative algebra. Their computation in rather complicated, and very performance-sensitive. We present here various low-level implementations of Gröbner basis computation algorithms; please see the previous section of the manual for usage. .. automodule:: diofant.polys.groebnertools :members: Algebraic number fields ======================= .. currentmodule:: diofant.polys.numberfields .. autofunction:: minpoly_groebner Factorization over algebraic number fields ========================================== .. automodule:: diofant.polys.factorization_alg_field :members: :private-members: Modular GCD =========== .. automodule:: diofant.polys.modulargcd :members: :private-members: Heuristic GCD ============= .. automethod:: diofant.polys.euclidtools._GCD._zz_heu_gcd Further tools ============= .. automodule:: diofant.polys.rootisolation :members: .. automodule:: diofant.polys.sqfreetools :members: Undocumented ============ Many parts of the polys module are still undocumented, and even where there is documentation it is scarce. Please contribute! .. automodule:: diofant.polys.polyoptions :members: Order .. automodule:: diofant.polys.polyerrors :members: .. currentmodule:: diofant.polys.fields .. autoclass:: FracElement :members: