§1. UFD  Lecture II  Page 1

Lecture II The GCD

Next to the four arithmetic operations, the greatest common denominator (GCD) is perhaps the most basic operation in algebraic computing. The proper setting for discussing GCD’s is in a unique factorization domain (UFD). For most common UFDs, that venerable algorithm of Euclid is available. In the domains and F [X], an efficient method for implementing Euclid’s algorithm is available. It is the so-called half-GCD approach, originating in ideas of Lehmer, Knuth and Schönhage.

The student will not go far amiss if she interprets all references to rings as either the integers or a polynomial ring F [X] over a field F (even taking F = ).

§1. Unique Factorization Domain

Let D be a commutative ring with unity 1 where 0 ≠ 1. For a, b ∈ D, we say b divides a, and write b|a, if there is a c ∈ D such that a = bc. We also call b a divisor of a, and a a multiple of b. Thus every element divides 0 but 0 does not divide any non-zero element. A zero-divisor is an element b such that bc = 0 for some non-zero c. We also call an element regular if it is not a zero-divisor. A domain D is a commutative ring with unity whose only zero-divisor is 0. A unit is an element that divides 1. (Alternatively, units are the invertible elements.) Thus the unity element 1 is always a unit and the zero element is never a unit. In a field, all non-zero elements are units. Two elements, a and b, are associates if a = ub for some unit u. Clearly the relation of being associates is an equivalence relation. So the elements of D are partitioned into equivalence classes of associates.

Exercise 1.1:

a) The set of units and the set of zero-divisors are disjoint.

b) a|b and b|a iff a, b are associates.

Convention. For each equivalence class of associates, we assume that a distinguished member is chosen. The unity element 1 is always distinguished. We follow the standard convention for our two main rings:

(i) In , the units are +1 and −1 and the equivalence classes are {−n, +n} for each n ∈ . The non-negative elements will be distinguished in . (ii) In the polynomial ring F [X] over a field F, the units are all non-zero constants and the distinguished elements are the monic elements, i.e., polynomials with leading coefficient 1.

A proper divisor of b is any divisor that is neither a unit nor an associate of b. An element is irreducible if it has no proper divisors; otherwise it is reducible. Since any divisor of a unit is a unit, it follows that units are irreducible. Furthermore, the zero element is irreducible if and only if D is a domain.

Definition: A unique factorization domain (abbreviated, UFD) is a domain D in which every non-unit b can be written as a product of irreducible non-units,

b = b_{1} b_{2} ⋯ b_{n} (n ≥ 1).

Moreover, these irreducible elements are unique up to reordering and associates. UFD’s are also called factorial domains.

The importance of UFD’s is that its elements are made up of “fundamental building blocks”, namely the irreducible elements. Note that is a UFD, by the Fundamental Theorem of arithmetic. In fact, a UFD can be said to be a domain that has an analogue to the Fundamental Theorem of arithmetic! The non-zero irreducible elements of are called primes. But in general, we define the prime elements of a ring R to be those non-units p ∈ R such that p ≠ 0 and if p divides any product a · b then p divides either a or b.

One sees that prime elements are irreducible but the converse is not generally true. For example (see [3, page 173]), in [X, Y, Z]/(Z^{2} − XY), Z is irreducible but not prime because Z divides XY without dividing X or Y. It is easy to see that in a UFD, every irreducible element is also a prime. Hence this is an example of a non-UFD.

Theorem 1 D is a UFD iff D[X] is a UFD.

It is clear that D[X] is not a UFD if D is not a UFD. The proof of the other direction is due to Gauss and is deferred to the next lecture. Trivially, a field F is a UFD. Hence, by induction on d ≥ 1, this theorem shows that F [X_{1},…, X_{d}] is a UFD. This is an important fact we will assume.

Greatest common divisor. Let D be a UFD and S ⊆ D be a finite non-empty set. We write a|S (read, a divides S) to mean a|b for all b ∈ S. An element d ∈ D is a greatest common divisor (abbreviated, GCD) of S if

1) d|S,

2) if c|S then c|d.

Exercise 1.2: Prove that S has a greatest common divisor, and this is determined up to associates.

We can therefore define the function GCD(S) by choosing the distinguished greatest common divisor of S. If S = {a_{1}, a_{2},…, a_{m}}, we write GCD(a_{1}, a_{2},…, a_{m}) for GCD(S). Unless otherwise noted, this lecture will assume that S has at most two elements S = {a, b}. In this case, the GCD function may be regarded as a two argument function, GCD(a, b). It is called the simple GCD function, as opposed to the general multiple GCD function for sets. Of course, if S has m ≥ 2 elements, we can compute GCD(S) using m − 1 simple GCD computations.

The following is easy.

GCD(1, b)  = 1
GCD(0, b) = b̂ where b̂ is the distinguished associate of b
GCD(a, b)  = GCD(b, a)
GCD(a + b, b) = GCD(a, b)
GCD(ua, b)  = GCD(a, b) where u is a unit

Say a, b are relatively prime if GCD(a, b) = 1.

For instance, GCD(123, 234) = 3 and GCD(3X^{4} − 6X^{3} + 13X^{2} − 8X + 12, 6x^{5} + 17X^{3} − 3X^{2} + 12X − 4) = 3X^{2} + 4.

GCD for ideals. Although we began with UFD’s such as and [X], our Fundamental Problems force us to consider more general domains such as number rings (§VI.3). These rings need not be UFD’s (exercise below). This led Kummer, Dedekind and Kronecker to develop ideal theory for algebraic numbers^{1}. To regain the UFD property, we generalize numbers to ideals and introduce the concept of prime ideals. The ideal theoretic analogue of UFD’s is this: a Dedekind domain is a one in which every ideal is a product of prime ideals. It can be proved that such prime ideal factorizations are unique (e.g., [12, p. 273]). Number rings are Dedekind domains.

We do not define the concept of ideal divisibility via ideal products. Instead, if I, J ⊆ D are ideals, we define I to be a divisor of J, and say I divides J, to mean I ⊇ J.

Perhaps this definition is unexpected, but it is a stroke of inspiration from Dedekind (1871). Consider ideals in : they have the form (n) where n ∈ since is a principal ideal domain (§3). Hence we can identify ideals of with numbers. Then m, n ∈ has the property that m|n iff (m) ⊇ (n), “agreeing” with our definition. In general, the relationship between ideal quotient and divisor property is only uni-directional: for ideals I, J ⊆ D, we have that I ⊇ IJ and so I divides IJ.

It follows that the GCD of a set S of ideals is defined to be the smallest ideal that divides each I ∈ S, i.e.,

GCD(S) = ⋂_{I∈S} I.

It is easy to check that the intersection of any set of ideals is an ideal, so GCD is well defined. In particular, for I = (a_{1},…, a_{m}) and J = (b_{1},…, b_{n}) then

GCD(I, J) = I + J = (a_{1},…, a_{m}, b_{1},…, b_{n}).  (1)

So the GCD problem for ideals is trivial unless we require some other conditions on the ideal generators. For instance, for the ideals of , the GCD of (a) and (b) is the ideal (a, b). But since is a principal ideal domain, we know that (a, b) = (d) for some d ∈ . We then interpret the GCD problem in to mean the computation of d from a, b. It is not hard to prove that d is what we have defined to be a greatest common divisor of a, b. Thus, the common notation ‘(a, b)’ for GCD(a, b) is consistent with the ideal theoretic notation! In general, for a, b in a UFD, one should not expect Ideal(a, b) to be generated by the GCD(a, b). For instance, [X] is a UFD, GCD(2, X) = 1 but Ideal(2, X) ≠ Ideal(1).

1The other historical root for ideal theory is rational function fields in one variable.

§2. Euclid’s Algorithm  Lecture II  Page 4

Exercise 1.3:

a) Is the set of ideals of a domain D under the ideal sum and ideal product operations a ring? The obvious candidates for the zero and unity elements are (0) and (1) = D.

b) Verify equation (1).

Exercise 1.4: Say a domain D is factorable if every non-unit of D is a finite product of irreducible elements. Prove that a factorable domain D is a UFD iff irreducible elements are prime.

Exercise 1.5: We will prove that the number ring [\sqrt{−5}] = \{x + y\sqrt{−5} : x, y ∈ \} (cf. §VI.3) is not a UFD. The norm of an element x + y\sqrt{−5} is N(x + y\sqrt{−5}) := x^{2} + 5y^{2}.

(a) Show that N(ab) = N(a)N(b).

(b) Show that N(a) = 1 iff a is a unit. Determine the units of [\sqrt{−5}].

(c) Show that if N(a) is a prime integer then a is irreducible in [\sqrt{−5}].

(d) Show that the numbers 2, 3, 1 + \sqrt{−5}, 1 − \sqrt{−5} are irreducible and not associates of each other. Since 6 = 2 · 3 = (1 + \sqrt{−5}) · (1 − \sqrt{−5}), conclude that [\sqrt{−5}] is not a UFD.

Exercise 1.6:

(a) In a principal ideal domain, the property “I ⊇ J” is equivalent to “there exists an ideal K such that IK = J”.

(b) In [X, Y], there does not exist an ideal K such that (X, Y) · K = (X^{2}, Y^{2}).

Exercise 1.7: (Lucas 1876) The GCD of two Fibonacci numbers is Fibonacci.

Exercise 1.8: (Kaplansky) Define a GCD-domain to be a domain in which any two elements has a GCD.

a) Show that if D is such a domain, then so is D[X].

b) Show that if for any two elements u, v ∈ D, either u|v or v|u (D is a valuation domain) then D is a GCD-domain.

§2. Euclid’s Algorithm

Let us recall Euclid’s algorithm for computing the GCD of two positive integers

m_{0} > m_{1} > 0.

We construct the sequence of remainders,

m_{0}, m_{1}, m_{2},…, m_{k},  (k ≥ 1)  (2)

where

m_{i+1} = m_{i−1} mod m_{i}  (i = 1,…, k − 1)
0 = m_{k−1} mod m_{k}.

Recall that a mod b is the remainder function that returns an integer in the range [0, b). But this is not the only possibility (next section).

Let us prove that m_{k} equals GCD(m_{0}, m_{1}). We use the observation that if any number d divides m_{i} and m_{i+1}, then d divides m_{i−1} (provided i ≥ 1) and d divides m_{i+2} (provided i ≤ k−2). Note that m_{k} divides m_{k} and m_{k−1}. So by repeated application of the observation, m_{k} divides both m_{0} and m_{1}. Next suppose d is any number that divides m_{0} and m_{1}. Then repeated application of the observation implies d divides m_{k}. Thus we conclude that m_{k} = GCD(m_{0}, m_{1}).

Certain data related to the GCD(m_{0}, m_{1}) is critical to some applications. Namely, there exist integers s, t such that

GCD(m_{0}, m_{1}) = sm_{0} + tm_{1}.  (3)

