Framework of algorithm

In this paper, we discuss non-singular case only, i.e., 𝐹 (𝑥, 0) • 𝐺(𝑥, 0) ̸ = 0 and 𝑓 𝑚 (0) • 𝑔 𝑛 (0) ̸ = 0. For non-singular case, every polynomial 𝑃 (𝑥, 𝑡) is transform to 𝑃 (𝑥, 𝑇, 𝑡) = 𝑃 (𝑥, 𝑇 𝑡 1 , . . . , 𝑇 𝑡 ℓ ), with 𝑇 is the total-degree variable. Every polynomial 𝑃 (𝑥, 𝑇, 𝑡) is represented as the sum of homogeneous polynomials w.r.t. the total-degree variable 𝑇 ; 𝑃 (𝑥, 𝑇, 𝑡) = 𝑃 (0) (𝑥) + 𝑇 • 𝛿𝑃 (1) (𝑥, 𝑡) + . . . + 𝑇 𝑤 • 𝛿𝑃 (𝑤) (𝑥, 𝑡) + . . . , 𝑃 (𝑤) (𝑥, 𝑇, 𝑡) = 𝑃 (0) (𝑥) + 𝑇 • 𝛿𝑃 (1) (𝑥, 𝑡) + . . . + 𝑇 𝑤 • 𝛿𝑃 (𝑤) (𝑥, 𝑡).

In non-singular case, the following two conditions exist: deg(gcd(𝐹, 𝐺)) ≤ deg(gcd(𝐹 (0) , 𝐺 (0) )) and gcd(𝐹, 𝐺)|gcd(𝐹 (0) , 𝐺 (0) ). In such situations, lifting algorithms can be applied. The proposed algorithm is discussed in the next section.

Computing cofactors of 𝐹 and 𝐺 via lifting method

Let 𝒮 𝑖 (𝐹, 𝐺) = 𝒮 𝑖 ∈ F[𝑡, 𝑇 ] (𝑚+𝑛−2𝑖)×(𝑚+𝑛−2𝑖) be an 𝑖th-subresultant matrix of 𝐹 and 𝐺 w.r.t. 𝑥, and be represented as

𝒮 𝑖 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝑓 𝑚 𝑔 𝑛 . . . . . . . . . . . . 𝑓 𝑚−𝑛−𝑘 𝑓 𝑚 𝑔 𝑛−𝑚−𝑘 𝑔 𝑛 . . . . . . . . . . . . . . . . . . 𝑓 𝑚−𝑛−𝑘 𝑔 𝑛−𝑚−𝑘 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = 𝒮 (0) 𝑖 + 𝑇 • 𝛿𝒮(1)

computation of cofactors for univariate part: SVD

Cofactors of 𝐹 (0) and 𝐺 (0) can be obtained from the null space of 𝒮 (0) 𝑘−1 . In this paper, we compute the null space of 𝒮 (0) 𝑘−1 using the SVD [1]. Using the SVD of 𝒮 (0) 𝑘−1 , we obtain the following decomposition:

𝒮 (0) 𝑘−1 = 𝑈 Σ𝑉 𝑇 = (︀ 𝑢 1 • • • 𝑢 𝐾 )︀ ⎛ ⎜ ⎝ 𝜎 1 . . . 𝜎 𝐾 ⎞ ⎟ ⎠ ⎛ ⎜ ⎝ 𝑣 𝑇 1 . . . 𝑣 𝑇 𝐾 ⎞ ⎟ ⎠ ,

where 𝐾 = 𝑚 − (𝑘 − 1) + 𝑛 − (𝑘 − 1), 𝑈 and 𝑉 are orthogonal matrices, and 𝜎 𝑖 are singular vectors with

