### March 19, 2013

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

### One Response to “Quadratic Sieve”

1. […] Quadratic Sieve (programmingpraxis.com) […]