We call (s, t) a pair of co-factors of (m_{0}, m_{1}). By definition, the extended Euclidean algorithm computes the co-factors of its input pair on numbers. It is easy to obtain the GCD from the co-factors although most extended Euclidean algorithms also produces the GCD. The existence of co-factors will be proved by our construction of the extended Euclidean algorithm next.

We proceed as follows. Suppose q_{i} is the quotient of the ith remaindering step in (2):

m_{i+1} = m_{i−1} − q_{i}m_{i}  (i = 2,…, k − 1)  (4)

We compute two auxiliary sequences

(5) (s_{0}, s_{1},…, s_{k}) and (t_{0}, t_{1},…, t_{k})

so that they satisfy the property

m_{i} = s_{i}m_{0} + t_{i}m_{1},  (i = 0,…, k).  (6)

Note that when i = k, this property is our desired equation (3). The auxiliary sequences are obtained by mirroring the remaindering step (4),

\begin{matrix}s_{i+1} = s_{i−1} − q_{i}s_{i} \\ t_{i+1} = t_{i−1} − q_{i}t_{i} \end{matrix} \}  (i = 2,…, k − 1)  (7)

To initialize the values of s_{0}, s_{1} and t_{0}, t_{1}, observe that

m_{0} = 1 · m_{0} + 0 · m_{1}

and

m_{1} = 0 · m_{0} + 1 · m_{1}.

Thus (6) is satisfied for i = 0, 1 if we set

(s_{0}, t_{0}) :=(1, 0), (s_{1}, t_{1}) :=(0, 1).

Inductively, (6) is satisfied because

m_{i+1} = m_{i−1} − q_{i}m_{i}
= (s_{i−1}m_{0} + t_{i−1}m_{1}) − q_{i}(s_{i}m_{0} + t_{i}m_{1})
= (s_{i−1} − q_{i}s_{i})m_{0} + (t_{i−1} − q_{i}t_{i})m_{1}
= s_{i+1}m_{0} + t_{i+1}m_{1}.

This completes the description and proof of correctness of the extended Euclidean algorithm.

§3. Euclidean Ring  Lecture II  Page 6

Application. Suppose we want to compute multiplicative inverses modulo an integer m_{0}. An element m_{1} has a multiplicative inverse modulo m_{0} if and only if GCD(m_{0}, m_{1}) = 1. Applying the extended algorithm to m_{0}, m_{1}, we obtain s, t such that

1 = GCD(m_{0}, m_{1}) = sm_{0} + tm_{1}.

But this implies

1 ≡ tm_{1}(mod m_{0}),

i.e., t is the inverse of m_{1} modulo m_{0}. Similarly s is the inverse of m_{0} modulo m_{1}.

Exercise 2.1: Show that every two steps of the Euclidean algorithm reduces the (bit) size of the larger integer by at least one. Conclude that the bit complexity of the Euclidean algorithm is Θ(nM_{B}(n)) = Θ(n^{2}L(n)) where M_{B}(n) is the bit complexity of integer multiplication.

Exercise 2.2: Consider the extended Euclidean algorithm.

(a) Show that for i ≥ 2, we have s_{i}t_{i} < 0 and s_{i} > 0 iff i is even.

(b) Show that the co-factors s, t computed by the algorithm satisfy max{|s|, |t|} < m_{0}.

Exercise 2.3: (Blankinship) The following is a simple basis for an extended multiple GCD algorithm. Let N = (n_{1},…, n_{k})^{T} be a k-column of integers and A the k × (k + 1) matrix whose first column is N, and the remaining columns forms an identity matrix. Now perform any sequence of row operations on A of the form “subtract an integer multiple of one row from another”. It is clear that we can a sequence of such operations so that the first column eventually contains only one non-zero entry d where d = GCD(n_{1},…, n_{k}). If the row containing d is (d, s_{1},…, s_{k}), prove that

d = ∑_{i=1}^{k} s_{i}n_{i}.

Exercise 2.4:

a) Let n_{1} > n_{2} > ⋯ > n_{k} > 1 (k ≥ 1) be integers. Let S = (s_{1},…, s_{k}) ∈ ^{k} be called a syzygy of N = (n_{1},…, n_{k}) if ∑_{i=1}^{k} s_{i}n_{i} = 0. Prove that the set of syzygies of N forms a -module. For instance, let s_{ij} (1 ≤ i < j ≤ n) be the k-vector (0,…, 0, n_{j}, 0,…, 0, −n_{i}, 0,…, 0) (where the only non-zero entries are at positions i and j as indicated). Clearly s_{ij} is a syzygy. It can be shown (see the Lecture on Hilbert basis theorem) that this module has a finite basis. Construct such a basis.

b) Two k-vector S, S′ are equivalent if S − S′ is a syzygy of N. Show that every S is equivalent to some S′ where each component c of S′ satisfies |c| < n_{0}.

§3. Euclidean Ring

We define the abstract properties that make Euclid’s algorithm work. A ring R is Euclidean if there is a function

φ : R → {−∞} ∪

such that i) b ≠ 0 and a|b implies φ(a) ≤ φ(b);

ii) for all r ∈ , the set {φ(a) : a ∈ R, φ(a) < r} is finite;

iii) for all a, b ∈ R (b ≠ 0), there exists q, r ∈ R such that

a = bq + r and φ(r) < φ(b)

.

A Euclidean domain is an Euclidean ring that is also a domain. We say that φ is an Euclidean value function for R, and call q and r (respectively) a quotient and remainder of a, b relative to φ. Property iii) is called the division property (relative to φ). We introduce the remainder rem(a, b) and quotient quo(a, b) functions that pick out some definite pair of remainder and quotient of a, b that simultaneously satisfies property iii). Note that these functions are only defined when b ≠ 0. Often it is convenient to write these two functions using infix operators mod and div:

rem(a, b) = a mod b, quo(a, b) = a div b.  (8)

Exercise 3.1:

a) rem(a, b) = 0 if and only if b|a (in particular, rem(0, b) = 0).

b) φ(a) = φ(b) when a and b are associates.

c) φ(0) < φ(b) for all non-zero b.

Example: Our two standard domains, and F [X], are Euclidean:

(A) is seen to be an Euclidean domain by letting φ(n) = |n|, the absolute value of n. There are two choices for rem(m, n) unless n|m, one positive and one negative. For instance, rem(8, 5) can be taken to be 3 or −2. There are two standard ways to make rem(m, n) functional. In these lectures, we always choose the non-negative remainder. The corresponding function rem(m, n) ≥ 0 is called the non-negative remainder function. An alternative is to choose the remainder that minimizes the absolute value (choosing the positive one in case of a tie); this corresponds to the symmetric remainder function. The function quo(a, b) is uniquely determined once rem(a, b) is fixed. (B) If F is any field, the following division property for polynomials holds: for A, B ∈ F [X] where B ≠ 0, there exists Q, R ∈ F [X] such that

A = BQ + R, deg(R) < deg(B).

This can be proved by the synthetic division algorithm which one learns in high school. It follows that the polynomial ring F [X] is an Euclidean domain, as witnessed by the choice φ(P) = deg P, for P ∈ F [X]. In fact, the synthetic division algorithm shows that rem(P, Q) and quo(P, Q) are uniquely determined. Despite property ii) in the definition of φ, there may be infinitely many a ∈ R with φ(a) < r. This is the case if R = F [X] with F infinite.

Lemma 2 If a is a proper divisor of b then φ(a) < φ(b).

Proof. By the division property, a = bq + r where φ(r) < φ(b). Now r ≠ 0 since otherwise b divides a, which contradicts the assumption that a properly divides b. Since a properly divides b, let b = ac for some c. Then r = a−bq = a(1−cq). Then property i) implies φ(a) ≤ φ(r) < φ(b). Q.E.D.

Theorem 3 An Euclidean ring is a principal ideal ring. Indeed, if b ∈ I \ {0} is such that φ(b) is minimum then I = Ideal(b).

Proof. Let I be any ideal. By property ii), there exists a b ∈ I \ {0} such that φ(b) is minimum. To show I = Ideal(b), it suffices to show that b divides any c ∈ I \{0}. By the division property, c = bq + r where φ(r) < φ(b). If r ≠ 0 then we have found an element r = c − bq ∈ I \ {0} with φ(r) < φ(b), contradicting our choice of b. If r = 0 then b|c. Q.E.D.

Lemma 4 In an Euclidean ring R, the non-zero irreducible elements are prime.

Proof. Let p ∈ R \ {0} be irreducible. If p divides the product bc, we must prove that p divides b or p divides c. By the division property, c = pq + r where φ(r) < φ(p). If r = 0 then p|c and we are done. Assume otherwise. Note that r ∉ Ideal(p) (otherwise p|r and φ(p) ≤ φ(r), contradiction). So Ideal(p, r) properly contains Ideal(p). Let Ideal(p, r) = Ideal(u) for some u ∈ R (since R is a principal ideal ring). We claim that u is a unit. If not, u|p and since p is irreducible, u and p are associates. This implies Ideal(u) = Ideal(p), contradiction. Hence 1 ∈ Ideal(p, r) or there exists α, β such that

1 = αp + βr
= αp + β(b − pq)
= (α − βq)p + βb.
c = (α − βq)pc + βbc
= (α − βq)pc + βpp′,

for some p′ (since p|bc). Since p divides the right hand side, p|c. Q.E.D.

Theorem 5 In a principal ideal ring, the factorization of a non-unit into irreducible non-units are unique, up to reordering and associates.

Proof. Suppose b ∈ R is a non-unit with two factorization into irreducible non-units:

b = p_{1}p_{2} ⋯ p_{m} = q_{1}q_{2} ⋯ q_{n}, 1 ≤ m ≤ n.

We use induction on m. If m = 1 then clearly n = 1 and q_{1} = p_{1}. Assume m > 1. Since p_{1} is a prime, it must divide some q_{i}, and we might as well assume p_{1}|q_{1}. But q_{1} is also a prime and so it must be an associate of p_{1}. Dividing by p_{1} on both sides of the expression, it follows that p_{2} ⋯ p_{m} = q′_{2}q_{3} ⋯ q_{n} where q′_{2} is an associate of q_{2}. By induction, m = n and the two factorization is unique up to reordering and associates. This implies our theorem. Q.E.D.

Corollary 6 An Euclidean domain is a UFD.

Remainder sequences. Relative to the remainder and quotient functions, we define a remainder sequence for any pair a, b ∈ R to be a sequence

a_{0}, a_{1},…, a_{k}  (k ≥ 1)  (9)

such that a_{0} = a, a_{1} = b and for i = 1,…, k − 1, a_{i+1} is an associate of rem(a_{i−1}, a_{i}), and rem(a_{k−1}, a_{k}) = 0. Note that termination of this sequence is guaranteed by property ii). The remainder sequence is strict if a_{i+1} is any remainder of a_{i−1}, a_{i} for all i; it is Euclidean if a_{i+1} = rem(a_{i−1}, a_{i}). In subsequent lectures we see non-strict remainder sequences as exemplified by Sturm sequences and subresultant sequences.

