Toward the best algorithm for approximate GCD of univariate polynomials

Abstract

Approximate polynomial GCD (greatest common divisor) of polynomials with a priori errors on their coefficients, is one of interesting problems in Symbolic-Numeric Computations. In fact, there are many known algorithms: QRGCD, UVGCD, STLN based methods, Fastgcd and so on. The fundamental question of this paper is “which is the best?” from the practical point of view, and subsequently “is there any better way?” by any small extension, any effect by pivoting, and any combination of sub-routines along the algorithms. In this paper, we consider a framework that covers those algorithms and their sub-routines, and makes their sub-routines being interchangeable between the algorithms (i.e. disassembling the algorithms and reassembling their parts). By this framework along with/without small new extensions and a newly adapted refinement sub-routine, we have done many performance tests and found the current answer. In summary, 1) UVGCD is the best way to get smaller tolerance, 2) modified Fastgcd is better for GCD that has one or more clusters of zeros with large multiplicity, and 3) modified ExQRGCD is better for GCD that has no cluster of zeros.

Introduction

Finding the greatest common divisor (GCD) of the given two polynomials is one of essential computations in computer algebra. In fact, most of other algebraic computations (e.g. reduction of fraction, factorization and so on) rely on the polynomial GCD. In general, finding the GCD can be done by the Euclidean algorithm or other fast algorithms (e.g. halfGCD and its variants (Roy and Sedjelmaci, 2013)). However, in some empirical or practical situations, the coefficients of the given polynomials may have a priori errors, and any conventional algorithm can not compute their expected GCD since their perturbations arisen from the errors make the given polynomials being coprime even if they had some non-trivial common factor in nature. For example, for $f(x)=0.999x^2+1.999x+1.001\approx {(x+1)}^2$ and $g(x)=1.001x^2-0.999\approx (x+1)(x-1)$, the Euclidean algorithm can not compute any non-trivial GCD, even if we would like to have their GCD like $d(x)=1.00063x+0.999375$.

To overcome this problem, for the couple of decades, several algorithms have been studies and the resulting quasi GCD is called “approximate polynomial GCD”. In the early stage of literatures, the first sight of the related problem is appeared in Dunaway (1974) where the Euclidean algorithm with numerical computations is proposed. In this stage, the other algorithms (Schönhage, 1985; Noda and Sasaki, 1991; Sasaki and Noda, 1989) follows the same approach hence their results can be considered as numerical variants of exact GCD. The turning point tending toward the modern definition (e.g. a non-trivial GCD of polynomials in the neighborhoods of the given polynomials over a number field without any numerical fashion) is given by Emiris et al., 1997, Emiris et al., 1996, and this definition or similar definitions are widely used after that.

In the present day, there are many algorithms for approximate GCD of univariate polynomials: QRGCD (Corless et al., 2004), UVGCD (Zeng, 2011), STLN based methods (Kaltofen et al., 2006, Kaltofen et al., 2007), Fastgcd (Bini and Boito, 2007), ExQRGCD (Nagasaka and Masui, 2013) and GPGCD (Terui, 2013) that are focused in this paper, and the others (Fazzi et al. (2019), Usevich and Markovsky (2017), Bourne et al. (2017), Christou et al. (2010) and so on). Therefore, one may have the very fundamental question “which is the best?” especially from the practical point of view. One of our contributions is answering to this question. One may think that each research paper has some comparison with known methods hence the question must exist in the papers. This is partially correct. However, in general, those comparisons are based on different programming languages, computer environments, implementors and code optimization levels, hence any well-considered comparison under an uniform treatment is required in addition to those comparisons.

Moreover, it is notable and important that most of these algorithms rely on matrix factorizations of some structured matrices and optimization techniques (for theoretical fundamentals and most of literatures can be found in Boito (2011)). This means that we can consider a framework that covers these algorithms and their sub-routines and makes their sub-routines being interchangeable between the algorithms. By this framework, we can disassemble the algorithms and reassemble their parts as new algorithms. Therefore, we have some answers to the subsequent open question “is there any better way?” along with/without small new extensions (including any effect by pivoting) and a newly adapted sub-routine (based on the approach in Giesbrecht et al. (2017)). This is also one of contributions of this paper and our result is given in Sections 4 and 5. In summary, 1) UVGCD is the best way to get smaller tolerance (a kind of radius of neighborhoods), 2) modified Fastgcd is better for GCD that has one or more clusters of zeros with large multiplicity, and 3) modified ExQRGCD is better for GCD that has no cluster of zeros.

After the preliminary subsection including some definitions and notations below, we start with the several brief reviews of known algorithms, and introduce some small new extensions and a newly adapted sub-routine in Section 2. Our framework is given in Section 3, and our result of some hard performance test is given in Section 4. Finally, in Section 5, we give our recommendation among the algorithms.

We assume that polynomials $f(x)={\sum }_{i=0}^m{f}_i{x}^i\in K[x]$ and $g(x)={\sum }_{i=0}^n{g}_i{x}^i\in K[x]$ are given with the assumption $m\ge n$. The domain $K[x]$ that $f(x)$ and $g(x)$ belong to is $\mathbb{R}[x]$ or $\mathbb{C}[x]$ depending on each problem setting. Each of these given polynomials has the unit Euclidean norm (2-norm) (i.e. ${\Vert f(x)\Vert }_2={\Vert g(x)\Vert }_2=1$). We scale the polynomials if they do not have the unit norm at the input. Although there are several variants of definitions, we use the following two definitions of approximate GCD.

Definition 1 approximate polynomial GCD with the given tolerance

For the tolerance $\epsilon \in {\mathbb{R}}_{\ge 0}$, we compute the polynomial $d(x)\in K[x]$ called “approximate polynomial GCD” of tolerance ε, which has the maximum degree and satisfies

$$f(x)+{\mathrm{\Delta }}_f(x)=f_1(x)d(x),g(x)+{\mathrm{\Delta }}_g(x)=g_1(x)d(x)$$

for some polynomials ${\mathrm{\Delta }}_f(x),{\mathrm{\Delta }}_g(x),f_1(x),g_1(x)\in K[x]$ such that

$$\mathrm{deg}({\mathrm{\Delta }}_f)\le \mathrm{deg}(f),\mathrm{deg}({\mathrm{\Delta }}_g)\le \mathrm{deg}(g),{\Vert {\mathrm{\Delta }}_f\Vert }_2<\epsilon {\Vert f\Vert }_2,{\Vert {\mathrm{\Delta }}_g\Vert }_2<\epsilon {\Vert g\Vert }_2.$$

Definition 2 approximate polynomial GCD with the given degree

For the degree $k\in \mathbb{N}$, we compute the polynomial $d(x)\in K[x]$ called “approximate polynomial GCD” of degree k, which minimizes ${\Vert {\mathrm{\Delta }}_f\Vert }_2^2+{\Vert {\mathrm{\Delta }}_g\Vert }_2^2$ (called perturbation) and satisfies

$$f(x)+{\mathrm{\Delta }}_f(x)=f_1(x)d(x),g(x)+{\mathrm{\Delta }}_g(x)=g_1(x)d(x),\mathrm{deg}(d)=k$$

for some polynomials ${\mathrm{\Delta }}_f(x),{\mathrm{\Delta }}_g(x),f_1(x),g_1(x)\in K[x]$ such that

$$\mathrm{deg}({\mathrm{\Delta }}_f)\le \mathrm{deg}(f),\mathrm{deg}({\mathrm{\Delta }}_g)\le \mathrm{deg}(g).$$

Remark 1

Most of algorithms for Definition 1 do not guarantee the maximum degree but try to find a factor of nearly maximum degree by floating-point calculations (see Emiris et al. (1997) for a maximum degree condition). Also for Definition 2, most of algorithms do not guarantee the minimum perturbation but try to find a factor with a nearly minimum perturbation by floating-point calculations, and there are similar but different objective functions (e.g. $\max ({\Vert {\mathrm{\Delta }}_f\Vert }_2,{\Vert {\mathrm{\Delta }}_g\Vert }_2)$). Moreover, there are algorithms for polynomials in $\mathbb{Z}[x]$ (hence by exact calculations), proposed by von zur Gathen et al. (2010) and Nagasaka (2011) for example.

Remark 2

In this paper, we seek the best algorithm among the algorithms we are interested in. The criteria of our evaluation is very simple. A larger detecting degree is better (and a smaller resulting tolerance is better if the same degree) for approximate polynomial GCD with the given tolerance (Definition 1), and a smaller resulting perturbation is better for approximate polynomial GCD with the given degree (Definition 2). However, too much elapsed time may not be a good behavior from the practical point of view. Therefore, we consider their elapsed time also.

For vectors and matrices, ${\Vert \cdot \Vert }_2$ and ${\Vert \cdot \Vert }_F$ denote the Euclidean norm (2-norm) and the Frobenius norm, respectively. $A^T$, $A^H$ and $A^{-1}$ are the transpose, the conjugate transpose, and the inverse of matrix A, respectively. The coefficient vector (in the descending term order) of $p(x)$ of degree k is denoted by $\overrightarrow{p}\in K^{k+1}$. The Sylvester matrix of $f(x)$ and $g(x)$ is denoted by $\mathrm{\text{Syl}}(f,g)$ defined as

$$\mathrm{\text{Syl}}(f,g)=\left(\begin{array}{ccccc}\hfill f_m\hfill & \hfill f_{m-1}\hfill & \hfill \cdots \hfill & \hfill f_1\hfill & \hfill f_0\hfill \\ \hfill \hfill & \hfill f_m\hfill & \hfill f_{m-1}\hfill & \hfill \cdots \hfill & \hfill f_1\hfill & \hfill f_0\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \cdots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill f_m\hfill & \hfill f_{m-1}\hfill & \hfill \cdots \hfill & \hfill f_1\hfill & \hfill f_0\hfill \\ \hfill g_n\hfill & \hfill g_{n-1}\hfill & \hfill \cdots \hfill & \hfill g_1\hfill & \hfill g_0\hfill & \hfill \hfill \\ \hfill \hfill & \hfill g_n\hfill & \hfill g_{n-1}\hfill & \hfill \cdots \hfill & \hfill g_1\hfill & \hfill g_0\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \cdots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill g_n\hfill & \hfill g_{n-1}\hfill & \hfill \cdots \hfill & \hfill g_1\hfill & \hfill g_0\hfill \end{array}\right)\in K^{(m+n)\times (m+n)}.$$

The subresultant matrix of order k of $f(x)$ and $g(x)$ ($n-k$ columns for $\overrightarrow{f}$ and $m-k$ columns for $\overrightarrow{g}$) is denoted by ${\mathrm{\text{S}}}_k(f,g)$ defined as

$${\mathrm{\text{S}}}_k(f,g)=\left(\begin{array}{ccccccc}\hfill f_m\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill g_n\hfill & \hfill \hfill & \hfill \hfill \\ \hfill f_{m-1}\hfill & \hfill f_m\hfill & \hfill \hfill & \hfill \hfill & \hfill g_{n-1}\hfill & \hfill g_n\hfill & \hfill \hfill \\ \hfill ⋮\hfill & \hfill f_{m-1}\hfill & \hfill \ddots \hfill & \hfill \hfill & \hfill ⋮\hfill & \hfill g_{n-1}\hfill & \hfill \ddots \hfill & \hfill \hfill \\ \hfill f_1\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill \hfill & \hfill g_1\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill \hfill \\ \hfill f_0\hfill & \hfill f_1\hfill & \hfill ⋮\hfill & \hfill \hfill & \hfill g_0\hfill & \hfill g_1\hfill & \hfill ⋮\hfill & \hfill \hfill \\ \hfill \hfill & \hfill f_0\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill & \hfill \hfill & \hfill g_0\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill & \hfill f_1\hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill & \hfill g_1\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill f_0\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill g_0\hfill \end{array}\right)\in K^{(m+n-k)\times (m+n-2k)}.$$

Throughout this paper, the capital letters D, L, U, P, Q and R denote diagonal, lower triangular, upper triangular, permutation, unitary (orthogonal) and upper triangular matrices, respectively.