𝜎 1 ≥ 𝜎 2 ≥ . . . ≥ 𝜎 𝐾−1 ≫ 𝜎 𝐾 ≥ 0, respectively. Then, 𝑣 𝐾 ∈ Ker(𝒮(0)

𝑘−1 ), and it is one of the solutions of 𝒮 (0

𝑘−1 𝑧 = 0 is 𝑧 = 𝑟 (0) = 𝑣 𝐾 ; 𝑣 𝐾 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝑔 ˜(0) 𝑛−𝑘 𝑔 ˜(0) 0 −𝑓 ˜(0) 𝑚−𝑘 −𝑓 ˜(0) 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , with ||𝑣 𝐾 || 2 = 1.

Computation of cofactors for multivariate part: lifting method

Suppose we have 𝑟 (𝑤) = 𝑟 (𝑤−1) +𝛿𝑟 (𝑤) . Here, 𝛿𝑟 (𝑤) is a vector generated by homogenious polynomials with total-degree 𝑤 w.r.t. 𝑇 . Then, 𝛿𝑟 (𝑤+1) is generated as follows. Note that the following consists.

𝒮 𝑘−1 𝑟 (𝑤+1) ≡ 0 (mod 𝑇 𝑤+2 ) 𝒮 (0) 𝑘−1 𝛿𝑟 (𝑤+1) = − 𝑤+1 ∑︁ 𝑗=1 𝒮 (𝑗)

𝑘−1 𝛿𝑟 (𝑤−𝑗) = 𝛿𝑝 (𝑤) . Now, 𝛿𝑟 (𝑤+1) and 𝛿𝑝 (𝑤+1) are transformed bases from 𝑒 1 , . . . , 𝑒 𝐾 to 𝑣 1 , . . . , 𝑣 𝐾 and 𝑢 1 , . . . , 𝑢 𝐾 , respectively, as follows:

𝛿𝑧 (𝑤) = ⎛ ⎜ ⎝ 𝛿𝑧 (𝑤) 1 . . . 𝛿𝑧 (𝑤) 𝐾 ⎞ ⎟ ⎠ = 𝛿𝑧 ˆ(𝑤) 1 𝑣 1 + . . . + 𝛿𝑧 ˆ(𝑤) 𝐾 𝑣 𝐾 , 𝛿𝑝 (𝑤) = ⎛ ⎜ ⎝ 𝛿𝑝 (𝑤) 1 . . . 𝛿𝑝 (𝑤) 𝐾 ⎞ ⎟ ⎠ = 𝛿𝑝 ˆ(𝑤) 1 𝑢 1 + . . . + 𝛿𝑝 ˆ(𝑤) 𝑁 𝑢 𝐾 .

Then, we obtain 𝛿𝑧 ˆ(𝑤+1)

𝑖 = 𝛿𝑝 ˆ(𝑤+1) 𝑖 /𝜎 𝑖 (𝑖 = 1, . . . , 𝐾 − 1)

. Therefore, 𝛿𝑟 (𝑤+1) is constructed, as follows.

𝛿𝑟 (𝑤+1) = 𝛿𝑝 ˆ(1) 1 /𝜎 1 𝑣 1 + . . . + 𝛿𝑝 ˆ(1) 𝐾−1 /𝜎 𝐾−1 𝑣 𝐾−1 + F[𝑇, 𝑡] 𝑤+1 • 𝑣 𝐾 ,

where F[𝑇, 𝑡] 𝑤+1 is homogeneous polynomial set with total-degree 𝑤 + 1 w.r.t. 𝑇 , and we have the following as a candidate solution.

𝑟 (𝑤) = 𝑡 𝐾 + 𝑤 ∑︁ 𝑗=1 𝛿𝑝 ˆ(𝑤) 1 /𝜎 1 𝑣 1 + . . . + 𝑤 ∑︁ 𝑗=1 𝛿𝑝 ˆ(𝑤) 𝐾−1 /𝜎 𝐾−1 𝑣 𝐾−1 + F[𝑇, 𝑡] [1,𝑤] • 𝑣 𝐾 ,

where

F[𝑇, 𝑡] [1,𝑤] = ∪ 𝑤 𝑗=1 F[𝑇, 𝑡] 𝑗 .

To compute the approximate GCD of 𝐹 and 𝐺, we need to determine 𝑞 (𝑤) = 𝛿𝑞 (1)

+ . . . + 𝛿𝑞 (𝑤) ∈ F[𝑇, 𝑡] [1,𝑤] s.t. 𝑟 (𝑤) (𝑞) = 𝑡 𝐾 + 𝑤 ∑︁ 𝑗=1 𝛿𝑝 ˆ(𝑤) 1 /𝜎 1 𝑣 1 + . . . + 𝑤 ∑︁ 𝑗=1 𝛿𝑝 ˆ(𝑤) 𝐾−1 /𝜎 𝐾−1 𝑣 𝐾−1 + 𝑞 (𝑤) • 𝑣 𝐾 ,

To determine the approximate GCD, we must determine 𝑞 (𝑖) . The following example shows one approach to determining each undetermined element 𝛿𝑞 (𝑖) for 1 ≤ 𝑖 ≤ 𝑤..

Example1

Polynomials 𝐹 (𝑥, 𝑡 1 , 𝑡 2 ) and 𝐺(𝑥, 𝑡 1 , 𝑡 2 ) having an approximate GCD 𝐶(𝑥, 𝑡 1 , 𝑡 2

) = 𝑥 3 + (1 + 𝑡 2 − 2𝑡 1 + 𝑡 1 2 )𝑥 + 3 are expressed as 𝐹 (𝑥, 𝑡 1 , 𝑡 2 ) = (𝑥 3 + (𝑡 2 2 + 𝑡 1 + 𝑡 1 − 2)𝑥 2 − 1) × 𝐶(𝑥, 𝑡 1 , 𝑡 2 ) + 𝑀 𝑒𝑝𝑠 , 𝐺(𝑥, 𝑡 1 , 𝑡 2 ) = (𝑥 3 + (2𝑡 2 2 − 𝑡 1 + 3)𝑥 2 − 1) × 𝐶(𝑥, 𝑡 1 , 𝑡 2 ) + 𝑀 𝑒𝑝𝑠 ,

where 𝑀 𝑒𝑝𝑠 is the machine epsilon.

In this example, 𝑘 = 2 is already known (𝐾 = 8). Then, one solution of

𝒮 (0) 1 𝑧 = 0 is 𝑧 = 𝑟 (0) = 𝑣 8 ; 𝑟 (0) = 𝑣 8 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ −0.242535625036333 −0.727606875108999 −2.24840273230668 × 10 −15 0.242535625036333 0.242535625036333 −0.485071250072665 −1.32375311946987 × 10 −15 −0.242535625036333 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . A candidate of 𝛿𝑟 (1) | 𝛿𝑞 (1) =0 = 𝛿𝑝 ˆ(1) 1 /𝜎 1 𝑣 1 + . . . + 𝛿𝑝 ˆ(1)

7 /𝜎 7 𝑣 7 + 𝛿𝑞 (1) × 𝑣 8 is

𝛿𝑧 (1) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0.0713340073 • • • 𝑡 1 + 0.0285336029 • • • 𝑡 2 −0.0285336029 • • • 𝑡 1 + 0.0856008088 • • • 𝑡 2 3.65419500 • • • × 10 −15 𝑡 1 − 4.96824803 • • • × 10 −15 𝑡 2 −0.0713340073 • • • 𝑡 1 − 0.0285336029 • • • 𝑡 2 −0.0713340073 • • • 𝑡 1 − 0.0285336029 • • • 𝑡 2 −0.0998676103 • • • 𝑡 1 − 0.185468419 • • • 𝑡 2 4.77577504 • • • × 10 −16 𝑡 1 − 1.10469359 • • • × 10 −15 𝑡 2 0.0713340073 • • • 𝑡 1 + 0.0285336029 • • • 𝑡 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ + 𝛿𝑞 (1) • 𝑣 8 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝛿𝑔 ˜(1) 𝑛−𝑘 𝛿𝑔 ˜(1) 𝑛−𝑘−1 . . . 𝛿𝑔 ˜(1) 0 −𝛿𝑓 (1) 𝑚−𝑘 −𝛿𝑓 (1) 𝑚−𝑘−1 . . . −𝛿𝑓 (1) 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .

Generally, it is difficult to determine 𝛿𝑞 (1) properly. However, assuming that cofactors are also not dense or the approximate GCD is monic, several coefficients will be zero. In this example, assume the 1st element is zero, 𝛿𝑞 (1) is 𝛿𝑞

(1) 1 = (0.0713340073636269 𝑡 1 + 0.0285336029454512 𝑡 2 )/0.242535625036333 and 𝛿𝑟 (1) becomes

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0.242535625036331𝑡 1 0 0 0 −0.242535625036331𝑡 1 − 0.242535625036334𝑡 2 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

.

It is unlikely that many factors will be close to zero simultaneously, and this can only happen if the result is correct. Unlike the EZ-GCD method, it is more efficient because it can extract each undetermined coefficient at each lifting step. So that, "check zeros" is very efficiency.

If the coefficients are dense, lc(lc(𝐹 ), lc(𝐺)) or lc(lc(gcd(𝐹, 𝐺))) should be calculated in advance so that the elements can be determined uniquely.

Computing approximate GCD

After obtaining cofactors, the approximate GCD is computed by solving and 𝐹

= (𝑓 𝑚 , 𝑓 𝑚−1 , . . . , 𝑓 0 ) 𝑇 ∈ F[𝑡] 𝑚+1 . ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝑓 ˜𝑚−𝑘 . . . . . . 𝑓 ˜0 𝑓 ˜𝑚−𝑘 . . . . . . 𝑓 ˜0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎝ 𝑐 𝑘 𝑐 𝑘−1 . . . 𝑐 0 ⎞ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎝ 𝑓 𝑚 𝑓 𝑚−1 . . . 𝑓 0 ⎞ ⎟ ⎟ ⎟ ⎠

This linear equation is solve as following step.

1. Solve 𝒞

(0) 𝑚+1,𝑘+1 (𝐹 ˜) • 𝑐 (0) = 𝐹 (0)

. Actually, we utilize the SVD as in the former case. 2. Lifting step: solve

𝒞 (0) 𝑚+1,𝑘+1 (𝐹 ˜) • 𝛿𝑐 (𝑤) = 𝛿𝐹 (𝑤) − ∑︀ 𝑖 𝐶 (𝑖)

𝑚+1,𝑘+1 (𝐹 ˜) × 𝛿𝑐 (𝑤−𝑖) . This step can also be solved using SVD. 3. Return 𝑐 𝑘 𝑥 𝑘 + • • • + 𝑐 0 as an approximate GCD.

Solve in ill-conditioned cases

In this section, we demonstrate that our method is stable for ill-conditioned cases [8,5]. On the other word, we deals with cases where the initial factor is an approximate common factor. In this case, the EZ-GCD method is unstable since large cancellation errors occur [6].

Example 2 (initial factors have approximate common factor)

Compute the approximate GCD of 𝐹 and 𝐺, where both polynomials are monic.

𝐹 (𝑥, 𝑡 1 , 𝑡 2 ) = (𝑥 3 + (𝑡 2 2 + 𝑡 1 + 𝑡 2 − 2)𝑥 2 − 1)(𝑥 − 1.0003 + 2𝑡 2 − 𝑡 1 2 )𝐶 + 𝑀 𝑒𝑝𝑠 , 𝐺(𝑥, 𝑡 1 , 𝑡 2 ) = (𝑥 3 + (2𝑡 2 2 − 𝑡 1 + 3)𝑥 2 − 1)(𝑥 − 1.0005 + 𝑡 1 + 𝑡 2 + 𝑡 1 𝑡 2 )𝐶 + 𝑀 𝑒𝑝𝑠 , 𝐶(𝑥, 𝑡 1 , 𝑡 2 ) = 𝑥 3 + (1 + 𝑡 2 − 2𝑡 1 + 𝑡 1 2 )𝑥 + 3.

Initial factors 𝐹 (0) and 𝐺 (0) have an approximate common factor (𝑥 − 1.0002) with tolerance 𝑂(10 −5 ). Sigular values of 𝒮

(0) 3 (𝐹, 𝐺) are 19.8 > 18.3 > 14.5 > 12.8 > 8.2 > 4.4 > 1.1 > 0.6 ≫ 1.5×10 −5 ≫ 1.1×10 −16 .

Because give polynomials are monic, the leading coefficient of cofactors and the approximate GCD are also monic, respectively.

When 𝑤 = 1, Adjusting the 1st element of 𝛿𝑟 (1) by 𝑧 𝐾 only, we obtained the following.

𝛿𝑟 (1) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0. −7.25814180424500 × 10 −8 𝑡 1 + 0.176741414005228𝑡 2 0.707053518826616𝑡 1 + 0.530224242015704𝑡 2 −6.96526170074208 × 10 −14 𝑡 1 + 6.29774010718620 × 10 −14 𝑡 2 −0.176741268890038𝑡 1 − 0.176741414005196𝑡 2 1.42941214420489 × 10 −15 𝑡 1 + 6.38378239159465 × 10 −16 𝑡 2 −0.176741268889989𝑡 1 − 0.530224096948041𝑡 2 0.176794218725546𝑡 1 + 0.883759874841642𝑡 2 −1.18932641512970 × 10 −14 𝑡 1 − 7.57727214306669 × 10 −15 𝑡 2 −7.25814328778052 × 10 −8 𝑡 1 + 0.353482755476628𝑡 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

The perturbation is ||𝐹 ˜(1) 𝐺 (1) − 𝐹 ˜(1) 𝐺 (1) || ≈ 𝑂(10 −8 ) ≈ 𝜎 𝐾 /𝜎 𝐾−1 . On the other hand, by adjusting 𝑧 𝐾−1 , we obtained the following, It can be confirmed that the solution is not accurate; ||𝐹 ˜(1) 𝐺 (1) − 𝐹 ˜(1) 𝐺 (1) || ≈ 𝜎 𝐾 .

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0. −0.1988511825793107𝑡 1 − 0.14357803747786974𝑡 2 0.11055165805122452𝑡 1 − 0.43065120320683287𝑡 2 0.0002976768351589665𝑡 1 + 0.00047951294108034004𝑡 2 0.022318043342197766𝑡 1 + 0.14391342019466513𝑡 2 −0.7183155936049679 × 10 −5 𝑡 1 − 0.000011570991832604571𝑡 2 0.02212670015656562𝑡 1 − 0.20987748805445672𝑡 2 −0.22096310056080598𝑡 1 + 0.243032215144183𝑡 2 0.00007010876634662433𝑡 1 + 0.00011293475597320968𝑡 2 −0.19880976332099576𝑡 1 + 0.03323002423516758𝑡 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

When 𝑤 = 2, by lifting step and adjusting the 1st element of 𝛿𝑟 (2) by 𝑧 𝐾 , we obtained the following vector. 𝑘 we have the following, and it obtains the expected solution one . In this case, perturbation becomes ||𝐹 ˜(2) 𝐺 (2) − 𝐹 ˜(2) 𝐺 (2) || ≈ 𝜎 𝐾 , it is better.

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0. 0.353120013467623𝑡 2 2 + 0.177466845429488𝑡 1 𝑡 2 − 0.000362979823316206𝑡 1 2 −0.354747432749536𝑡 2 2 + 0.355659122269873𝑡 1 𝑡 2 − 0.177830208344029𝑡 1 2 8.84819995 • • • × 10 −13 𝑡 2 2 − 4.59634934 • • • × 10 −12 𝑡 1 𝑡 2 + 1.03052288 • • • × 10 −12 𝑡 1 2 0.000362669479650586𝑡 2 2 − 0.177466845422696𝑡 1 𝑡 2 + 0.000362979818887450𝑡 1 2 −9.08967345 • • • × 10 −13 𝑡 2 2 − 3.63698827 • • • × 10 −12 𝑡 1 𝑡 2 − 1.36436695 • • • × 10 −

The SVD is stable even if the matrix is irregular. Thus, the SVD of 𝒮 (0) is also stable even if initial factors have an approximate common factor. On the other hand, a lifting method using the Bezout matrix is unstable since initial matrix is assumed to be regular [4]. Hence, our method is more stable and efficient than existing methods.