Example: In , (13, 8, 5, 3, 2, 1), (13, 8, −3, 2, ±1) and (13, 8, 3, −1) are all strict remainder sequences for (13, 8). A non-strict remainder sequence for (13, 8) is (13, 8, −5, 2, 1).

Associated to each remainder sequence (9) is another sequence

q_{1}, q_{2},…, q_{k}  (10)

where a_{i−1} = a_{i}q_{i} + u_{i}a_{i+1} (i = 1,…, k − 1, u_{i} is a unit) and a_{k−1} = a_{k}q_{k}. We call (10) the quotient sequence associated to the remainder sequence.

Norms. In some books, the function φ is restricted to the range . This restriction does not materially affect the concept of Euclidean domains, and has the advantage that property ii) is automatic. Our formulation makes it more convenient to formulate functions φ that have other desirable properties. For instance, we often find the properties:

φ(ab) = φ(a)φ(b)

and

φ(a + b) ≤ φ(a) + φ(b).

In this case, we call φ a multiplicative norm (or valuation), and might as well (why?) assume φ(0) = 0 and φ(1) = 1. Similarly, if φ(ab) = φ(a)+φ(b) and φ(a+b) ≤ O(1)+max{φ(a), φ(b)} then we call φ an additive norm, and might as well assume φ(0) = −∞ and φ(1) = 0. Clearly φ is multiplicative implies log φ is additive.

Remarks. The number rings _{α} (§VI.3) have properties similar to the integers. In particular, they support the concepts of divisibility and factorization. Gauss pointed out that such rings may not be a UFD (class number 1). Even when _{α} is a UFD, it may not be Euclidean; the “simplest” example is \sqrt{−19} (see [8]). An obvious candidate for the Euclidean value function φ is the norm of algebraic numbers, but other functions are conceivable. Turning now to the quadratic number rings (i.e., _{α} where α = \sqrt{d} and d is squarefree), Kurt Heegner [Diophantische analysis und modulfunktionen, Mathematische Zeitschrift, vol. 56, 1952, pp.227– 253] was the first ^{2} to prove that there are exactly nine such UFD’s in which d < 0, viz., d = −1, −2, −3, −7, −11, −19, −43, −67, −163. In contrast, it is conjectured that there are infinitely many UFD’s among real (i.e., d > 0) quadratic number fields. It is known that there are precisely 21 real quadratic domains that support the Euclidean algorithm. Currently, the most general GCD algorithms are from Kaltofen and Rolletschek [5] who presented polynomialtime GCD algorithms for each quadratic number ring \sqrt{d} that is a UFD, not necessarily Euclidean.

Exercise 3.2: Justify the above remarks about multiplicative and additive norms.

2This result is often attributed to Stark and Baker who independently proved this in 1966. See Buell [4].

§4. HGCD  Lecture II  Page 10

Exercise 3.3: Verify that the Euclidean algorithm computes the GCD in an Euclidean domain, relative to the function rem(a, b).

Exercise 3.4: An Euclidean domain is a UFD.

Exercise 3.5: Show that the number ring _{ı} (= [ı] = {a+ıb : a, b ∈ }) of Gaussian integers forms an Euclidean domain with respect to the multiplicative norm φ(a + ıb) = a^{2} + b^{2}. Use the identity of Fibonacci,

φ(xy) = φ(x)φ(y), x, y ∈ [ı].

What are the possible choices for defining the remainder and quotient functions here?

Exercise 3.6: Show that \sqrt{−19} is not an Euclidean domain with respect to the standard norm function. Note: \sqrt{−19} = \{m + nω : m, n ∈ \} where ω = \frac{1+\sqrt{−19}}{2}. The norm of x + y\sqrt{−19} ∈ (\sqrt{−19}) is by x^{2} + 19y^{2}. HINT: What is the remainder of 5 divided by

\sqrt{−19}?

§4. The Half-GCD problem

An exercise in §2 shows that that Euclid’s algorithm for integers has bit complexity Θ(n^{2}L(n)). Knuth [6] is the first to obtain a subquadratic complexity for this problem. In 1971, Schönhage [9] improved it to the current record of O(log nM_{B}(n)) = nL^{2}(n).

Since F [X] is an Euclidean domain, Euclid’s algorithm can be applied to polynomials as well. Given P_{0}, P_{1} ∈ F [X] with n = deg P_{0} > deg P_{1} ≥ 0, consider its Euclidean remainder sequence

P_{0}, P_{1}, P_{2},…, P_{h}  (h ≥ 1).  (11)

Call the sequence normal if deg P_{i−1} = 1 + deg P_{i} for i = 2,…, h. A random choice for P_{0}, P_{1} gives rise to a normal sequence with high probability (this is because non-normal sequences arise from the vanishing of certain determinants involving the the coefficients of P_{0}, P_{1}, as shown in the next lecture). The algebraic complexity of this Euclidean algorithm is therefore

O(nM_{A}(n))  (12)

where M_{A}(n) = O(n log n) is the algebraic complexity of polynomial multiplication (Lecture 1). Moenck [7] improves (12) to O(log nM_{A}(n)) in case the remainder sequence is normal. Aho-Hopcroft-Ullman [1] incorrectly proved that the Moenck algorithm works in general. Brent-Gustavson-Yun [2] presented a corrected version, but without a proof. Independently, Thull and Yap [11] discovered the error and presented a fairly simple correctness proof, reproduced below. In this lectures, we follow the development in [11] who gave a unified framework for both the polynomial and integer GCD algorithms.

Let us motivate the approach of Schönhage and Moenck. These ideas are best illustrated in the case of the polynomials. If the sequence (11) is normal with n = h and deg P_{i} = n − i (i = 0,…, n) then

∑_{i=0}^{h} deg P_{i} = n(n − 1)/2.

So any algorithm that explicitly computes each member of the remainder sequence has at least quadratic complexity. On the other hand, if

Q_{1}, Q_{2},…, Q_{h}

is the quotient sequence associated to (11), then it is not hard to show that

∑_{i=1}^{h} deg Q_{i} = n.  (13)

Indeed, we can quickly (in O(n log^{2} n) time, see Exercise below) obtain any member of the remainder sequence from the Q_{i}’s. This suggests that we redirect attention to the quotient sequence.

Matrix Terminology. To facilitate description of our algorithms, we resort to the language of matrices and vectors. In this lecture, all matrices will be 2×2 matrices and all vectors will be column 2-vectors. The identity matrix is denoted by E. The Euclidean algorithm as embodied in (11) can be viewed as a sequence of transformations of 2-vectors:

( \begin{matrix}P_{0} \\ P_{1} \end{matrix} ) ⟶^{M_{1}} ( \begin{matrix}P_{1} \\ P_{2} \end{matrix} ) ⟶^{M_{2}} ⋯ ⟶^{M_{h−1}} ( \begin{matrix}P_{h−1} \\ P_{h} \end{matrix} ) ⟶^{M_{h}} ( \begin{matrix}P_{h} \\ 0 \end{matrix} ) .  (14)

Precisely, if U, V are vectors and M a matrix, we write

U ⟶^{M} V

to mean that U = MV. Hence equation (14) can be correctly interpreted if we define

M_{i} = ( \begin{matrix}Q_{i} & 1 \\ 1 & 0 \end{matrix} ) .

In general, an elementary matrix refers to a matrix of the form M = ( \begin{matrix}Q & 1 \\ 1 & 0 \end{matrix} ) where Q is a

polynomial with positive degree. We call Q the partial quotient in M. A regular matrix M is a product of zero or more elementary matrices,

M = M_{1}M_{2} ⋯ M_{k}  (k ≥ 0).  (15)

When k = 0, M is interpreted to be E. The sequence Q_{1},…, Q_{k} of partial quotients associated with the elementary matrices M_{1},…, M_{k} in equation (15) is called the sequence of partial quotients of M. Also, Q_{k} is called its last partial quotient. Note that regular matrices have determinant ±1 and so are invertible. Regular matrices arise because

U ⟶^{M} V and V ⟶^{M′} W implies U ⟶^{MM′} W.

Our terminology here is motivated by the connection to continued fractions (for instance, regular matrices are related to regular continued fractions).

We are ready to define the Half-GCD (or, HGCD) problem for a polynomial ring F [X]:

Given P_{0}, P_{1} ∈ F [X] where n = deg P_{0} > deg P_{1}, compute the regular matrix M such that if

( \begin{matrix}P_{0} \\ P_{1} \end{matrix} ) ⟶^{M} ( \begin{matrix}P_{2} \\ P_{3} \end{matrix} )

then

deg P_{2} ≥ n/2 > deg P_{3}.  (16)

In general, we say two numbers a, b straddle a third number c if a ≥ c > d. Thus deg P_{2}, deg P_{3} straddle n/2 in equation (16).

Now we show how to compute the GCD using a HGCD subroutine. In fact we present the “extended version” of the GCD algorithm:

Algorithm Extended Polynomial GCD(P_{0}, P_{1}): Input: A pair of polynomials with deg P_{0} > deg P_{1} Output: A regular matrix M such that

( \begin{matrix}P_{0} \\ P_{1} \end{matrix} ) ⟶^{M} ( \begin{matrix}GCD(P_{0}, P_{1}) \\ 0 \end{matrix} )
[1] Compute M_{0} ← Hgcd(P_{0}, P_{1}).
[2] Recover P_{2}, P_{3} via

( \begin{matrix}P_{2} \\ P_{3} \end{matrix} ) ← M_{0}^{−1} · ( \begin{matrix}P_{0} \\ P_{1} \end{matrix} ) .

[3] if P_{3} = 0 then return(M_{0}).

else, perform one step of the Euclidean algorithm,

( \begin{matrix}P_{2} \\ P_{3} \end{matrix} ) ⟶^{M_{1}} ( \begin{matrix}P_{3} \\ P_{4} \end{matrix} )

where M_{1} is an elementary matrix.

[4] if P_{4} = 0 then return(M_{0}M_{1}). else, recursively compute M_{2} ← EGCD(P_{3}, P_{4}) return(M_{0}M_{1}M_{2}).

The correctness of this algorithm is clear. The reason for step [3] is to ensure that in our recursive call, the degree of the polynomials is less than n/2. The algebraic complexity T (n) of this GCD algorithm satisfies

T (n) = T′(n) + O(M_{A}(n)) + T (n/2)

where T′(n) is the complexity of the HGCD algorithm. Let us assume that

M_{A}(n) = O(T′(n)), T′(αn) ≤ αT′(n).

For instance, the first relation holds if T′(n) = Ω(M_{A}(n)); the second relation holds if T′(n) is bounded by a polynomial. In particular, they hold for T′(n) = Θ(M(n) log n), which is what we will later prove. Then

T (n) = O(T′(n) + T′(n/2) + T′(n/4) + ⋯) = O(T′(n)).

