Multiple Polynomial Quadratic Sieve
June 21, 2013
We’ll only hit the highlights here; see the next page for the full code.
We begin with the sieve function, which changes from the prior versions. We are back to a straight sieve, rather than a segmented sieve, both for clarity of code and also because the sieve is very much smaller with mpqs than with the basic qs. The primary difference is in the way that a relation is built, which we packaged in a local function:
(define (sieve n fb ts ls e f m a b q-inv) ; => (values rels parts)
(define (make-rel x ys)
(cons (modulo (* (+ (* a x) b) q-inv) n) ys))
; a relation, whether full or partial, has (ax+b)/q in its car and a
; list of factors of g(x)=ax^2+2bx+c, in descending order, in its cdr;
; a large prime, if it exists, is at the cadr of the relation
(let* ((c (/ (- (* b b) n) a)) (rels (list)) (parts (list))
(sieve (make-vector (+ m m 1) (+ e e))))
(do ((fb (cdr fb) (cdr fb)) (ts ts (cdr ts)) (ls ls (cdr ls))
(invs (map (lambda (f) (inverse a f)) (cdr fb)) (cdr invs)))
((null? fb))
(let ((f (car fb)) (t (car ts)) (l (car ls)) (v (car invs)))
(do ((i (modulo (* v (- t b)) f) (+ i f))) ((<= (+ m m 1) i))
(vector-set! sieve i (+ (vector-ref sieve i) l)))
(do ((i (modulo (* v (- (- t) b)) f) (+ i f))) ((<= (+ m m 1) i))
(vector-set! sieve i (+ (vector-ref sieve i) l)))))
(do ((i 0 (+ i 1))) ((= (+ m m 1) i))
(let* ((x (- i m)) (g (+ (* a x x) (* 2 b x) c)))
(if (< (ilog 2 g) (vector-ref sieve i))
(let* ((ys (smooth g fb)) (rel (make-rel x ys)))
(if (pair? ys)
(if (<= (car ys) f)
(set! rels (cons rel rels))
(set! parts (cons rel parts))))))))
(values rels parts)))
If you look closely, you will see the calculation of the modular inverses of all the primes in the factor base, repeated for each polynomial. This calculation is a major part of the cost of sieving; in a future exercise we will see a way to nearly eliminate that cost.
The other big change is the qs function, which now has to compute q, a, b and c. Variable q counts up from the square root of n.
(define (qs n f m)
(define (make-odd q) (if (odd? q) q (+ q 1)))
(call-with-values
(lambda () (factor-base n f))
(lambda (fb ts ls e len-fb)
(when verbose? (display "Factor base of ")
(display len-fb) (display " primes.") (newline))
(let loop ((q (make-odd (isqrt (quotient (isqrt (+ n n)) m)))) (rels (list)) (parts (list)))
(if (not (and (prime? q) (= (jacobi n q) 1))) (loop (+ q 2) rels parts)
(let* ((a (* q q)) (b (apply min (lift-root n q))) (q-inv (inverse q n)))
(when verbose? (display "q=") (display q) (display ": "))
(call-with-values
(lambda () (sieve n fb ts ls e f m a b q-inv))
(lambda (rs ps)
(let* ((rels (append rs rels)) (parts (append ps parts))
(matches (match parts)) (len-rels (length rels))
(len-parts (length parts)) (len-matches (length matches)))
(when verbose? (display len-rels) (display " smooths, ")
(display len-parts) (display " partials, ")
(display len-matches) (display " matches.") (newline))
(if (< (+ len-rels len-matches -10) len-fb) (loop (+ q 2) rels parts)
(let ((fact (solve n f fb (append rels matches))))
(if fact fact (loop (+ q 2) rels parts)))))))))))))
Here is an example:
> (qs 294729242679158229936006281 1500 100000)
Factor base of 122 primes.
q=15601: 11 smooths, 504 partials, 5 matches.
q=15607: 24 smooths, 1006 partials, 30 matches.
q=15619: 34 smooths, 1477 partials, 55 matches.
q=15629: 46 smooths, 1942 partials, 97 matches.
2971215073
You can clearly see the power of the multiple polynomial quadratic sieve as compared to the basic quadratic sieve. In the previous exercise, we used f = 2000; here, we reduce f to 1500. We also reduced m from 3000000 to 100000. Quite an improvement.
You can see the entire program assembled at the next page, or run the program at http://programmingpraxis.codepad.org/j3nstX1Y.
Programming Praxis 