codepad

Language:

; quadratic sieve
;     f: maximum prime in factor base
;     2*m: number of entries in sieve
;     relations with x in car and ys in cdr
; examples
;   (qs 87463 30 30) => 587
;   (qs 13290059 150 300) => 4261
;   (qs (* 2971215073 99194853094755497) 2000 3000000) => 99194853094755497

(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 (square x) (* x x))

(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 (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 (msqrt a p) ; principal (smaller) value of modular square root
  (apply min (mod-sqrt a p)))

(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 (factor-base n f)
  (let loop ((ps (cdr (primes f))) (fs (list 2)))
    (cond ((null? ps) (reverse fs))
          ((= (jacobi n (car ps)) 1)
            (loop (cdr ps) (cons (car ps) fs)))
          (else (loop (cdr ps) fs)))))

(define (smooth n fb) ; list of factors with -1, or null if not smooth
  (let ((sign (if (negative? n) -1 1)) (n (abs n)))
    (let loop ((n n) (fb fb) (fs (list)))
      (cond ((null? fb) (list))
            ((= n 1) (if (negative? sign) (cons -1 (reverse fs)) (reverse fs)))
            ((zero? (modulo n (car fb)))
              (loop (/ n (car fb)) fb (cons (car fb) fs)))
            (else (loop n (cdr fb) fs))))))

(define (qs n f m)
  (let* ((e 10) ; fudge factor 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))))))

(define (solve n fb rels)
  (display rels) (newline))

(define (make-expo-vector fb rel)
  (define (add-1bit x) (if (zero? x) 1 0))
  (let loop ((fb fb) (rel rel) (prev -2) (es (list)))
    (cond ((null? fb) (list->vector (reverse es)))
          ((null? rel) (loop (cdr fb) rel prev (cons 0 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 fb) (cons 1 es)))
          (else (loop (cdr fb) rel (car fb) (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 (right-most-one vec r)
  (let ((row (vector-ref vec r)))
    (let loop ((i (- (vector-length row) 1)))
      (if (negative? i) i
        (if (= (vector-ref row i) 1) i
          (loop (- i 1)))))))

(define (pivot-row expo c)
  (let ((max-r (vector-length expo)))
    (let loop ((r 0))
      (if (= r max-r) r
        (if (= (right-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)
  (let* ((h (vector-ref hist r))
         (h-len (vector-length h)))
    (let loop ((i 0) (x 1) (y2 1))
      (if (= i h-len)
          (let ((g (gcd (- x (isqrt y2)) n)))
            (if (< 1 g n) g #f))
          (if (= (vector-ref h i) 1)
              (loop (+ i 1)
                    (* x (car (vector-ref rels i)))
                    (apply * y2 (cdr (vector-ref rels i))))
              (loop (+ i 1) x y2))))))

(define (solve n fb rels)
  (let* ((fb (cons -1 fb)) (fb-len (length fb)) (rel-len (length rels))
         (expo (list->vector (map (lambda (rel) (make-expo-vector fb (cdr rel))) rels)))
         (hist (make-identity-matrix rel-len))
         (rels (list->vector rels)))
    (do ((c (- fb-len 1) (- c 1))) ((negative? c))
      (let ((p (pivot-row expo c)))
        (do ((r (+ p 1) (+ r 1))) ((<= rel-len r))
          (when (= (right-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)))))))

(display (qs 13290059 150 300))

Private [?]

Run code