### June 21, 2013

`; multiple polynomial quadratic sieve`

`(define verbose? #t)`

```(define (primes n)     (let ((bits (make-vector (+ n 1) #t)))       (let loop ((p 2) (ps '()))         (cond ((< n p) (reverse ps))               ((vector-ref bits p)                 (do ((i (* p p) (+ i p))) ((< n i))                   (vector-set! bits i #f))                 (loop (+ p 1) (cons p ps)))               (else (loop (+ p 1) ps))))))```

```(define prime?   (let ((seed 3141592654))     (lambda (n)       (define (rand)         (set! seed (modulo (+ (* 69069 seed) 1234567) 4294967296))         (+ (quotient (* seed (- n 2)) 4294967296) 2))       (define (expm b e m)         (define (times x y) (modulo (* x y) m))         (let loop ((b b) (e e) (r 1))           (if (zero? e) r             (loop (times b b) (quotient e 2)                   (if (odd? e) (times b r) r)))))       (define (spsp? n a)         (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))             ((odd? d)               (let ((t (expm a d n)))                 (if (or (= t 1) (= t (- n 1))) #t                   (do ((s (- s 1) (- s 1))                        (t (expm t 2 n) (expm t 2 n)))                       ((or (zero? s) (= t (- n 1)))                         (positive? s))))))))       (if (not (integer? n))           (error 'prime? "must be integer")           (if (< n 2) #f             (do ((a (rand) (rand)) (k 25 (- k 1)))                 ((or (zero? k) (not (spsp? n a)))                   (zero? k))))))))```

```(define (isqrt n)   (if (not (and (positive? n) (integer? n)))       (error 'isqrt "must be positive integer")       (let loop ((x n))         (let ((y (quotient (+ x (quotient n x)) 2)))           (if (< y x) (loop y) x)))))```

```(define (ilog b n)   (let loop1 ((lo 0) (b^lo 1) (hi 1) (b^hi b))     (if (< b^hi n) (loop1 hi b^hi (* hi 2) (* b^hi b^hi))       (let loop2 ((lo lo) (b^lo b^lo) (hi hi) (b^hi b^hi))         (if (<= (- hi lo) 1) (if (= b^hi n) hi lo)           (let* ((mid (quotient (+ lo hi) 2))                  (b^mid (* b^lo (expt b (- mid lo)))))             (cond ((< n b^mid) (loop2 lo b^lo mid b^mid))                   ((< b^mid n) (loop2 mid b^mid hi b^hi))                   (else mid))))))))```

```(define (expm b e m)   (define (m* x y) (modulo (* x y) m))   (cond ((zero? e) 1)         ((even? e) (expm (m* b b) (/ e 2) m))         (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))```

```(define (mod-sqrt a p)   (define (both n) (list n (- p n)))   (cond ((not (and (odd? p) (prime? p)))           (error 'mod-sqrt "modulus must be an odd prime"))         ((not (= (jacobi a p) 1))           (error 'mod-sqrt "must be a quadratic residual"))         (else (let ((a (modulo a p)))                 (case (modulo p 8)                   ((3 7) (both (expm a (/ (+ p 1) 4) p)))                   ((5) (let* ((x (expm a (/ (+ p 3) 8) p))                               (c (expm x 2 p)))                          (if (= a c) (both x)                            (both (modulo (* x (expm 2 (/ (- p 1) 4) p)) p)))))                   (else (let* ((d (let loop ((d 2))                                     (if (= (jacobi d p) -1) d                                       (loop (+ d 1)))))                                (s (let loop ((p (- p 1)) (s 0))                                     (if (odd? p) s                                       (loop (quotient p 2) (+ s 1)))))                                (t (quotient (- p 1) (expt 2 s)))                                (big-a (expm a t p))                                (big-d (expm d t p))                                (m (let loop ((i 0) (m 0))                                     (cond ((= i s) m)                                           ((= (- p 1)                                               (expm (* big-a (expm big-d m p))                                                     (expt 2 (- s 1 i)) p))                                             (loop (+ i 1) (+ m (expt 2 i))))                                           (else (loop (+ i 1) m))))))                           (both (modulo (* (expm a (/ (+ t 1) 2) p)                                            (expm big-d (/ m 2) p)) p)))))))))```

```(define (lift-root n p) ; hensel's lemma   (let* ((r (apply min (mod-sqrt n p)))          (s (/ (- (modulo n (* p p)) (* r r)) p))          (t (modulo (* (inverse (+ r r) p) s) p))          (u (+ r (* t p))))     (list u (- (* p p) u))))```

