Algorithms for algebraic computations¶

Arithmetics of polynomials¶

Arithmetics, i.e. addition, subtraction, multiplication, exponentiation and division, of polynomials form a basis for all other polynomials manipulation algorithms. In SymPy we currently implement only, so called, classical algorithms for this purpose, i.e. repeated squaring algorithm for exponentiation, and O(n^2) algorithms for multiplication and division, where n is the maximal number of terms in the set of input polynomials. Over finite fields we use algorithms of [Monagan1993inplace], which slightly improve speed of computations over this very specific domain.

An alternative to classical algorithms are, so called, fast algorithms, which are sub–quadratic time algorithms for doing arithmetics polynomials [Moenck1976practical], [Bernstein2008fast]. Fast algorithms are usually limited to specific domains of computation, like integers or rationals. The family includes Karatsuba’s and FFT (Fast Fourier Transform) algorithms. The decision was made to use classical algorithms at this point, because it is not a trivial task to make fast algorithms really advantageous, especially for small or ill conditioned polynomials. In future, when the module will stabilize, we will consider implementing fast algorithms, as a companion to classical algorithms, and use them where it makes sense.

There are other ideas to improve arithmetics of polynomials, especially over integers and rationals. An interesting example of such optimisation is algorithm of [Fateman2005encoding], where polynomials with integer coefficients are encoded as sufficiently large integers and arithmetics are done using long integer arithmetics, and later results are converted back to polynomials. This is possible by using an isomorphism (a reversible transformation) between polynomials and integers. To make this algorithm efficient, we need a very fast integer arithmetics library, like gmpy, which is not always available on the system. One has to also be aware of the fact that transformations back and forth may be very costly, especially for large polynomials.

Evaluation of polynomials¶

In SymPy we employ Horner scheme [Geddes1992algorithms] for evaluation of univariate polynomials over arbitrary domains. The algorithm of Horner was proved to be the optimal method for evaluation of polynomials, which takes minimum number of additions and multiplications necessary to compute a result. In the multivariate case, Horner scheme is non–unique, thus there are many different schemes possible, which lead to different evaluation times. Currently we use a natural scheme, which is a consequence of using by default recursive polynomial representation. In the past we experimented with greedy algorithms for optimizing multivariate Horner scheme [Ceberio2004greedy], however, without much success. This was because, although evaluation was considerably faster comparing to the standard scheme, however optimisation times were often comparable to computation times. In future we may reconsider using some sort of optimisation of polynomial evaluation algorithm in the multivariate case.

Horner scheme is a very general algorithm, which is also used in SymPy for computing compositions and rational transformations of polynomials, and many other, which require some sort of efficient evaluation of polynomials.

The Greatest Common Divisor¶

Another fundamental algorithm of polynomials manipulation is the GCD (Greatest Common Divisor) algorithm, which allows us to compute common factors of two or more polynomials. This is very useful on its own and as a component of other algorithms, e.g. square–free decomposition or simplification of rational expressions.

We implement three algorithms for computing GCDs. The most general algorithm is based on subresultants (a special case of polynomial remainder sequences) and can be used regardless of the ground domain and polynomial representation. This is also the slowest algorithm, but very useful if other algorithms fail or are not applicable. Where possible, we use heuristic GCD algorithm [Liao1995heuristic], which transforms a problem of computing GCD of polynomials to integer GCD problem. Although the algorithm is heuristic, with current parametrization it never failed in SymPy. Daily practice shows that this approach is superior to the algorithm based on subresultants, however heuristic GCD is only limited to integers and rationals (by clearing denominators). It also requires very efficient integer GCD algorithm, so it is beneficial to use gmpy library for this purpose. We also implement an algorithm that uses Gröbner bases for computing GCD of multivariate polynomials [Cox1997ideals]. This was historically the first implementation of multivariate polynomial GCD algorithm in SymPy (see section Special case 2: Euclid’s algorithm for details).

In future we plan to implement EEZ–GCD algorithm of Wang [Wang1980eezgcd], [Moses1973ezgcd], which we hope will improve computation of GCDs of sparse multivariate polynomials over integers and rationals. Implementation of EEZ–GCD algorithm should be rather straightforward, because we already have EEZ polynomial factorization algorithm implemented, and both algorithms share a common core (variable–by–variable Hensel lifting algorithm), and they differ mainly in the initialization phase. Having a fast algorithm for sparse polynomials is very important, because most polynomials that we encounter in real life problems are sparse.

An alternative would be to implement sparse modular algorithm (SPMOD) of Zippel, which is also optimized for sparse multivariate case. This algorithm is even more interesting in the light of recent developments [Monagan2004algebraic], [Javadi2007spmod], where the algorithm was successfully employed over algebraic number and function fields. Although there is a simple idea standing behind SPMOD, this algorithm is considered, in the literature, to be very hard to implement, because there are many special cases, which have to properly worked out. Thus, at least for optimizing GCD computations over integers and rationals, we see EEZ algorithm more beneficial at the moment.

Some parts of modular algorithms can be relatively easily parallelized on multiple processors, because it often happens that several computations over a smaller domain have to be performed to compute the GCD in the original domain. For example this is the case in EEZ–GCD algorithm where, as one of initialization steps, we need to compute several univariate GCDs to compute the original multivariate GCD. Those univariate GCDs can be very costly, but we have to do several such computations (at least three) to guarantee correctness of EEZ algorithm (otherwise the algorithm may fail and would have to be restarted and another sequence of univariate polynomial GCDs would have to be computed). Those univariate GCDs can be computed in parallel, greatly improving speed of multivariate GCD algorithm.

Square–free decomposition¶

Given a polynomial f, square–free decomposition (factorization) of f gives a list of polynomials (factors) f_1, f_2, \ldots, f_n, such that all pairs of polynomials (f_i, f_j), for i \not= j, are co–prime, and f = f_1 f_2^2 \ldots f_n^n. Thus each f_i has no repeated roots. Note that square–free decomposition does not give a true factorization into irreducibles, although is a very important step in any factorization algorithm (which we will describe in the following section).

In SymPy we implement the fast algorithm of Yun [Yun1976squarefree] for computing square–free decompositions in domains of characteristic zero. The cost of computing square–free decomposition is equivalent to the computation of the greatest common divisor of f and its derivative. Over finite fields we currently use less efficient algorithm, due to odd but well known behaviour of derivatives over domains with finite number of elements, where f' might vanish even if f is non–constant polynomial (e.g. f = x^k over \F_k, for any k >= 2), which leads to complications in the algorithm. In future we should implement Yun’s algorithm also in this case as well.

Factorization of polynomials¶

Algorithms for computing factorizations of polynomials into irreducibles over various domains are the landmark of symbolic mathematics. The work in this area started early, in ninetieth century, and algorithms for factoring of univariate and multivariate polynomials over rationals were invented by Kronecker. Those algorithms had exponential time complexity and were impractical for any real–life computations. Kronecker’s algorithms were the first polynomial factorization algorithms that were implemented in SymPy.

Over the years mathematicians tried to invent a more efficient method for factoring large polynomials, focusing their research on univariate polynomials with integer coefficients. A partial breakthrough came with Hensel’s lemma (Hensel lifting algorithm), which allowed to perform computations with integer valued polynomials over finite fields. Thus, an integer polynomial problem can be transformed into finite field polynomial problem, then computations can be done in a much smaller (finite) domain and results can be transformed (lifted) back to the integer polynomial domain.

In 1967, Berlekamp gave first complete factorization algorithm over finite fields. Moreover, this algorithm had polynomial time complexity. Soon after, Zassenhaus combined Berlekamp’s algorithm and Hensel’s lemma, giving first efficient algorithm for polynomial factorization over integers. Although, the algorithm of Zassenhaus had still exponential time worst–case complexity, on average it behaved as polynomial time algorithm, allowing to compute with large polynomials with integer coefficients. This breakthrough stimulated the society and, in the following twenty years after Zassenhaus algorithm was invented, many other algorithms were invented, covering wider range of coefficient domains and introducing modular techniques (Hensel’s lemma) to the case of multivariate polynomials.

Gröbner bases¶

The method of Gröbner bases is a powerful technique for solving problems in commutative algebra (polynomial ideal theory, algebraic geometry). In chapter Gröbner bases and their applications we will describe Gröbner bases in very detail, so we will skip any further discussion on this topic in this section, to avoid redundancy.

Root isolation¶

Polynomials are solvable by radicals only up to degree 4 (inclusive). This is an unfortunate but well known consequence of Abel–Ruffini theorem. SymPy implements heuristic algorithms for solving polynomials in terms of radicals in roots() function. In some cases it is possible to find roots of higher degree polynomials, by taking advantage of polynomial factorization and decomposition algorithms, and pattern matching.

This is an obviously limited approach and there is a need, in various areas of symbolic mathematics, e.g. when solving of systems of polynomial equations [Strzebonski1997computing], to compute values of roots of a polynomial to a desired precision. This could be done by using numerical root finding algorithms, like Durand–Kerner’s, which has its implementation in mpmath library and is exposed to the top–level via nroots() function. However, in pathological cases, numerical algorithms may fail to compute correct values of polynomials’ roots.

To tackle this problem, when the user needs guaranteed error bounds of the computed roots, symbolic root isolation algorithms should be used. SymPy can isolate roots of polynomials with rational coefficients over real and complex domains, taking advantage of most recent algorithmic developments in the field. Symbolic root isolation is not that efficient as numerical root finding, but it is always successful for arbitrary polynomials, giving, as the result, isolation intervals of the roots of a polynomial, in the real case, or isolation rectangles, in the complex case.

Conclusions¶

In this chapter we gave a brief description to the most important algorithms of polynomials manipulation module. There are other algorithms that were implemented in the module, which also deserve attention and a few words of explanation. Hopefully in some foreseeable future we will be able to write a more capacious volume, in which we will describe all of them. We also gave references to the most influential literature that was used to implement those algorithms or seems promising for further developments in near future. We consider this a good starting point for people who would be interested in picking up some development tasks to improve the module.