## Quadratic Sieve

### 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.

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