Programming Praxis

Living in the age of computers, I’ve always been fascinated by how much arithmetic it was possible to do before the advent of computers. The Egyptians knew the peasant method of multiplication when the built the pyramids. Archimedes knew the value of pi two centuries before Christ, with an approximation we still use today. Heron had an algorithm for calculating square roots two thousand years ago. And so on.

By the time of Fermat, Euler and Gauss, arithmetic was done with simple algorithms, rules of thumb, and large tables. Sadly, the tables were frequently wrong; Charles Babbage invented his Difference Engine for the purpose of automatically constructing accurate tables.

We studied Fermat’s algorithm for factoring integers in a previous exercise:

function fermat(n)
    x := floor(sqrt(n))
    r := x * x - n
    t := 2 * x + 1
    while r is not a square
        r := r + t
        t := t + 2
    x := (t - 1) / 2
    y := sqrt(r)
    return x - y, x + y

But how did Fermat actually work his algorithm? In addition to the basic arithmetic, he used some simple rules for identifying squares and a table of the squares of numbers similar to that on the next page, with the largest square at least as large as n.

Let’s consider the factorization of n = 13290059 as an example. First, Fermat removed small factors by trial division, say for factors less than 30 or 50. Then, Fermat looked up n in the table. If n appeared in the body of the table, then it is a perfect square, and the square root is its factor. But in the general case, Fermat finds x as the square root of the next smaller entry in the table. For 13290059, x = 3645, because 3645² < 13290059 < 3646². Then t = 7291 and r = -4034, which is not a perfect square (because it’s negative).

Now the iteration begins. Fermat calculated r and t in 45 steps as follows, first r, then t below it, then a sum line, then the sum r = r + t, then t incremented by 2 from the previous t, then a sum line, then a new sum r, at each step checking if r is a square:

0:  r =  -4034 = 3645 * 3645 - 13290059
    t =   7291 = 2 * 3645 + 1
        ------
1:  r =   3257 ends in 7; not a square
    t =   7293 = 7291 + 2
        ------
2:  r =  10550 ends in 50; not a square
    t =   7295 = 7293 + 2
        ------
3:  r =  17845 ends in 45; not a square
           ...
44: r = 318662 ends in 2; not a square
    t =   7381 = 7379 + 2
        ------
45: r = 326041 in table: 571 * 571

At that point, Fermat could calculate the factors of n: x = (7381 – 1) / 2 = 3690, y = 571, and the factors of n are x − y = 3690 − 571 = 3119 and x + y = 3690 + 571 = 4261.

The key to Fermat’s calculation was identifying the case where r is a square. Fermat used three rules: First, the last two digits of the number had to be 00, e1, e4, 25, o6 or e9, where e is any even digit and o is any odd digit. Second, the digital root of the number had to be 1, 4, 7 or 9: sum the digits of n; if the result is less than 10, then return the sum of the digits, otherwise return the digital root of the sum of the digits. In practice, this goes much more quickly than it appears by “casting out nines” wherever they appear in the digits of n or the result. Third, if he had not already determined that the number was not a square, he looked it up in his table of squares. The whole calculation goes much more quickly than I can describe it here.

And that’s how Fermat factored integers. Of course, sometimes he failed, when it took too many iterations to find a square, or when his table of squares was too small, or when he made an error in the arithmetic. Fortunately, it was easy enough for him to confirm that any factorization he found was correct by multiplying the factors.

Your task is to write a program that factors integers the same way Fermat factored integers. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the program in the comments below.

Pages: 1 2 3