## How Fermat Factored Integers

### July 18, 2014

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^{2} < 13290059 < 3646^{2}. 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, *e*1, *e*4, 25, *o*6 or *e*9, 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.

## How Many Distinct Products In A Times Table?

### July 15, 2014

We all learned our times tables in elementary school. For instance, here is the standard 9 by 9 times table:

`1 2 3 4 5 6 7 8 9`

2 4 6 8 10 12 14 16 18

3 6 9 12 15 18 21 24 27

4 8 12 16 20 24 28 32 36

5 10 15 20 25 30 35 40 45

6 12 18 24 30 36 42 48 54

7 14 21 28 35 42 49 56 63

8 16 24 32 40 48 56 64 72

9 18 27 36 45 54 63 72 81

That table has 36 distinct products: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 24, 25, 27, 28, 30, 32, 35, 36, 40, 42, 45, 48, 49, 54, 56, 63, 64, 72 and 81.

Your task is to write a program that computes the number of distinct products in an *n* by *n* times table. There is an obvious O(*n*^{2}) algorithm; can you do better? 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.

## FSM Generator

### July 11, 2014

A few weeks ago we had an exercise that solved the problem of removing singleton occurrences of a given character from a string while leaving multiple occurrences of the given character intact. Our solution used a finite state machine that iterated over the input string writing output as it went.

The finite state machine used in that solution was hand-crafted to solve the given problem. But it is possible to write a program that takes a definition of a finite state machine and creates a function that takes an input string and applies the finite state machine to it. Such a program isn’t hard to create and provides a useful utility for small parsing problems.

Your task is to write a program that generates a finite state machine from a simple description; the format and capabilities of the machine are up to you, but it should be at least powerful enough to solve the “remove singleton” problem. 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.

## SQUFOF, Continued Fraction Version

### July 8, 2014

Square Forms Factorization, abbreviated SQUFOF, is a special-purpose method for factoring integers. It is commonly used to split small semi-primes in the double-large-prime variants of the quadratic sieve and number field sieve factoring algorithms. We studied the “binary quadratic forms” variant of SQUFOF in a previous exercise; today we study the “continued fraction” variant of SQUFOF.

Although the mathematical basis of SQUFOF is the quadratic forms of Gauss, SQUFOF was discovered by Daniel Shanks while he was investigating the failures of the continued fraction factoring algorithm; you should read his paper, which is as good a demonstration of the process of mathematical discovery as anything I have ever read. Shanks’ algorithm expands the continued fraction of the square root of *n*, the number being factored, searching among the convergents for a quadratic form that is a perfect square; in this, it differs from the CFRAC algorithm, which combines multiple convergents to assemble a perfect square. We copy the algorithm from the paper by Jason Gower and Sam Wagstaff that provides the standard description of SQUFOF:

[ Note: In the description that follows, be careful to distinguish the symbols

Q̂andQ; the first one has a circumflex over theQ, the second one doesn’t. Some browsers render the “combining circumflex accent” readably, but some don’t. On the browsers available to me, only one out of four rendered the symbol readably, sorta kinda, and only when I expanded the font to be quite large. If in doubt, look at the program on the next page, where the first symbol is named`qhat`

and the second symbol is named`q`

. ]The variable

Nis the integer to factor,Sremembersq_{0},qholds the currentq,_{i}PandP’hold two consecutive values ofP,_{i}Q̂andQhold two consecutive values ofQ, and_{i}tis a temporary variable used in updatingQ. Finally,Bis an upper bound on the number of forms tested for being square forms before the algorithm gives up.

1. Initialize:Read the odd positive integerNto be factored. IfNis the square of an integer, output the square root and stop. IfN≡ 1 mod 4, then setD← 2N; otherwise, setD←N. In any case, setS← ⌊√D⌋,Q̂ ← 1,P←S,Q←D−P·P,L← ⌊2√2√D⌋,B← 2 ·L, andi← 0.

2. Cycle forward to find a proper square form:Steps 2a through 2e are repeated fori= 1, 2, 3, ….

2a:Setq← ⌊(S+P) /Q⌋ andP’←q·Q−P.

2b:IfQ≤l, then:

IfQis even, put the pair (Q/ 2,Pmod (Q/ 2) onto the QUEUE; otherwise, ifQ≤L/ 2, then put the pair (Q,PmodQ) onto the QUEUE.

2c:Sett←Q̂ +q· (P−P’),Q̂ ←t, andP←P’.

2d:Ifiis odd, go to Step 2e. IfQis not the square of an integer, go to Step 2e. Otherwise, setr← √Q, a positive integer. If there is no pair (r,t) in the QUEUE for whichrdividesP−t, then go to Step 3. Ifr> 1 and there is a pair (r,t) in the QUEUE for whichrdividesP−t, then remove all pairs from the beginning of the QUEUE up to and including this pair and go to Step 2e. Ifr= 1 and there is a pair (1,t) in the QUEUE, then the algorithm fails because we have traversed the entire principal period of quadratic forms of discriminant 4Nwithout finding a proper square form.

2e:Leti←i+ 1. Ifi>B, give up. Otherwise, go to Step 2a.

3. Compute an inverse square root of the square form:SetQ̂ ←r,P←P+r· ⌊(S−P) /r⌋, andQ← (D−P·P) /Q̂.

4. Cycle in the reverse direction to find a factor ofN:

4a:Setq← ⌊(S+P) /Q⌋ andP’←q·Q−P.

4b:IfP=P’, then go to Step 5.

4c:Sett←Q̂ +q· (P−P’,Q̂ ←Q,Q←t, andP←P’and go to Step 4a.

5. Print the factor ofIfN:Qis even, setQ←Q/ 2. Output the factorQofN.

Later, in Section 5.2 of their paper, Gowers and Wagstaff give instructions for using a multiplier *m* to help find a factor of *N*, changing the algorithm given above as follows:

Change Step 1 to this: Read the odd positive integer

Nto be factored. IfNis the square of an integer, output the square root and stop. IfmN≡ 1 mod 4, then setD← 2mN; otherwise, setD←mN. In any case, setS← ⌊√D⌋,Q̂ ← 1,P←S,Q←D−P·P,L← ⌊2√2√D⌋,B← 2 ·L, andi← 0.Change Step 2b to: Let

g←Q/ gcd(Q, 2m). Ifg≤L, put (g,Pmodg) onto the QUEUE. [ Note: If you’re reading the linked version of the Gowers and Wagstaff paper, this formula is erroneously given as (g,Bmodg); that error cost me several hours of pain. ]In Step 5, replace

QbyQ/ gcd(Q, 2m) before the factorQis output.

To find a factor, first use multiplier *m* = 1. If that fails, try another multiplier from the set 3, 5, 7, 11, 3 · 5 = 15, 3 · 7 = 21, 3 · 11 = 33, 5 · 7 = 35, 5 · 11 = 55, 7 · 11 = 77, 3 · 5 · 7 = 105, 3 · 5 · 11 = 165, 3 · 7 · 11 = 231, 5 · 7 · 11 = 385, 3 · 5 · 7 · 11 = 1155, or any other squarefree product of small odd primes.

Your task is to implement SQUFOF with multipliers. 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.

## Minimum And Maximum

### July 4, 2014

We have today a two-part interview question:

First, write a program that takes an array of integers and returns the minimum and maximum elements of the array. Use as few comparisons as possible.

Second, write a program that takes an array of integers and returns the second minimum and maximum elements of the array; that is, the second-smallest element and the second-largest element. Again, use as few comparisons as possible.

Your task is to write the two programs described above, and state the number of comparisons each makes in the general case; be sure they are proof against malicious input. 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.

## Moonrise, Moonset

### July 1, 2014

We have previously studied the times of sunrise and sunset and the phases of the moon. Today we look at the times of moonrise and moonset.

The standard source for astronomical calculations is the Naval Observatory; they point to a 1989 article in *Sky & Telescope* for the calculation. The article isn’t online, but the code is: a horribly ugly BASIC program, reproduced on the next page.

Your task is to write a program that calculates the times of moonrise and moonset. 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.

## Raindrops

### June 27, 2014

We’re going to try something different today. We often have interview questions here, but always of the type that require writing a program. Today, we will have one of the brain-teaser type of interview question:

How many raindrops fall on the planet every year?

Your task is to estimate the answer to the question posed above. When you are finished, you are welcome to read a suggested solution or to discuss your solution in the comments below.

## Busy Beaver

### June 24, 2014

Yesterday was Alan Turing’s birthday, so today we will write a program for a turing machine.

The Busy Beaver is a turing machine that performs the maximum work for a given configuration of machine; the concept was invented by Tibor Radó in his 1962 paper *On Non-Computable Functions*. We won’t look at the underlying theory, although it is fascinating if you have the time. Instead, we’ll be content to implement the first few busy beavers. Here are the two-symbol busy beavers for one, two, three and four states, from Wikipedia; the two symbols are 0 and 1, the states are letters A through D plus the halting state H, and the moves L and R are left and right:

A | |

0 | 1RH |

1 | unused |

A | B | |

0 | 1RB | 1LA |

1 | 1LB | 1RH |

A | B | C | |

0 | 1RB | 0RC | 1LC |

1 | 1RH | 1RB | 1LA |

A | B | C | D | |

0 | 1RB | 1LA | 1RH | 1RD |

1 | 1LB | 0RC | 1LD | 0RA |

Your task is to implement the four busy beavers. 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.

## Birthday Paradox

### June 20, 2014

The birthday paradox, which we studied in a previous exercise, states that in any group of 23 people there is a 50% chance that two of them share a birthday. The BBC recently published an article that shows 16 of the 32 World Cup teams, each consisting of 23 players, have shared birthdays, thus demonstrating the paradox precisely. Today’s exercise asks you to recreate their calculation.

You can obtain the same listing of player birthdays that the BBC used from FIFA. Another source is the player rosters at WikiPedia.

Your task is to demonstrate that the 2014 World Cup rosters honor the birthday paradox. 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.

## Collinearity

### June 17, 2014

Beware: today’s exercise sounds simple but is actually quite complex if you don’t look at it properly.

Your task is to write a function that takes three points in the *x*,*y* plane and determines if they are collinear; be sure to handle vertical lines and horizontal lines properly. 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.