Programming Praxis

In the seventeenth century, the French jurist and amateur mathematician Pierre de Fermat investigated numbers of the form 2²ⁿ+1, now known as Fermat numbers: 3, 5, 17, 257, 65537, 4294967297, 18446744073709551617, 340282366920938463463374607431768211457, … (Sloane’s A000215). Fermat noted that the first five numbers in the sequence, F₀ to F₄, were prime, and conjectured that all the remaining numbers in the sequence were prime. He was wrong, as Leonhard Euler proved in 1732 with the factorization of F₅ = 641 · 6700417; then Thomas Clausen factored F₆ = 274177 · 67280421310721 in 1855. Those first five numbers are the only members of the sequence that are known to be prime, and it now conjectured that all the remaining Fermat numbers are composite. There is something of a cottage industry among mathematicians and programmers to find new factors of Fermat numbers; you can find the current status of that friendly competition at http://www.prothsearch.net/fermat.html.

On September 13th, 1970, Michael A. Morrison and John Brillhart factored the seventh Fermat number: F₇ = 2²⁷+1 = 2¹²⁸+1 = 340282366920938463463374607431768211457 = 59649589127497217 · 5704689200685129054721. Their feat inaugurated the modern culture of factorization by computer, and their now-deprecated method was the direct predecessor of the quadratic sieve that is today the most powerful factorization method available for personal computers, able to factor numbers up to about a hundred digits (the most powerful factoring algorithm known is the number field sieve, which can factor numbers up to about 150 digits, but requires a large cluster of computers). On the fortieth anniversary of their achievement we recreate their computation in this exercise and the next. We are following their description of their work in “A Method of Factoring and the Factorization of F₇” in Mathematics of Computation, Volume 29, Number 129, January 1975, Pages 183-205, available on-line at http://www.ams.org/journals/mcom/1975-29-129/S0025-5718-1975-0371800-5/S0025-5718-1975-0371800-5.pdf.

The math behind the method dates to Fermat, who noted that a factor of N could be found when x² – N was a perfect square; thus, Fermat’s factorization method started at the square root of N and worked x backwards until it reached a solution, as we did in a previous exercise. Early in the twentieth century, Maurice Kraitchik, a Ukrainian emigré and recreational mathematician living in Paris, observed that Fermat’s formula could be written as x² ≡ y² (mod N). Then, in 1931, D. H. Lehmer and R. E. Powers discovered that the convergents of the continued fraction representation of the square root of N provided a suitable x and y, but their method was considered unusable because the calculations were too tedious and time-consuming for manual calculation, and the method frequently failed to find a factorization. In 1965, Brillhart realized that the tedious calculations were fit for computer calculation, and he and Morrison spent the summer of 1970 developing their method.

The continued fraction method works in two stages, a first stage that computes successive convergents of the continued fraction representing the square root of the number being factored, and a second stage that performs linear algebra on the factors of the convergents and computes a factor of the original number. Morrison and Brillhart, using an IBM 360/91 mainframe computer at UCLA, coding in PL/1 with assembly-language libraries for large-integer processing, divided the work into programs RESIDUE and ANSWER to make it fit into the machine. We will likewise divide the work into two exercises, looking at RESIDUE today and ANSWER in the next exercise.

We saw continued fractions previously, in the exercise on the golden ratio. Square roots can be computed by continued fractions, and are always periodic; we’ll refer to the article by Morrison and Brillhart for the math and skip directly to their algorithm to expand √N, or √kN for some suitable multiplier k ≥ 1, into a simple continued fraction:

Initialize A₋₂ = 0, A₋₁ = 1, Q₋₁ = kN, r₋₁ = g, P₀ = 0, Q₀ = 1, and g = ⌊√kN⌋.

For each n from 0 to a user-specified limit:

• Compute q_n and r_n using the formula g + P_n = q_nQ_n + r_n, where 0 ≤ r_n < Q_n.
• Compute A_n (mod N) using the formula A_n ≡ q_nA_n−1 + A_n−2 (mod N).
• Compute g + P_n+1 using the formula g + P_n+1 = 2g − r_n.
• Compute Q_n+1 using the formula Q_n+1 = Q_n−1 + q_n(r_n −r_n−1).

Some of the Q_n must be factored. That “some” sounds odd. The idea is to factor those Q_n that have all their prime factors less than some pre-specified limit, said primes having a Jacobi symbol of 0 or 1 indicating that they are quadratic residues of kN and thus potential factors; we saw the Jacobi symbol in the exercise on modular arithmetic. The procedure is to choose a limit, compute a factor base of those primes that are potential factors, then use trial division to determine if the Q_n should be added to the out-going list; Morrison and Brillhart call such a Q_n an A − Q pair.

Morrison and Brillhart, based on Lehmer and Powers before them, give this example: let N = 13290059 and k = 1; then g = 3645. Then selected results from the expansion of √kN are given below:

The output from RESIDUE is a list of the following items for each Q_n that could be factored over the factor base: n, A_n−1, Q_n, and a list of the odd-power primes dividing Q_n (thus, factors like 5² of Q₃₁ are omitted); forty years ago, the output was given on punched cards! In the example above, the rows for Q₅, Q₁₀, Q₂₂, Q₂₃, Q₂₆, Q₃₁, and Q₄₀ should appear in the output, based on a factor base containing 2, 5, 31, 41, 43, 53, and 113. Note that 5 is part of the factor base even though it doesn’t appear as an odd power in any of the Q_n.

Your task is to write the RESIDUE function that computes a factor base and returns a list of factored Q_n; we will complete the factorization of F₇ in the next exercise. Apply your function to F₇ using a factor base of the first 2700 suitable primes less than 60000 and a multiplier of 257. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2