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 = i2 (mod len) and set the kth 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)2n (mod len) 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 #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 b2 = a2n and, if b is an integer, report factors ab 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.

About these ads

Pages: 1 2

5 Responses to “Factoring With Bicycle Chains”

  1. Lucas A. Brown said

    Typo in the second paragraph: I believe you mean n = x² − y²?

  2. programmingpraxis said

    @Lucas: Fixed. Thank you.

  3. Lucas A. Brown said

    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)
    
  4. Lucas A. Brown said

    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)
    
  5. Lucas A. Brown said

    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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 609 other followers

%d bloggers like this: