COMPUTING AURIFEUILLIAN FACTORS

RICHARD P. BRENT

Abstract. For odd square-free n > 1, the cyclotomic polynomial Φ_{n}(x) satisfies an identity Φ_{n}(x) = C_{n}(x)^{2} ± nxD_{n}(x)^{2} of Aurifeuille, Le Lasseur and Lucas. Here C_{n}(x) and D_{n}(x) are monic polynomials with integer coefficients. These coefficients can be computed by simple algorithms which require O(n^{2}) arithmetic operations over the integers. Also, there are explicit formulas and generating functions for C_{n}(x) and D_{n}(x). This paper is a preliminary report which states the results for the case n = 1 mod 4, and gives some numerical examples. The proofs, generalisations to other square-free n, and similar results for the identities of Gauss and Dirichlet, will appear elsewhere.

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

and the Möbius inversion formula [9] gives

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

Equation (1) is useful for theoretical purposes, but (2) 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

(3) Φ_{n}(x) = \begin{cases}x − 1 & if n = 1, \\ Φ_{n/p}(x^{p})/Φ_{n/p}(x) & if p is prime and p|n, \end{cases}

give another convenient recursion for computing Φ_{n}(x).

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

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

To appear in Proceedings of a Conference on Computational Algebra and Number Theory (held at the University of Sydney, November 1992). This paper is a preliminary report and introduction to [1].

Copyright \begin{matrix}c \\ ◯ \end{matrix} 1992, R. P. Brent.

rpb127 typeset using AMS-L^{A}TEX

In this preliminary report we omit proofs, and assume from now on that

(4) n > 1 is square-free and n = 1 mod 4.

The results can be generalized to other square-free n, and similar results hold for the identities of Gauss and Dirichlet. The interested reader is referred to [1] for details.

Φ_{n}(x) satisfies an identity

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

of Aurifeuille, Le Lasseur and Lucas^{1}. For a proof, see Lucas [15] or Schinzel [17]. Here C_{n}(x) and D_{n}(x) are symmetric, monic polynomials with integer coefficients. For example, if n = 5, we have

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

so

C_{5}(x) = x^{2} + 3x + 1 and D_{5}(x) = x + 1. (6)

In Section 1.1 we summarize our notation. Then, in Section 2, we outline the theoretical basis for our algorithm for computing C_{n}(x) and D_{n}(x). The algorithm (Algorithm L) is presented in Section 3. Algorithm L appears to be new, although the key idea (using Newton’s identities to evaluate polynomial coefficients) is due to Dirichlet [8]. A different algorithm, due to Stevenhagen [18], is discussed in Section 3.1.

In Section 4 we give explicit formulas for C_{n}(x), D_{n}(x) etc. These may be regarded as generating functions if x is an indeterminate, or may be used to compute C_{n}(x) and D_{n}(x) for given argument x. In the special case x = 1 there is an interesting connection with Dirichlet L-functions and the theory of class numbers of quadratic fields.

One application of cyclotomic polynomials is to the factorization of integers of the form a^{n}±b^{n}: see for example [3, 4, 5, 6, 11, 12, 16]. If x = m^{2}n for any integer m, then (5) is a difference of squares, giving integer factors C_{n}(x) ± mnD_{n}(x) of x^{n} ± 1. Examples are given in Sections 3–4.

1.1. Notation. For consistency we follow the notation of [1] where possible, although there are simplifications due to our assumption (4).

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 [9]. Note that μ(1) = ϕ(1) = 1.

