Programming Praxis

We begin with make-row, which returns the exponent vector for a single Q_n:

(define (make-row rs fb len)
  (let ((vec (make-vector (+ len 1) 0)))
    (vector-set! vec 0 (if (odd? (car rs)) 1 0))
    (let loop ((i 1) (rs (cdddr rs)) (fb fb))
      (cond ((or (null? rs) (null? fb)) vec)
            ((= (car fb) (car rs))
              (vector-set! vec i 1)
              (loop (+ i 1) (cdr rs) (cdr fb)))
            (else (loop (+ i 1) rs (cdr fb)))))))

Make-ev and make-hv build the exponent and history matrices, respectively:

(define (make-ev fb len residues)
  (let loop ((rss residues) (es '()) (as '()) (qs '()))
    (if (null? rss)
        (values (list->vector (reverse es))
                (list->vector (reverse as))
                (list->vector (reverse qs)))
        (loop (cdr rss)
              (cons (make-row (car rss) fb len) es)
              (cons (caddar rss) as)
              (cons (cadar rss) qs)))))

(define (make-hv ev-len)
  (let ((hv (make-vector ev-len)))
    (do ((i 0 (+ i 1))) ((= i ev-len) hv)
      (let ((h (make-vector ev-len 0)))
        (vector-set! h i 1)
        (vector-set! hv i h)))))

The calculation of the rightmost 1 is shown below. Later, in the answer function, we will keep track of the rightmost 1 in each row rather than recomputing it each time it is needed, just as Morrison and Brillhart did:

(define (right-most-one rv start)
  (let loop ((i start))
    (cond ((negative? i) i)
          ((= (vector-ref rv i) 1) i)
          (else (loop (- i 1))))))

Finally, we give the code for ANSWER, which initializes the exponent vector and history vector, creates a pointer vector for the rightmost 1 in each column, performs the reduction procedure, and then examines S-sets until it finds one that gives a factorization:

(define (answer n fb residues)
  (let-values (((ev av qv) (make-ev fb (length fb) residues)))
    (let* ((fb-len (length fb)) (ev-len (vector-length ev))
           (hv (make-hv ev-len)) (ptr (make-vector ev-len)))
      (do ((i 0 (+ i 1))) ((= i ev-len))
        (vector-set! ptr i (right-most-one (vector-ref ev i) fb-len)))
      ; reduce ev and hv
      (do ((j fb-len (- j 1))) ((negative? j))
        (let loop ((i 0))
          (when (< i ev-len)
            (if (= (vector-ref ptr i) j)
                (do ((m (+ i 1) (+ m 1))) ((= m ev-len))
                  (when (= (vector-ref ptr m) j)
                    (vector-set! ev m (vector-add (vector-ref ev i) (vector-ref ev m)))
                    (vector-set! hv m (vector-add (vector-ref hv i) (vector-ref hv m)))
                    (vector-set! ptr m (right-most-one (vector-ref ev m) j))))
                (loop (+ i 1))))))
      ; examine s-sets
      (let loop ((i 0))
        (cond ((= i ev-len) #f)
              ((zero? (vector-sum (vector-ref ev i)))
                (let ((d (gcd (a-minus-q n (vector-ref hv i) ev-len av qv) n)))
                  (if (< 1 d n) d (loop (+ i 1)))))
              (else (loop (+ i 1))))))))

Answer uses isqrt from the Standard Prelude, plus the utility functions shown below:

(define (vector-add v1 v2)
  (let* ((len (vector-length v1))
         (v3 (make-vector len)))
    (do ((i 0 (+ i 1))) ((= i len) v3)
      (vector-set! v3 i
        (if (= (vector-ref v1 i)
               (vector-ref v2 i)) 0 1)))))

(define (vector-sum vec)
  (do ((i 0 (+ i 1)) (sum 0 (+ sum (vector-ref vec i))))
      ((= i (vector-length vec)) sum)))

Answer also uses a-minus-q, which takes an S-set and returns a divisor of N; it is up to answer to determine if that factor is useful or trivial; note that the calculation of the square root of q-prod is exact:

(define (a-minus-q n vec ev-len av qv)
  (let loop ((i 0) (a-prod 1) (q-prod 1))
    (cond ((= i ev-len) (modulo (- a-prod (isqrt q-prod)) n))
          ((= (vector-ref vec i) 1)
            (loop (+ i 1)
                  (modulo (* a-prod (vector-ref av i)) n)
                  (* q-prod (vector-ref qv i))))
          (else (loop (+ i 1) a-prod q-prod)))))

Using the factor base and residues calculated in the previous exercise, the factorization of F₇ is shown below:

> (answer f7 fb residues)
59649589127497217
> (/ f7 59649589127497217)
5704689200685129054721

The S-set that completed the factorization contained 466 A – Q pairs, which are listed on the next page. You can see the entire program at http://programmingpraxis.codepad.org/JGMHnA0v.

Pages: 1 2 3