Programming Praxis

The previous version of the quadratic sieve used three functions: factor-base calculated the factors in the factor base, smooth determined if a number was smooth with respect to the factor base and returned its factors if it was, and qs was the main calling function that organized the computation. All three of those functions will change. In particular, qs will become more of a driver function with less in-line calculation, and we introduce two new functions: sieve performs the actual sieving and match combines partials into full relations. The supporting functions like jacobi and mod-sqrt are unchanged. The linear algebra changes slightly; make-expo-vector needs to ignore large prime factors, and the calculation of the square root is improved.

We begin with factor-base, combining the search for quadratic residues, calculation of the starting offsets, and the calculation of the approximate logarithm into a single loop. Quadratic residues are calculated as they were previously. Logarithms are simple if you keep your wits about you: any factor base prime between 2^x−1/2 and 2^x+1/2 will have an approximate logarithm of x. That means we can keep track of the next place where the logarithm changes, and recalculate the next change point each time we breach the last change point. The error factor e is set at double the maximum logarithm; this generally gives us a larger range of possible smooth numbers than the original version of the program, but still keeps us in a relatively small range of potential smooth numbers.

(define (factor-base n f) ; => (values fb ts ls e)
  (let loop ((ps (cdr (primes f))) (fb (list 2)) (ts (list)) (ls (list)) (x 2) (limit 5))
    (when (and (pair? ps) (< limit (car ps)))
      (set! x (+ x 1)) (set! limit (isqrt (expt 2 (+ x x 1)))))
    (cond ((null? ps) (values (reverse fb) (reverse ts) (reverse ls) (max 5 x)))
          ((= (jacobi n (car ps)) 1)
            (loop (cdr ps) (cons (car ps) fb)
                  (cons (msqrt n (car ps)) ts) (cons x ls) x limit))
          (else (loop (cdr ps) fb ts ls x limit)))))

Next is smooth, which has to change to recognize large primes. Fortunately, the change is small; if we ever get to the point where f × f > n, then n is prime and we return the smooth relation, whether or not n > f.

(define (smooth n fb) ; list of smooth factors of n in descending order, or null
  (let loop ((n (abs n)) (fb fb) (fs (if (negative? n) (list -1) (list))))
    (cond ((null? fb) (list))
          ((< n (* (car fb) (car fb))) (cons n fs))
          ((zero? (modulo n (car fb)))
            (loop (/ n (car fb)) fb (cons (car fb) fs)))
          (else (loop n (cdr fb) fs)))))

The new function sieve implements the segmented sieve; we extract that code to a separate function because it is rather more complicated than we want to have in the driver. Function sieve takes as inputs the fb, ts, ls and e computed by factor-base, as well as n and m, and returns two values: a list of relations rels and a list of partial relations parts.

(define (sieve n fb ts ls e f m) ; => (values rels parts)
  ; a relation, whether full or partial, has x in its car and a list
  ; of smooth factors of g(x)=x^2-n, in descending order, in its cdr;
  ; a large prime, if it exists, is at the cadr of the relation
  (let* ((sqrt-n (isqrt n)) (b (max (- sqrt-n m) 2)) (rels (list)) (parts (list))
         (lo b) (hi (+ b m m)) (delta 100000) (sieve (make-vector delta)))
    (do ((lo lo (+ lo delta))) ((< hi lo))
      (do ((i 0 (+ i 1))) ((= delta i)) (vector-set! sieve i (+ e e)))
      (do ((fb (cdr fb) (cdr fb)) (ts ts (cdr ts)) (ls ls (cdr ls))) ((null? fb))
        (do ((i (modulo (- (car ts) lo) (car fb)) (+ i (car fb)))) ((<= delta i))
          (vector-set! sieve i (+ (vector-ref sieve i) (car ls))))
        (do ((i (modulo (- (- (car ts)) lo) (car fb)) (+ i (car fb)))) ((<= delta i))
          (vector-set! sieve i (+ (vector-ref sieve i) (car ls)))))
      (do ((i 0 (+ i 1))) ((= delta i))
        (let* ((x (+ i lo)) (g (- (* x x) n)))
          (if (< (ilog 2 g) (vector-ref sieve i))
            (let* ((ys (smooth g fb)) (rel (cons x ys)))
              (if (pair? ys) (if (<= (car ys) f)
                                 (set! rels (cons rel rels))
                                 (set! parts (cons rel parts)))))))))

The other new function match takes the list of partial relations and returns as many full relations as it can find; it first sorts on the large primes, then scans, writing a full relation (including the two copies of the large prime) whenever it finds a large prime the same as its predecessor.

(define (match parts)
  (define (lt? a b) (< (cadr a) (cadr b)))
  (let loop ((parts (sort lt? parts)) (prev (list 0 0)) (zs (list)))
    (cond ((null? parts) zs)
          ((= (cadar parts) (cadr prev))
            (loop (cdr parts) prev
                  (cons (cons (* (caar parts) (car prev))
                              (merge > (cdar parts) (cdr prev))) zs)))
          (else (loop (cdr parts) (car parts) zs)))))

The last changed function is qs, which connects the rest of the functions as needed; it first computes the factor base, then sieves, computes matches, and finally solves.

(define (qs n f m)
  (call-with-values
    (lambda () (factor-base n f))
    (lambda (fb ts ls e)
      (when verbose? (display "Factor base of ")
        (display (length fb)) (display " primes.") (newline))
      (call-with-values
        (lambda () (sieve n fb ts ls e f m))
        (lambda (rels parts)
          (let ((matches (match parts)))
            (when verbose? (display "Found ") (display (length rels))
              (display " smooth relations,") (newline)
              (display (length parts)) (display " partial relations,")
              (newline) (display "and ") (display (length matches))
              (display " matches.") (newline))
            (solve n f fb (append rels matches))))))))

Here are some examples.

> (qs 87463 30 30)
Factor base of 6 primes.
Found 8 smooth relations,
101 partial relations,
and 47 matches.
587
> (qs 13290059 150 300)
Factor base of 18 primes.
Found 125 smooth relations,
1695 partial relations,
and 1071 matches.
4261
> (qs 294729242679158229936006281 2000 3000000)
Factor base of 149 primes.
Found 160 smooth relations,
6118 partial relations,
and 695 matches.
2971215073

We used the same examples as last time, including the same f and m. If you compare, you will see that with the same effort we found more relations; there are more smooth relations because of the change in calculating the error adjustment, and the matches are completely new. By reducing f or m or both, we could have found sufficient relations with less effort. And the effect grows as n gets larger. The improvements are well worth the effort we made and the additional complexity they bring to the code.

You can run the complete program at http://programmingpraxis.codepad.org/ohwDTFhz. We’re not finished with the quadratic sieve just yet, so keep watching for additional improvements in the future.

Pages: 1 2