(m|n) denotes the Jacobi symbol^{2} except that (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. As usual, m|n without parentheses means that m divides n.

n denotes a positive integer satisfying (4), which implies that (−1|n) = 1. It is convenient to write g_{k} for (k, n).

For given n, we define s′ = (2|n). In view of (4), the following are equivalent:

s′ = (2|n) = (−1)^{(n^{2}−1)/8} = (−1)^{(n−1)/4} = \begin{cases}+1 & if n = 1 mod 8, \\ −1 & if n = 5 mod 8. \end{cases}

1Lucas [13, page 276] states “Les formules et les conséquences précédentes sont dues à la collaboration de M. Aurifeuille, ancien Professeur au lycée de Toulouse, actuellement décédé, et de M. Le Lasseur, de Nantes”. See also [14, page 785].

2See, for example, Riesel [16]. To avoid ambiguity, we never write the Jacobi symbol as (\frac{m}{n}).

The Aurifeuillian factors of Φ_{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).

From (5) we have Φ_{n}(x) = F_{n}^{−}(x)F_{n}^{+}(x).

1.2. Acknowledgements. Thanks are due to Emma Lehmer, Brendan McKay, Hans Riesel, Sam Wagstaff, and Hugh Williams for their comments and assistance.

2. Theoretical results

In this section we summarise some theoretical results which form the basis for Algorithm L. 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 (9). Consider the polynomial

L_{n}(x) = ∑_{{0<j<n \atop (j|n)=(−1)^{j}}} (x − ζ^{j})(x − ζ^{−j}),

which we may write as

L_{n}(x) = ∑_{{0<j<n \atop (j|n)=(−1)^{j}}} (x^{2} − 2 (cos \frac{πj}{n} ) x + 1) . (7)

L_{n}(x) has degree ϕ(n). Also, from (7), L_{n}(x) is symmetric and has real coefficients. Schinzel [17] shows that

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

where C_{n}(x) and D_{n}(x) are the polynomials of (5), and s′ = (2|n) as usual. Clearly F_{n}^{−}(x) = L_{n}(s′\sqrt{x}) and F_{n}^{+}(x) = L_{n}(−s′\sqrt{x}).

For example, suppose n = 5. Then (7) gives

L_{5}(x) = (x^{2} − 2 (cos \frac{3π}{5} ) x + 1) (x^{2} − 2 (cos \frac{4π}{5} ) x + 1) ,

but cos 3π/5 = (1 − \sqrt{5})/4 and cos 4π/5 = −(1 + \sqrt{5})/4, so it is easily verified that

L_{5}(x) = x^{4} + \sqrt{5}x^{3} + 3x^{2} + \sqrt{5}x + 1 = C_{5}(x^{2}) + x\sqrt{5}D_{5}(x^{2}),

where C_{5}(x) and D_{5}(x) are as in (6).

Let g_{k} = (k, n). It may be shown that the Gaussian sums p_{k} of k-th powers of roots of L_{n}(x) are

(9) p_{k} = \begin{cases}(n|k)s′\sqrt{n} & if k is odd, \\ μ(n/g_{k})ϕ(g_{k}) & if k is even. \end{cases}

3. An algorithm for computing C_{n} and D_{n}

In this section we consider the computation of C_{n} and D_{n}. Define d = ϕ(n)/2. Thus deg L_{n} = 2d, deg C_{n} = d, and deg D_{n} = d − 1. From (8) it is enough to compute the coefficients a_{k} of L_{n}(x). Using (9), the coefficients of L_{n}(x), and hence of C_{n}(x) and D_{n}(x), may be evaluated from Newton’s identities. In order to work over the integers, we define

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

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

(10) q_{k} = \begin{cases}(n|k) & if k is odd, \\ μ(n/g_{k})ϕ(g_{k}) & if k is even. \end{cases}

If

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

then, from (8),

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

In particular, γ_{0} = δ_{0} = 1. Using Newton’s identities, we obtain the recurrences

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

and

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

for k > 0.

We can use the fact that C_{n}(x) and D_{n}(x) are symmetric to reduce the number of times the recurrences (11)–(12) need to be applied. An algorithm which incorporates this refinement is:

Algorithm L 1. Evaluate q_{k} for k = 1, …, d using the definition (10).

2. Set γ_{0} ← 1 and δ_{0} ← 1.

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

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

Examples.

1. Consider the case n = 5. We have s′ = (2|5) = −1, d = ϕ(5)/2 = 2. Thus

q_{1} = (5|1) = 1 and q_{2} = μ(5)ϕ(1) = −1.

The initial conditions are γ_{0} = δ_{0} = 1. The recurrence (11) gives

γ_{1} = (5q_{1}δ_{0} − q_{2}γ_{0})/2 = 3.

Using symmetry we obtain γ_{2} = γ_{0} = 1 and δ_{1} = δ_{0} = 1. Thus

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

and it is easy to verify that Φ_{5}(x)^{2} = C_{5}(x)^{2} − 5xD_{5}(x)^{2}, as expected from (5).

2. Now consider n = 33. We have s′ = (2|33) = 1, d = ϕ(33)/2 = 10. Thus

q_{1} = (33|1) = 1,  q_{2} = μ(33)ϕ(1) = 1,
q_{3} = (33|3) = 0,  q_{4} = μ(33)ϕ(1) = 1,
q_{5} = (33|5)= −1,  q_{6} = μ(11)ϕ(3) = −2,
q_{7} = (33|7)= −1,  q_{8} = μ(33)ϕ(1) = 1,
q_{9} = (33|9) = 0,  q_{10}= μ(33)ϕ(1) = 1.

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

γ_{1} = (33q_{1}δ_{0} − q_{2}γ_{0})/2 = 16,
δ_{1} = (γ_{1} + q_{3}γ_{0} − q_{2}δ_{0})/3 = 5,
γ_{2} = (33q_{3}δ_{0} − q_{4}γ_{0} + 33q_{1}δ_{1} − q_{2}γ_{1})/4 = 37,
δ_{2} = (γ_{2} + q_{5}γ_{0} − q_{4}δ_{0} + q_{3}γ_{1} − q_{2}δ_{1})/5 = 6,
γ_{3} = (33q_{5}δ_{0} − q_{6}γ_{0} + 33q_{3}δ_{1} − q_{4}γ_{1} + 33q_{1}δ_{2} − q_{2}γ_{2})/6 = 19,
δ_{3} = (γ_{3} + q_{7}γ_{0} − q_{6}δ_{0} + q_{5}γ_{1} − q_{4}δ_{1} + q_{3}γ_{2} − q_{2}δ_{2})/7 = −1,
γ_{4} = (33q_{7}δ_{0} − q_{8}γ_{0} + ⋯ + 33q_{1}δ_{3} − q_{2}γ_{3})/8 = −32,
δ_{4} = (γ_{4} + q_{9}γ_{0} − q_{8}δ_{0} + ⋯ + q_{3}γ_{3} − q_{2}δ_{3})/9 = −9,
γ_{5} = (33q_{9}δ_{0} − q_{10}γ_{0} + ⋯ + 33q_{1}δ_{4} − q_{2}γ_{4})/10 = −59.

Using symmetry, we obtain

C_{33}(x) = x^{10} + 16x^{9} + 37x^{8} + 19x^{7} − 32x^{6} − 59x^{5} − 32x^{4} + 19x^{3} + 37x^{2} + 16x + 1

and

D_{33}(x) = x^{9} + 5x^{8} + 6x^{7} − x^{6} − 9x^{5} − 9x^{4} − x^{3} + 6x^{2} + 5x + 1.

From the recurrence (3),

Φ_{33}(x) = Φ_{3}(x^{11})/Φ_{3}(x) = \frac{x^{22} + x^{11} + 1}{x^{2} + x + 1},

and it is straightforward to verify that Φ_{33}(x) = C_{33}(x)^{2} − 33xD_{33}(x)^{2}.

3.1. Stevenhagen’s algorithm. Stevenhagen [18] gives a different 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 at most ϕ(n)/2 is generated by the Euclidean algorithm. Thus, the algorithm requires O(n^{2}) arithmetic operations, the same order^{3} as our Algorithm L.

Unfortunately, Stevenhagen’s algorithm suffers from a well-known problem of the Euclidean algorithm [10] – although the initial and final polynomials have small integer coefficients, the intermediate results grow exponentially large. When implemented in 32-bit integer arithmetic, Stevenhagen’s algorithm fails due to integer overflow for n ≥ 35.

Algorithm L does not suffer from this problem. It is easy to see from the recurrences (11)–(12) 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 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 (11) or by 2k + 1 in (12) gives a non-integer result.

3The complexity of both algorithms can be reduced to O(n(log n)^{2}) arithmetic operations by standard “divide and conquer” techniques.

4. Explicit expressions for C_{n} and D_{n}

In this section we give generating functions for the coefficients of C_{n} and D_{n}. These generating functions seem to be new. They can be used to evaluate the coefficients of C_{n}(x) and D_{n}(x) in O(n log n) arithmetic operations, via the fast power series algorithms of [2, Section 5]. Also, where they converge, they give explicit formulas which can be used to compute C_{n}(x) and D_{n}(x) at particular arguments x. However, it may be more efficient to compute the coefficients of the polynomials using Algorithm L, and then evaluate the polynomials by Horner’s rule.

The generating functions may be written in terms of an analytic function g_{n}, which we now define. We continue to assume that (4) holds.

4.1. The analytic functions f_{n} and g_{n}. In this subsection x is a complex variable. For x inside the unit circle, and on the boundary |x| = 1 where the series converge, define

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

and

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

Observe that g_{n}(x) is an odd function, so g_{n}(−x) = −g_{n}(x). Our assumption (4) implies that (n|2j + 1) = (2j + 1|n), so

2g_{n}(x) = f_{n}(x) − f_{n}(−x). (15)

Also, since (2j|n) = (2|n)(j|n) = s′(j|n), it is easy to see that

f_{n}(x) + f_{n}(−x) = s′f_{n}(x^{2}). (16)

In [1] it is shown that analytic continuations of f_{n}(x) and g_{n}(x) outside the unit circle are given by the simple functional equations

f_{n}(x) = f_{n}(1/x), g_{n}(x) = g_{n}(1/x).

f_{n}(1) is related to the class number h(n) of the quadratic field Q[\sqrt{n}] with discriminant n. In the notation of Davenport [7], f_{n}(1) = L_{−1}(1) = L(1) = L(1, χ), where χ(j) = (j|n) is the real, nonprincipal Dirichlet character appearing in (13). Thus, using well-known results,

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

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 we have f_{5}(1) = (ln ε)/\sqrt{5} = 0.4304 …

Using (15)–(16), we obtain a simple relation between g_{n}(1) and f_{n}(1):

g_{n}(1) = (1 − \frac{s′}{2} ) f_{n}(1).

Thus, in our example, g_{5}(1) = 3f_{5}(1)/2.

4.2. Generating functions. In [1] it is shown that

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

This leads to the following theorem. As usual, we continue to assume that n satisfies (4).

Theorem 1. The Aurifeuillian factors F_{n}^{±}(x) = C_{n}(x) ± \sqrt{nx}D_{n}(x) of Φ_{n}(x) are given by

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

Also,

C_{n}(x) = \sqrt{Φ_{n}(x)} cosh (\sqrt{n}g_{n}(\sqrt{x}))

and

D_{n}(x) = \sqrt{\frac{Φ_{n}(x)}{nx}} sinh (\sqrt{n}g_{n}(\sqrt{x})) .

4.3. Application to integer factorization. In this section we illustrate how the results of Sections 3 and 4.2 can be used to obtain factors of integers of the form a^{n} ± b^{n}. Other examples can be found in [1, 3, 4].

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 Φ_{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.

Before giving numerical examples, we state explicitly how Theorem 1 can be used to compute F_{n}^{±}(m^{2}n) with a finite number of arithmetic operations. The following theorem shows how many terms have to be taken in the infinite series (14) defining g_{n}. Because there is a little slack in the proof of the theorem, there is no practical difficulty in evaluating the exponential and square root to sufficient accuracy.

Theorem 2. Let m, n be positive integers, n > 1 square-free, n = 1 mod 4, x = m^{2}n, and λ = ϕ(n)/2. Then the Aurifeuillian factors of Φ_{n}(x) are

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

and

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

where

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

Examples.

3. Consider n = 5, m = 3, so x = m^{2}n = 45 and λ = ϕ(5)/2 = 2. Thus

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

F̂ = \sqrt{Φ_{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 Φ_{5}(x). By division we obtain the other factor 2851. Thus

45^{5} − 1 = 44Φ_{5}(x) = 2^{2} · 11 · 1471 · 2851.

In this example the Aurifeuillian factors are prime.

Table 1. Some Aurifeuillian factorisations

a^{n}  a^{n} − 1
21^{189} 2^{2} · 5 · 43 · 109 · 127 · 163 · 379 · 463 · 631 · 757 · 3319 · 4789 ·

6427 · 51787 · 4779433 · 85775383 · 227633407 · 4167831781 · 22125429901 · 7429452749713 · 27186384126763 · 100595851688887003 · 559529226207687925351 · 592823611828574163154462624637481670158792334981 · P_{60}

33^{99} 2^{5} · 37 · 67 · 199 · 991 · 1123 · 2113 · 19009 · 90619 ·

34905511 · 91402147 · 747487377451 · 4098986195943739 · 987839961952536875400662210432222899 · P_{46}

33^{165} 2^{5} · 31 · 67 · 331 · 1123 · 1321 · 2113 · 4951 · 8581 · 9241 · 39451 ·

90619 · 9540301 · 91402147 · 204970261 · 275465191 · 10125617371 · 47284185301 · 180115639771 · 747487377451 · 4098986195943739 · 11193560623980192151 · 1076141944549238849546221 · 142336076865537701527905793791583051 · P_{44}

77^{77} 2^{2} · 19 · 23 · 617 · 757 · 25411 · 52344007 · 278949511 · 6165802127 ·

12416123247268023977 · 18845698508450782105492211746760179 · P_{53}

97^{97} 2^{5} · 3 · 389 · 363751 · 684640163 · 11943728733741294764390602153 ·

549180361199324724418373466271912931710271534073773 · P_{95}

101^{101} 2^{2} · 5^{2} · 607 · 1213 · 5657 · 157561 ·

9931988588681 · 102208068907493 · 393101595766008847 · 12602965626536109872384216297085760308823294522746017 · P_{89}

105^{105} 2^{3} · 13 · 151 · 211 · 421 · 631 · 1009 · 1201 · 1621 · 2731 · 11131 · 102181 ·

485689 · 18416161 · 1340912959 · 59785910251 · 3662332210521480889 · 23965462949313970771 · 49743995480142943374722277091 · 5384579552746854831338204156683983031 · P_{43}

4. Consider n = 5, m = 40, so x = m^{2}n = 8000 and x^{n} − 1 = 20^{15} − 1. We have

F̂ = \sqrt{Φ_{5}(x)} exp (−\frac{1}{m} + \frac{1}{3m^{3}n})
= 64004000.37 … × 0.9753109279 … = 62423800.99 …

and rounding to the nearest integer gives an Aurifeuillian factor F_{5}^{−} = 62423801 of Φ_{5}(x). By division we obtain the other Aurifeuillian factor F_{5}^{+} = 65624201. Alternatively, we can find the same factors from (6) by evaluating C_{5}(x) and D_{5}(x). Neither of the Aurifeuillian factors is prime, but 20^{15} − 1 also has “algebraic” factors 20^{3} − 1 = 19 · 421 and 20^{5} − 1 = 11 · 19 · 61 · 251. Thus, it is easy to find that F_{5}^{−} = 11 · 19 · 61 · 3001, F_{5}^{+} = 251 · 261451, and

20^{15} − 1 = 11 · 19 · 31 · 61 · 251 · 421 · 3001 · 261451.

5. For an example with larger λ, consider m = 1, n = 13, so x = m^{2}n = 13, x^{n} − 1 = 13^{13} − 1, and λ = ϕ(13)/2 = 6. Theorem 2 gives

F̂ = \sqrt{Φ_{13}(13)} exp (− ∑_{j=0}^{5} \frac{(13|2j + 1)}{(2j + 1)13^{j}} )

= \sqrt{\frac{13^{13} − 1}{13 − 1}} exp (−1 − \frac{1}{3 · 13} + \frac{1}{5 · 13^{2}} + \frac{1}{7 · 13^{3}} − \frac{1}{9 · 13^{4}} + \frac{1}{11 · 13^{5}} )

= 5023902.0906 … × 0.3590131665 … = 1803646.998 …,

and rounding to the nearest integer gives an Aurifeuillian factor F_{13}^{−} = 1803647 of Φ_{13}(13). The same factor could have been found from the polynomials

C_{13}(x) = x^{6} + 7x^{5} + 15x^{4} + 19x^{3} + 15x^{2} + 7x + 1

and

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

It is easy to deduce that

13^{13} − 1 = 2^{2} · 3 · 53 · 264031 · 1803647.

6. An illustrative sample of other factorisations which can be obtained from Algorithm L or Theorem 2, and would have been difficult to obtain in any other way, is given in Table 1. The factors given explicitly in Table 1 are prime. As usual, large k-digit primes are written as P_{k} if they can be found by division.

References

1. R. P. Brent, “On computing factors of cyclotomic polynomials”, Mathematics of Computation, D. H. Lehmer memorial issue, 1993, to appear. Preprint available by anonymous ftp from

dcssoft.anu.edu.au (rpb135.dvi.Z in directory pub/Brent).

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

3. 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. Available by anonymous ftp from dcssoft.anu.edu.au (rpb134.*.Z in directory pub/Brent).

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

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

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

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

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

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

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

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

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

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

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

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

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

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

18. Peter Stevenhagen, “On Aurifeuillian factorizations”, Nederl. Akad. Wetensch. Indag. Math. 49 (1987), 451-468.

Computer Sciences Lab, Australian National University, Canberra, ACT 0200 E-mail address: rpb@cslab.anu.edu.au