[ create a new paste ] login | about

Link: http://codepad.org/j3nstX1Y    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Jun 20:
; multiple polynomial quadratic sieve

(define verbose? #t)

(define sort #f)
(define merge #f)
(let ()
  (define dosort
    (lambda (pred? ls n)
      (if (= n 1)
          (list (car ls))
          (let ((i (quotient n 2)))
            (domerge pred?
                     (dosort pred? ls i)
                     (dosort pred? (list-tail ls i) (- n i)))))))
  (define domerge
    (lambda (pred? l1 l2)
      (cond
        ((null? l1) l2)
        ((null? l2) l1)
        ((pred? (car l2) (car l1))
         (cons (car l2) (domerge pred? l1 (cdr l2))))
        (else (cons (car l1) (domerge pred? (cdr l1) l2))))))
  (set! sort
    (lambda (pred? l)
      (if (null? l) l (dosort pred? l (length l)))))
  (set! merge
    (lambda (pred? l1 l2)
      (domerge pred? l1 l2))))

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

(display (qs 5624788913579353957 800 4000))


Output:
1
2
3
4
5
6
Factor base of 76 primes.
q=919: 16 smooths, 243 partials, 0 matches.
q=929: 30 smooths, 475 partials, 18 matches.
q=937: 41 smooths, 756 partials, 40 matches.
q=953: 57 smooths, 995 partials, 72 matches.
2572067131


Create a new paste based on this one


Comments: