[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/ohwDTFhz    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Mar 31:
; 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 (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) ; => (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)))))

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

(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)
  (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 (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 87463 30 32)) (newline)
(display (qs 13290059 150 300)) (newline)
(display (qs 294729242679158229936006281 1500 1000000)) (newline)


Output:
1
2
3
4
5
6
7
8
9
10
11
12
13
Factor base of 6 primes.
Found 8 smooth relations,
101 partial relations,
and 47 matches.
587
Factor base of 18 primes.
Found 125 smooth relations,
1695 partial relations,
and 1071 matches.
4261
Factor base of 122 primes.

Timeout


Create a new paste based on this one


Comments: