Computing GCRDs of Approximate Differential Polynomials
CComputing GCRDs of ApproximateDifferential Polynomials
Mark Giesbrecht and Joseph HaraldsonCheriton School of Computer ScienceUniversity of WaterlooWaterloo, Ontario { mwg,jharalds } @uwaterloo.caOctober 29, 2018 Abstract
Differential (Ore) type polynomials with approximate polynomialcoefficients are introduced. These provide a useful representation ofapproximate differential operators with a strong algebraic structure,which has been used successfully in the exact, symbolic, setting. Wethen present an algorithm for the approximate Greatest CommonRight Divisor (GCRD) of two approximate differential polynomials,which intuitively is the differential operator whose solutions are thosecommon to the two inputs operators. More formally, given approxi-mate differential polynomials f and g , we show how to find “nearby”polynomials (cid:101) f and (cid:101) g which have a non-trivial GCRD. Here “nearby” isunder a suitably defined norm. The algorithm is a generalization of theSVD-based method of Corless et al. (1995) for the approximate GCDof regular polynomials. We work on an appropriately “linearized”differential Sylvester matrix, to which we apply a block SVD. Thealgorithm has been implemented in Maple and a demonstration of itsrobustness is presented. a r X i v : . [ c s . S C ] J u l Introduction
The ring of differential (Ore) polynomials R ( t )[ ∂ ; (cid:48) ] over the real numbers R provides a (non-commutative) polynomial ring structure to the linear ordi-nary differential operators. Differential polynomials have found great utilityin symbolic computation, as they allow us to apply algebraic tools to thesimplification and solution of linear differential equations; see Bronstein andPetkovˇsek (1994) for a nice introduction to the mathematical and compu-tational aspects. The ring of differential polynomials R ( t )[ ∂ ; (cid:48) ] is defined asthe usual polynomials in R ( t )[ ∂ ] (i.e., polynomials in ∂ over the commutativefield of rational functions R ( t )), under the usual polynomial addition and thenon-commutative multiplication rule ∂f ( t ) = f ( t ) ∂ + f (cid:48) ( t ) for f ( t ) ∈ R ( t ) , where f (cid:48) ( t ) is the usual formal derivative of f ( t ) with respect to t . Thisdefinition of R ( t )[ ∂ ; (cid:48) ] is useful because there is a natural action of R ( t )[ ∂ ; (cid:48) ]on the space C ∞ of infinitely differentiable functions y ( t ) : R → R . Inparticular, for any y ( t ) ∈ C ∞ , f ( ∂ ) = (cid:88) ≤ i ≤ k f i ( t ) ∂ i acts on y ( t ) as (cid:88) ≤ i ≤ k f i ( t ) d i dt i y ( t ) . A primary benefit of viewing differential operators in this way is thatthey have the structure of a left (and right) Euclidean domain. In particular,for any two polynomials f, g ∈ R ( t )[ ∂ ; (cid:48) ], there is a unique polynomial h ∈ R ( t )[ ∂ ; (cid:48) ] of maximal degree in ∂ such that f = uh and g = vh for u, v ∈ R ( t )[ ∂ ; (cid:48) ] (i.e., h divides f and g exactly on the right). This polynomial h iscalled the GCRD of f and g and it is unique up to a scalar multiplication by anon-zero element of R ( t ) (we could make this GCRD have leading coefficient1, but that would introduce denominators from R [ t ], as well as potentialnumerical instability, as we shall see).The important geometric interpretation is that there is an algorithm todetermine the differential polynomial whose solution spaces is the intersectionof the solution space of f and g ; this is precisely h = gcrd( f, g ).The goal of this paper is to devise an efficient, numerically robust algo-rithm to compute the GCRD when the coefficients in R are given approxi-mately. Specifically, given f, g ∈ R ( t )[ ∂ ; (cid:48) ], we wish to find (cid:101) f , (cid:101) g ∈ R ( t )[ ∂ ; (cid:48) ],where (cid:101) f is “near” f and (cid:101) g is “near” g such that deg ∂ gcrd( (cid:101) f , (cid:101) g ) ≥
1. That2s, (cid:101) f and (cid:101) g have a non-trivial GCRD. The precise definition of nearness isgiven below. The approach we will take is similar to the one developed inCorless et al. (1995) for regular polynomials, which is to reduce the problemto a singular value decomposition or total least squares problem.The problem of computing the GCRD in a symbolic and exact settingdates back to Ore (1933), who presents a Euclidean like algorithm. SeeBronstein and Petkovˇsek (1994) for an elaboration of this approach. Li andNemes (1997) introduce a differential-resultant-based algorithm which makescomputation of the GCRD very efficient using modular arithmetic. We gen-eralize and adapt their approach to a numerical setting here.The analogous approximate GCD problem for approximate regular (com-mutative) polynomials has been a key topic of research in symbolic-numericcomputing since its inception. A full survey is not possible here, but wenote the deep connection between our current work and that of Corless et al.(1995); see also Karmarkar and Lakshman (1996), Sasaki and Sasaki (1997),and Zeng and Dayton (2004). Also important to this current work is the useof so-called structured (numerical) matrix methods for approximate GCD,such as structured total least squares (STLS) and structured total least norm(STLN); see Botting et al. (2005) and Kaltofen et al. (2005). A structuredapproach to relative primality is taken in Beckermann and Labahn (1998).More directly employed later in this paper is the multiple polynomial approx-imate GCD method of Kaltofen et al. (2006). This latter paper also providesa nice survey of recent developments.
The ring R ( t )[ ∂ ; (cid:48) ] is a non-commutative principal left (and right) ideal do-main. For f, g ∈ R ( t )[ ∂ ; (cid:48) ], with deg ∂ f = n and deg ∂ g = m , we have thefollowing properties (Ore, 1933).1. deg ∂ ( f g ) = deg ∂ f + deg ∂ g ,deg ∂ ( f + g ) ≤ max { deg ∂ f, deg ∂ g } .2. There exist q, r ∈ R ( t )[ ∂ ; (cid:48) ] with deg ∂ r < deg ∂ g such that f = qg + r (Division with Remainder).3. There exists h ∈ R ( t )[ ∂, (cid:48) ] of maximal degree in ∂ with f = w h and g = w h . h is called the GCRD (Greatest Common Right Divisor) of f and g . 3. There exist w , w ∈ R ( t )[ ∂, (cid:48) ] such that w f = w g = h for h ofminimal degree. h is called the LCLM (Least Common Left Multiple)of f and g .5. deg ∂ lclm( f, g ) = deg ∂ f + deg ∂ g − deg ∂ gcrd( f, g ).These immediately imply the following characterization of a non-trivialGCRD. Lemma 1.1.
Suppose f, g ∈ R ( t )[ ∂ ; (cid:48) ] , with deg ∂ f = m and deg ∂ g = n .Then deg ∂ gcrd( f, g ) ≥ if and only if there exists u, v ∈ R ( t )[ ∂ ; (cid:48) ] such that deg ∂ u < n , deg ∂ v < m , and uf + vg = 0 . This lemma will enable us to set up a resultant-like linear system over R for the GCRD, which will lead to the desired algorithms.Because of the non-commutative property of R ( t )[ ∂ ; (cid:48) ], it will be importantto maintain a canonical form for any f ∈ R ( t )[ ∂ ; (cid:48) ]. We will always write f = 1 f − ( t ) m (cid:88) i =0 f i ( t ) ∂ i , for polynomials f − , f , . . . , f m ∈ R [ t ], with coefficients in R ( t ) always writtento the left of powers of ∂ . Moreover, for f as above, and (cid:96) > deg f , we defineΨ (cid:96) ( f ) = 1 f − · ( f , f , . . . , f m , , . . . , ∈ R ( t ) (cid:96) . I.e., Ψ (cid:96) maps polynomials in R ( t )[ ∂ ; (cid:48) ] of degree (in ∂ ) less than (cid:96) into R ( t ) (cid:96) .It will also be useful to ensure that our differential polynomials are prim-itive . Definition 1.2 (Primitive Differential Polynomial) . Let f ∈ R [ t ][ ∂ ; (cid:48) ] where deg ∂ f = m is in standard form. The content of f is given by cont( f ) =gcd( f , f , . . . , f m ) . If cont( f ) = 1 , we say that f is primitive. For our primary problem of computing GCRDs of differential polynomials f, g ∈ R ( t )[ ∂ ; (cid:48) ], we will assume both that the coefficients of f, g are poly-nomials in R [ t ], and that f and g are primitive. In the exact setting this isclearly without loss of generality, since non-zero elements of R ( t ) are units(and hence we can multiply and divide on the left by them). For approximatedifferential polynomials we must compute the approximate GCD of a numberof polynomials in R [ t ]. This is in itself an important research problem, buthas been considered deeply in Kaltofen et al. (2006), which also contains auseful and current survey of related approximate GCD results.4 .2 Norms of differential polynomials To provide our notion of approximate differential polynomials a formal mean-ing we need a proper definition of the norm of a differential polynomial. Forthis we will use the coefficient 2-norm as follows.
Definition 1.3. (i) For (a regular polynomial) p = (cid:80) ≤ i ≤ k p i t i ∈ R [ t ] , define (cid:107) p (cid:107) = (cid:107) p (cid:107) =( (cid:80) ≤ i ≤ k p i ) / .(ii) For f = f + f ∂ + · · · + f m ∂ m ∈ R [ t ][ ∂ ; (cid:48) ] define (cid:107) f (cid:107) = (cid:107) f (cid:107) = (cid:0)(cid:80) ≤ i ≤ m (cid:107) f i (cid:107) (cid:1) / . Note that we are assuming that our coefficients are polynomials in R [ t ],and not rational functions. One could extend the definition of norm to en-compass coefficients in R ( t ), but it will not be necessary in this paper. We can now formally state the main problem under consideration in thispaper.
Problem 1.4.
Given f, g ∈ R [ t ][ ∂ ; (cid:48) ] , find a small ε > and (cid:101) f , (cid:101) g ∈ R [ t ][ ∂ ; (cid:48) ] with (cid:107) f − (cid:101) f (cid:107) < ε (cid:107) f (cid:107) and (cid:107) g − (cid:101) g (cid:107) < ε (cid:107) g (cid:107) such that deg ∂ gcrd( (cid:101) f , (cid:101) g ) ≥ . That is, we are looking for “nearby” differential polynomials, in the co-efficient 2-norm, which possess a non-trivial GCRD.
In this section we demonstrate how to reduce the computation of the GCRDto that of linear algebra over R ( t ), and then over R itself. This approach hasbeen used in the exact computation of GCRDs (see (Li and Nemes, 1997))and Hermite forms (Giesbrecht and Kim, 2013), and has the benefit of re-ducing differential, and more general Ore problems, to a system of equationsover a commutative field. Here we will show that it makes our approximateversion of the GCRD problem amenable to numerical techniques.5 .1 Reduction to linear algebra over R ( t ) Let f, g ∈ R ( t )[ ∂ ; (cid:48) ] have degrees in ∂ of m and n respectively. Then byLemma 1.1, deg ∂ gcrd( f, g ) ≥ u, v ∈ R ( t )[ ∂ ; (cid:48) ]such that deg ∂ u < n and deg ∂ v < m and uf + vg = 0. We can encode theexistence of u, v as an ( m + n ) × ( m + n ) matrix over R ( t ) as follows. Forconvenience define the matrix V = V ( f, g ) = Ψ m + n ( f )Ψ m + n ( ∂f )...Ψ m + n ( ∂ n − f )Ψ m + n ( g )Ψ m + n ( ∂g )...Ψ m + n ( ∂ m − g ) ∈ R ( t ) ( m + n ) × ( m + n ) , the differential Sylvester matrix of f and g (see (Li and Nemes, 1997)),analogous to the Sylvester matrix for usual polynomials (see, e.g., (von zurGathen, 2003, Chapter 6)).The utility of this comes in the following observation. Let u = (cid:88) ≤ i ≤ n − u i ∂ i , v = (cid:88) ≤ i ≤ m − v i ∂ i ∈ R ( t )[ ∂ ; (cid:48) ] , and w = ( u , u , . . . , u n − , v , v , . . . , v m − ) ∈ R ( t ) × ( m + n ) . Then uf + vg = 0 implies wV = 0. This means that w is a non-trivial vectorin the nullspace of V , and in particular, V is singular. Clearing denominatorsof f and g we may assume that u, v ∈ R [ t ][ ∂ ; (cid:48) ], i.e., they have polynomialcoefficients, which implies that V ∈ R [ t ] ( m + n ) × ( m + n ) . Moreover, if f, g ∈ R [ t ][ ∂ ; (cid:48) ] have degrees in t at most d then deg t V ij ≤ d . Lemma 2.1.
Suppose f, g ∈ R [ t ][ ∂ ; (cid:48) ] , where deg ∂ f = m , deg ∂ g = n and deg t f ≤ d and deg t g ≤ d . • V = V ( f, g ) is singular if and only deg ∂ gcrd( f, g ) ≥ . • deg ∂ gcrd( f, g ) = dim null (cid:96) ( V ) , where null (cid:96) ( V ) is the left nullspace of V . For any w = ( u , . . . , u n − , v , . . . , v m − ) ∈ R ( t ) × ( m + n ) such that wV =0 , we have uf + vg = 0 , where u = (cid:80) ≤ i Example 2.2. Let f = ∂ + (0 . t + 1 . ∂ + 0 . t + 0 . t + 0 . and g = ∂ + (cid:0) . t + 1 . . t (cid:1) ∂ + 0 . . t + 0 . t . The corresponding differential Sylvester matrix V is given by . t + 0 . t + 0 . . t + 1 . . . t . . t + 0 . t . t + 1 . . . t + 0 . t . t + 1 . . t . t + 0 . t . t + 0 . t + 1 . t + 0 . . t + 1 . . t .V has rank with the (left) null space vector − . t + 9 . t + 60 . t − . − . t + 10 . t . t − . t − . t + 10 . . t − . t T . Definition 2.3. For any matrix V ∈ R ( t )[ ∂ ; (cid:48) ] , we define the Frobenius norm (cid:107) V (cid:107) F by (cid:107) V (cid:107) F = (cid:88) ij (cid:107) V ij (cid:107) . Lemma 2.4. Let f, g ∈ R [ t ][ ∂ ; (cid:48) ] have deg ∂ f ≤ m , deg ∂ g ≤ n , and both havedegree in t at most d . Let V = V ( f, g ) be the differential Sylvester matrix of f and g . Then (cid:107) V (cid:107) F ≤ d n (cid:107) f (cid:107) + d m (cid:107) g (cid:107) . roof. First note that (cid:107) ∂f (cid:107) ≤ ( d + 1) (cid:107) f (cid:107) , and hence (cid:107) ∂ k f (cid:107) ≤ ( d +1) k (cid:107) f (cid:107) . Thus (cid:107) V (cid:107) F = (cid:88) ≤ i This follows directly from the definition of Γ and Lemma 2.1. Example 2.6. Consider f = (0 . t + 0 . ∂ +0 . t +0 . and g = 0 . ∂ +0 . t . Then the matrix (cid:98) V ( f, g ) is given by . 42 0 . 11 0 0 0 0 0 . 45 0 . 84 0 0 0 00 0 . 42 0 . 11 0 0 0 0 0 . 45 0 . 84 0 0 00 0 0 . 42 0 . 11 0 0 0 0 0 . 45 0 . 84 0 00 0 0 0 . 42 0 . 11 0 0 0 0 0 . 45 0 . 84 00 0 0 0 0 . 42 0 . 11 0 0 0 0 0 . 45 0 . 840 0 . 92 0 0 0 0 0 . 66 0 0 0 0 00 0 0 . 92 0 0 0 0 0 . 66 0 0 0 00 0 0 0 . 92 0 0 0 0 0 . 66 0 0 00 0 0 0 0 . 92 0 0 0 0 0 . 66 0 00 0 0 0 0 0 . 92 0 0 0 0 0 . 66 0 We can now bound the norm of the inflated differential Sylvester matrix. Lemma 2.7. Let f, g ∈ R [ t ][ ∂ ; (cid:48) ] have deg ∂ f ≤ m , deg ∂ g ≤ n and bothhave degree degree at most d in t . Let (cid:98) V ∈ R ( m + n )( µ +1) × ( m + n )( µ + d +1) be theinflated differential Sylvester matrix of f and g , where µ = 2( m + n ) d . Then (cid:107) (cid:98) V (cid:107) ≤ µ · ( d n (cid:107) f (cid:107) + d m (cid:107) g (cid:107) ) .Proof. Each row of (cid:98) V consists precisely of entries of V = V ( f, g ), shifted inposition with respect to the previous row. Thus (cid:107) (cid:98) V (cid:107) F ≤ µ · (cid:107) V (cid:107) F ≤ µ · ( d n (cid:107) f (cid:107) + d m (cid:107) g (cid:107) ) , and (cid:107) (cid:98) V (cid:107) ≤ (cid:107) (cid:98) V (cid:107) F . See (Golub and van Loan, 2013, § Computing an approximate GCRD We have now formulated the problem of determining the existence of GCRD’sof differential polynomials in R [ t ][ ∂ ; (cid:48) ] as one of computing left null vectorsof the inflated differential Sylvester matrix over R . We can now adapt theapproach of Corless et al. (1995) of using the SVD to find the nearest singularmatrix. While this will not be perfect, in that the nearest singular matrixwill not generally have the same structure as the inflated differential Sylvestermatrix, if our input differential polynomials are “nearby” polynomials witha non-trivial GCRD we will generally recover them.For convenience we will generally assume throughout this section that ourinput differential polynomials are normalized, that is have coefficient 2-norm1 under the definition of Section 1.2. This can, of course, be enforced by asimple a priori renormalization, i.e., dividing through by the actual norm,and does not affect the generality or quality of the results. It is well understood how to find the nearest singular unstructured matrix toa given matrix via the singular value decomposition (SVD); see (Golub andvan Loan, 2013, § V ∈ R [ t ] N × N is thedifferential Sylvester matrix from Subsection 2.1, of differential polynomials f, g ∈ R [ t ][ ∂ ; (cid:48) ] of degrees (in ∂ ) of m and n respectively, with N = m + n . From this we construct the inflated differential Sylvester matrix (cid:98) V ∈ R N ( µ +1) × N ( µ + d +1) as in Subsection 2.2. Using the SVD we can find the matrix∆ (cid:98) V of minimal 2-norm such that (cid:98) V + ∆ (cid:98) V has a prescribed rank. First, wecompute the SVD of (cid:98) V as (cid:98) V = P Σ Q, where P ∈ R N ( µ +1) × N ( µ +1) , and Q ∈ R N ( µ + d +1) × N ( µ + d +1) are orthogonal andΣ = diag( σ , . . . , σ N ( µ +1) ) ∈ R N ( µ +1) × N ( µ + d +1) , satisfies σ ≥ σ . . . ≥ σ ≥ σ N ( µ +1) . Note that Σ is not square (it has morecolumns than rows), and we simply pad it with zeros to obtain the desiredshape. 10ow, by Lemma 2.5, we want to find a nearby matrix whose left nullspacehas reduced dimension by multiples of ( µ + d + 1), that is P Σ Q = (cid:98) V + ∆ (cid:98) V , where Σ = diag( σ , σ , . . . , σ ( N − (cid:37) )( µ + d +1) , , . . . , ∈ R N ( µ +1) × N ( µ + d +1) , where (cid:37) = dim null (cid:96) ( (cid:98) V +∆ (cid:98) V ) µ + d +1 . Then (cid:98) V will be by the singular matrix, (cid:98) V + ∆ (cid:98) V of prescribed rank. Of course, (cid:98) V + ∆ (cid:98) V is probably an unstructured matrix,and in particular, not an inflated differential Sylvester matrix.Next we show that a matrix of the desired rank deficiency and (inflateddifferential) structure exists within a relatively small radius of (cid:98) V . Supposethere is an (cid:101) f , (cid:101) g ∈ R [ t ][ ∂ ; (cid:48) ], with (cid:107) (cid:101) f − f (cid:107) ≤ ε and (cid:107) (cid:101) g − g (cid:107) ≤ ε , such thatdeg ∂ gcrd( (cid:101) f , (cid:101) g ) = (cid:37) ≥ 1. Let ∆ f = f − (cid:101) f and ∆ g = g − (cid:101) g , so (cid:107) ∆ f (cid:107) , (cid:107) ∆ g (cid:107) < ε .Moreover, the differential resultant matrix W ∈ R [ t ][ ∂ ; (cid:48) ] N × N formed from∆ f , and ∆ g has (cid:107) W (cid:107) F < ( d n + d m ) · ε by Lemma 2.4. Thus, the inflated differential resultant matrix (cid:99) W ∈ R N ( µ +1) × N ( µ + d +1) has (cid:107) (cid:99) W (cid:107) < µ · ( d n + d m ) · ε by Lemma 2.7. Moreover, dim null (cid:96) ( (cid:98) V + (cid:99) W ) = (cid:37) ( µ + d + 1). Thus, forsufficiently small ε there exists a perturbation (cid:99) W such that (cid:107) (cid:99) W (cid:107) is small (atleast assuming small n, m ) and (cid:98) V + (cid:99) W has appropriate rank and structure.Due to the unstructured nature of (cid:98) V + ∆ (cid:98) V , one must take care in workingwith a reasonable approximation for f and g . However, there is considerableredundancy of the coefficients of f and g in their inflated differential Sylvestermatrix, if only because each entry of f and g shows up multiple times underthe map Γ; see Section 2.2. There is, in fact, even more redundancy becauseof the different derivatives in rows of the differential Sylvester matrix, butwe will not capitalize on this. 11 .2 Computing V + ∆ V and reconstructing f + ∆ f and g + ∆ g To form V +∆ V we take a weighted average of the t -shifted blocks of (cid:98) V +∆ (cid:98) V that correspond to f and g . This involves identifying the “blocks” of (cid:98) V + ∆ (cid:98) V that correspond to f and g and re-constructing them entrywise ensuring theentries in degree t of f and g do not increase. The f block consists of rows1 through µ + 1 and the g block consists of rows deg ∂ g ( µ + 1) + 1 throughdeg ∂ g ( µ + 1) + µ + 2. The columns in both cases are the µ + d + 1 columnsfor each block entry.The reason that this reconstruction is often satisfactory is we have thatif (cid:37) > (cid:98) V has rank ( N − (cid:37) )( µ + d + 1), and if σ ( N − (cid:37) )( µ + d +1)+1 < ε , then (cid:107) Σ − Σ (cid:107) F = ( µ +1) N (cid:88) i =( N − (cid:37) )( µ + d +1)+1 | σ i | ≤ ε [( µ + 1) N − ( N − (cid:37) )( µ + d + 1)] ≤ ε (cid:37) ( µ + d + 1) . We have that (cid:107) (cid:99) W (cid:107) F = (cid:107) Σ − Σ (cid:107) F (Golub and van Loan, 2013, Corollory2.4.3) because the singular values of (cid:99) W = P (Σ − Σ) Q are a permutation ofthe entries along the main diagonal of Σ − Σ.In our construction we require that deg t (cid:101) f i ≤ deg t f i for 0 ≤ i ≤ deg ∂ f and a similar condition on g in order to preserve the structure of V + ∆ V .Furthermore, if the perturbation from adjusting the singular values is small,then the non-zero entries are “small” and can usually be ignored withoutlosing too much information.We should now have a matrix that is numerically singular, V + ∆ V andperturbations ∆ f and ∆ g such that NumericGCRD ( f + ∆ f, g + ∆ g ) is nontrivial and satisfies the conditions deg t ∆ f i ≤ deg t f i for 0 ≤ i ≤ m anddeg t ∆ g j ≤ deg t g j for 0 ≤ j ≤ n . Let f, g ∈ R [ t ][ ∂ ; (cid:48) ] have degrees m and n respectively. Let G = gcrd( f, g )and deg ∂ ( G ) = D . Then one may obtain an R ( t ) multiple of G by solving wV = (cid:0) ∗ ∗ · · · ∗ D · · · (cid:1) , ∗ . Solving this system will give us amultiple of G , which we may assume is in R [ t ] by clearing fractions from thedenominator. Computing an Approximate Primitive GCRD When computing the GCRD numerically we obtain a result that is an R [ t ]-multiple of a primitive GCRD upon clearing fractions. In some applicationsit is desirable to remove this content. However the coefficients are not exactlyknown so taking an exact GCD of the coefficients will yield unsatisfactoryanswers.Consider the case of our NumericGCRD algorithm, Algorithm 2. Our solu-tion will have approximate content even if we use rational arithmetic because f + ∆ f and g + ∆ g have an approximate GCRD but may not have a GCRDalgebraically due to round off errors in their recovery. If a primitive solutionis desired, then techniques to remove the content are required. Leading Coefficient Known in Advance Since leading coefficients of GCRDs are propagated through multiplication,we often know the leading coefficient of a GCRD in advance. In general,the leading coefficient of the GCRD should be approximately one (i.e., aconstant) save a few special cases where the leading coefficients of f and g have a non-trivial approximate GCD.As an observation, given f, g in R [ t ][ ∂ ; (cid:48) ] where deg ∂ f = m and deg ∂ g = n and gcd( f m , g n ) = 1 then for a suitable approximate GCD algorithm (seeCorless et al. (2004) or Zeng and Dayton (2004)) we have that a primitivenumeric GCRD of f and g satisfies lcoeff( G ) = 1. The reason that weneed the notion of approximate GCD is that it is possible that we may havedifferent algorithms returning different answers. Consider t and t + 2 − s where s is large. Algebraically both polynomials are co-prime but if s issufficiently large then some GCD algorithms will return a non trivial GCD.To justify this, let deg ∂ G = D . G is a GCRD so we have that G divides both f and g on the right. If G is primitive, it follows that G D | f m and G D | g n so G D | gcd( f m , g n ) = 1. This occurs if and only if G D = 1.If we are given a candidate GCRD (cid:101) G that is not primitive, where deg ∂ (cid:101) G = D , and we know the primitive GCRD has leading coefficient 1, then lcoeff( (cid:101) G ) ≈ cont (cid:101) G . It follows that (cid:101) G D | (cid:101) G i for 0 ≤ i < D . This means that the remainders13re numerically trivial, so we can assume they are zero. Using a method ofapproximate division we can recover a primitive approximate GCRD.If the primitive GCRD is known to be 1, it is not sufficient to solve thesystem wV = (cid:0) ∗ ∗ . . . ∗ D − . . . (cid:1) T because we are performing numerical linear algebra and we will often obtaina solution over R ( t )[ ∂ ; (cid:48) ]. One particular method of approximate division isby interpolation. This particular method of approximate division yields an-swers that one would expect with prior knowledge of a GCRD and a uniformdistribution of noise. As expected with an interpolation based method, itwill break down if we are unable to accurately compute the degrees of termsdue to artifacts and round off errors. The method of Bini and Pan (1986) forapproximate division via a Fast Fourier Transform (FFT) is used for contentremoval in this case, and proves both fast and numerically robust in practice. Algorithm 1 Content Removal via Approximate Division Input: (cid:101) G ∈ R [ t ][ ∂ ; (cid:48) ] of degree D in ∂ with appropriate degree of entries(leading coefficient has minimal degree). Output: G ≈ (cid:101) G/ cont( (cid:101) G ). Vectorize G as ( G , G , . . . , G D ). for i from 0 to D − do Remove left over artifacts from our linearalgebra below a given threshold. Compute FFT ( (cid:101) G i ) and FFT ( (cid:101) G D ) using a deg t ( (cid:101) G i ) + 1 root of unity. Compute the element wise quotient (cid:101) G i / (cid:101) G D for each entry. Compute G i = InverseFFT ( (cid:101) G i / (cid:101) G D ) Set the last deg ∂ terms to 0 and remove other artifacts from the FFTbelow a given threshold. end for set G D = 1 Devectorize ( G , G , . . . , G D ) into G . return G Example 3.1 (Numeric GCRD) . Consider f = − . ∂ − . t∂ − . t − . and g = ∂ + ( t + 0 . ∂ + (2 . . t ) ∂ + 0 . 66 + 0 . t . numeric GCRD of f and g is given by G = (cid:0) . t + 0 . t − . t − . (cid:1) ∂ +0 . t + 0 . t − . t − . t − . . Given the low degrees and leading coefficients, a primitive GCRD is probablya unit in R . Removing content with an FFT gives us a primitive numericGCRD of ∂ + 1 . t . Leading Coefficient of G Unknown It is not possible to determine the leading coefficient of a primitive GCRDin advance when the leading coefficients of f and g have a non-trivial GCD,or if they share a nearby common solution. In order work around this onewould need to approximate cont( G ) = gcd( G , G , . . . , G D ) numerically, thenperform an approximate polynomial division. We used the method of Corlesset al. (2004) to compute pair-wise GCDs and obtained somewhat mixedresults. In some instances the content had a degree that was too small basedon our construction or our division algorithm did not provide an answerconsistent with the GCRD we constructed. We might hope to overcomesome of these problems using a more specific method for the GCD of multiplepolynomials, as developed in Kaltofen et al. (2005).15 .4 Algorithms This section provides a high-level description of the primary algorithms weare using, in pseudocode.In practice, we will demonstrate that our algorithms work well on low de-gree differential polynomials as input. However, in general our algorithms arenot provably guaranteed to give close polynomials with a non-trivial GCRD.In practice, if (cid:107) f (cid:107) = (cid:107) g (cid:107) = 1 and (cid:107) ∆ f (cid:107) = (cid:107) ∆ g (cid:107) < . f, g ). Algorithm 2 : NumericGCRDInput: • f, g ∈ R [ t ][ ∂ ; (cid:48) ] non-zero with (cid:107) f (cid:107) = (cid:107) g (cid:107) = 1; • A search radius ε > Output: G ≈ gcrd( f, g ) ∈ R [ t ][ ∂ ; (cid:48) ] with deg ∂ G ≥ 1, or an indication that f and g are co-prime within search radius ε . m ← deg ∂ f , n ← deg ∂ g , d ← max { deg t f, deg t g } and µ ← m + n ) d . Form the differential Sylvester matrix V ( f, g ) ∈ R [ t ] ( m + n ) × ( m + n ) . Form the inflated differential Sylvester matrix (cid:98) V = (cid:98) V ( f, g ) ∈ R ( m + n )( µ +1) × ( m + n )( µ + d +1) of V . Compute the numerical rank r of V using Algorithm 4 on (cid:98) V with searchradius ε . If r > D = m + n − r . Otherwise indicate that f and g areco-prime with respect to ε and return. Solve for w from wV = (cid:0) ∗ ∗ . . . ∗ D . . . (cid:1) T . Set G = wV Optionally remove the content from G numerically.As an observation, for f, g ∈ R [ t ][ ∂ ; (cid:48) ] of degrees m and in n in ∂ respec-tively, we can use some heuristics to detect failures of our algorithms. Let G = gcrd( f, g ) and deg ∂ G = D . It is clear that if D > min { n, m } then G cannot be a GCRD. Simmilarily if deg t ( G i ) < deg t ( G D ) for some 0 ≤ i < D then G cannot be a GCRD of f and g over R [ t ] if gcd(lcoeff( f ) , lcoeff( g )) = 1 . lgorithm 3 : Nearest With GCRDInput: • f, g ∈ R [ t ][ ∂ ; (cid:48) ] with (cid:107) f (cid:107) = (cid:107) g (cid:107) = 1; • A search radius ε > Output: f + ∆ f , g + ∆ g where deg t ∆ f ≤ deg t f , deg ∂ ∆ f ≤ deg ∂ f ,deg t ∆ g ≤ deg t g , deg ∂ ∆ g ≤ deg ∂ g and G ≈ gcrd( f + ∆ f, g + ∆ g ) ∈ R [ t ][ ∂ ; (cid:48) ] with deg ∂ G ≥ f and g are co-prime within search radius ε . m ← deg ∂ f , n ← deg ∂ g , d ← max { deg t f, deg t g } and µ ← m + n ) d . Form the differential Sylvester matrix V ( f, g ) ∈ R [ t ] ( m + n ) × ( m + n ) . Form the inflated differential Sylvester matrix (cid:98) V = (cid:98) V ( f, g ) ∈ R ( m + n )( µ +1) × ( m + n )( µ + d +1) of V . Compute the SVD of (cid:98) V , (cid:98) V = P Σ Q with P, Σ and Q as discussed in § Compute the numerical rank r of V using Algorithm 4 on (cid:98) V with searchradius ε . If r > r ( µ + d + 1) singular values to 0 and compute Σas discussed in § f and g are co-prime withrespect to ε and return. Compute (cid:98) V + ∆ (cid:98) V = P Σ Q . Recover f + ∆ f and g + ∆ g from (cid:98) V + ∆ (cid:98) V as discussed in § Compute G = NumericGCRD ( f + ∆ f, g + ∆ g ) using Algorithm 2 with ε used to validate the degree of our approximate GCRD. Optionally remove the content from G numerically. Example 3.2 (Nearest With GCRD) . Consider f = (1 . . t ) ∂ + (3 . t − . ∂ + 2 . t + 1 . and g = t ∂ + (cid:0) − . t + t + 0 . (cid:1) ∂ + t with a given search radius ε = . · − . Applying Algorithm 4 on the singularvalues of (cid:98) V we obtain a rank of where the expected rank is , which iswithin reason. From this we conclude that the degree of our GCRD is and lgorithm 4 : Deflated RankInput: • An inflated differential Sylvester matrix (cid:98) V ∈ R ( m + n )( µ +1) × ( m + n )( µ + d +1) with m, n, d and µ defined as in Algorithm 2or Algorithm 3. • A search radius ε > Output: The the numeric rank of the (non-inflated) differential Sylvestermatrix V from Algorithm 2 or 3. Find the maximum k such that σ k > ε √ ( m + n )(2 µ + d +2) µ + d +1 and σ k +1 < ε . if σ k > ε for all k then (cid:98) V has full rank. If there is no significant change between σ k and σ k +1 for all k as deter-mined by step 2 then return failure. Set r = (cid:108) kµ + d +1 (cid:109) , the scaled rank. we compute (cid:101) f =0 . . t + ( − . . t ) ∂ + (0 . . t ) ∂ and (cid:101) g = − . − . t + 1 . t + (cid:0) . − . t + 0 . t + 0 . t (cid:1) ∂ + (cid:0) − . − . t + 0 . t (cid:1) ∂ . Furthermore, we obtain NumericGCRD ( ˆ f , ˆ g ) ≈ − . . t + ∂ ,after removing content. We have that the size of the perturbations are (cid:107) f − (cid:101) f (cid:107) = 1 . × − and (cid:107) g − (cid:101) g (cid:107) = 0 . . In this example the largestsingular value we removed was known to three decimal places so we can expectthe accuracy of our answer to reflect this. As a remark, we constructed (cid:101) f and (cid:101) g from their first occurrence in (cid:98) V for illustration purposes. In general weighted approaches work well and aredemonstrated in Section 5. 18 n the Rank Algorithm Algorithm 4 follows similarly to the method for determination of numericalrank in the SVD GCD of Corless et al. (1995). In our circumstances thesingular values are in (column) “blocks” of size µ + d + 1. We choose todeclare the scaled rank as (cid:100) kµ + d +1 (cid:101) because it will tend to ignore spurioussingular values provided there are fewer spurious singular values than theblock size. Since we know that non-spurious singular values come in theblock size, this allows us to accurately compute the rank of V even if wedon’t correctly compute the rank of (cid:98) V in the case of low degree differentialpolynomials.Another important reason for this decision is if we set a singular value to0 that is reasonably far away from the next singular value then the resultingmatrix (cid:98) V + ∆ (cid:98) V becomes highly unstructured and our algorithm fails toproduce meaningful results.It is possible that the singular values are not clearly separated. Thishappens frequently with dense large degree differential polynomials. In thisscenario it is possible to under estimate the rank of V . If f and g are knownto have a GCRD, then it may be possible to eliminate bad GCRD candidatesusing heuristics. In this section we investigate the computational complexity in terms of theinput size and some of the numerical stability of the algorithms. We assumeas usual that for f, g ∈ R [ t ][ ∂ ; (cid:48) ] that deg ∂ f = m , deg ∂ g = n , V = V ( f, g ) ∈ R [ t ] ( m + n ) × ( m + n ) and (cid:98) V ∈ R ( m + n )( µ +1) × ( m + n )( µ + d +1) with deg V = d . Weassume that all operations over R can be performed in a constant amountof time, i.e., are floating point operations or flops . Our algorithms will beanalyzed in a bottom-up approach, reflecting their dependencies. Analysis of Algorithm 4: Deflated Rank The cost of the algorithm is dependent on computing the singular values of (cid:98) V .The cost of the computing the singular values of (cid:98) V is O (( m + n ) d ) operationsover R or flops, using the standard method of matrix multiplication.19 nalysis of Algorithm 2: NumericGCRD The cost of the algorithm is dominated by that of computing the singularvalues of (cid:98) V to determine the rank of V . The cost of performing the linearalgebra on V is O ( m + n ) d ) operations over R . The cost of the SVD on (cid:98) V is O (( m + n ) d ) operations and hence dominates.In general the numerical stability of the algorithms depends on the meth-ods used to solve the problem and the conditioning of V . It is difficult tosay much else without making further assumptions. We hope to explore thisfurther in subsequent work. Analysis of Algorithm 3: Nearest With GCRD The cost of the algorithm is dominated by that of computing the singularvalues of (cid:98) V to determine the rank of V . This requires O ((( m + n ) d ) oper-ations over R or flops, with the usual method of matrix multiplication. Thisalgorithm will call Algorithm 2, but again the cost of computing the singularvalues of (cid:98) V is the dominating cost.The accuracy of our answer is subject to the singular values we set to zeroor remove. In particular, the larger the singular values we remove, the lessaccurate our answer typically becomes. In our experiments this became no-ticeable with singular values around 10 − . Under the assumption that noiseis distributed uniformly over f and g then the singular values of (cid:98) V will gen-erally remain small, which means Algorithm 3 typically produce a reasonablenumeric GCRD. Again, we hope to bolster these heuristic observations withmore careful analysis in subsequent work. In order to verify the robustness of the algorithm and whether it is able tocompute (cid:101) f and (cid:101) g reliably, we performed 100 random trials with differentsearch neighborhoods and adding random amounts of noise. In our experi-ments we are adding a noise factor of size δ and working with a search radiusof size ρ . The goal of our experiments is to demonstrate that given a pair ofrelatively prime polynomials (cid:98) f , (cid:98) g ∈ R [ t ][ ∂ ; (cid:48) ] where (cid:107) f − (cid:98) f (cid:107) = (cid:107) g − (cid:98) g (cid:107) = δ for some f and g such that gcrd( f, g ) is non trivial, we can recover a pairof polynomials such that Algorithm 3 returns a non trivial answer in agiven search radius of size ρ . More precisely, the perturbations (cid:107) (cid:98) f − (cid:101) f (cid:107) (cid:107) (cid:98) g − (cid:101) g (cid:107) are minimal approximations to the 2-norm we are using in ourleast squares setting. We can think of ε from Problem 1.4 as being quantifiedby max {(cid:107) (cid:98) f − (cid:101) f (cid:107) , (cid:107) (cid:98) g − (cid:101) g (cid:107)} in these tests.We perform tests on two sets of examples where the input size wasbounded. We found the algorithm worked quite well in practice on examplestaken uniformly at random, however the lack normalization or structure madeit difficult to obtain comparable data. In the data tables we provide statisticsfor two different reconstructions. The first approach is reconstructing (cid:101) f and (cid:101) g from the first row they appear in (cid:98) V + ∆ (cid:98) V . The other reconstruction ap-proach is using the weighted approach over the entire block to recover (cid:101) f w and (cid:101) g w respectively. On average the weighted approach tends to smooth valuesover at the cost of structure where as taking the first row can preserve theunderlying structure of f and g , especially if they are sparse.We note that as the noise decreases so does the size of the perturbation.We will continue to get smaller perturbations until the roundoff error fromwriting (cid:98) V ≈ P Σ Q dominates the perturbation sizes. If the noise is sufficientlysmall then we are executing NumericGCRD (Algorithm 2) as the perturbationswill become indistinguishable from roundoff errors. Despite the noise in sometests appearing large, it is distributed uniformly into the coefficients resultingin each coefficient being perturbed by a fraction of the total amount of noiseadded. Although the worst case perturbations in (cid:101) f and (cid:101) g were relativelylarge, such perturbations were uncommon in our experiments and we stillmanaged to obtain valid candidates. in the context of Problem 1.4. Wejustify adding noise uniformly since in practice we would expect this workto be applied on data which suffers from round off errors which tend to beuniformly distributed across the data. In this section we perform our tests on f and g whose coefficients are boundedand somewhat structured. The following steps detail the construction of ourexamples.1. Generate h , h , h ∈ R [ t ][ ∂ ; (cid:48) ] where deg t h i ≤ ≤ deg ∂ h i ≤ h ij = h ij / (cid:107) h ij (cid:107) for i = 1 , , ≤ j ≤ deg h i . More precisely, (cid:107) h i (cid:107) = deg t h i + 1. 21. Compute f = h · h + δ f and g = h · h + δ g , where the noise isdistributed uniformly.4. Run Algorithm 3, the Nearest With GCRD algorithm on f and g .5. If Nearest With GCRD ( f, g ) = 1 then ignore the result in the statistics.The construction of these examples seems peculiar and counterintuitive,but it provides a suitable set of tests. The justification for the set of testsis to get an idea how the algorithm performs on random data that is notnormalized, but is bounded with a random structure.The table in Figure 1 details the relevant statistics obtained from runningAlgorithm 3. The table in Figure 2 provides the number of trivial GCRDsthat occurred for a given ρ and δ .Figure 1: Perturbation statistics for bounded coefficients( ρ, δ, Reconstructed ) Max Average Standard Deviation ( . , . , (cid:101) f ) 0.269443 0.022653 0.040056( . , . , (cid:101) g ) 0.233732 0.025051 0.037960( . , . , (cid:101) f w ) 0.269443 0.023507 0.038163( . , . , (cid:101) g w ) 0.233732 0.023625 0.032776( . , . , (cid:101) f ) 0.047198 0.010184 0.011120( . , . , (cid:101) g ) 0.061742 0.011490 0.013434( . , . , (cid:101) f w ) 0.045457 0.009629 0.010101( . , . , (cid:101) g w ) 0.049726 0.010082 0.010634( . , . , (cid:101) f ) 0.233401 0.006653 0.030810( . , . , (cid:101) g ) 0.166854 0.005394 0.021614( . , . , (cid:101) f w ) 0.233401 0.005647 0.027110( . , . , (cid:101) g w ) 0.166854 0.005647 0.018869( . , . , (cid:101) f ) 0.283611 0.006625 0.037541( . , . , (cid:101) g ) 0.182687 0.004582 0.025346( . , . , (cid:101) f w ) 0.283611 0.006277 0.037003( . , . , (cid:101) g w ) 0.182687 0.004317 0.02493322igure 2: Trivial Numeric GCRDs for Bounded Coefficients ρ δ Trivial GCRD.5 .5 9.5 .1 3.5 .01 3.5 .001 0.05 .5 95.05 .1 42.05 .01 0.05 .001 2 In this section we perform our tests on f and g that were normalized andadded noise. The following steps detail the construction of our examples.1. Generate h , h , h ∈ R [ t ][ ∂ ; (cid:48) ] where deg t h i ≤ ≤ deg ∂ h i ≤ f = h · h (cid:107) h · h (cid:107) + δ f and g = h · h (cid:107) h · h (cid:107) + δ g , where the noise isdistributed uniformly.3. Run Algorithm 3, the Nearest With GCRD algorithm on f and g .4. If Nearest With GCRD ( f, g ) = 1 then ignore the result in the statistics.The table in Figure 3 details the relevant statistics obtained from runningAlgorithm 3. The table in Figure 4 provides the number of trivial GCRDsthat occurred for a given ρ and δ . 23igure 3: Perturbation statistics for normalized f and g ( ρ, δ, Reconstructed ) Max Average Standard Deviation ( . , . , (cid:101) f ) 0.108768 0.015880 0.022180( . , . , (cid:101) g ) 0.158056 0.016359 0.024371( . , . , (cid:101) f w ) 0.052060 0.010584 0.010404( . , . , (cid:101) g w ) 0.097296 0.011577 0.014039( . , . , (cid:101) f ) 0.058105 0.017043 0.014834( . , . , (cid:101) g ) 0.063637 0.025312 0.017561( . , . , (cid:101) f w ) 0.057251 0.012524 0.013178( . , . , (cid:101) g w ) 0.046134 0.019487 0.012614( . , . , (cid:101) f ) 0.053014 0.005987 0.009114( . , . , (cid:101) g ) 0.057982 0.006851 0.009757( . , . , (cid:101) f w ) 0.031045 0.004263 0.006710( . , . , (cid:101) g w ) 0.038133 0.005996 0.007523( . , . , (cid:101) f ) 0.066271 0.003778 0.009018( . , . , (cid:101) g ) 0.093175 0.004416 0.011236( . , . , (cid:101) f w ) 0.067082 0.003556 0.009915( . , . , (cid:101) g w ) 0.038098 0.003213 0.006429Figure 4: Trivial Numeric GCRDs for Normalized f and gρ δ Trivial GCRD.5 .5 6.5 .1 0.5 .01 2.5 .001 1.05 .5 86.05 .1 22.05 .01 2.05 .001 124 Conclusions and Future Work We have developed a framework for approximate differential polynomials anddemonstrated an algorithm for computing the greatest common right divisorof two approximate differential polynomials. This corresponds to finding arepresentation of the common solutions of two approximate linear differentialoperators.Our algorithm makes use of the (unstructured) SVD approach introducedin Corless et al. (1995), and works well when differential polynomials with anon-trivial GCRD are nearby. Unfortunately, it also suffers some of the samedrawbacks, especially when the nearest polynomials with a non-trivial GCRDare relatively far away. Even more so, the lack of a geometric root space (asfor conventional polynomials) makes the analysis substantially more difficult.One possible remedy we are exploring is to this is to use a structured matrixapproach, as explored in Beckermann and Labahn (1998), Botting et al.(2005), Kaltofen et al. (2006). In particular, Riemannian SVD-like methodswould seem relatively easy to adapt in this case (though again, analysis willnot be easy).Another problem is the factorial-like scaling introduced by multiple dif-ferentiations when constructing the (inflated) differential Sylvester matrix.This leads to numerical instability when the degree in ∂ gets even modestlylarge. Some row scaling may alleviate this somewhat, but an alternative dif-ferential resultant formulation would seem a better approach overall, and isa path we are investigating.The algorithms described in this paper have been implemented in Mapleand are primarily based on the LinearAlgebra package. This allows flexibilityin determining the method used to solve the linear systems and dealing withother problems such as content removal and numerical stability.In general the approximate GCRD algorithm performs quite well for lowdegree GCRDs and low noise. If the degree of the GCRD increases relativeto the degrees of f and g then (cid:98) V + ∆ (cid:98) V becomes more unstructured andthe perturbations become larger which leads to unsuitable numeric GCRDs.In the case of degree one or two GCRDs we were still able to reconstructmeaningful answers even if the noise was outside of our search radius becausethe matrix still retained a Sylvester-like structure.Despite the speculative nature of our algorithm we find that it worksreasonably well on the test data when we have a priori bound on the noise andchoose a search radius accordingly. Approximating GCRDs of higher order25ecomes more difficult with the unstructured approach as the perturbationsin the data will become larger. If other assumptions are made about thestructure of f and g and the distribution of noise then it would be possible toobtain a better approximate GCRD by exploiting this underlying structure. References B. Beckermann and G. Labahn. When are two numerical polynomials rela-tively prime? Journal of Symbolic Computation , 26:677–689, 1998.D. Bini and V. Pan. Polynomial division and its computational complexity. Journal of Complexity , 2(3):179–203, 1986.B. Botting, M. Giesbrecht, and J.P. May. Using the Riemannian SVD forproblems in approximate algebra. In Proc. Workshop on Symbolic-NumericComputation (SNC’05) , pages 209–219, 2005.M. Bronstein and M. Petkovˇsek. On Ore rings, linear operators and factori-sation. Programmirovanie , 20:27–45, 1994.R. M. Corless, P. M. Gianni, B. M. Trager, and S. M. Watt. The singu-lar value decomposition for polynomial systems. In Proc. InternationalSymposium on Symbolic and Algebraic Computation (ISSAC’95) , pages189–205, 1995.R. M. Corless, S. M. Watt, and L. Zhi. QR factoring to compute the GCDof univariate approximate polynomials. IEEE Transactions on Signal Pro-cessing , 52(12), 2004.J. von zur Gathen, J. Gerhard. Modern Computer Algebra . CambridgeUniversity Press, New York, NY, USA, 2 edition, 2003.M. Giesbrecht and M. Kim. Computing the Hermite form of a matrix of Orepolynomials. Journal of Algebra , 376:341–362, 2013.G. Golub and C. van Loan. Matrix Computations . Johns Hopkins UniversityPress, Baltimore, USA, 4th edition, 2013.E. Kaltofen, Z. Yang, and L. Zhi. Structured low rank approximation of aSylvester matrix. In Proc. Workshop on Symbolic-Numeric Computation(SNC’05) , pages 69–83, 2005. 26. Kaltofen, Z. Yang, and L. Zhi. Approximate greatest common divisorsof several polynomials with linearly constrained coefficients and singularpolynomials. In Proc. International Symposium on Symbolic and AlgebraicComputation (ISSAC’06) , pages 169–176, 2006.N. Karmarkar and Y. N. Lakshman. Approximate polynomial greatest com-mon divisors and nearest singular polynomials. In Proc. InternationalSymposium on Symbolic and Algebraic Computation (ISSAC’96) , pages35–39, 1996.Z. Li and I. Nemes. A modular algorithm for computing greatest commonright divisors of Ore polynomials. In Proc. International Symposium onSymbolic and Algebraic Computation (ISSAC’97) , pages 282–289, 1997.O. Ore. Theory of non-commutative polynomials. Annals of Mathematics.Second Series , 34:480–508, 1933.T. Sasaki and M. Sasaki. Polynomial remainder sequence and approximateGCD. ACM SIGSAM Bulletin , 31:4–10, 1997.Z. Zeng and B. H. Dayton. The approximate GCD of inexact polynomials.In