In conclusion, the complexity of the GCD problem is the same order of the complexity as the HGCD problem. Henceforth, we focus on the HGCD problem.

Remarks: the above complexity counts ring operations from F. If we count operations in F [X], the complexity is O(n log n). This counting is more general because it applies also to the case of integer HGCD to be discussed. Strassen [10] has proved that this complexity is optimal: O(n log n) is both necessary and sufficient.

Exercise 4.1: Recall the auxilliary sequences (s_{0}, s_{1},…, s_{k}) and (t_{0}, t_{1},…, t_{k}) computed in the Extended Euclidean algorithm (§2) for the GCD of pair of integers a_{0} > a_{1} > 1. Show

that the appropriate elementary matrices have the form ( \begin{matrix}s_{i} & t_{i} \\ s_{i+1} & t_{i+1} \end{matrix} )^{−1}.

§5. Norm  Lecture II  Page 13

Exercise 4.2:

a) Show equation (13).

b) Show that in O(n log^{2} n) time, we can reconstruct the polynomials S, T from the quotient sequence Q_{1},…, Q_{k} where SP_{0} + T P_{1} = GCD(P_{0}, P_{1}). HINT: note that

( \begin{matrix}P_{i} \\ P_{i+1} \end{matrix} ) = ( \begin{matrix}0 & 1 \\ 1 & −Q_{i} \end{matrix} ) ( \begin{matrix}P_{i−1} \\ P_{i} \end{matrix} ), and use a balanced binary tree scheme.

§5. Properties of the Norm

For the rest of this Lecture, the domain D refers to either or F [X]. Define the (additive) norm ∥a∥ of a ∈ D thus:

∥a∥:= \{ \begin{matrix}log_{2} |a| & if a ∈ \\ deg(a) & if a ∈ F [X] \end{matrix}

This ∥ · ∥ notation is special to the present Lecture.

The previous section describes the polynomial HGCD problem. A similar, albeit more complicated, development can be carried out for integers. We now describe a common framework for both the integer and polynomial HGCD algorithms.

The following properties are easy to check:

a) ∥a∥ ∈ {−∞} ∪ ^{∗} where ^{∗} is the set of non-negative real numbers.

b) ∥a∥ = −∞ ⟺ a = 0

c) ∥a∥ = 0 ⟺ a is a unit.

d) ∥ − a∥ = ∥a∥

e) ∥ab∥ = ∥a∥ + ∥b∥

f) ∥a + b∥ ≤ 1 + max{∥a∥, ∥b∥}.

The last two properties imply that the norm is additive (§3). However, polynomials satisfy the stronger non-Archimedean property:

∥a + b∥ ≤ max{∥a∥, ∥b∥}.

It is the non-Archimedean property that makes polynomials relatively easier than integers. This property implies

∥a + b∥ = max{∥a∥, ∥b∥}  if ∥a∥ ≠ ∥b∥.  (17)

Exercise 5.1: Prove this.

This norm function serves as the Euclidean value function for D. In particular, the division property relative to the norm holds: for any a, b ∈ D, b ≠ 0, there exists q, r ∈ D such that

a = qb + r, ∥r∥ < ∥b∥.

The remainder and quotient functions, rem(a, b) and quo(a, b), can be defined as before. Recall that these functions are uniquely defined for polynomials, but for integers, we choose rem(a, b) to be the non-negative remainder function. Note that in the polynomial case,

∥a mod X^{m}∥ ≤ min{∥a∥, m − 1}  (18)
∥a div X^{m}∥ = \{ \begin{matrix}∥a∥ − m & if ∥a∥ ≥ m \\ −∞ & else. \end{matrix}  (19)

Matrices and vectors. The previous section (§4) introduced the matrix concepts we needed. Those definitions extend in the obvious way to our present setting, except for one place, where

we need special care: A matrix of the form M = ( \begin{matrix}q & 1 \\ 1 & 0 \end{matrix} ) (where q ∈ D) is denoted ⟨q⟩. A

matrix is elementary if it has the form ⟨q⟩ where ∥q∥ > 0 in the case of polynomials (as before), and q > 0 in the case of integers. A finite product ⟨q_{1}⟩⟨q_{2}⟩ ⋯ ⟨q_{k}⟩ (k ≥ 0) of elementary matrices is again called regular and may be denoted ⟨q_{1},…, q_{k}⟩. When k = 2, the careful reader will note the clash with our notation for scalar products, but this ambiguity should never confuse.

A regular matrix M = ( \begin{matrix}p & q \\ r & s \end{matrix} ) satisfies the following ordering property:

M ≠ E ⇒ ∥p∥ ≥ max{∥q∥, ∥r∥} ≥ min{∥q∥, ∥r∥} ≥ ∥s∥, ∥p∥ > ∥s∥.  (20)

Exercise 5.2:

a) Prove the ordering property.

b) If all the inequalities in the definition of the ordering property are in fact strict and ∥s∥ ≥ 0, we say M satisfies the strict ordering property. Show that the product of three or more elementary matrices has the strict ordering property.

c) Bound the norms of the entries of the matrix ⟨q_{1},…, q_{k}⟩ in terms of the individual norms ∥q_{i}∥.

For vectors U, V and matrix M, we write

U ⟶^{M} V

(or simply, U ⟶ V) if U = MV. We say M reduces U to V if M is a regular matrix. If, in

addition, U = ( \begin{matrix}a \\ b \end{matrix} ) , V = ( \begin{matrix}a′ \\ b′ \end{matrix} ) such that ∥a∥ > ∥b∥ and ∥a′∥ > ∥b′∥, then we say this is an

Euclidean reduction.

A matrix is unimodular if^{3} it has determinant ±1. Clearly regular matrices are unimodular.

Thus their inverses are easily obtained: if M = ( \begin{matrix}p & q \\ r & s \end{matrix} ) is regular with determinant det M = δ then

3Note that in some literature, “unimodular” refers to determinant +1.

M^{−1} = δ ( \begin{matrix}s & −q \\ −r & p \end{matrix} ) = ± ( \begin{matrix}s & −q \\ −r & p \end{matrix} ) .  (21)

If U = ( \begin{matrix}a \\ b \end{matrix} ) then we write GCD(U) for the GCD of a and b. We say U, V are equivalent if

U = MV for some unimodular matrix M.

Exercise 5.3: Prove that this binary relation between vectors is indeed a mathematical equivalence relation.

Lemma 7 U and V are equivalent if and only if GCD(U) = GCD(V).

Proof. It is easy to check that if the two vectors are equivalent then they must have the same

GCD. Conversely, by Euclid’s algorithm, they are both equivalent to the vector ( \begin{matrix}g \\ 0 \end{matrix} ) where

g is their common GCD. Q.E.D.

The following is a key property of Euclidean remainder sequences:

Lemma 8 Given a, b, a′, b′ such that ∥a∥ > ∥b∥ ≥ 0. The following are equivalent:

(i) a′, b′ are consecutive elements in the Euclidean remainder sequence of a, b. (ii) There is a regular matrix M such that

( \begin{matrix}a \\ b \end{matrix} ) ⟶^{M} ( \begin{matrix}a′ \\ b′ \end{matrix} )  (22)

and either ∥a′∥ > ∥b′∥ ≥ 0 (polynomial case) or a′ > b′ > 0 (integer case).

Proof. If (i) holds then we can (by Euclid’s algorithm) find some regular matrix M satisfying (ii). Conversely assume (ii). We show (i) by induction on the number of elementary matrices in the product M. The result is immediate if M = E. If M is elementary, then (i) follows from the division property for our particular choices for D. Otherwise let M = M′′M′ where M′ is elementary and M′′ is regular. Then for some a′′, b′′,

( \begin{matrix}a \\ b \end{matrix} ) ⟶^{M′′} ( \begin{matrix}a′′ \\ b′′ \end{matrix} ) ⟶^{M′} ( \begin{matrix}a′ \\ b′ \end{matrix} ) .

But a′′ = a′q′ + b′ and b′′ = a′ where q′ is the partial quotient of M′. We verify that this means ∥a′′∥ > ∥b′′∥. By induction, a′′, b′′ are consecutive elements in a strict remainder sequence of a, b. Then (i) follows. Q.E.D.

Exercise 5.4: In exercise 2.1, we upper bound the length of the integer Euclidean remainder sequence of a > b > 0 by 2 log_{2} a. We now give a simple improvement.

a) Prove that for k ≥ 1,

\underbrace{⟨1,…, 1⟩}_{k} = ⟨1⟩^{k} = ( \begin{matrix}F_{k+1} & F_{k} \\ F_{k} & F_{k−1} \end{matrix} )

§6. Polynomial HGCD  Lecture II  Page 16

where {F_{i}}_{i≥0} is the Fibonacci sequence defined by: F_{0} = 0, F_{1} = 1 and F_{i+1} = F_{i} + F_{i−1} (i ≥ 1).

b) Let ϕ be the positive root of the equation X^{2}−X −1 = 0 (so ϕ = (1+\sqrt{5})/2 = 1.618...). Prove inductively that

(1.5)^{k−1} ≤ F_{k+1} ≤ ϕ^{k}.

c) Say (q_{1}, q_{2},…, q_{k}) is the quotient sequence associated to the remainder sequence of

a > b > 0. If ⟨q_{1},…, q_{k}⟩ = ( \begin{matrix}p & q \\ r & s \end{matrix} ) prove that

∥p∥ ≤ ∥a∥ − ∥b∥.

d) Conclude that k < 1 + log_{1.5} a.

§6. Polynomial HGCD

We describe the polynomial HGCD algorithm and prove its correctness.

Parallel reduction. The idea we exploit in the HGCD algorithm might be called “parallel reduction”. Suppose we want to compute the HGCD of the pair of polynomials A, B ∈ F [X] where deg A = 2m. First we truncate these polynomials to define

( \begin{matrix}A_{0} \\ B_{0} \end{matrix} ):= ( \begin{matrix}A div X^{m} \\ B div X^{m} \end{matrix} ) .  (23)

Suppose that R is the matrix returned by HGCD(A_{0}, B_{0}); so

( \begin{matrix}A_{0} \\ B_{0} \end{matrix} ) ⟶^{R} ( \begin{matrix}A′_{0} \\ B′_{0} \end{matrix} )  (24)

for some A′_{0}, B′_{0}. Then we can define A′, B′ via

( \begin{matrix}A \\ B \end{matrix} ) ⟶^{R} ( \begin{matrix}A′ \\ B′ \end{matrix} ) .  (25)

Two reductions by the same matrix are said to be parallel. Thus (24) and (25) are parallel reductions. If A′, B′ turn out to be two consecutive terms in the remainder sequence of A, B, then we may have gained something! This is because we had computed R without looking at the lower order coefficients of A, B. But we need another property for R to be useful. We want the degrees of A′, B′ to straddle a sufficiently small value below 2m. By definition of HGCD, the degrees of A′_{0}, B′_{0} straddle m/2. A reasonable expectation is that the degrees of A′, B′ straddle 3m/2. This would be the case if we could, for instance, prove that

