On computing factors of cyclotomic polynomials^{∗}^{†}

Richard P. Brent Computer Sciences Laboratory Australian National University Canberra, ACT 0200

                      rpb@cslab.anu.edu.au

Report TR-CS-92-13 September 1992

In memory of Derrick H. Lehmer 1905–1991

Abstract

For odd square-free n > 1 the cyclotomic polynomial Φ_{n}(x) satisfies the identity of Gauss

4Φ_{n}(x) = A_{n}^{2} − (−1)^{(n−1)/2}nB_{n}^{2}.

A similar identity of Aurifeuille, Le Lasseur and Lucas is

Φ_{n}((−1)^{(n−1)/2}x) = C_{n}^{2} − nxD_{n}^{2}

or, in the case that n is even and square-free,

±Φ_{n/2}(−x^{2}) = C_{n}^{2} − nxD_{n}^{2},

Here A_{n}(x), …, D_{n}(x) are polynomials with integer coefficients. We show how these coefficients can be computed by simple algorithms which require O(n^{2}) arithmetic operations and work over the integers. We also give explicit formulae and generating functions for A_{n}(x), …, D_{n}(x), and illustrate the application to integer factorization with some numerical examples.

1991 Mathematics Subject Classification. Primary 11-04, 05A15; Secondary 11T06, 11T22, 11T24, 11Y16, 12-04, 12E10, 12Y05,

Key words and phrases. Aurifeuillian factorization, class number, cyclotomic field, cyclotomic polynomial, Dirichlet series, exact computation, Gauss’s identities, generating functions, integer factorization, Lucas’s identities, Newton’s identities.

∗To appear in Mathematics of Computation D. H. Lehmer memorial issue, 1993.

†Copyright \begin{matrix}c \\ ◯ \end{matrix} 1992, R. P. Brent. rpb135 typeset using LaTEX

1 Introduction

For integer n > 0 let Φ_{n}(x) denote the cyclotomic polynomial

Φ_{n}(x) = ∏_{{0<j≤n \atop (j,n)=1}}(x − ζ^{j}),  (1)

where ζ is a primitive n-th root of unity. Clearly

x^{n} − 1 = ∏_{d|n} Φ_{d}(x),  (2)

and the Möbius inversion formula [13] gives

