Programming Praxis

Here is our version of the quadratic sieve:

(define (qs n f m)
  (let* ((e 10) ; error permitted on sum of logarithms
         (sqrt-n (isqrt n)) (b (- sqrt-n m)) (fb (factor-base n f))
         (sieve (make-vector (+ m m) (- e (inexact->exact (round (log (* 2 sqrt-n)))))))
         (ts (map (lambda (f) (msqrt n f)) (cdr fb))) ; exclude 2
         (ls (map (lambda (f) (inexact->exact (round (log f)))) (cdr fb))))
    (when verbose? (display "Factor base of ") (display (length fb))
                   (display " primes") (newline))
    (do ((fb (cdr fb) (cdr fb)) (ts ts (cdr ts)) (ls ls (cdr ls))) ((null? fb))
      (do ((i (modulo (- (car ts) b) (car fb)) (+ i (car fb)))) ((<= (+ m m) i))
        (vector-set! sieve i (+ (vector-ref sieve i) (car ls))))
      (do ((i (modulo (- (- (car ts)) b) (car fb)) (+ i (car fb)))) ((<= (+ m m) i))
        (vector-set! sieve i (+ (vector-ref sieve i) (car ls)))))
    (let loop ((i 0) (rels (list)))
      (if (= i (+ m m))
          (begin
            (when verbose? (display "Found ") (display (length rels))
                           (display " smooth relations") (newline))
            (solve n fb rels))
          (if (positive? (vector-ref sieve i))
              (let ((ys (smooth (- (square (+ i b)) n) fb)))
                (if (pair? ys)
                    (loop (+ i 1) (cons (cons (+ i b) ys) rels))
                    (loop (+ i 1) rels)))
              (loop (+ i 1) rels))))))

We made two changes from the previous description. First, we introduced an error term e on the sum of the logarithms, which partially covers the case where we don’t sieve on prime powers. Second, we decided not to sieve on the factor base prime 2; it requires additional code to find the starting point, it takes much time to access m sieve locations, and it contributes little to the sum of the logarithms.

Sieving is performed in the three do loops. The outer do iterates over the factor base primes, excluding 2. The first inner do iterates over soln1 and the second inner do iterates over soln2. The named-let scans the sieve, adding smooth relations to the rels list, and calls solve to perform the linear algebra after the scan is complete.

The rest of the code is straight forward. The factor-base function iterates over the primes p less than f, returning those for which (jacobi n p) equals 1. The smooth function performs trial division and returns a list of prime factors (possibly including -1) if the input is smooth over the factor base or a null list if it is not smooth. We steal much from previous exercises, which you can see assembled at http://programmingpraxis.codepad.org/k6dqWasV.

Here are some examples:

> (set! verbose? #t)
> (qs 87463 30 30)
Factor base of 6 primes
Found 6 smooth relations
587
> (qs 13290059 150 300)
Factor base of 18 primes
Found 23 smooth relations
4261
> (qs 294729242679158229936006281 2000 3000000)
Factor base of 149 primes
Found 156 smooth relations
2971215073

There are several improvements that can be made to this basic version of the quadratic sieve, which we will discuss in future exercises.

Pages: 1 2 3