## Factoring With Bicycle Chains

### April 4, 2014

The algorithm that determines where to place the pegs in a chain of length *len* consists of two steps:

1) Compute quadratic residues: Create a vector *qs* of length *len*, initialize each element of the vector to `#f`

. Then for each *i* from 0 to *len*−1, calculate *k* = *i*^{2} (mod *len*) and set the *k*th element of the vector to `#t`

, because *k* is a quadratic residue mod *n*.

2) Compute pegs: Create a vector *ps* of length *len*, initialize each element of the vector to `#f`

. Initialize *a* = ⌈ √*n* ⌉ where *n* is the number being factored. For each *i* from 0 to *len*−1, calculate *k* = (*a*+*i*)^{2} − *n* (mod *len*) and set the *i*th element of the vector *ps* to the value of the *k*th 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 `#t`

.

Here is our version of `chain`

:

`(define (chain len)`

(let ((a (+ (isqrt n) 1))

(qs (make-vector len #f)) ; quadratic residue?

(cs (make-vector len #f))) ; chain link has peg?

(do ((i 0 (+ i 1))) ((= i len))

(vector-set! qs (modulo (square i) len) #t))

(do ((i 0 (+ i 1)))

((= i len) (cycle (vector->list cs)))

(vector-set! cs i

(vector-ref qs (modulo (- (square (+ a i)) n) len))))))

To perform the sieving: Set *a* = ⌈ √*n* ⌉ and set all chains at zero. Then loop forever: when all chains are `#t`

, 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 `#f`

, or if all chains are `#t`

but *b* is not an integer, increment *a* by 1, advance all the chains by one link, and loop. Here is our version of the machine:

`(define (bicycle n chains)`

(let loop ((a (+ (isqrt n) 1))

(cs (map (lambda (len) (chain n len)) chains)))

(if (not (every? (map car cs))) (loop (+ a 1) (map cdr cs))

(let* ((b2 (- (square a) n)) (b (isqrt b2)))

(if (= (square b) b2) (list (- a b) (+ a b))

(loop (+ a 1) (map cdr cs)))))))

We use the same set of chains as Lehmer:

`(define chains ; lehmer's 19-chain machine`

'(64 27 25 49 22 26 17 19 23 29 31 37 41 43 47 53 59 61 67))

Here’s an example:

`> (time (bicycle 2019210335106439 chains))`

(time (bicycle 2019210335106439 ...))

518 collections

2996 ms elapsed cpu time, including 15 ms collecting

3002 ms elapsed real time, including 47 ms collecting

4370972848 bytes allocated, including 4365536448 bytes reclaimed

(25709599 78539161)

We used `isqrt`

, `square`

, `every?`

, `last-pair`

and `cycle`

from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/twQwAWN5.

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.

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:

And here it is as a single function:

When executed with CPython, this finishes in about 37 seconds on my computer while PyPy only needs 4.6 seconds.