Factoring With Bicycle Chains
April 4, 2014
As penance for teasing you on April Fool’s Day, I have another exercise on factoring. Today we will simulate a machine invented by Derrick H Lehmer (the son of the famous father-son-daughter-in-law trio of mathematicians named Lehmer) that uses bicycle chains to factor integers. Unlike the previous exercise, this factoring method is real; in fact, a variant of this method was the best available method to factor integers from the mid-1930s to the mid-1960s, and Lehmer’s machine found many valuable factorizations for the Cunningham Project. We’ll discuss a simplified version of the machine today, as described by Uli Meyer, and the full version of the machine in a future exercise. I am unable to find a copyright-free picture of Lehmer’s machine, but if you’re interested you could visit the Computer History Museum, or see Lehmer himself describe the machine.
The machine is based on Fermat’s method of factoring by a difference of squares: if n = x2 − y2 then the factors of n are (x − y) and (x + y). Fermat’s method starts with y = ⌈ √n ⌉, computes y2 − n, and stops when the result is a square.
Lehmer’s machine identifies likely squares by looking at quadratic residues. An integer a is a quadratic residue modulo n if there exists some x such that x2 ≡ a (mod n). For a number y2 − n to be a square in Fermat’s algorithm, it must be a quadratic residue to many small bases. Lehmer’s machine represents quadratic residues as “open” points on a bicycle chain, then spins many bicycle chains of different lengths along a common axis and stops when all are quadratic residues; in many but not all cases, such a number will will indicate a factor of n.
Uli Meyer gives a better description, shows a useful diagram, and provides photos of his Lego sieve machine.
Your task is to write a program that simulates Lehmer’s machine to factor integers. 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.
Programming Praxis 
Typo in the second paragraph: I believe you mean n = x² − y²?
@Lucas: Fixed. Thank you.
Here’s my work so far… It doesn’t seem to work, but I can’t figure out what’s wrong.
from itertools import cycle def chain(n, l): # 0: The input is the number to be factored and the length of the chain to construct. # 1: Compute the quadratic residues modulo l: "Create a vector qs of length l, initialize each element of the vector to # False. Then for each i from 0 to l-1, calculate k = i**2 (mod l) and set the kth element of the vector to True, # because k is a quadratic residue mod n. qs = [False] * l k = 0 for i in xrange(l): qs[k] = True k += 2*i + 1 k %= l # 2: Compute pegs: "Create a vector ps of length l, initialize eac helement of the vector to False. Initialize # a = ceil(sqrt(n)) where n is the number being factored. For each i from 0 to l-1 calculate k = (a+i)**2 - n mod l and # set the ith element of the vector ps to the value of the kth element of the qs vector. Then form the ps vector into a # cyclical list (the chain) where the pegs are at the places where the vector is True." ps = [False] * l a = isqrt(n-1) + 1 # ceil(sqrt(n)) for i in xrange(l): k = ( (a+i)**2 - n ) % l # TODO: rewrite this in incremental fashion as in the previous loop. + beats * in sieves. ps[i] = qs[k] return cycle(ps) def lehmer_bikechain(n, chains=[64,27,25,49,22,26,17,19,23,29,31,37,41,43,47,53,59,61,67]): # 0: The input is the number to be factored and the list of chain lengths to be used in the sieve. # 1: Replace the chain lengths with the cycled, pegged chain of that length. chains = [chain(n, l) for l in chains] # 2: Sift. "Set a = ceil(sqrt(n)) and set all chains at zero. Then loop forever: when all chains are True, calculate # b**2 = a**2 - n and, if b is an integer, report factors a-b and a+b and halt. If any chain is False, or if all chains # are True but b is not an integer, increment a by 1, advance all the chains by one link, and loop." a = isqrt(n-1) + 1 # ceil(sqrt(n)) while True: if all(c.next() for c in chains): b2 = a**2 - n b = isqrt(b2) print a, b, b2 - b**2, a+b, a-b if b**2 == b2: return a-b a += 1 n = 25709599 * 78539161 print lehmer_bikechain(n)Figured out what I did wrong. Using “all()” at line 38 doesn’t advance every chain; it just advances them until it reaches a chain that reads “False” and then stops, completely ignoring the rest of the chains. Here’s the fixed code:
from itertools import cycle def chain(n, l): qs = [False] * l k = 0 for i in xrange(l): qs[k] = True k += 2*i + 1 k %= l ps = [False] * l a = isqrt(n-1) + 1 for i in xrange(l): k = ( (a+i)**2 - n ) % l ps[i] = qs[k] return cycle(ps) def lehmer_bikechain(n, chains=[64,27,25,49,22,26,17,19,23,29,31,37,41,43,47,53,59,61,67]): chains = [chain(n, l) for l in chains] a = isqrt(n-1) + 1 # ceil(sqrt(n)) while True: if reduce(lambda a, b: a and b, [c.next() for c in chains], True): b2 = a**2 - n b = isqrt(b2) print a, b, b2 - b**2, a+b, a-b if b**2 == b2: return gcd(a-b, n) a += 1 n = 25709599 * 78539161 print lehmer_bikechain(n)And here it is as a single function:
from itertools import cycle def lehmer_bikechain(n, chains=[64,27,25,49,22,26,17,19,23,29,31,37,41,43,47,53,59,61,67]): #chains = [chain(n, l) for l in chains] for j in xrange(len(chains)): l = chains[j] qs = [False] * l k = 0 for i in xrange(l): qs[k] = True k += 2*i + 1 k %= l ps = [False] * l a = isqrt(n-1) + 1 for i in xrange(l): k = ( (a+i)**2 - n ) % l # TODO: rewrite in incremental fashion as above. + beats * in sieves. ps[i] = qs[k] chains[j] = cycle(ps) a = isqrt(n-1) + 1 # ceil(sqrt(n)) while True: if reduce(lambda a, b: a and b, [c.next() for c in chains], True): b2 = a**2 - n b = isqrt(b2) if b**2 == b2: return gcd(a-b, n) a += 1 n = 25709599 * 78539161 print lehmer_bikechain(n)When executed with CPython, this finishes in about 37 seconds on my computer while PyPy only needs 4.6 seconds.