deg(A′) = m + deg(A′_{0}), deg(B′) = m + deg(B′_{0}).

This is not quite correct, as we will see. But it will serve to motivate the following outline of the HGCD algorithm.

Outline. Given A, B with deg(A) = 2m > deg(B) > 0, we recursively compute R ← HGCD(A_{0}, B_{0}) as above. Now use R to carry out the reduction of (A, B) to (A′, B′). Note that although the degrees of A′, B′ straddle 3m/2, we have no upper bound on the degree of A′. Hence we perform one step of the Euclidean reduction:

( \begin{matrix}A′ \\ B′ \end{matrix} ) ⟶^{⟨Q⟩} ( \begin{matrix}C \\ D \end{matrix} ) .

Now the degree of C = B′ is less than 3m/2. We can again truncate the polynomials C, D via

( \begin{matrix}C_{0} \\ D_{0} \end{matrix} ):= ( \begin{matrix}C div X^{k} \\ D div X^{k} \end{matrix} ) .

for a certain k ≥ 0. Intuitively, we would like to pick k = m/2. Then we make our second recursive call to compute S ← HGCD(C_{0}, D_{0}). Using M_{1} we reduce (C, D) to (C′, D′). Hopefully, the degrees of C′ and D′ straddle m, which would imply that our output matrix is

R · ⟨Q⟩ · S

The tricky part is that k cannot be simply taken to be m/2. This choice is correct only if the remainder sequence is normal, as Moenck assumed. Subject to a suitable choice of k, we have described the HGCD algorithm.

We are ready to present the actual algorithm. We now switch back to the norm notation, ∥A∥ instead of deg(A), to conform to the general framework.

Algorithm Polynomial HGCD(A, B): Input: A, B are univariate polynomials with ∥A∥ > ∥B∥ ≥ 0. Output: a regular matrix M which reduces (A, B) to (C′, D′)

where ∥C′∥, ∥D′∥ straddle ∥A∥/2.

[1] m ← ⌈\frac{∥A∥}{2}⌉; \{This is the magic threshold\} if ∥B∥ < m then return(E);

[2] ( \begin{matrix}A_{0} \\ B_{0} \end{matrix} ) ← ( \begin{matrix}A div X^{m} \\ B div X^{m} \end{matrix} ).

{now ∥A_{0}∥ = m′ where m + m′ = ∥A∥} R ← Hgcd(A_{0}, B_{0});

\{⌈\frac{m′}{2} ⌉ is the magic threshold for this recursive call\}

( \begin{matrix}A′ \\ B′ \end{matrix} ) ← R^{−1} ( \begin{matrix}A \\ B \end{matrix} );

[3] if ∥B′∥ < m then return(R);

[4] Q ← A′ div B′; ( \begin{matrix}C \\ D \end{matrix} ) ← ( \begin{matrix}B′ \\ A′ mod B′ \end{matrix} );

[5] l ← ∥C∥; k ← 2m − l; \{now l − m < ⌈\frac{m′}{2} ⌉\}

[6] C_{0} ← C div X^{k}; D_{0} ← D div X^{k}; {now ∥C_{0}∥ = 2(l − m)}

S ← Hgcd(C_{0}, D_{0});

{l − m is magic threshold for this recursive call.}

[7] M ← R · ⟨Q⟩ · S; return(M);

The programming variables in this algorithm are illustrated in the following figure.

Figure 1. Variables in the polynomial HGCD algorithm

To prove the correctness of this lemma, we must show that the output matrix M satisfies

( \begin{matrix}A \\ B \end{matrix} ) ⟶^{M} ( \begin{matrix}C′ \\ D′ \end{matrix} ) ,  ∥C′∥ ≥ ⌈\frac{∥A∥}{2}⌉ > ∥D′∥.  (26)

The Basic Setup. Let A, B ∈ F [X], ∥A∥ > ∥B∥ ≥ 0 and m ≥ 1 be given. Define A_{0}, B_{0} as in equation (23). This determines A_{1}, B_{1} via the equation

( \begin{matrix}A \\ B \end{matrix} ) = ( \begin{matrix}A_{0}X^{m} + A_{1} \\ B_{0}X^{m} + B_{1} \end{matrix} ) = ( \begin{matrix}A_{0} & A_{1} \\ B_{0} & B_{1} \end{matrix} ) ( \begin{matrix}X^{m} \\ 1 \end{matrix} ) .  (27)

Now let M be any given regular matrix. This determines A′_{0}, B′_{0}, A′_{1}, B′_{1} via

( \begin{matrix}A′_{0} & A′_{1} \\ B′_{0} & B′_{1} \end{matrix} ):= M^{−1} ( \begin{matrix}A_{0} & A_{1} \\ B_{0} & B_{1} \end{matrix} ) .  (28)

Finally, define A′, B′ via

( \begin{matrix}A′ \\ B′ \end{matrix} ):= ( \begin{matrix}A′_{0} & A′_{1} \\ B′_{0} & B′_{1} \end{matrix} ) ( \begin{matrix}X^{m} \\ 1 \end{matrix} ) .  (29)

Note the parallel reductions,

( \begin{matrix}A_{0} \\ B_{0} \end{matrix} ) ⟶^{M} ( \begin{matrix}A′_{0} \\ B′_{0} \end{matrix} ) , ( \begin{matrix}A \\ B \end{matrix} ) ⟶^{M} ( \begin{matrix}A′ \\ B′ \end{matrix} ) .

Lemma 9 (Correctness Criteria) Let A, B, m, M be given as in the Basic Setup, and define the remaining notations A_{i}, B_{i}, A′_{i}, B′_{i}, A′, B′ (i = 0, 1) as indicated. If

∥A′_{0}∥ > ∥B′_{0}∥,  (30)
∥A_{0}∥ ≤ 2∥A′_{0}∥  (31)

then

∥A′∥ = m + ∥A′_{0}∥, ∥B′∥ ≤ m + max{∥B′_{0}∥, ∥A_{0}∥ − ∥A′_{0}∥ − 1}.

In particular,

∥A′∥ > ∥B′∥.

Proof. Let M = ( \begin{matrix}P & Q \\ R & S \end{matrix} ). First observe that ∥A′_{0}∥ > ∥B′_{0}∥ and A_{0} = A′_{0}P + B′_{0}Q implies

∥A_{0}∥ = ∥A′_{0}∥ + ∥P ∥. Hence (31) is equivalent to

∥P ∥ ≤ ∥A′_{0}∥.  (32)

Since M^{−1} = ± ( \begin{matrix}S & −Q \\ −R & P \end{matrix} ) and A′_{1} = ±(A_{1}S − B_{1}Q),

∥A′_{1}∥ ≤ max{∥A_{1}S∥, ∥B_{1}Q∥} < m + ∥P ∥ ≤ m + ∥A′_{0}∥

Hence A′ = A′_{0}X^{m} + A′_{1} implies ∥A′∥ = m + ∥A′_{0}∥, as desired.

From B′_{1} = ±(−A_{1}R + B_{1}P) we get ∥B′_{1}∥ ≤ m − 1 + ∥P ∥ = m − 1 + ∥A_{0}∥ − ∥A′_{0}∥. From B′ = B′_{0}X^{m} + B′_{1} we get the desired inequality ∥B′∥ ≤ m + max{∥B′_{0}∥, ∥A_{0}∥ − ∥A′_{0}∥ − 1}. Q.E.D.

We call the requirement (31) the (lower) “threshold” for ∥A′_{0}∥. This threshold is the reason for the lower bound on ∥C′∥ in the HGCD output specification (26).

Finally we prove the correctness of the HGCD algorithm.

Lemma 10 (HGCD Correctness) Algorithm HGCD is correct: with input polynomials A, B where ∥A∥ > ∥B∥ ≥ 0, it returns a regular matrix M satisfying (26).

Proof. To keep track of the proof, the following sequence of reductions recalls the notations of the algorithm:

( \begin{matrix}A \\ B \end{matrix} ) ⟶^{R} ( \begin{matrix}A′ \\ B′ \end{matrix} ) ⟶^{⟨Q⟩} ( \begin{matrix}C \\ D \end{matrix} ) ⟶^{S} ( \begin{matrix}C′ \\ D′ \end{matrix} ) .  (33)

The algorithm returns a matrix in steps [1], [3] or [7]. Only when the algorithm reaches step [7] does the full sequence (33) take effect. It is clear that the returned matrix is always regular. So it remains to check the straddling condition of equation (26). In step [1], the result is clearly correct.

Consider the matrix R returned in step [3]: the notations m, A_{0}, B_{0}, A′, B′ in the algorithm conforms to those in Correctness Criteria (lemma 9), after substituting R for M. By induction hypothesis, the matrix R returned by the first recursive call (step [2]) satisfies

∥A′_{0}∥ ≥ ⌈m′/2⌉ > ∥B′_{0}∥, (m′ = ∥A_{0}∥)

where ( \begin{matrix}A_{0} \\ B_{0} \end{matrix} ) ⟶^{R} ( \begin{matrix}A′_{0} \\ B′_{0} \end{matrix} ). Then lemma 9 implies ∥A′∥ = m + ∥A′_{0}∥ ≥ m. Since m > ∥B′∥

is a condition for exit at step [3], it follows that the straddling condition (26) is satisfied at this exit.

Finally consider the matrix M returned in step [7]. Since we did not exit in step [3], we have m ≤ ∥B′∥. In step [4] we form the quotient Q and remainder D of A′ divided by B′. Also we renamed B′ to C. Hence m ≤ l where l = ∥C∥. To see that C_{0} is properly computed, let us verify

l ≥ k ≥ 0.  (34)

The first inequality in (34) follows from l ≥ m ≥ m + (m − l) = k. To see the second, l = ∥B′∥ ≤ m+max{∥B′_{0}∥, ∥A_{0}∥−∥A′_{0}∥+1} (Correctness Criteria) and so l ≤ m+max{⌈m′/2⌉− 1, ⌊m′/2⌋ + 1} ≤ m + ⌊m′/2⌋ + 1. Thus l − m ≤ ⌊m′/2⌋ + 1 ≤ m. Hence k = m − (l − m) ≥ 0, proving (34). In the second recursive call, HGCD(C_{0}, D_{0}) returns S. By induction,

∥C′_{0}∥ ≥ ⌈∥C_{0}∥/2⌉ > ∥D′_{0}∥,  where ( \begin{matrix}C_{0} \\ D_{0} \end{matrix} ) ⟶^{S} ( \begin{matrix}C′_{0} \\ D′_{0} \end{matrix} ) .  (35)

But ∥C_{0}∥ = l − k = 2(l − m) so (35) becomes

∥C′_{0}∥ ≥ l − m > ∥D′_{0}∥.

Now let ( \begin{matrix}C \\ D \end{matrix} ) ⟶^{S} ( \begin{matrix}C′ \\ D′ \end{matrix} ). Another application of lemma 9 shows that

