{'fraction_non_alphanumeric': 0.06409855543743356, 'fraction_numerical': 0.01838534175473704, 'mean_word_length': 4.063964534515516, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
GPGCD, an Iterative Method for Calculating Approximate GCD of Univariate Polynomials, with the Complex Coefficients

We present an extension of our GPGCD method, an iterative method for calculating approximate greatest common divisor (GCD) of univariate polynomials, to polynomials with the complex coefficients. For a given pair of polynomials and a degree, our algorithm finds a pair of polynomials which has a GCD of the given degree and whose coefficients are perturbed from those in the original inputs, making the perturbations as small as possible, along with the GCD. In our GPGCD method, the problem of approximate GCD is transfered to a constrained minimization problem, then solved with a so-called modified Newton method, which is a generalization of the gradient-projection method, by searching the solution iteratively. While our original method is designed for polynomials with the real coefficients, we extend it to accept polynomials with the complex coefficients in this paper.


Introduction
For algebraic computations on polynomials and matrices, approximate algebraic algorithms are attracting broad range of attentions recently. These algorithms take inputs with some "noise" such as polynomials with floating-point number coefficients with rounding errors, or more practical errors such as measurement errors, then, with minimal changes on the inputs, seek a meaningful answer that reflect desired property of the input, such as a common factor of a given degree. By this characteristic, approximate algebraic algorithms are expected to be applicable to more wide range of problems, especially those to which exact algebraic algorithms were not applicable.
Among methods in the above, we focus our attention on optimization strategies. Already proposed algorithms utilize iterative methods including the Levenberg-Marquardt method [2], the Gauss-Newton method [19] and the structured total least norm (STLN) method ( [7], [8]). Among them, STLN-based methods have shown good performance calculating approximate GCD with sufficiently small perturbations efficiently.
In this paper, we discuss an extension of the GPGCD method, proposed by the present author [17], an iterative method with transferring the original approximate GCD problem into a constrained optimization problem, then solving it by a so-called modified Newton method [16], which is a generalization of the gradient-projection method [12]. In the previous paper [17], we have shown that our method calculates approximate GCD with perturbations as small as those calculated by the STLN-based methods and with significantly better efficiency than theirs. While our original method accepts polynomials with the real coefficients as inputs and outputs in the previous paper, we extend it to handle polynomials with the complex coefficients in more generalized settings in this paper.
The rest part of the paper is organized as follows. In Section 2, we transform the approximate GCD problem into a constrained minimization problem for the case with the complex coefficients. In Section 3, we show details for calculating the approximate GCD, with discussing issues in minimizations. In Section 4, we demonstrate performance of our algorithm with experiments.

Formulation of the Approximate GCD Problem
Let F (x) and G(x) be univariate polynomials of degree m and n, respectively, with the complex coefficients and m ≥ n > 0. We permit F and G be relatively prime in general. For a given integer d satisfying n ≥ d > 0, let us calculate a deformation of F (x) and G(x) in the form of where ∆F (x) and ∆G(x) are polynomials with the complex coefficients, whose degrees do not exceed those of F (x) and G(x), respectively, H(x) is a polynomial of degree d, andF (x) andḠ(x) are pairwise relatively prime. In this situation, H(x) is an approximate GCD of F (x) and G(x). For a given d, we try to minimize ∆F (x) 2 2 + ∆G(x) 2 2 the norm of the deformations.
In the caseF (x) andG(x) have a GCD of degree d, then the theory of subresultant tells us that the (d−1)-th subresultant ofF andG becomes zero, namely we have S d−1 (F ,G) = 0, where S k (F ,G) denotes the k-th subresultant ofF andG. Then, the (d− 1)-th subresultant matrix N d−1 (F, G), where the k-th subresultant matrix N k (F, G) is a submatrix of the Sylvester matrix N (F, G) by taking the left n − k columns of coefficients of F and the left m − k columns of coefficients of G, has a kernel of dimension equal to 1. Thus, there exist polynomials with deg(A) < n−d and deg(B) < m−d and A(x) and B(x) are relatively prime. Therefore, for the given F (x), G(x) and d, our problem is to find ∆F (x), ∆G(x), A(x) and B(x) satisfying Eq. (1) with making ∆F 2 2 + ∆G 2 2 as small as possible. Let us assume that F (x) and G(x) are represented as G(x) = (g n,1 + g n,2 i )x n + · · · + (g 0,1 + g 0, where f j,1 , g j,1 , f j,2 , g j,2 are the real numbers and i is the imaginary unit, and F Re (x) and G Re (x) represent the real part of F (x) and G(x), respectively, while F Im (x) and G Im (x) represent the imaginary part of F (x) and G(x), respectively. Furthermore, we represent F (x),G(x), A(x) and B(x) with the complex coefficients as A(x) = (a n−d,1 + a n−d,2 i )x n−d + · · · + (a 0,1 + a 0,2 i ) respectively, wheref j,1 ,f j,2 ,g j,1 ,g j,2 , a j,1 , a j,2 , b j,1 , b j,2 are the real numbers, and, as in For the objective function, For the constraint, Eq. (1) becomes as By expressing the subresultant matrix and the column vector in (4) separated into the real and the complex parts, respectively, we express (4) as We can expand the left-hand-side of Eq. (5) as ( (5) is equivalent to a system of equations: Furthermore, we add another constraint for the coefficient of A(x) and B(x) as where Eq. (8) has been put on the top of Eq. (7). Note that, in Eq. (9), we have total of 2(m + n − d + 1) + 1 equations in the coefficients of polynomials in (2) as a constraint, with the j-th row of which is expressed as q j = 0.

The Algorithm for Approximate GCD
We calculate an approximate GCD by solving the constrained minimization problem (11), (12) with the gradient projection method by Rosen [12] (whose initials become the name of our GPGCD method) or a modified Newton method by Tanabe [16] (for review, see the author's previous paper [17]). Our preceding experiments [17, Section 5.1] have shown that a modified Newton method was more efficient than the original gradient projection method while the both methods have shown almost the same convergence property, thus we adopt a modified Newton method in this paper.
In applying a modified Newton method to the approximate GCD problem, we discuss issues in the construction of the algorithm in detail, such as • Representation of the Jacobian matrix J g (x) and certifying that J g (x) has full rank (Section 3.1), • Setting the initial values (Section 3.2), • Regarding the minimization problem as the minimum distance problem (Section 3.3), • Calculating the actual GCD and correcting the coefficients ofF andG (Section 3.4), as follows.

Representation and the rank of the Jacobian Matrix
For a polynomial P (x) ∈ C[x] represented as P (x) = p n x n + · · · + p 0 x 0 , let C k (P ) be a complex (n + k, k + 1) matrix defined as For co-factors A(x) and B(x) in (2), define matrices A 1 and A 2 as (Note that A 1 and A 2 are matrices of the real numbers of m + n − d + 1 rows and m + n + 2 columns.) Then, by the definition of the constraint (12), we have the Jacobian matrix J g (x) (with the original notation of variables (10) for x) as with A 1 and A 2 as in (13) and N 1 , N 2 , v 1 and v 2 as in (6), respectively, which can be easily constructed in every iteration.
In executing iterations, we need to keep that J g (x) has full rank: otherwise, we are unable to decide proper search direction. For this requirement, we have the following observations. Proposition 1. Let x * ∈ V g be any feasible point satisfying Eq. (12). Then, if the corresponding polynomials do not have a GCD whose degree exceeds d, then J g (x * ) has full rank.
Proof. Let x * = (f m,1 , . . . ,f 0,1 ,g n,1 , . . . ,g 0,1 ,f m,2 , . . . ,f 0,2 ,g n,2 , . . . ,g 0,2 , a n−d,1 , . . . , a 0,1 , b m−d,1 , . . . , b 0,1 , a n−d,2 , . . . , a 0,2 , b m−d,2 , . . . , b 0,2 ) with its polynomial representation expressed as in (2) (note that this assumption permits the polynomialsF (x) andG(x) to be relatively prime in general). To verify our claim, we show that we have rank(J g (x * )) = 2(m + n − d + 1) + 1. Let us express J g (x * ) = J L | J R , where J L and J R are column blocks expressed as respectively. Then, we have the following lemma.  N (A, B), the Sylvester matrix of A(x) and B(x). By the assumption, polynomials A(x) and B(x) are relatively prime, and there exist no nonzero elements inC except for the bottom m + n − 2d rows, thus we have rank(C) = m + n − 2d. By the structure ofC and the lower triangular structure of C m (A) 1 and C n (B) 1 , we can take the left d + 1 columns of C m (A) 1 or C n (B) 1 satisfying linear independence along withC, which implies that there exist a nonsingular square matrix T of order m + n + 2 satisfying where R is a lower triangular matrix, thus we have rank(A 1 ) = rank(R) = m + n − d + 1.
Furthermore, by using T and R in (14), we have followed by a suitable transformation on columns on the matrix in the right-hand-side of (15), we can make A 2 T to zero matrix, which implies that This proves the lemma.
Proof of Proposition 1 (continued). By the assumptions, we have at least one nonzero coordinate in the top row in J R , while we have no nonzero coordinate in the top row in J L , thus we have rank(J g (x)) = 2(m + n − d + 1) + 1, which proves the proposition.
Proposition 1 says that, so long as the search direction in the minimization problem satisfies that corresponding polynomials have a GCD of degree not exceeding d, then J g (x) has full rank, thus we can safely calculate the next search direction for approximate GCD.
However, we can regard our problem as finding a point x ∈ V g which has the minimum distance to the initial point x 0 with respect to the (x 1 , . . . , x 2(m+n+2) )-coordinates which correspond to the coefficients in F (x) and G(x). Therefore, as in the real case (see the authors previous paper [17]), we change the objective function asf (x) = 1 2 f (x), then solve the minimization problem off (x), subject to q(x) = 0.

Calculating the Actual GCD and Correcting the Deformed Polynomials
After successful end of the iterations, we obtain the coefficients ofF (x),G(x), A(x) and B(x) satisfying (1) with A(x) and B(x) are relatively prime. Then, we need to compute the actual GCD H(x) ofF (x) andG(x). Although H can be calculated as the quotient of F divided by B orG divided by A, naive polynomial division may cause numerical errors in the coefficient. Thus, we calculate the coefficients of H by the so-called least squares division [19], followed by correcting the coefficients inF andG by using the calculated H, as follows. For polynomialsF ,G, A and B represented as in (2) and H represented as H(x) = (h d,1 + h d,2 i )x d + · · · + (h 0,1 + h 0,2 i )x 0 , solve the equations HB =F and HA =G with respect to H as solving the least squares problems of linear systems . . , h 0,1 + h 0,2 i ) = t (g n,1 +g n,2 i , . . . ,g 0,1 +g 0,2 i ), respectively. Then, we transfer the linear systems (16) and (17), as follows. For (17), let us express the matrices and vectors as the sum of the real and the imaginary part of which, respectively, as . . ,f 0,1 +f 0,2 i ) = f 1 + i f 2 . Then, (17) is expressed as By equating the real and the imaginary parts in Eq. (18), respectively, we have ( Thus, we can calculate the coefficients of H(x) by solving the real least squares problem (19). We can solve (16) similarly.
be the candidates for the GCD whose coefficients are calculated as the least squares solutions of (16) and (17), respectively. Then, for i = 1, 2, calculate the norms of the residues as r i = F − H i B 2 2 + G − H i A 2 2 , respectively, and set the GCD H(x) be H i (x) giving the minimum value of r i .

Experiments
We have implemented the GPGCD algorithm for polynomials with the complex coefficients on the computer algebra system Maple and compared its performance with a method based on the structured total least norm (STLN) method [7] for randomly generated polynomials with approximate GCD. The tests have been carried out on Intel Core2 Duo Mobile Processor T7400 (in Apple MacBook "Mid-2007" model) at 2.16 GHz with RAM 2GB, under MacOS X 10.5.
In the tests, we have generated random polynomials with GCD then added noise, as follows. First, we have generated a pair of monic polynomials F 0 (x) and G 0 (x) of degrees m and n, respectively, with the GCD of degree d. The GCD and the prime parts of degrees m − d and n − d are generated as monic polynomials and with random coefficients c ∈ [−10, 10] of floating-point numbers. For noise, we have generated a pair of polynomials F N (x) and G N (x) of degrees m − 1 and n − 1, respectively, with random coefficients as the same as for F 0 (x) and G 0 (x). Then, we have defined a pair of test polynomials F (x) and G(x) as respectively, scaling the noise such that the 2-norm of the noise for F and G is equal to e F and e G , respectively. In the present test, we set e F = e G = 0.1.
In this test, we have compared our implementation against a method based on the structured total least norm (STLN) method [7], using their implementation (see Acknowledgments). In their STLN-based method, we have used the procedure C_con_mulpoly which calculates the approximate GCD of several polynomials in C[x]. The tests have been carried out on Maple 12 with Digits=15 executing hardware floating-point arithmetic. For every example, we have generated 100 random test polynomials as in the above. In executing a modified Newton method, we set a threshold of the 2-norm of the "update" vector in each iteration ε = 1.0 × 10 −8 ; in C_con_mulpoly, we set the tolerance e = 1.0 × 10 −8 . Table 1 shows the results of the test: m and n denotes the degree of a pair F and G, respectively, and d denotes the degree of approximate GCD. The columns with "STLN" are the data for the STLN-based method, while those with "GPGCD" are the data for the GPGCD method. "Error" is the perturbation F − F 2 2 + G − G 2 2 , where "ae−b" denotes a × 10 −b ; "#Iterations" is the number of iterations; "Time" is computing time in seconds.
We see that, in the most of tests, both methods calculate approximate GCD with almost the same amount of perturbations, while the GPGCD method runs much faster than STLNbased method by approximately from 10 to 30 times. Note that, in contrast to the real coefficient case [17], both methods have converged in all the test cases with the number of iterations and sufficiently small amount of perturbations as approximately equal to those shown as in Table 1.

Concluding Remarks
Based on our previous research [17], we have extended our GPGCD method for polynomials with the complex coefficients.
Our experiments have shown that, as in the real coefficients case [17], our algorithm calculates approximate GCD with perturbations as small as those calculated by methods based on the structured total least norm (STLN) method, while our method has shown significantly better performance over the STLN-based methods in its speed, by approximately up to 30 times, which seems to be sufficiently practical for inputs of low or moderate degrees. This result shows that, in contrast to their structure preserving method, our simple method can achieve accurate and efficient computation as or more than theirs in calculating approximate GCDs. Our future research includes theoretical investigation of convergence properties, investigation for efficient computation in solving a linear system in each iteration by analysis of the structure of matrices, generalization of our method to several input polynomials, and so on.