On computing factors of cyclotomic polynomials

For integer n > 0 let Φ n (x) denote the cyclotomic polynomial where ζ is a primitive n-th root of unity. Clearly and the Möbius inversion formula [13] gives Equation ( 1) is useful for theoretical purposes, but (3) is more convenient for computation as it leads to a simple algorithm for computing the coefficients of Φ n (x) or evaluating Φ n (x) at integer arguments using only integer arithmetic. If n is square-free the relations Φ n (x) = x -1, if n = 1; Φ n/p (x p )/Φ n/p (x), if p|n, p prime; (4) give another convenient recursion for computing Φ n (x). Although Φ n (x) is irreducible over Z (see for example [29]), Φ n (x) may be reducible over certain quadratic fields. For example, 4Φ 5 (x) = (2x 2 + x + 2) 2 -5x 2 , (5) so Φ 5 (x) has factors x 2 + 1± √ 5 2 x + 1 whose coefficients are algebraic integers in Q[ √ 5]. For odd square-free n > 1 the cyclotomic polynomial Φ n (x) satisfies the identity Gauss [12] proved (6) for odd prime n; the generalization to other odd square-free n is due to Dirichlet [11]. Related identities of Aurifeuille and Le Lasseur [2] are for odd square-free n, and for even square-free n > 2. For a proof, see Lucas [22] or Schinzel [26]. In (6)(7)(8), A n (x), . . . , D n (x) are polynomials with integer coefficients, and without loss of generality we can assume that A n (x)/2, B n (x), C n (x) and D n (x) are monic. In Section 3 we show how the coefficients of A n , . . . , D n can be computed by simple algorithms which require O(n 2 ) arithmetic operations and work entirely over the integers. In Section 1.1 we summarize our notation for future reference. Some numerical examples are given in Sections 1.2-1.3, and Newton's identities are discussed in Section 1.4. Then, in Section 2, we discuss the theoretical basis for the algorithms. The results for A n and B n are known (though perhaps forgotten) -they may be found in Dirichlet [11]. We present them in Section 2.2 for the sake of completeness and to aid the reader in understanding the results for C n and D n . The algorithms are presented in Section 3. The algorithm (Algorithm D) for computing A n and B n is essentially due to Dirichlet [11], who illustrated it with some numerical examples but did not state it in general terms. The algorithm (Algorithm L) for computing C n and D n appears to be new. In Section 3.3 we comment briefly on Stevenhagen's algorithm [27] and compare it with Algorithm L. Finally, in Section 4 we give some explicit formulas for A n (x), . . . , D n (x). These may be regarded as generating functions if x is an indeterminate, or may be used to compute A n (x), . . . , D n (x) for given argument x. In the special case x = 1 the results for A n (1), B n (1) reduce to known formulas involving the class number of the quadratic field Q[ √ ±n]. One application of cyclotomic polynomials is to the factorization of integers of the form a n ± b n : see for example [6,7,8,9,15,16,24,26,27]. If x = m2 n for any integer m, then (7-8) are differences of squares, giving rational integer factors of x n ± 1. Examples may be found in Section 4.4. For the reader interested in integer factorization, our most significant results are Algorithm L of Section 3.2 and Theorem 3 of Section 4.4. Unless qualified by "algebraic", the term "integer" means a rational integer. x usually denotes an indeterminate, occasionally a real or complex variable. µ(n) denotes the Möbius function, φ(n) denotes Euler's totient function, and (m, n) denotes the greatest common divisor of m and n. For definitions and properties of these functions, see for example [13]. Note that µ(1) = φ(1) = 1. (m|n) denotes the Jacobi symbol1 except that, as is usual for the Kronecker symbol 2 , (m|n) is defined as 0 if (m, n) > 1. Thus, when specifying a condition such as (m|n) = 1 we may omit the condition (m, n) = 1. n denotes a positive integer (square-free from Section 2.2 on). For given n, we define integers n ′ , s and s ′ as follows: It is convenient to write g k for (k, n) and g ′ k for (k, n ′ ). Define Thus we can write (7-8) as The factor (-1) φ(n/2) in the definition of F n is only relevant if n = 2, and ensures that (10) is valid for n = 2 (with C 2 (x) = x + 1, D 2 (x) = 1). The Aurifeuillian factors of F n (x) are and From (10) we have F n (x) = F - n (x)F + n (x). We may write F ± n for one of F + n , F - n . We sometimes need to specify a particular complex square root. If m < 0 then √ m means i |m|. d is usually the degree of a polynomial, while D is the discriminant of a quadratic form. For odd square-free n we always have D = sn, so D = 1 mod 4. Some additional notation is introduced in Section 1.4. Taking n = 15, we have and the reader may easily verify that ( 6) and ( 7) are satisfied. As an example of (8), for n = 14 we have and D 14 (x) = x 5 + 2x 4 -x 3 -x 2 + 2x + 1. Taking n = 5 in (7) we obtain and by division (taking the factors with opposite signs of √ 5) we obtain factors This is not of the form (10) because it gives a factorization of Φ 15 (±x) over (11) is an example of the more general identities of Beeger [3] and Schinzel [26]. These identities can all be obtained in a similar manner from (10), so in the application to integer factorization they do not give any factors which could not be found from several applications of (10) and some greatest common divisor calculations. For this reason we have restricted our attention to identities of the form ( 6) and (10). be a polynomial of degree d with arbitrary roots ξ j and coefficients a 0 = 1, a 1 , . . . , a d . For k > 0, define Newton (1707) 3 showed how to express the elementary symmetric functions a 1 , a 2 , . . . in terms of the sums of powers p 1 , p 2 , . . . We may find a 1 , . . . , a d by solving a lower triangular linear system of special form [28]. Writing the solution explicitly in the form of a linear recurrence, we have for k = 1, . . . , d. An alternative expression for a k as a determinant may be obtained by applying Cramer's rule to the lower triangular system. However, for computational purposes ( 13) is more convenient. In Section 4 we use the following generating function [25] for (a 0 , a 1 , . . .): Differentiating both sides of ( 14) and equating coefficients shows that ( 13) and ( 14) are formally equivalent. An independent proof of ( 14) is the following: for sufficiently small x we have ln( In all our applications of ( 14) the p j are bounded, so the infinite series converges for |x| < 1. In the following, p j and a j are not fixed, but depend on the particular polynomial under consideration at the time. This should not cause any confusion. Our idea is to compute sums of powers of certain roots of the polynomials occurring on the left side of ( 6) and (10), and then use Newton's identities in the form (13) to compute the coefficients of A n , . . . , D n . First consider the computation of the coefficients of the cyclotomic polynomial Φ n (x) for n > 1. This is presented to illustrate a simple case of the technique; in practice it is more efficient to compute Φ n (x) from (3). Let ζ be a primitive n-th root of unity. To apply Newton's identities we need to evaluate for k = 1, 2, . . . , φ(n). This problem is well-known 4 . If n is prime the problem is easy: from then the same argument as before shows that p k = p 1 . If g k > 1 then the sum (15) defining p k consists of φ(n)/φ(n/g k ) copies of a sum of primitive (n/g k )-th roots of unity. Thus, the result is Using ( 16) the coefficients a 1 , . . . , a φ(n) of Φ n (x) may be evaluated from the recurrence (13). As an application of ( 14) and ( 16) we prove two Lemmas which give upper bounds on |Φ n (x)| and |F n (x)| for complex x outside the unit circle. Here F n (x) is the modified cyclotomic polynomial defined by ( 9). Lemma 2 is used in Section 4.4. From ( 14) with x replaced by 1/x, we have but from ( 16) we have This completes the proof. ✷ in both cases. The bound on |F n (x)| follows from Lemma 1 applied to Φ n (±x) if n is odd, and from Lemma 1 applied to Φ n/2 (-x 2 ) if n is even. ✷ From now on we assume that n > 1 is square-free. In this subsection we also assume that n is odd. Consider the polynomial of degree φ(n)/2, where ζ = e 2πi/n is a primitive n-th root of unity. The particular choice of primitive root is only significant for the sign of the square root of sn appearing in the equations below. From Dirichlet [11], where A n and B n are as in (6), and s is defined in Section 1.1. Since n is odd, we have s = (-1|n) = (-1) so Gauss's identity (6) may be written as The sums of k-th powers of roots of G n (x) are Let The result (22) is essentially due to Dirichlet [11], but we sketch a proof. If g k = 1, then (22) follows from the discussion in Section 2.1 (where p k has a different meaning!) and the classical result that the Gaussian sum 0 1, observe that (g k , n/g k ) = 1 because n is square-free. Thus, we can write the summation index j in (21) in the form j = j 0 g k + j 1 (n/g k ), and (j|n) = (j 1 |g k )(j 0 |(n/g k )). Since jk = j 0 g k k mod n, ζ jk is independent of j 1 , and it follows that the sum (21) defining p k consists of φ(g k )/2 copies of a complete sum of primitive (n/g k )-th roots of unity. Thus, (22) follows as in the proof of (16). Although (22) has been written with two cases for the sake of clarity, our convention that (k|n) = 0 if (k, n) > 1 implies that the expression is valid in both cases. Similarly for Gn (x), with the sign of √ sn in ( 23) reversed. Observe that Using (22), the coefficients of G n (x), and hence of A n (x) and B n (x), may be evaluated from the recurrence (13). Moreover, it is possible to perform the computation using only integer arithmetic. Details are given in Section 3. Here we assume that n > 1 is square-free, but not necessarily odd. Recall the definitions of n ′ , s and s ′ from Section 1.1. Let ζ = e πi/n ′ be a primitive 2n ′ -th root of unity. The particular choice of primitive root is only significant for the sign of the square root in (29). Consider the polynomial where Observe that L n (x) has degree φ(n ′ ) = φ(2n). Also, j ∈ S n iff 2n ′ -j ∈ S n , so the coefficients of L n (x) are real. In fact, from (29) below, they are in Q[ √ n]. We later use the fact that L n (x) is symmetric. Schinzel [26] essentially shows (with a different notation) that where C n (x) and D n (x) are the polynomials of (10). Define so after a change of variable (10) may be written as Clearly Observe that the cosine in ( 29) is 0 or ±1, and depends only on n mod 4 and k/2 mod 4. The proof of ( 29) is similar to that of (22), but tedious because of the number of cases to be considered. Thus, we omit the details. Using (29), the coefficients of L n (x), and hence of C n (x) and D n (x), may be evaluated from the recurrence (13). Details are given in Section 3. In this section we use the analytic results of Section 2 to derive efficient algorithms for computing the coefficients of the polynomials A n , . . . , D n . Consider the computation of A n and B n for odd square-free n. Our notation is the same as in Section 2.2. Write Recall the definition (21) of p k . For k > 0 we have, from (23), where q k and r k are integers given by and Using ( 13) and (18), we obtain the recurrences and for k = 1, 2, . . . , d. The algorithm is now clear: Algorithm D (for Dirichlet) 1. Evaluate q k and r k for k = 1, . . . , d using (30-31). 2. Set α 0 ← 2 and β 0 ← 0. 3. Evaluate α k and β k for k = 1, . . . , d using (32-33). 1. (32-33) should give exact integer results; in practice a sum not divisible by 2k is a symptom of integer overflow. 2. The operation count can be reduced by a factor of close to four if advantage is taken of the following properties of A n and B n : A n is anti-symmetric if its degree d = φ(n)/2 is odd, otherwise A n is symmetric (except for the trivial case A 3 (x) = 2x + 1). Thus, we may use B n /x is antisymmetric if n is composite and n = 3 mod 4, otherwise B n /x is symmetric. Thus, we may use Using these properties, the recurrences (32-33) need only be applied for k ≤ max(1, ⌊d/2⌋). Consider the case n = 15 as in Section 1.2. We have s = -1, d = φ(15)/2 = 4. Thus q 3 , q 4 , r 3 , and r 4 are not required if we use symmetry. The initial conditions are α 0 = 2 and β 0 = 0. The recurrences (32-33) give Using symmetry of A n (x) and anti-symmetry of B n (x)/x, or continuing with the recurrences (32-33), we obtain 26) it is enough to compute the coefficients a k of L n (x). In order to work over the integers, we define where p k is the sum of k-th powers of roots of L n (x). Thus, from (29), then, from (26), In particular, γ 0 = δ 0 = 1. Using ( 13) we obtain and We may use the fact that C n (x) and D n (x) are symmetric to reduce the number of times the recurrences (35-36) need to be applied. An algorithm which incorporates this refinement is: Algorithm L (for Lucas) 1. Evaluate q k for k = 1, . . . , d using (34). 2. Set γ 0 ← 1 and δ 0 ← 1. Consider the case n = 15 as in Section 1.2. We have n ′ = 2n = 30, s ′ = 1, d = φ(30)/2 = 4. Thus q 1 = (15|1) = 1, The initial conditions are γ 0 = δ 0 = 1. The recurrences (35-36) give Using symmetry we obtain γ 3 = γ 1 = 8, γ 4 = γ 0 = 1, δ 2 = δ 1 = 3, and δ 3 = δ 0 = 1. Thus C 15 (x) = x 4 + 8x 3 + 13x 2 + 8x + 1 and D 15 (x) = x 3 + 3x 2 + 3x + 1, as expected. Stevenhagen [27] gives an algorithm for computing the polynomials C n (x) and D n (x). His algorithm depends on the application of the Euclidean algorithm to two polynomials with integer coefficients and degree O(n). C n (x) and D n (x) may be computed as soon as a polynomial of degree ≤ φ(n)/2 is generated by the Euclidean algorithm. Thus, the algorithm requires O(n 2 ) arithmetic operations, the same order5 as our Algorithm L. Unfortunately, Stevenhagen's algorithm suffers from a well-known problem of the Euclidean algorithm [14] -although the initial and final polynomials have small integer coefficients, the intermediate results grow exponentially large. When implemented in 32-bit integer arithmetic we found that Stevenhagen We now use (14) to give generating functions for the coefficients of A n , . . . , D n . The generating functions can be used to evaluate the coefficients of A n (x), . . . , D n (x) in O(n log n) arithmetic operations, via the fast power series algorithms of Section 5 of Brent and Kung [5]. Also, where the generating functions converge, they give explicit formulas which can be used to compute A n (x), . . . , D n (x) at particular arguments x. However, it is often more efficient to compute the coefficients of the polynomials by the algorithms of Section 3 and then evaluate the polynomials by Horner's rule. The generating functions may be written in terms of certain analytic functions f n and g n , which we now define. For odd square-free n > 1 and |x| ≤ 1, define Similarly, for square-free n > 1 and |x| ≤ 1, define Observe that g n (x) is an odd function, so g n (-x) = -g n (x). It follows from ( 46) and (57) below that exp( √ snf n (x)) and exp(2 √ ng n (x)) are rational functions with zeros and poles at certain roots of unity. From these representations it follows that the analytic continuations outside the unit circle are given by (39) and The functions f n (x) and g n (x) are closely related. For example, taking the odd terms in the sum (37) and using the law of quadratic reciprocity, we obtain Such identities are a consequence of relationships between the polynomials G n (x) and L n (x). f n (1) is related to the class number h(D) of the quadratic field Q[ √ D] with discriminant D = sn. In the notation of Davenport [10], f n (1) = L -1 (1) = L(1) = L(1, χ), where χ(j) = (j|n) is the real, nonprincipal Dirichlet character appearing in (37). Known results [10,18,19,30] in the case n = 3 mod 4 (so Here In the case n = 1 mod 4, we have D = n, and Here ε is the "fundamental unit", i. Let n > 3 be odd and square-free. We exclude n = 3 to avoid the special case in (42), but the results apply with minor modifications when n = 3. Let s, G n , Gn be as in Section 2.2, and From ( 14) and ( 23), we have Also, from (44), where If n = 1 mod 4 then by grouping the terms for j and n -j (j < n/2) in (49) we have n|σ. If n = 3 mod 4 then, from (41), we have σ = -nh(-n) so again n|σ. Thus, in both cases ζ σ = 1, and from (48) we have We see from ( 46) or (51) that, as claimed above, exp( √ snf n (x)) is a rational function. It has zeros at ζ j , (j|n) = -s; and poles at ζ j , (j|n) = +s. From ( 20) and ( 51), taking a square root, we obtain If ( 52) is interpreted as a generating function for G n (x) then Φ n (x) and f n (x) should be interpreted as power series in x, and the correct sign of the square root is positive. On the other hand, if (52) is regarded as an exact expression for G n (x), then the sign of the square root is positive for real x, because G n (x) and Φ n (x) have no real roots, and the exponential never vanishes, so a change in sign would contradict the continuity of G n (x). An extension of this argument shows that the same branch of the square root must be taken in any simply-connected, closed region which does not contain any of the zeros of Φ n (x). (We omit similar comments below.) From (52) we easily deduce the corresponding expressions for A n (x) = G n (x) + Gn (x) and B n (x) = ( Gn (x) -G n (x))/ √ sn. We state the results as a Theorem: Theorem 1 For odd, square-free n > 3, the polynomials A n (x) and B n (x) occurring in Gauss's identity (6) are and If n = 3 mod 4, so s = -1, then it is natural to replace cosh(iz) by cos(z) in (53) and sinh(iz) by i sin(z) in (54), giving and Consider the case n = 15. We expand the right side of (55) as a power series in x, keeping enough terms to find A 15 (x) without using symmetry. From Section 1.2 we have ) as expected. We can ignore the O(x 5 ) term (which in fact vanishes) since we know that deg A 15 (x) = φ(15)/2 = 4. A reader who attempts similar computations for larger n will soon be convinced that Algorithm D of Section 3.1 is more convenient, if only because all intermediate results are integers and there is an easy check on the accuracy of the inner product accumulations. Using (41-43) and the fact (an immediate consequence of ( 4)) that We now consider the analogues of (53-54) for the Lucas polynomials C n and D n . The argument is similar to that of Section 4.2, but simpler because the polynomial L n (x) is symmetric, which leads to the simple functional equation (40) for g n (x). Assume that n > 1 is square-free, and adopt the notation of Section 2.3. Using (29) in the generating function (14), we have This shows that exp(2 we deduce the functional equation (40), which gives the analytic continuation of g n (x) outside the unit circle. We may also write (57) more simply as From ( 28) and (58), taking a square root, we obtain Theorem 2 Let n > 1 be square-free. The Aurifeuillian factors of F n (x) are given by and Recall that the Aurifeuillian factors . Thus (60) follows from (59). Since we easily deduce (61-62). ✷ In this section we illustrate how the results of Sections 3.2 and 4.3 can be used to obtain factors of integers of the form a n ± b n . Our examples are for illustrative purposes, so are small enough to be verified by hand. Many larger examples can be found in [6,7]. As usual, n > 1 is a square-free integer. Recall the definition of F n (x) from Section 1.1. Note that the polynomial F n (x) is a factor of x n ± 1 (where the sign is "-" if n = 1 mod 4, and "+" otherwise). If x has the form m 2 n, where m is a positive integer, then √ nx = mn is an integer, and the Aurifeuillian factors F ± n (x) = C n (x) ± mnD n (x) give integer factors of F n (x), and hence of x n ± 1 = m 2n n n ± 1. For example, if m = n k , we obtain factors of n (2k+1)n ± 1. More generally, if m = p/q is rational, we obtain rational factors of x n ± 1 = p 2n q -2n n n ± 1, and thus integer factors of p 2n n n ± q 2n . We consider one example later, but for the moment we continue to assume that m is an integer. Before giving numerical examples, we state explicitly how the results of Section 4.3 can be used to compute F ± n (m 2 n). Theorem 3 Let m, n be positive integers, n > 1 square-free, x = m 2 n, and λ = φ(2n)/2 . Then the Aurifeuillian factors F ± n (x) of F n (x) are given by To conclude, we give an example where m = p/q is rational but not an integer. Consider n = 7, p = 2, q = 5, so x = p 2 n/q 2 = 28/25 and √ nx = pn/q = 14/5. We have Theorem 3 is not applicable. Because x is close to 1, the series for g 7 (1/ √ x) converges rather slowly, and we need to take at least 35 terms to obtain sufficient accuracy. However, using Algorithm L we easily find that See, for example, Riesel[24]. To avoid ambiguity, we never write the Jacobi symbol as m n . Note that m|n without parentheses means that m divides n. See, for example, Landau[18]. See Turnbull[28], where the notation s k is used in place of our p k . It would be confusing to use Turnbull's notation because we have used s and s ′ for other purposes. Note that our p k are not generally prime numbers. If ζ = e 2πi/n then our p k is "Ramanujan's sum" cn(k), in the notation of Ramanujan[23] or Chapter 26 of Davenport[10]. The complexity of both algorithms can be reduced to O(n(log n) 2 ) arithmetic operations by standard "divide and conquer" techniques[1,4], but this is not of practical significance.