∥C′∥ = k + ∥C′_{0}∥ ≥ k + l − m = m

and

∥D′∥ ≤ k + max{∥D′_{0}∥, ∥C_{0}∥ − ∥C′_{0}∥ − 1}
≤ k + max{l − m − 1, l − m − 1}
= k + l − m − 1 = m − 1.

This shows ∥C′∥ ≥ m > ∥D′∥ and hence (26). Q.E.D.

Remark: The proof shows we could have used k ← 2m − l − 1 as well. Furthermore, we could modify our algorithm so that after step [4], we return R · ⟨Q⟩ in case ∥D∥ < m. This may be slightly more efficient.

Complexity analysis. The HGCD algorithm makes two recursive calls to itself, Hgcd(A_{0}, B_{0}) and Hgcd(C_{0}, D_{0}). We check that ∥A_{0}∥ and ∥C_{0}∥ are both bounded by n/2. The work in each call to the algorithm, exclusive of recursion, is O(M_{B}(n)) = O(n log n). Hence the algebraic complexity T (n) of this HGCD algorithm satisfies

T (n) = 2T (n/2) + O(n log n).

This yields T (n) = O(n log^{2} n).

Exercise 6.1: Generalize the HGCD problem to the following: the function FGCD(A, B, f) whose arguments are polynomials A, B as in the HGCD problem, and f is a rational number between 0 and 1. FGCD(A, B, f) returns a matrix M that reduces the pair (A, B) to (A′, B′) such that ∥A′∥, ∥B′∥ straddle f∥A∥. Thus FGCD(A, B, 1/2) = Hgcd(A, B). Show that FGCD can be computed in the same complexity as Hgcd by using Hgcd as a subroutine.

Exercise 6.2: Modify the polynomial HGCD algorithm so that in step [5], the variable k is set to ⌈m/2⌉. This is essentially the algorithm of Moenck-Aho-Hopcroft-Ullman [1]. We want to construct inputs to make the algorithm return wrong answers. Note that since the modified algorithm works for inputs with normal remainder sequence (see §3), we are unlikely to find such inputs by generating random polynomials. Suppose the output M reduces the input (A, B) to (A′, B′). The matrix M may be wrong for several reasons:

§A. APPENDIX: Integer HGCD  Lecture II  Page 21

(i) ∥B′∥ > ⌈∥A∥/2⌉. (ii) ∥A′∥ ≤ ⌈∥A∥/2⌉. (iii) A′, B′ are not consecutive entries of the Euclidean remainder sequence of A, B. Construct inputs to induce each of these possibilities. (The possibilities (i) and (ii) are known to occur.)

§A. APPENDIX: Integer HGCD

We develop an integer version of the HGCD algorithm. We note two simple tricks. The first is to recover the non-Archimedean property thus: for a, b ∈ ,

ab ≤ 0  ⟹  ∥a + b∥ ≤ max{∥a∥, ∥b∥}.

One consequence of the non-Archimedean property we exploited was that if ∥a∥ ≠ ∥b∥ then ∥a + b∥ = max{∥a∥, ∥b∥}. Here is an integer analogue:

∥a∥ − ∥b∥ ≥ 1  ⟹  ∥a + b∥ = ∥a∥ ± ϵ, (0 ≤ ϵ ≤ 1).

To carry out a parallel reduction, the expected integer analogue would be to call HGCD on a div 2^{m}, b div 2^{m} for suitable m. Instead, the second trick will call HGCD on

a_{0} := 1 + (a div 2^{m}), b_{0} := b div 2^{m}.  (36)

The Basic Setup. We begin by trying to prove the analogue of the Correctness Criteria (lemma 9). The following notations will be fixed for the next two lemmas.

Assume that we are given a > b > 0 and m ≥ 1 where a ≥ 2^{m}. This determines the non-negative values a_{0}, a_{1}, b_{0}, b_{1} via

( \begin{matrix}a \\ b \end{matrix} ) = ( \begin{matrix}a_{0} & −a_{1} \\ b_{0} & b_{1} \end{matrix} ) ( \begin{matrix}2^{m} \\ 1 \end{matrix} ) , 0 < a_{1} ≤ 2^{m}, 0 ≤ b_{1} < 2^{m}.  (37)

Note that both tricks are incorporated in (37). Defining a_{0}, b_{0} as in (36) is the same as choosing a_{1} := 2^{m} − (a mod 2^{m}) and b_{1} := b mod 2^{m}. This choice ensures a_{0} > b_{0}, as we assume in the recursive call to the algorithm on a_{0}, b_{0}.

We are also given a regular matrix M. This determines the values a′_{0}, b′_{0}, a′_{1}, b′_{1}, a′, b′ via

( \begin{matrix}a′_{0} & a′_{1} \\ b′_{0} & b′_{1} \end{matrix} ):= M^{−1} ( \begin{matrix}a_{0} & −a_{1} \\ b_{0} & b_{1} \end{matrix} )  (38)

and

( \begin{matrix}a′ \\ b′ \end{matrix} ):= ( \begin{matrix}a′_{0} & a′_{1} \\ b′_{0} & b′_{1} \end{matrix} ) ( \begin{matrix}2^{m} \\ 1 \end{matrix} ) .  (39)

Hence we have the parallel reductions

( \begin{matrix}a_{0} \\ b_{0} \end{matrix} ) ⟶^{M} ( \begin{matrix}a′_{0} \\ b′_{0} \end{matrix} ) , ( \begin{matrix}a \\ b \end{matrix} ) ⟶^{M} ( \begin{matrix}a′ \\ b′ \end{matrix} ) .

Finally, we assume two key inequalities:

a′_{0} > b′_{0} ≥ 0  (40)
2∥a′_{0}∥ − 2 ≥ ∥a_{0}∥  (41)

Now write M as

M = ( \begin{matrix}p & q \\ r & s \end{matrix} ) ,  M^{−1} = δ ( \begin{matrix}s & −q \\ −r & p \end{matrix} )  (42)

where δ = det M = ±1. From (38) we obtain

a′_{1} = −δ(sa_{1} + qb_{1}), b′_{1} = δ(ra_{1} + pb_{1}).  (43)

The proof below uses (43) to predict the signs of a′_{1}, b′_{1}, assuming the sign of δ. This is possible thanks to the second trick.

The following is the integer analogue of lemma 9:

Lemma 11 (Partial Correctness Criteria) With the Basic Setup:

(–) Suppose det M = −1. (–a) ∥a′∥ = m + ∥a′_{0}∥ + ϵ_{1}, (0 ≤ ϵ_{1} < 1).

(–b) ∥b′∥ ≤ m + max{∥b′_{0}∥, ∥a_{0}∥ − ∥a′_{0}∥ + 1}.
Moreover, ∥a′∥ > ∥b′∥.

(+) Suppose det M = +1. (+a) ∥a′∥ = m + ∥a′_{0}∥ − ϵ_{2}, (0 ≤ ϵ_{2} < 1). (+b) ∥b′∥ ≤ 1 + m + max{∥b′_{0}∥, ∥a_{0}∥ − ∥a′_{0}∥ + 1}. Furthermore b′ ≥ 0.

In both cases (–) and (+), a′ > 0.

Proof. Since a_{0} = pa′_{0} + qb′_{0}, the ordering property (20) and (40) yields

∥a_{0}∥ = ∥p∥ + ∥a′_{0}∥ + ϵ_{3} (0 ≤ ϵ_{3} < 1).

Hence (41) is equivalent to

∥p∥ + ϵ_{3} ≤ ∥a′_{0}∥ − 2.

We now prove cases (–) and (+) in parallel.

Part (a). From equation (43),

∥a′_{1}∥ ≤ max{∥sa_{1}∥, ∥qb_{1}∥} + 1
< ∥p∥ + m + 1  (by (20), ∥a_{1}∥ ≤ m, ∥b_{1}∥ < m)
≤ ∥a′_{0}∥ + m − 1.

Hence ∥a′_{0}2^{m}∥ > 1 + ∥a′_{1}∥ and so a′ = a′_{0}2^{m} + a′_{1} > 0. This proves the desired a′ > 0. If δ = −1 then a′_{1} ≥ 0 (by equation (43)) and hence ∥a′∥ = m + ∥a′_{0}∥ + ϵ_{1} as required by subcase (–a). On the other hand, if δ = +1 then a′_{1} ≤ 0 and a′ = a′_{0}2^{m} + a′_{1} > a′_{0}2^{m−1} and hence subcase (+a) follows.

Part (b). Again from (43),

∥b′_{1}∥ ≤ max{∥ra_{1}∥, ∥pb_{1}∥} + 1
≤ ∥p∥ + m + 1
≤ ∥a_{0}∥ − ∥a′_{0}∥ + m + 1.

In case δ = +1, b′_{1} ≥ 0 and hence b′ = b′_{0}2^{m} + b′_{1} ≥ 0, as desired. Also subcase (+b) easily follows. In case δ = −1, b′_{1} ≤ 0 and b′_{0}b′_{1} ≤ 0. This gives the non-Archimedean inequality:

∥b′∥ ≤ max{∥b′_{0}2^{m}∥, ∥b′_{1}∥},

which proves subcase (–b).

Finally, we must show that δ = −1 implies ∥a′∥ > ∥b′∥: this follows immediately from (40),

(41) and subcases (–a) and (–b). Q.E.D.

To see inadequacies in the Partial Correctness Criteria, we state our precise algorithmic goal.

Integer HGCD Criteria: On input a > b ≥ 0, the integer HGCD algorithm outputs a regular matrix M such that

a ≤ 3 ⇒ M = E  (44)
a ≥ 4 ⇒ ∥a′∥ ≥ 1 + ⌈\frac{∥a∥}{2}⌉ > ∥b′∥, a′ > b′ ≥ 0  (45)

where ( \begin{matrix}a \\ b \end{matrix} ) ⟶^{M} ( \begin{matrix}a′ \\ b′ \end{matrix} ) .

Note that if a ≥ 4 then ∥a∥ ≥ 1 + ⌈∥a∥/2⌉ and so the desired matrix M exists.

Discussion. We see that the pair a′, b′ obtained in the Partial Correctness Criteria lemma may fail two properties needed for our HGCD algorithm: Case det M = −1: b′ can be negative. Case det M = +1: the inversion b′ ≥ a′ may occur. Clearly these two failures are mutually exclusive. On deeper analysis, it turns out that we only have to modify M slightly to obtain some regular matrix M^{∗} such that

( \begin{matrix}a \\ b \end{matrix} ) ⟶^{M^{∗}} ( \begin{matrix}a^{∗} \\ b^{∗} \end{matrix} )

and a^{∗}, b^{∗} satisfy the correctness criteria, ∥a^{∗}∥ ≥ m > ∥b^{∗}∥, a^{∗} > b^{∗} ≥ 0. The “Fixing Up lemma” below shows how to do this. The fixing up is based on three simple transformations of regular matrices: advancing, backing up and toggling.