Φ_{n}(x) = ∏_{d|n}(x^{d} − 1)^{μ(n/d)}.  (3)

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) = \{ \begin{matrix}x − 1, & if n = 1; \\ Φ_{n/p}(x^{p})/Φ_{n/p}(x), & if p|n, p prime; \end{matrix}  (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} + (\frac{1±\sqrt{5}}{2}) x + 1 whose coefficients are algebraic integers in Q[\sqrt{5}]. For odd square-free n > 1 the cyclotomic polynomial Φ_{n}(x) satisfies the identity

4Φ_{n}(x) = A_{n}^{2} − (−1)^{(n−1)/2}nB_{n}^{2}.  (6)

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

Φ_{n}((−1)^{(n−1)/2}x) = C_{n}^{2} − nxD_{n}^{2}  (7)

for odd square-free n, and

Φ_{n/2}(−x^{2}) = C_{n}^{2} − nxD_{n}^{2}  (8)

for even square-free n > 2. For a proof, see Lucas [22] or Schinzel [26].

In (6–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[\sqrt{±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 = m^{2}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.

1.1 Notation

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 symbol^{1} 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:

n′ = \{ \begin{matrix}n, & if n = 1 mod 4; \\ 2n, & otherwise. \end{matrix}

s = \{ \begin{matrix}−1, & if n = 3 mod 4; \\ +1, & otherwise. \end{matrix}

s′ = \{ \begin{matrix}−1, & if n = 5 mod 8; \\ +1, & otherwise. \end{matrix}

It is convenient to write g_{k} for (k, n) and g′_{k} for (k, n′).

Define

F_{n}(x) = \{ \begin{matrix}Φ_{n}(sx), & if n is odd; \\ (−1)^{ϕ(n/2)}Φ_{n/2}(−x^{2}), & if n is even. \end{matrix}  (9)

Thus we can write (7–8) as

F_{n}(x) = C_{n}^{2} − nxD_{n}^{2}.  (10)

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

F_{n}^{+}(x) = C_{n}(x) + \sqrt{nx}D_{n}(x)

and

F_{n}^{−}(x) = C_{n}(x) − \sqrt{nx}D_{n}(x).

1See, for example, Riesel [24]. To avoid ambiguity, we never write the Jacobi symbol as (\frac{m}{n}). Note that m|n without parentheses means that m divides n.

2See, for example, Landau [18].

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 \sqrt{m} means i\sqrt{|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.

1.2 Examples

Taking n = 15, we have

Φ_{15}(x) = \frac{(x^{15} − 1)(x − 1)}{(x^{5} − 1)(x^{3} − 1)} = x^{8} − x^{7} + x^{5} − x^{4} + x^{3} − x + 1,

A_{15}(x) = 2x^{4} − x^{3} − 4x^{2} − x + 2,
B_{15}(x) = x^{3} − x,
C_{15}(x) = x^{4} + 8x^{3} + 13x^{2} + 8x + 1,
D_{15}(x) = x^{3} + 3x^{2} + 3x + 1,

and the reader may easily verify that (6) and (7) are satisfied. As an example of (8), for n = 14 we have

F_{14}(x) = Φ_{7}(−x^{2}) = \frac{x^{14} + 1}{x^{2} + 1} = x^{12} − x^{10} + x^{8} − x^{6} + x^{4} − x^{2} + 1,

C_{14}(x) = x^{6} + 7x^{5} + 3x^{4} − 7x^{3} + 3x^{2} + 7x + 1,

and

D_{14}(x) = x^{5} + 2x^{4} − x^{3} − x^{2} + 2x + 1.

1.3 The identities of Beeger and Schinzel

Taking n = 5 in (7) we obtain

Φ_{5}(x) = (x^{2} + 3x + 1)^{2} − 5x(x + 1)^{2},

so Φ_{5}(x^{2}) has factors x^{4} + 3x^{2} + 1 ± \sqrt{5}(x^{3} + x) in Q[\sqrt{5}]. Replacing x by x^{3} we obtain factors of Φ_{5}(x^{6}). Now

Φ_{15}(x^{2}) = Φ_{5}(x^{6})/Φ_{5}(x^{2}),

and by division (taking the factors with opposite signs of \sqrt{5}) we obtain factors

x^{8} + 2x^{6} + 3x^{4} + 2x^{2} + 1 ± \sqrt{5}(x^{7} + x^{5} + x^{3} + x)

of Φ_{15}(x^{2}). Thus

Φ_{15}(x) = (x^{4} + 2x^{3} + 3x^{2} + 2x + 1)^{2} − 5x(x^{3} + x^{2} + x + 1)^{2}.  (11)

This is not of the form (10) because it gives a factorization of Φ_{15}(±x) over Q[\sqrt{±5}] instead of Q[\sqrt{∓15}]. Instead, (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).

1.4 Newton’s identities

Let

P (x) = ∏_{j=1}^{d}(x − ξ_{j}) = ∑_{j=0}^{d} a_{j}x^{d−j}

be a polynomial of degree d with arbitrary roots ξ_{j} and coefficients a_{0} = 1, a_{1}, …, a_{d}. For k > 0, define

p_{k} = ∑_{j=1}^{d} ξ_{j}^{k}.  (12)

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

ka_{k} = − ∑_{j=0}^{k−1} p_{k−j}a_{j}  (13)

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}, …):

x^{d}P (1/x) = ∑_{j=0}^{d} a_{j}x^{j} = exp (− ∑_{j=1}^{∞} p_{j}x^{j}/j)  (14)

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(x^{d}P (1/x)) = ∑_{k=1}^{d} ln(1 − ξ_{k}x) = − ∑_{k=1}^{d} ∑_{j=1}^{∞} ξ_{k}^{j}x^{j}/j = − ∑_{j=1}^{∞} (∑_{k=1}^{d} ξ_{k}^{j}) x^{j}/j = − ∑_{j=1}^{∞} p_{j}x^{j}/j.

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.

2 Theoretical Basis for the Algorithms

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}.

2.1 Cyclotomic polynomials

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).

3See 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.

Let ζ be a primitive n-th root of unity. To apply Newton’s identities we need to evaluate

p_{k} = ∑_{{0<j<n \atop (j,n)=1}} ζ^{jk}  (15)

for k = 1, 2, …, ϕ(n). This problem is well-known^{4}.

If n is prime the problem is easy: from

1 + ζ + ζ^{2} + ⋯ + ζ^{n−1} = \frac{1 − ζ^{n}}{1 − ζ} = 0,

we have p_{1} = −1. Moreover, for any k with (k, n) = 1, the map z ↦ z^{k} merely permutes {ζ, …, ζ^{n−1}}, so p_{k} = p_{1}.

Now consider the general case, n not necessarily prime. From (3) it is clear that p_{1} = μ(n). Let g_{k} = (k, n). If g_{k} = 1 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

p_{k} = \frac{μ(n/g_{k})ϕ(n)}{ϕ(n/g_{k})}.  (16)

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.

Lemma 1 If |x| ≥ R > 1 then

|Φ_{n}(x)| < R^{ϕ(n)} exp (\frac{1}{R − 1}) .

Proof

Let d = ϕ(n). From (14) with x replaced by 1/x, we have

Φ_{n}(x)/x^{d} = exp (− ∑_{j=1}^{∞} p_{j}x^{−j}/j) ,

but from (16) we have

|p_{j}| ≤ g_{j} ≤ min(j, n),

so

|Φ_{n}(x)/x^{d}| < exp (∑_{j=1}^{∞} R^{−j}) = exp (\frac{1}{R − 1}) .

This completes the proof. □

Lemma 2 If n > 1 is square-free and |x| ≥ R > 1 then

|F_{n}(x)| < R^{ϕ(2n)} exp (\frac{1}{R − 1}) .

4If ζ = e^{2πi/n} then our p_{k} is “Ramanujan’s sum” c_{n}(k), in the notation of Ramanujan [23] or Chapter 26 of Davenport [10].

Proof

From the definition (9) of F_{n}(x),

deg F_{n} = \{ \begin{matrix}ϕ(n), & if n is odd; \\ 2ϕ(n/2), & if n is even; \end{matrix}

so it is easy to see that

deg F_{n} = ϕ(2n)

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. □

2.2 The identity of Gauss

From now on we assume that n > 1 is square-free. In this subsection we also assume that n is odd. Consider the polynomial

G_{n}(x) = ∏_{{0<j<n \atop (j|n)=1}}(x − ζ^{j})  (17)

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],

2G_{n}(x) = A_{n}(x) − \sqrt{sn}B_{n}(x),  (18)

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)^{(n−1)/2}.

Define

2G̃_{n}(x) = A_{n}(x) + \sqrt{sn}B_{n}(x),  (19)

so Gauss’s identity (6) may be written as

Φ_{n}(x) = G_{n}(x)G̃_{n}(x).  (20)

The sums of k-th powers of roots of G_{n}(x) are

p_{k} = ∑_{{0<j<n \atop (j|n)=1}} ζ^{jk}.  (21)

Let g_{k} = (k, n). Then

2p_{k} = \{ \begin{matrix}μ(n) + (k|n)\sqrt{sn}, & if g_{k} = 1; \\ μ(n/g_{k})ϕ(g_{k}), & otherwise. \end{matrix}  (22)

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<j<n}(j|n)ζ^{j}

is \sqrt{sn}. On the other hand, if g_{k} > 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

2p_{k} = μ(n/g_{k})ϕ(g_{k}) + (k|n)\sqrt{sn}  (23)

is valid in both cases. Similarly for G̃_{n}(x), with the sign of \sqrt{sn} in (23) reversed.

Observe that p_{k} ∈ Q[\sqrt{n}] is real if n = 1 mod 4, but p_{k} ∈ Q[\sqrt{−n}] is complex if g_{k} = 1 and n = 3 mod 4.

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.

2.3 The identities of Aurifeuille, Le Lasseur and Lucas

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

L_{n}(x) = ∏_{j∈S_{n}}(x − ζ^{j}),  (24)

where

S_{n} = \{ \begin{matrix}\{j | 0 < j < 2n′, (j, n′) = 1, (j|n) = (−1)^{j}\}, & if n = 1 mod 4; \\ \{j | 0 < j < 2n′, (j, n′) = 1, (n|j) = 1\}, & otherwise. \end{matrix}  (25)

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[\sqrt{n}]. We later use the fact that L_{n}(x) is symmetric.

Schinzel [26] essentially shows (with a different notation) that

L_{n}(x) = C_{n}(x^{2}) − s′x\sqrt{n}D_{n}(x^{2})  (26)

where C_{n}(x) and D_{n}(x) are the polynomials of (10). Define

L̃_{n}(x) = L_{n}(−x) = C_{n}(x^{2}) + s′x\sqrt{n}D_{n}(x^{2}),  (27)

so after a change of variable (10) may be written as

F_{n}(x^{2}) = L_{n}(x)L̃_{n}(x).  (28)

Clearly F_{n}^{−}(x) = L_{n}(s′\sqrt{x}) and F_{n}^{+}(x) = L̃_{n}(s′\sqrt{x}).

Let g′_{k} = (k, n′). The sums p_{k} of k-th powers of roots of L_{n}(x) are

p_{k} = \{ \begin{matrix}(n|k)s′\sqrt{n}, & if k is odd; \\ μ(n′/g′_{k})ϕ(g′_{k}) cos ((n − 1)kπ/4) , & if k is even. \end{matrix}  (29)

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.

3 Algorithms

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}.

3.1 An algorithm for computing A_{n} and B_{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

A_{n}(x) = ∑_{j=0}^{d} α_{j}x^{d−j},

B_{n}(x) = ∑_{j=0}^{d} β_{j}x^{d−j},

where d = ϕ(n)/2, α_{0} = 2, β_{0} = 0, β_{1} = 1.

Recall the definition (21) of p_{k}. For k > 0 we have, from (23),

2p_{k} = q_{k} + r_{k}\sqrt{sn},

where q_{k} and r_{k} are integers given by

q_{k} = μ(n/g_{k})ϕ(g_{k})  (30)

and

r_{k} = (k|n).  (31)

Using (13) and (18), we obtain the recurrences

α_{k} = \frac{1}{2k} ∑_{j=0}^{k−1} (snr_{k−j}β_{j} − q_{k−j}α_{j})  (32)

and

β_{k} = \frac{1}{2k} ∑_{j=0}^{k−1} (r_{k−j}α_{j} − q_{k−j}β_{j})  (33)

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).

Comments on Algorithm D

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

α_{k} = (−1)^{d}α_{d−k}

if 2k > d and n > 3.

B_{n}/x is antisymmetric if n is composite and n = 3 mod 4, otherwise B_{n}/x is symmetric. Thus, we may use

β_{k} = \{ \begin{matrix}β_{d−k}, & in the symmetric case; \\ −β_{d−k}, & in the anti-symmetric case. \end{matrix}

Using these properties, the recurrences (32–33) need only be applied for k ≤ max(1, ⌊d/2⌋).

Example

Consider the case n = 15 as in Section 1.2. We have s = −1, d = ϕ(15)/2 = 4. Thus

q_{1} = q_{2} = q_{4} = μ(15)ϕ(1) = 1, q_{3} = μ(5)ϕ(3) = −2,
r_{1} = (1|15) = 1, r_{2} = (2|15) = (2|3)(2|5) = 1, r_{3} = (3|15) = 0, r_{4} = (4|15) = 1.

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

α_{1} = (−15r_{1}β_{0} − q_{1}α_{0})/2 = −1,
β_{1} = (r_{1}α_{0} − q_{1}β_{0})/2 = 1,
α_{2} = (−15r_{2}β_{0} − 15r_{1}β_{1} − q_{2}α_{0} − q_{1}α_{1})/4 = −4,
β_{2} = (r_{2}α_{0} + r_{1}α_{1} − q_{2}β_{0} − q_{1}β_{1})/4 = 0.

Using symmetry of A_{n}(x) and anti-symmetry of B_{n}(x)/x, or continuing with the recurrences (32–33), we obtain α_{3} = α_{1} = −1, β_{3} = −β_{1} = −1, α_{4} = α_{0} = 2, β_{4} = −β_{0} = 0. Thus A_{15}(x) = 2x^{4} − x^{3} − 4x^{2} − x + 2 and B_{15}(x) = x^{3} − x, as expected.

3.2 An algorithm for computing C_{n} and D_{n}

Consider the computation of C_{n} and D_{n} for square-free n > 1. Define n′, s, s′ and L_{n} as in Section 2.3, and d = ϕ(n′)/2. Thus deg L_{n} = 2d, deg C_{n} = d, and deg D_{n} = d − 1. From (26) it is enough to compute the coefficients a_{k} of L_{n}(x). In order to work over the integers, we define

q_{k} = \{ \begin{matrix}s′p_{k}/\sqrt{n}, & if k is odd; \\ p_{k}, & if k is even; \end{matrix}

where p_{k} is the sum of k-th powers of roots of L_{n}(x). Thus, from (29),

q_{k} = \{ \begin{matrix}(n|k), & if k is odd; \\ μ(n′/g′_{k})ϕ(g′_{k}) cos ((n − 1)kπ/4) , & otherwise. \end{matrix}  (34)

If

C_{n}(x) = ∑_{j=0}^{d} γ_{j}x^{d−j}

and

D_{n}(x) = ∑_{j=0}^{d−1} δ_{j}x^{d−1−j},

then, from (26),

γ_{k} = a_{2k}

and

δ_{k} = −s′a_{2k+1}/\sqrt{n}.

In particular, γ_{0} = δ_{0} = 1. Using (13) we obtain

γ_{k} = \frac{1}{2k} ∑_{j=0}^{k−1} (nq_{2k−2j−1}δ_{j} − q_{2k−2j}γ_{j})  (35)

and

δ_{k} = \frac{1}{2k + 1} (γ_{k} + ∑_{j=0}^{k−1} (q_{2k+1−2j}γ_{j} − q_{2k−2j}δ_{j}))  (36)

for k = 1, 2, …

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.

3. Evaluate γ_{k} for k = 1, …, ⌊d/2⌋ and δ_{k} for k = 1, …, ⌊(d − 1)/2⌋ using (35–36).

4. Evaluate γ_{k} for k = ⌊d/2⌋ + 1, …, d using γ_{k} = γ_{d−k}.

5. Evaluate δ_{k} for k = ⌊(d + 1)/2⌋, …, d − 1 using δ_{k} = δ_{d−1−k}.

Example

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,
q_{2} = μ(15)ϕ(2) cos(7π) = −1,
q_{3} = (15|3) = 0,
q_{4} = μ(15)ϕ(2) cos(14π) = 1.

The initial conditions are γ_{0} = δ_{0} = 1. The recurrences (35–36) give

γ_{1} = (15q_{1}δ_{0} − q_{2}γ_{0})/2 = 8,
δ_{1} = (γ_{1} + q_{3}γ_{0} − q_{2}δ_{0})/3 = 3,
γ_{2} = (15q_{3}δ_{0} + 15q_{1}δ_{1} − q_{4}γ_{0} − q_{2}γ_{1})/4 = 13.

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.

3.3 Stevenhagen’s algorithm

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 order^{5} 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’s algorithm failed due to integer overflow for n = 35.

Algorithm L does not suffer from this problem. It is easy to see from the recurrences (35– 36) that intermediate results can grow only slightly larger than the final coefficients γ_{k} and δ_{k}. A straightforward implementation of Algorithm L can compute C_{n} and D_{n} for all square-free n < 180 without encountering integer overflow in 32-bit arithmetic. When it does eventually occur, overflow is easily detected because the division by 2k in (35) or by 2k + 1 in (36) gives a non-integer result.

4 Explicit expressions for A_{n}, …, D_{n}

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.

4.1 The analytic functions f_{n} and g_{n}

For odd square-free n > 1 and |x| ≤ 1, define

f_{n}(x) = ∑_{j=1}^{∞} (j|n)\frac{x^{j}}{j} .  (37)

Similarly, for square-free n > 1 and |x| ≤ 1, define

g_{n}(x) = ∑_{j=0}^{∞} (n | 2j + 1)\frac{x^{2j+1}}{2j + 1}.  (38)

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(\sqrt{sn}f_{n}(x)) and exp(2\sqrt{n}g_{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

\{ \begin{matrix}f_{n}(x) = f_{n}(1/x), & if n = 1 mod 4; \\ f_{n}(x) + f_{n}(1/x) = 2f_{n}(1), & if n = 3 mod 4; \end{matrix}  (39)

5The 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.

and

g_{n}(x) = g_{n}(1/x).  (40)

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

f_{n}(x) − f_{n}(−x) = 2g_{n}(x\sqrt{s})/\sqrt{s}.

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[\sqrt{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 D = −n) are

f_{n}(1) = \frac{−π}{n^{3/2}} ∑_{j=1}^{n−1}(j|n)j = \frac{π}{(2 − (2|n))\sqrt{n}} ∑_{j=1}^{(n−1)/2}(j|n) = \frac{2π}{w\sqrt{n}}h(−n) > 0.  (41)

Here

w = \{ \begin{matrix}6, & if n = 3; \\ 2, & if n = 3 mod 4, n > 3 \end{matrix}

is the number of roots of unity in Q[\sqrt{D}]. Since h(−n) is an integer and h(−3) = 1, we have

exp(2i\sqrt{n}f_{n}(1)) = \{ \begin{matrix}(−1 + \sqrt{−3})/2, & if n = 3; \\ 1, & if n = 3 mod 4, n > 3. \end{matrix}  (42)

In the case n = 1 mod 4, we have D = n, and

f_{n}(1) = \frac{ln ε}{\sqrt{n}}h(n).  (43)

Here ε is the “fundamental unit”, i.e. ε = (|u| + \sqrt{n}|v|)/2, where (u, v) is a minimal nontrivial solution of u^{2} − nv^{2} = 4. For example, if n = 5 then ε = (3 + \sqrt{5})/2, h(5) = 1 and f_{5}(1) = (ln ε)/\sqrt{5} = 0.4304 …

4.2 The polynomials A_{n} and B_{n}

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}, G̃_{n} be as in Section 2.2, and d = ϕ(n)/2. Recall that

G_{n}(x) = ∏_{{0<j<n \atop (j|n)=1}}(x − ζ^{j})  (44)

and

G̃_{n}(x) = ∏_{{0<j<n \atop (j|n)=−1}}(x − ζ^{j}).  (45)

From (14) and (23), we have

G̃_{n}(1/x)/G_{n}(1/x) = exp(\sqrt{sn}f_{n}(x)).  (46)

Also, from (44),

(−x)^{d}G_{n}(1/x) = ∏_{{0<j<n \atop (j|n)=1}}(ζ^{j}x − 1) = ∏_{{0<j<n \atop (j|n)=1}} ζ^{j}(x − ζ^{−j}),  (47)

so

G_{n}(1/x)/G̃_{n}(1/x) = ζ^{σ} ∏_{j=1}^{n−1}(x − ζ^{−j})^{(j|n)},  (48)

where

σ = ∑_{j=1}^{n−1}(j|n)j.  (49)

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

G_{n}(1/x)/G̃_{n}(1/x) = \{ \begin{matrix}G_{n}(x)/G̃_{n}(x), & if n = 1 mod 4; \\ G̃_{n}(x)/G_{n}(x), & if n = 3 mod 4. \end{matrix}  (50)

It follows from (46) that

G̃_{n}(x)/G_{n}(x) = exp(s\sqrt{sn}f_{n}(x)).  (51)

We see from (46) or (51) that, as claimed above, exp(\sqrt{sn}f_{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

G_{n}(x) = \sqrt{Φ_{n}(x)} exp (\frac{−s\sqrt{sn}}{2}f_{n}(x))  (52)

If (52) is interpreted as a generating function for G_{n}(x) then \sqrt{Φ_{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) + G̃_{n}(x) and B_{n}(x) = (G̃_{n}(x) − G_{n}(x))/\sqrt{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

A_{n}(x) = 2\sqrt{Φ_{n}(x)} cosh (\frac{\sqrt{sn}}{2}f_{n}(x))  (53)

and

B_{n}(x) = 2\sqrt{\frac{Φ_{n}(x)}{sn}} sinh (\frac{\sqrt{sn}}{2}f_{n}(x)) .  (54)

Remark

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

A_{n}(x) = 2\sqrt{Φ_{n}(x)} cos (\frac{\sqrt{n}}{2}f_{n}(x))  (55)

and

B_{n}(x) = 2\sqrt{\frac{Φ_{n}(x)}{n}} sin (\frac{\sqrt{n}}{2}f_{n}(x)) .  (56)

Example

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

Φ_{15}(x) = 1 − x + x^{3} − x^{4} + x^{5} − x^{7} + x^{8},

so

\sqrt{Φ_{15}(x)} = 1 − \frac{x}{2} − \frac{x^{2}}{8} + \frac{7x^{3}}{16} − \frac{37x^{4}}{128} + ⋯

Also,

f_{15}(x) = x + \frac{x^{2}}{2} + \frac{x^{4}}{4} − \frac{x^{7}}{7} + ⋯

so

cos (\frac{\sqrt{15}}{2}f_{15}(x)) = 1 − \frac{15x^{2}}{8} − \frac{15x^{3}}{8} + \frac{15x^{4}}{128} + ⋯

and

2\sqrt{Φ_{15}(x)} cos (\frac{\sqrt{15}}{2}f_{15}(x)) = 2 − x − 4x^{2} − x^{3} + 2x^{4} + ⋯,

which is A_{15}(x) + O(x^{5}) 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

Φ_{n}(1) = \{ \begin{matrix}n, & if n is prime; \\ 1, & for other square-free n > 1; \end{matrix}

we can verify that (46) and (51–56) give sensible results in the case x = 1. For example, we have h(−15) = 2, so f_{15}(1) = 2π/\sqrt{15} and cos(\frac{\sqrt{15}}{2}f_{15}(1)) = cos(π) = −1. Thus (55) gives the correct value −2 of A_{15}(1).

4.3 The polynomials C_{n} and D_{n}

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

L̃_{n}(1/x)/L_{n}(1/x) = exp (2s′\sqrt{n}g_{n}(x)) .  (57)

This shows that exp(2\sqrt{n}g_{n}(x)) is a rational function. From

L̃_{n}(1/x)/L_{n}(1/x) = L̃_{n}(x)/L_{n}(x)

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

L̃_{n}(x)/L_{n}(x) = exp (2s′\sqrt{n}g_{n}(x)) .  (58)

From (28) and (58), taking a square root, we obtain

L_{n}(x) = \sqrt{F_{n}(x^{2})} exp (−s′\sqrt{n}g_{n}(x)) .  (59)

Theorem 2 Let n > 1 be square-free. The Aurifeuillian factors F_{n}^{±}(x) = C_{n}(x) ± \sqrt{nx}D_{n}(x) of F_{n}(x) are given by

F_{n}^{±}(x) = \sqrt{F_{n}(x)} exp (±\sqrt{n}g_{n}(\sqrt{x})) .  (60)

Also,

C_{n}(x) = \sqrt{F_{n}(x)} cosh (\sqrt{n}g_{n}(\sqrt{x}))  (61)

and

D_{n}(x) = \sqrt{\frac{F_{n}(x)}{nx}} sinh (\sqrt{n}g_{n}(\sqrt{x})) .  (62)

Proof

Recall that the Aurifeuillian factors F_{n}^{±}(x) = C_{n}(x) ± \sqrt{nx}D_{n}(x) of F_{n}(x) are L_{n}(±\sqrt{x}). Thus (60) follows from (59). Since

L_{n}(x) + L_{n}(−x) = 2C_{n}(x^{2})

and

L_{n}(x) − L_{n}(−x) = −2s′x\sqrt{n}D_{n}(x^{2}),

we easily deduce (61–62). □

4.4 Application to integer factorization

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 \sqrt{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

F_{n}^{−}(x) = ⌊F̂ + 1/2⌋

and

F_{n}^{+}(x) = F_{n}(x)/F_{n}^{−}(x),

where

F̂ = \sqrt{F_{n}(x)} exp (− \frac{1}{m} ∑_{j=0}^{λ−1} \frac{(n | 2j + 1)}{(2j + 1)x^{j}} ) .

Proof

From (60), using the functional equation (40) and the power series (38) for g_{n}(1/\sqrt{x}), we have

F_{n}^{−}(x) = \sqrt{F_{n}(x)} exp (− \frac{1}{m} ∑_{j=0}^{∞} \frac{(n | 2j + 1)}{(2j + 1)x^{j}} ) .

Thus, we only have to show that the error incurred by truncating the power series after λ terms is less than 1/2 in absolute value, i.e. that |F̂ − F_{n}^{−}| < 1/2.

If x < 5 we must have m = 1, n = 2 or 3, x = n, λ = 1. In both cases F_{n}^{−}(n) = 1 and it is easy to verify that 1/2 < F̂ < 3/2. Thus, from now on we assume that x ≥ 5.

Let

t = \frac{1}{m} ∑_{j=λ}^{∞} \frac{(n | 2j + 1)}{(2j + 1)x^{j}}.

Since m ≥ 1, λ ≥ 1, x ≥ 5, and the Jacobi symbol is 0 or ±1, we have

|t| ≤ ∑_{j=λ}^{∞} \frac{x^{−j}}{(2j + 1)} ≤ (\frac{1}{3} + \frac{1}{5 · 5} + \frac{1}{7 · 5^{2}} + ⋯) x^{−λ}

so

|t| ≤ 5 (\sqrt{5} ln (\frac{1 + \sqrt{5}}{2}) − 1) x^{−λ} < 0.3802x^{−λ}.  (63)

In particular, |t| < 0.3802/x < 0.08, so

| exp(t) − 1| < \frac{|t|}{1 − |t|/2} < \frac{0.3802x^{−λ}}{0.96} < 0.4x^{−λ}.  (64)

Now F̂ = F_{n}^{−} exp(t), so

|F̂ − F_{n}^{−}| = F_{n}^{−}| exp(t) − 1|.  (65)

Applying Lemma 2 of Section 2.1 with R = 5 gives

F_{n}(x) < exp(1/4)x^{ϕ(2n)},

but ϕ(2n) = 2λ and F_{n}^{−} is the smaller factor of F_{n}, so

F_{n}^{−}(x) ≤ \sqrt{F_{n}(x)} < exp(1/8)x^{λ}.  (66)

From (64–66) we finally obtain the bound

|F̂ − F_{n}^{−}| < 0.4 exp(1/8) < 0.5,

which completes the proof. □

Since the argument of the exponential in Theorem 3 is −1/m+O(1/n) as n → ∞, we have the following result, which sheds some light on the ratio of the Aurifeuillian factors. (The Corollary strictly follows from Theorem 3 only if m is an integer, but from the proof of Theorem 3 it is clear that the Corollary is also valid for rational m.)

Corollary 1 Let x = m^{2}n where m > 0 is rational n > 1 is an integer. Consider m fixed as n → ∞ through square-free values. Then the Aurifeuillian factors F_{n}^{±}(x) of F_{n}(x) satisfy

F_{n}^{+}(x)/F_{n}^{−}(x) = exp(2/m) + O(1/n).

Examples

1. To start with a simple example, consider n = 2, m = 2, so x = m^{2}n = 8 and λ = 1. In Theorem 3 we have F_{2}(x) = x^{2} + 1 = 65,

F̂ = \sqrt{F_{2}(x)} exp(−1/m) = \sqrt{65} exp(−1/2) = 4.89 ⋯

and rounding to the nearest integer gives the factor 5 of 65.

2. A similar but less trivial example is n = 2, m = 2^{5}, x = m^{2}n = 2^{11}. Here F_{2}(x) = x^{2} + 1 = 2^{22} + 1 and

F̂ = \sqrt{F_{2}(x)} exp(−1/m) = \sqrt{2^{22} + 1} exp(−2^{−5}) = 1984.98 ⋯

which on rounding gives the factor F_{2}^{−} = 1985 of 2^{22} + 1. By division we find F_{2}^{+} = 2113, so the complete factorization is 2^{22} + 1 = 5 · 397 · 2113 .

3. Now consider n = 5, m = 3, so x = m^{2}n = 45 and λ = ϕ(10)/2 = 2. In this case

F_{5}(x) = Φ_{5}(x) = (x^{5} − 1)/(x − 1) = 4193821,

F̂ = \sqrt{F_{5}(x)} exp (−\frac{1}{m} + \frac{1}{3m^{3}n}) = \sqrt{4193821} exp(−134/405) = 1470.99924 ⋯

and rounding to the nearest integer gives the factor 1471 of F_{5}(x). By division we obtain the other factor 2851. In this example, but not in general, the Aurifeuillian factors are prime.

4. Now consider a composite n, say n = 15. To keep the arithmetic easy we take m = 1, so x = 15 and

F_{15}(x) = Φ_{15}(−x) = x^{8} + x^{7} − x^{5} − x^{4} − x^{3} + x + 1 = 2732936641.

The example in Section 3.2 shows how we can use Algorithm L to compute

C_{15}(x) = x^{4} + 8x^{3} + 13x^{2} + 8x + 1

and

D_{15}(x) = x^{3} + 3x^{2} + 3x + 1.

Evaluating the polynomials, we find

C_{15}(15) = 80671, D_{15}(15) = 4096,

so the Aurifeuillian factors of F_{15}(15) are 80671 ± 15 · 4096. This gives 19231 and 142111, which can easily be verified to be factors of 15^{15} + 1. The “algebraic” factors (15^{3} + 1)/16 and (15^{5} + 1)/16 allow us to complete the factorization:

15^{15} + 1 = 2^{4} · 31 · 211 · 1531 · 19231 · 142111

Alternatively, instead of evaluating C_{15}(15) and D_{15}(15), we can use Theorem 3. We have λ = ϕ(15)/2 = 4. Since (15|3) = (15|5) = 0 and (15|7) = 1, the computation gives

F̂ = \sqrt{F_{15}(x)} exp (−1 − \frac{1}{7n^{3}} ) = 19231.00217 ⋯

and rounding to the nearest integer gives the factor 19231 of F_{15}(15).

5. 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 \sqrt{nx} = pn/q = 14/5. We have

F_{7}(x) = Φ_{7}(−x) = \frac{x^{7} + 1}{x + 1} = \frac{28^{7} + 25^{7}}{53 · 25^{6}}

Theorem 3 is not applicable. Because x is close to 1, the series for g_{7}(1/\sqrt{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

C_{7}(x) = x^{3} + 3x^{2} + 3x + 1, D_{7}(x) = x^{2} + x + 1,

so Horner’s rule gives

C_{7}(x) = 148877/5^{6}, D_{7}(x) = 2109/5^{4},

and C_{7}(x) ± \sqrt{nx}D_{7}(x) gives 296507/5^{6} and 1247/5^{6}. Thus, we have obtained factors 296507 and 1247 of (25^{7} + 28^{7})/53. The larger factor is prime, so it is easy to deduce the complete factorization

25^{7} + 28^{7} = 29 · 43 · 53 · 296507

Acknowledgements

Thanks are due to Brendan McKay for suggesting the use of the generating function (14), to Emma Lehmer and an anonymous referee for helpful comments on the exposition, to Hans Riesel for his kind assistance with the solution of exercise A6.2 of [24], and to Sam Wagstaff and Hugh Williams for providing copies of several references which were difficult to find in Australia.

References

[1] A. V. Aho, J. E. Hopcroft and J. D. Ullman, The Design and Analysis of Computer Algorithms, Addison-Wesley, Menlo Park, Calif., 1974, Chapter 8.

[2] A. Aurifeuille and H. Le Lasseur, see [20], page 276, or [21], page 785.

[3] N. G. W. H. Beeger, “On a new quadratic form for certain cyclotomic polynomials”, Nieuw Arch. Wisk. (2), 23 (1951), 249-252.

[4] R. P. Brent, F. G. Gustavson and D. Y. Y. Yun, “Fast solution of Toeplitz systems of equations and computation of Padé approximants”, J. Algorithms 1 (1980), 259-295.

[5] R. P. Brent and H. T. Kung, “Fast algorithms for manipulating formal power series”, J. ACM 25 (1978), 581-595.

[6] R. P. Brent and H. J. J. te Riele, Factorizations of a^{n} ± 1, 13 ≤ a < 100, Report NM-R9212, Department of Numerical Mathematics, Centrum voor Wiskunde en Informatica, Amsterdam, June 1992.

[7] J. Brillhart, D. H. Lehmer, J. L. Selfridge, B. Tuckerman and S. S. Wagstaff, Jr., Factorizations of b^{n} ± 1, b = 2, 3, 5, 6, 7, 10, 11, 12 up to high powers, second edition, American Mathematical Society, Providence, Rhode Island, 1988.

[8] A. J. C. Cunningham, “Factorisation of N = y^{y} ∓ 1 and x^{xy} ∓ y^{xy}”, Messenger of Math. (2), 45 (1915), 49-75.

[9] A. J. C. Cunningham and H. J. Woodall, Factorisation of y^{n} ∓ 1, y = 2, 3, 5, 6, 7, 10, 11, 12 up to high powers (n), Hodgson, London, 1925.

[10] H. Davenport, Multiplicative Number Theory, second edition (revised by H. L. Montgomery), Springer-Verlag, New York, 1980.

[11] P. G. Lejeune Dirichlet, Vorlesungen über Zahlentheorie, fourth edition, Friedr. Vieweg & Sohn, Braunschweig, 1894, Chapter 5 and Supplement 7.

[12] C. F. Gauss, Disquisitiones Arithmeticæ, G. Fleischer, Leipzig, 1801, Art. 356-357. Reprinted in Carl Friedrich Gauss Werke, Band 1, Georg Olms Verlag, Hildesheim, 1981.

[13] G. H. Hardy and E. M. Wright, An Introduction to the Theory of Numbers, fifth edition, Clarendon Press, Oxford, 1984, Ch. 16.

[14] D. E. Knuth, The Art of Computer Programming, Volume 2: Seminumerical Algorithms (second edition), Addison-Wesley, Menlo Park, 1981, Chapter 3.

[15] M. Kraitchik, “Décomposition de a^{n} ± b^{n} en facteurs dans le cas où nab est un carré parfait avec une table des décompositions numériques pour toutes les valeurs de a et b inférieures à 100”, Gauthiers-Villars, Paris, 1922.

[16] M. Kraitchik, Recherches sur la Théorie des Nombres, Volume 1, Gauthiers-Villars, Paris, 1924.

[17] M. Kraitchik, Introduction à la Théorie des Nombres, Gauthiers-Villars, Paris, 1952.

[18] Edmund Landau, Vorlesungen über Zahlentheorie, Band 1(1): Aus der elementaren und additiven Zahlentheorie, Leipzig, 1927. English translation: Elementary Number Theory, Chelsea, New York, 1958.

[19] Serge Lang, Cyclotomic Fields I and II, combined second edition, Graduate Texts in Mathematics 126, Springer-Verlag, New York, 1990.

[20] E. Lucas, “Théorèmes d’arithmétique”, Atti. R. Acad. Sc. Torino 13 (1877-8), 271-284.

[21] E. Lucas, “Sur la série récurrente de Fermat”, Bull. Bibl. Storia Sc. Mat. e Fis. 11 (1878), 783-789.

[22] E. Lucas, “Sur les formules de Cauchy et de Lejeune-Dirichlet”, Ass. Française pour l’Avanc. des Sci., Comptes Rendus 7 (1878), 164-173.

[23] S. Ramanujan, “On certain trigonometrical sums and their applications in the theory of numbers”, Trans. Cambridge Philos. Soc. 22, 13 (1918), 259-276. Reprinted in Collected Papers of Srinivasa Ramanujan (edited by G. H. Hardy, P. V. Seshu Aiyar and B. M. Wilson), Cambridge Univ. Press, 1927.

[24] Hans Riesel, Prime Numbers and Computer Methods for Factorization, Birkhäuser, Boston, 1985.

[25] John Riordan, An Introduction to Combinatorial Analysis, Princeton University Press, New Jersey, 1978, Ch. 2, exercise 27.

[26] A. Schinzel, “On the primitive prime factors of a^{n} − b^{n}”, Proc. Cambridge Philos. Soc. 58 (1962), 555-562. MR 26#1280

[27] Peter Stevenhagen, “On Aurifeuillian factorizations”, Nederl. Akad. Wetensch. Indag. Math. 49 (1987), 451-468. MR 89a:11015 [28] H. W. Turnbull, Theory of Equations (fifth edition), Oliver and Boyd, Edinburgh, 1952,

Sec. 32.

[29] B. L. van der Waerden, Algebra, Vol. 1 (English translation by Fred Blum, fifth edition), Frederick Ungar, New York, 1953, Ch. 7.

[30] L. C. Washington, Introduction to Cyclotomic Fields, Graduate Texts in Mathematics 83, Springer-Verlag, New York, 1982.