```(define (jacobi a m)   (if (not (integer? a)) (error 'jacobi "must be integer")     (if (not (and (integer? m) (positive? m) (odd? m)))         (error 'jacobi "modulus must be odd positive integer")         (let loop1 ((a (modulo a m)) (m m) (t 1))           (if (zero? a) (if (= m 1) t 0)             (let ((z (if (member (modulo m 8) (list 3 5)) -1 1)))               (let loop2 ((a a) (t t))                 (if (even? a) (loop2 (/ a 2) (* t z))                   (loop1 (modulo m a) a                          (if (and (= (modulo a 4) 3)                                   (= (modulo m 4) 3))                              (- t) t))))))))))```

```(define (inverse x m)   (let loop ((x x) (a 0) (b m) (u 1))     (if (positive? x)         (let ((q (quotient b x)) (r (remainder b x)))           (loop (modulo b x) u x (- a (* q u))))         (if (= b 1) (modulo a m) (error 'inverse "must be coprime"))))) ```

```(define (factor-base n f) ; => (values fb ts ls e fb-len)   (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) (length fb)))           ((= (jacobi n (car ps)) 1)             (loop (cdr ps) (cons (car ps) fb)                   (cons (apply min (mod-sqrt n (car ps))) ts) (cons x ls) x limit))           (else (loop (cdr ps) fb ts ls x limit)))))```

```(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)))))```

```(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)))```

```(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)))))```

```(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)))))))))))))```

```(define (make-expo-vector f fb rel)   (define (add-1bit x) (if (zero? x) 1 0))   (let loop ((fb fb) (rel rel) (prev -2) (es (list)))     (cond ((and (null? fb) (null? rel)) (list->vector (reverse es)))           ((null? fb) (loop fb (cdr rel) prev (cons (add-1bit (car es)) (cdr es))))           ((null? rel) (loop (cdr fb) rel prev (cons 0 es)))           ((< f (car rel)) (loop fb (cdr rel) prev es))           ((= (car rel) prev) (loop fb (cdr rel) prev (cons (add-1bit (car es)) (cdr es))))           ((= (car rel) (car fb)) (loop (cdr fb) (cdr rel) (car rel) (cons 1 es)))           (else (loop (cdr fb) rel prev (cons 0 es))))))```

```(define (make-identity-matrix n)   (let ((id (make-vector n)))     (do ((i 0 (+ i 1))) ((= i n) id)       (let ((v (make-vector n 0)))         (vector-set! v i 1)         (vector-set! id i v)))))```

```(define (left-most-one vec r) ; column of left-most 1, or -1 if all zero   (let* ((row (vector-ref vec r)) (len (vector-length row)))     (let loop ((i 0))       (cond ((= i len) -1)             ((= (vector-ref row i) 1) i)             (else (loop (+ i 1)))))))```

```(define (pivot-row expo c)   (let ((max-r (vector-length expo)))     (let loop ((r 0))       (if (= r max-r) r         (if (= (left-most-one expo r) c) r           (loop (+ r 1)))))))```

```(define (add-rows matrix r1 r2)   (define (add a b) (if (= a b) 0 1))   (let ((row1 (vector-ref matrix r1)) (row2 (vector-ref matrix r2)))     (do ((i 0 (+ i 1))) ((= i (vector-length row1)) row2)       (vector-set! row2 i (add (vector-ref row1 i) (vector-ref row2 i))))))```

```(define (any-one? vec r)   (let* ((row (vector-ref vec r)) (r-len (vector-length row)))     (let loop ((i 0))       (if (= i r-len) #f         (if (positive? (vector-ref row i)) #t           (loop (+ i 1)))))))```

```(define (factor n hist rels r)   (define (root y2)     (let loop ((y2 (sort < y2)) (s 1))       (if (null? y2) s         (loop (cddr y2) (* s (car y2))))))   (let* ((h (vector-ref hist r)) (h-len (vector-length h)))     (let loop ((i 0) (x 1) (y2 (list)))       (cond ((= i h-len)               (let ((g (gcd (- x (root y2)) n)))                 (if (< 1 g n) g #f)))             ((= (vector-ref h i) 1)               (loop (+ i 1) (* x (car (vector-ref rels i)))                     (append (cdr (vector-ref rels i)) y2)))             (else (loop (+ i 1) x y2))))))```

```(define (solve n f fb rels)   (let* ((fb (reverse (cons -1 fb))) (fb-len (length fb)) (rel-len (length rels))          (expo (list->vector (map (lambda (rel) (make-expo-vector f fb (cdr rel))) rels)))          (hist (make-identity-matrix rel-len))          (rels (list->vector rels)))     (do ((c 0 (+ c 1))) ((= c fb-len))       (let ((p (pivot-row expo c)))         (do ((r (+ p 1) (+ r 1))) ((<= rel-len r))           (when (= (left-most-one expo r) c)             (vector-set! expo r (add-rows expo p r))             (vector-set! hist r (add-rows hist p r))))))     (let loop ((r 0))       (cond ((= r rel-len) #f)             ((any-one? expo r) (loop (+ r 1)))             ((factor n hist rels r) =>               (lambda (f) (if f f (loop (+ r 1)))))             (else (loop (+ r 1)))))))```