nt.number theory - How did Cole factor $2^{67}-1$ in 1903? - MathOverflow

How did Cole factor $2^{67}-1$ in 1903?

I just heard a This American Life episode which recounted the famous anecdote about Frank Nelson Cole factoring N:=267−1 as 193,707,721×761,838,257,287. There doesn't seem to be a historical record of how Cole achieved this; all we have I could find his statement that it took "three years of Sundays". Ira Glass's guest, Paul Hoffman, suggests that this was done by trial division.

But this is nuts, unless I am missing something. Three years of Sundays is 156 days. If he works 10 hours a day, that's 93,600 minutes. There are 10,749,692 primes up to 193,707,721. So that is more than 100 trial divisions a minute. Worse than that, existing prime tables didn't go high enough: According to Chapter XIII of Dickson's History of the Theory of Numbers, existing tables of primes only ran to something like 10,000,000 (664,579 primes), so for the vast majority of the trial divisions, he'd have to find the primes first. (Lehmer, in 1914, went up to 10,006,721.)

But I'm puzzled thinking what else Cole could have done. I skimmed Chapter XIV in Dickson. The methods which seem to have existed at the time are:

• Various ways to speed up trial division for the first 1000's of prime numbers. That only helps at the start.

• Writing N as x2−y2. But y would be 380,822,274,783, which is an even larger search.

• Since 2 is a square modulo N (namely, (234)2≡2modN), we know that all prime factors must be ±1mod8, which cuts the time in half. But that's only a factor of 2.

• Since N≡3mod4, there must be a prime factor which is 3mod4, so we could try only checking those primes. But this turns out to make things worse, since the SMALL factor is the one which is 1mod4.

• If we could write N as a sum of squares in two ways, we'd be done. But N isn't a sum of squares, since it is 3mod4. Generalizations to other positive definite quadratic forms were known at the time, but how would Cole know which quadratic form to try?

• A variant of the above would be to use the quadratic form 2x2−y2, since we already have one solution. Dickson doesn't mention any work using mixed signature forms, but it would work. And since Z[√2] is a PID, there must be a second way to write N as 2x2−y2, not related to the previous by units of Z[√2]. I'm not sure how large this second solution is.

So, my question is:

How could someone find the prime factors of N in 100,000 minutes of hand computation?