In the following, let a > b > 0 and M = ⟨q_{1},…, q_{k}⟩ be a regular matrix such that

( \begin{matrix}a \\ b \end{matrix} ) ⟶^{M} ( \begin{matrix}a′ \\ b′ \end{matrix} ) .

(I) Advancing: If q′ = a′ div b′, then we say that M has advanced by one step to the matrix ⟨q_{1},…, q_{k}, q′⟩. Note that this operation defines a regular matrix iff q ≥ 1, i.e., ∥a′∥ ≥ ∥b′∥. In general, we may speak of advancing M by more than one step.

(II) Backing Up: We call the matrix ⟨q_{1},…, q_{k−i}⟩ the backing up of M by i steps (0 ≤ i ≤ k); in case i = 1, we simply call it the backup of M. To do backing up, we need to recover

the last partial quotient x from a regular matrix M = ( \begin{matrix}p & q \\ r & s \end{matrix} ). Note that M = E if

and only if q = 0, but in this case x is undefined. Hence assume M ≠ E. Then M is elementary if and only if s = 0, and in this case x = p. So we next assume that M is not elementary. Write

M′ = ( \begin{matrix}p′ & q′ \\ r′ & s′ \end{matrix} ) , M = M′ · ( \begin{matrix}x & 1 \\ 1 & 0 \end{matrix} )

where M′ ≠ E and p = xp′ + q′, q = p′. There are two cases. Case of q = 1: Clearly p′ = 1. Since p′ ≥ q′ ≥ 1 (ordering property), we must have q′ = 1. Hence x equals p − 1. Case of q > 1: Then p′ > q′ (there are two possibilities to check, depending on whether M′ is elementary or not). This implies x = p div q. In summary, the last partial quotient of M is given by

x = \{ \begin{matrix}undefined & if q = 0 \\ p & if s = 0 \\ p − 1 & if q = 1 \\ p div q & otherwise \end{matrix}

(III) Toggling: We call T = ( \begin{matrix}1 & 1 \\ 0 & −1 \end{matrix} ) the toggling matrix, so-called because T is idempo-

tent (T^{2} = E). The matrix MT is the toggle of M. We observe that MT is equal to ⟨q_{1},…, q_{k−1}, q_{k} − 1, 1⟩ in case q_{k} > 1, and MT = ⟨q_{1},…, q_{k−2}, q_{k−1} + 1⟩ in case q_{k} = 1 and k > 1. However, if q_{k} = 1 and k = 1, MT is not a regular matrix. In any case, we have

( \begin{matrix}a \\ b \end{matrix} ) ⟶^{MT} ( \begin{matrix}a′ + b′ \\ −b′ \end{matrix} ) .

Exercise A.1: Verify the remarks on the toggling matrix T.

Lemma 12 (Fixing Up) With the notations from the Basic Setup, let t be any number (the “fixup threshold”) such that

∥a′_{0}∥ ≥ t > max{∥b′_{0}∥, ∥a_{0}∥ − ∥a′_{0}∥ + 1}.  (46)

Moreover, if we write M as ⟨q_{1},…, q_{k}⟩ and M^{∗} is as specified below, then

∥a^{∗}∥ ≥ m + t > ∥b^{∗}∥  (47)

and

b^{∗} ≥ 0,  (48)

where ( \begin{matrix}a \\ b \end{matrix} ) ⟶^{M^{∗}} ( \begin{matrix}a^{∗} \\ b^{∗} \end{matrix} ). Here M^{∗} is the regular matrix specified as follows:

(–) Suppose det M = −1. (–A) If b′ ≥ 0 then M^{∗}:= M. (–B) Else if ∥a′ + b′∥ ≥ m + t then M^{∗} is the toggle of M. (–C) Else if q_{k} ≥ 2 then M^{∗}:=⟨q_{1},…, q_{k−1}, q_{k} − 1⟩ is the backup of the toggle of M. (–D) Else M^{∗} is the backing up of M by two steps.

\begin{matrix}c \\ ◯ \end{matrix} Chee-Keng Yap September 8, 1993 (+) Suppose det M = +1.

(+A) If ∥a′∥ ≤ ∥b′∥ then M^{∗} is the advancement of ⟨q_{1},…, q_{k−1}⟩ by at most two steps. (+B) Else if ∥a′∥ < m + t then M^{∗} is the backing up of M by one or two steps. (+C) Else M^{∗} is the advancement of M by at most two steps.

Proof. The Partial Correctness Criteria lemma will be repeatedly exploited. First assume det M = −1.

Subcase (–A). In this subcase, (48) is automatic, and (47) follows from case (–) of the Partial Correctness Criteria lemma.

Subcase (–B). In this subcase, M^{∗} = MT reduces ( \begin{matrix}a \\ b \end{matrix} ) to ( \begin{matrix}a^{∗} \\ b^{∗} \end{matrix} ) = ( \begin{matrix}a′ + b′ \\ −b′ \end{matrix} ), as noted

earlier. Recall that MT is not regular in case k = 1 and q_{1} = 1. But if this were the case then

( \begin{matrix}a \\ b \end{matrix} ) = ( \begin{matrix}1 & 1 \\ 1 & 0 \end{matrix} ) ( \begin{matrix}a′ \\ b′ \end{matrix} ) .

This implies a′+b′ = a > b = a′ and so b′ > 0, contradicting the fact that subcase (–A) has been excluded. Again, (48) is immediate, and (47) follows from case (–) of the Partial Correctness Criteria Lemma (∥a^{∗}∥ = ∥a′ + b′∥ ≥ m + t by assumption).

Subcase (–C). In this subcase, M^{∗} can be written as M ( \begin{matrix}1 & 0 \\ −1 & 1 \end{matrix} ), and so ( \begin{matrix}a^{∗} \\ b^{∗} \end{matrix} ) =

( \begin{matrix}a′ \\ a′ + b′ \end{matrix} ). We see that (48) holds by virtue of a′ + b′ > 0 (since ∥a′∥ > ∥b′∥, a′ > 0).

Also (47) holds because the Partial Correctness Criteria lemma implies ∥a′∥ ≥ m + t and, since subcase (–B) fails, ∥a′ + b′∥ < m + t.

Subcase (–D). Now q_{k} = 1 and M^{∗} omits the last two partial quotients (q_{k−1}, q_{k}) = (x, 1) where we write x for q_{k−1}. We ought to show k ≥ 2, but this is the same argument as in

subcase (–B). Hence M^{∗} = M ( \begin{matrix}1 & −x \\ −1 & x + 1 \end{matrix} ) and ( \begin{matrix}a^{∗} \\ b^{∗} \end{matrix} ) = ( \begin{matrix}a′(x + 1) + b′x \\ a′ + b′ \end{matrix} ). Hence

∥a^{∗}∥ = ∥xb^{∗} + a′∥ > ∥b^{∗}∥ and so a^{∗}, b^{∗} are consecutive elements of the remainder sequence of a, b. Then (48) holds because a′ + b′ > 0. To see (47), it is clear that m + t > ∥b^{∗}∥ (otherwise subcase (–B) applies) and it remains to show ∥a^{∗}∥ ≥ m + t. But this follows from ∥a^{∗}∥ = ∥a′ + x(a′ + b′)∥ > ∥a′∥ ≥ m + t.

Now consider the case det M = +1.

Subcase (+A). So there is inversion, a′ ≤ b′. Let us back up M to ⟨q_{1},…, q_{k−1}⟩:

( \begin{matrix}a \\ b \end{matrix} ) ⟶^{⟨q_{1},...,q_{k−1}⟩} ( \begin{matrix}a′′ \\ a′ \end{matrix} ) ⟶^{⟨q_{k}⟩} ( \begin{matrix}a′ \\ b′ \end{matrix} ) .

Hence a′′ = a′q_{k} + b′ > a′. Thus a′′, a′ are consecutive members of the remainder sequence of a, b. Now ∥a′′∥ ≥ ∥2a′∥ = ∥a′∥ + 1 > m + ∥a′_{0}∥ ≥ m + t. Also ∥a′∥ ≤ ∥b′∥ < 1 + m + t (by Partial Correctness Criteria). Therefore, if we advance ⟨q_{1},…, q_{k−1}⟩ by at most two steps, we would reduce a′′, a′ to a^{∗}, b^{∗} where ∥b^{∗}∥ < m + t.

Subcase (+B). Now a′ > b′ ≥ 0 and ∥a′∥ < m + t. So a′, b′ are consecutive members of the remainder sequence of a, b. Consider the entry a′′ = a′q_{k} + b′ preceding a′ in the remainder sequence of a, b. If ∥a′′∥ ≥ m + t, we are done since ∥a′′∥, ∥a′∥ straddle m + t. Otherwise, consider the entry a′′′ = a′′q_{k−1} + a′ preceding a′′. We have

∥a′′′∥ = ∥(a′q_{k} + b′)q_{k−1} + a′∥ ≥ ∥a′∥ + 1 ≥ m + t.

Thus ∥a′′′∥, ∥a′′∥ straddle m + t.

Subcase (+C). Again a′, b′ are consecutive members of the remainder sequence with ∥a′∥ ≥ m+t. But ∥b′∥ − 1 < m + t implies that if we advance M by at most two steps, the pair a′, b′ would be reduced to a^{∗}, b^{∗} where ∥b^{∗}∥ < m + t.

Q.E.D.

It is interesting to note that in tests on randomly generated numbers of about 45 digits, subcase (–D) never arose.

The procedure Fixup. The Fixing Up lemma and its proof provides a specific procedure to convert the tentative output matrix M of the HGCD algorithm into a valid one. To be specific, let

Fixup(M, a, b, m, t)

denote the subroutine that returns M^{∗}, as specified in the Fixing Up lemma:

( \begin{matrix}a \\ b \end{matrix} ) ⟶^{M^{∗}} ( \begin{matrix}a^{∗} \\ b^{∗} \end{matrix} ) , ∥a^{∗}∥ ≥ m + t > ∥b^{∗}∥.

The correct behavior of the Fixup procedure depends on its input parameters fulfilling the conditions of the Fixing Up lemma. In particular, it must fulfil the conditions of the Basic Setup (mainly the inequalities (40) and (41)) and also the inequality (46).

In a typical invocation of Fixup(M, a, b, m, t), the values M, a, b, m are available as in the Basic Setup. To pick a value of t, we use the fact that the following typically holds:

∥a′_{0}∥ ≥ 1 + ⌈∥a_{0}∥/2⌉ > ∥b′_{0}∥  (49)

(cf. (45)). In this case, it is easily verified that the choice t = 1 + ⌈∥a_{0}∥/2⌉ will satisfy (46). Of course inequality (49) also implies inequality (41).

However, our Fixup procedure may also be called in a situation when the Fixing Up lemma does not hold, namely, when a_{0} = 1 + (a div 2^{m}) ≤ 3. In this case, no choice of t satisfying inequality (46) is possible. Note that b_{0} < a_{0} ≤ 3 implies b = b_{0}2^{m} + b_{1} < 3 · 2^{m}. It is easy to check that if we take at most three of the usual Euclidean remaindering steps, we reduce a, b to a^{∗}, b^{∗} where ∥a^{∗}∥ ≥ m > ∥b^{∗}∥. In such a situation, if we assume that the Fixup procedure is called with M = E and t = 0, the returned matrix M^{∗} is the advancement of E by at most three steps. More generally, if ∥a∥, ∥b∥ straddle m + i where i = 0, 1, 2, and we call Fixup with the arguments

Fixup(E, a, b, m, 0),

we say this is the “easy fixup” case, because M^{∗} is the advancement of E by at most 4 steps.

We present the integer HGCD algorithm.

Algorithm Integer HGCD(a, b): Input: integers a, b with a > b ≥ 0. Output: a regular matrix M satisfying the integer HGCD criteria (44) or (45).

[1] m ← 1 + ⌈\frac{∥a∥}{2}⌉; \{this is the magic threshold\}

if a ≤ 3 or ∥b∥ < m then return(E);

[2] a_{0} ← 1 + (a div 2^{m}); b_{0} ← b div 2^{m}; if a_{0} ≤ 3 then t ← 0 else t ← 1 + ⌈\frac{∥a_{0}∥}{2}⌉; R ← Fixup(Hgcd(a_{0}, b_{0}), a, b, m, t);

( \begin{matrix}a′ \\ b′ \end{matrix} ) ← R^{−1} ( \begin{matrix}a \\ b \end{matrix} );

[3] if ∥b′∥ < m then return(R);

[4] q ← a′ div b′; ( \begin{matrix}c \\ d \end{matrix} ) ← ( \begin{matrix}b′ \\ a′ mod b′ \end{matrix} );

if 1 + (c div 2^{m}) ≤ 3 then return(R · Fixup(E, c, d, m, 0));

[5] l ← ⌈∥c∥⌉; k ← 2m − l − 1; {Now ∥c∥ ≥ m + 1 ≥ 4. We claim ∥c∥ − 1 ≥ k ≥ 0}

[6] c_{0} ← 1 + (c div 2^{k}); d_{0} ← d div 2^{k};

if c_{0} ≤ 3 then t′ ← 0 else t′ ← 1 + ⌈\frac{∥c_{0}∥}{2}⌉;

S ← Fixup(Hgcd(c_{0}, d_{0}), c, d, k, t′); {We claim k + t′ = m + 1.}

[7] ( \begin{matrix}c′ \\ d′ \end{matrix} ) ← S^{−1} ( \begin{matrix}c \\ d \end{matrix} ); \{So ∥c′∥, ∥d′∥ straddle k + t′\}

T ← Fixup(E, c′, d′, m, 0);
M ← R · ⟨q⟩ · S · T; return(M);

Correctness. HGCD returns in four places (steps [1],[3],[4] and [7]) in the algorithm. We show that the matrix returned at each of these places is correct. Since these matrices are regular, we basically have to check the straddling property (45) when a ≥ 4. We will also need to check that each call to the Fixup procedure is proper.

a) In case the algorithm returns the identity matrix E in step [1], the correctness is trivial.

