In the seventeenth century, the French jurist and amateur mathematician Pierre de Fermat investigated numbers of the form 22n+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, F0 to F4, 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 F5 = 641 · 6700417; then Thomas Clausen factored F6 = 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: F7 = 227+1 = 2128+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 F7” 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 x2N 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 x2y2 (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−2 = 0, A−1 = 1, Q−1 = kN, r−1 = g, P0 = 0, Q0 = 1, and g = ⌊√kN⌋.

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

  • Compute qn and rn using the formula g + Pn = qnQn + rn, where 0 ≤ rn < Qn.
  • Compute An (mod N) using the formula AnqnAn−1 + An−2 (mod N).
  • Compute g + Pn+1 using the formula g + Pn+1 = 2grn.
  • Compute Qn+1 using the formula Qn+1 = Qn−1 + qn(rnrn−1).

Some of the Qn must be factored. That “some” sounds odd. The idea is to factor those Qn 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 Qn should be added to the out-going list; Morrison and Brillhart call such a Qn an AQ 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:

n g+Pn Qn qn rn An−1 (mod N) Qn factored
−1 13290059 3645 0
0 3645 1 3645 0 1
1 7290 4034 1 3256 3645 2 · 2017
2 4034 3257 1 777 3646 3257
3 6513 1555 4 293 7291 5 · 311
4 6997 1321 5 392 32810 1321
5 6898 2050 3 748 171341 2 · 52 · 41
10 6318 1333 4 986 6700527 31 · 43
22 4779 4633 1 146 5235158 41 · 113
23 7144 226 31 138 1914221 2 · 113
26 5622 3286 1 2336 11455708 2 · 31 · 53
31 6248 5650 1 598 1895246 2 · 52 · 53
40 6576 4558 1 2018 3213960 2 · 43 · 53
52 7273 25 290 23 2467124 52

The output from RESIDUE is a list of the following items for each Qn that could be factored over the factor base: n, An−1, Qn, and a list of the odd-power primes dividing Qn (thus, factors like 52 of Q31 are omitted); forty years ago, the output was given on punched cards! In the example above, the rows for Q5, Q10, Q22, Q23, Q26, Q31, and Q40 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 Qn.

Your task is to write the RESIDUE function that computes a factor base and returns a list of factored Qn; we will complete the factorization of F7 in the next exercise. Apply your function to F7 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

Follow

Get every new post delivered to your Inbox.

Join 574 other followers