b) In step [2], we must check that the proper conditions are fulfilled for calling Fixup. When a_{0} ≤ 3 we have the “easy fixup” case. Otherwise a_{0} ≥ 4 and the first recursive call to HGCD returns some regular matrix R′ which is fixed up as R by Fixup. The conditions of the Basic

Setup are fulfilled with a, b, m as usual and M = R′. If ( \begin{matrix}a_{0} \\ b_{0} \end{matrix} ) ⟶^{R′} ( \begin{matrix}a′_{0} \\ b′_{0} \end{matrix} ), then inductively,

the correctness of the HGCD procedure implies equation (49) (and hence (41)) holds. As discussed following equation (49), the choice t = 1 + ⌈∥a_{0}∥/2⌉ then satisfies (46),

c) Suppose the matrix R is returned at step [3]. This is correct since ∥a′∥, ∥b′∥ straddle m + t (by correctness of the Fixup procedure) and the condition for exit is ∥b′∥ < m.

d) In step [4], the call to Fixup is the “easy fixup” case since ∥c∥ ≥ m and ∥c∥ ≤ m + 2. The returned matrix is clearly correct.

e) Suppose we reach step [5]. We show the claim ∥c∥ − 1 ≥ k ≥ 0. We have ∥c∥ − 1 ≥ m − 1 ≥ (m − 1) + (m − l) = k. Next, k ≥ 0 is equivalent to 2m − 1 ≥ l. This follows from:

l = ⌈∥b′∥⌉ ≤ m + t  (from the first FIXUP)

= m + 1 + ⌈\frac{∥1 + (a div 2^{m})∥}{2}⌉
≤ m + 1 + ⌈\frac{1 + ⌊∥a div 2^{m}∥⌋}{2}⌉ (since ∥1 + x∥ ≤ 1 + ⌊∥x∥⌋)
= m + 1 + ⌈\frac{1 + ⌊∥a/2^{m}∥⌋}{2}⌉ (since ⌊∥x div y∥⌋ = ⌊∥x/y∥⌋)
= m + 1 + ⌈\frac{1 + ⌊∥a∥⌋ − m}{2}⌉
≤ m + ⌈(m − 1)/2⌉ (since m = 1 + ⌈∥a∥/2⌉)
≤ 2m − 1  (m ≥ 2)

f) The call to Fixup in step [6] fulfills the appropriate conditions. [Reasoning as before: note that ∥c∥−1 ≥ k implies that c_{0} ≥ 3. Hence, the “easy fixup” case occurs iff c_{0} = 3. Otherwise, the Basic Setup conditions prevail with a, b, m, M in the Basic Setup replaced by c, d, k, Hgcd(c_{0}, d_{0}).] Next we prove the claim k + t′ = m + 1:

t′ = 1 + ⌈∥c_{0}∥/2⌉
= 1 + ⌈\frac{⌈∥c_{0}∥⌉}{2}⌉
= 1 + ⌈\frac{⌈∥ϵ + (c/2^{k})∥⌉}{2}⌉ (0 < ϵ ≤ 1, c_{0} = 1 + ⌊c/2^{k}⌋)
= 1 + ⌈\frac{l − k + δ}{2} ⌉  (δ = 0 or 1, l = ⌈∥c∥⌉)
= 1 + ⌈\frac{2l − 2m + 1 + δ}{2} ⌉ (k = 2m − l)
= 2 + l − m.

Thus k + t′ = k + (2 + l − m) = m + 1, as desired.

g) We reach step [7]. By the correctness of the Fixup procedure, ∥c′∥, ∥d′∥ straddle k+t′ = m+1. Hence we have the right conditions for the “easy fixup” case. The final output is clearly correct.

This concludes our correctness proof. Q.E.D.

Computational details and analysis. The algorithm has to perform comparisons of the kind

∥a∥ : m,

and compute ceiling functions in the special forms

⌈∥a∥⌉, ⌈∥a∥/2⌉,

where a, m are positive integers. (In the algorithm a may be zero, but we can treat those as special cases.) Since ∥a∥ is generally not rational, we do not want to explicitly compute it. Instead we reduce these operations to integer comparisons, checking if a is a power of two, and to computing the function bit(a), which is defined to be the number of bits in the binary representation of a positive integer a. So bit(a) = 1 + ⌊log_{2} a⌋ and clearly this function is easily computed in linear time in the usual Turing machine model. Then we have

∥a∥ ≥ m ⇔ bit(a) − 1 ≥ m

and

∥a∥ > m ⇔ \{ \begin{matrix}bit(a) − 1 & > & m & if a is a power of 2 \\ bit(a) & > & m & else. \end{matrix}

and finally,

⌈\frac{∥a∥}{2}⌉ = \{ \begin{matrix}⌈\frac{bit(a)−1}{2}⌉ & if a is a power of 2 \\ ⌈\frac{bit(a)}{2}⌉ & else \end{matrix}

The global structure of the complexity analysis is similar to the polynomial HGCD case: with M_{B}(n) = nL(n) denoting as usual the bit complexity integer multiplication, it is not hard to see that Fixup takes time O(M_{B}(n)), under the conditions stipulated for its invocation. In the two recursive calls to HGCD, it is easy to check that the integers has bit size \frac{n}{2} + O(1). Hence, if T (n) is the bit complexity of our HGCD algorithm on inputs of size at most n, then

T (n) = O(M_{B}(n)) + 2T (\frac{n}{2} + O(1)).

This has solution T (n) = O(M_{B}(n) log n) = nL^{2}(n).

Exercise A.2:

a) Verify the remarks on reducing operations involving ∥a∥ to integer operations, the function bit(a) and testing if a is a power of 2.

b) Derive the time bound T (n) = nL^{2}(n) for the HGCD algorithm.

Exercise A.3: Try to simplify the integer HGCD algorithm by separating the truncation value t (as in a_{0} := a div 2^{t}) from the straddling value s (as in ∥a′∥ ≥ s > ∥b′∥). Currently

t = s = 1 + ⌈∥a∥/2⌉.

References

[1] A. V. Aho, J. E. Hopcroft, and J. D. Ullman. The design and analysis of computer algorithms. Addison-Wesley, 1974.

[2] R. P. Brent, F. G. Gustavson, and D. Y. Y. Yun. Fast solution of Toeplitz systems of equations and computation of Pade approximants. J. Algorithms, 1:259–295, 1980.

[3] E. Brieskorn and H. Knörrer. Plane Algebraic Curves. Birkhäuser Verlag, 1986.

[4] D. A. Buell. Binary Quadratic Forms: classical theory and modern computations. Springer-Verlag, 1989.

[5] E. Kaltofen and H. Rolletschek. Computing greatest common divisors and factorizations in quadratic number fields. Math. Comp., 52:697–720, 1989.

[6] D. E. Knuth. The analysis of algorithms. In Actes du Congrés International des Mathématiciens, pages 269–274, Nice, France, 1970. Gauthier-Villars.

[7] R. Moenck. Fast computations of GCD’s. Proc. 5th ACM Symposium on Theory of Computation, pages 142–171, 1973.

[8] P. Samuel. About Euclidean rings. J. Algebra, 19:282–301, 1971.

[9] Arnold Schönhage. Schnelle Berechnung von Kettenbruchentwicklungen. Acta Informatica, 1:139–144, 1971.

[10] Volker Strassen. The computational complexity of continued fractions. SIAM Journal Computing, 12:1–27, 1983.

[11] Klaus Thull and Chee Yap. A unified approach to fast GCD algorithms for polynomials and integers. Technical report, NYU-Courant Institute, Robotics Laboratory, 1992.

[12] Oscar Zariski and Pierre Samuel. Commutative Algebra, volume 1. Springer-Verlag, New York, 1975.