Improved Standard Continuation

November 8, 2011

(define (primes . args) ; (primes [lo] hi) inclusive at both ends
  (let* ((lo (if (null? (cdr args)) 0 (car args)))
         (hi (if (null? (cdr args)) (car args) (cadr args))))
    (cond ((and (<= lo 100000) (<= hi 1000000)) ; simple sieve
           (let* ((max-index (quotient (- hi 3) 2))
                  (v (make-vector (+ max-index 1) #t)))
             (let loop ((i 0) (ps (list 2)))
               (let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
                 (cond ((< hi (* p p))
                        (let loop ((j i) (ps ps))
                          (cond ((< max-index j)
                                 (let loop ((ps (reverse ps)))
                                   (if (<= lo (car ps)) ps
                                     (loop (cdr ps)))))
                                ((vector-ref v j)
                                 (loop (+ j 1) (cons (+ j j 3) ps)))
                                (else (loop (+ j 1) ps)))))
                       ((vector-ref v i)
                        (let loop ((j startj))
                          (when (<= j max-index)
                            (vector-set! v j #f) (loop (+ j p))))
                        (loop (+ i 1) (cons p ps)))
                       (else (loop (+ i 1) ps)))))))
          ((< lo (sqrt hi))
           (let* ((r (inexact->exact (ceiling (sqrt hi))))
                  (r (if (even? r) r (+ r 1))))
             (append (primes lo r) (primes r hi))))
          (else ; segmented sieve
           (let* ((lo (if (even? lo) lo (- lo 1)))
                  (b 50000) (bs (make-vector b #t))
                  (r (inexact->exact (ceiling (sqrt hi))))
                  (ps (cdr (primes r)))
                  (qs (map (lambda (p)
                             (modulo (* -1/2 (+ lo 1 p)) p)) ps))
                  (zs (list)) (z (lambda (p) (set! zs (cons p zs)))))
             (do ((t lo (+ t b b))
                  (qs qs (map (lambda (p q) (modulo (- q b) p))
                              ps qs)))
                 ((<= hi t)
                   (let loop ((zs zs))
                     (if (<= (car zs) hi) (reverse zs)
                       (loop (cdr zs)))))
               (do ((i 0 (+ i 1))) ((= i b)) (vector-set! bs i #t))
               (do ((ps ps (cdr ps)) (qs qs (cdr qs))) ((null? qs))
                 (do ((j (car qs) (+ j (car ps)))) ((<= b j))
                   (vector-set! bs j #f)))
               (do ((j 0 (+ j 1))) ((= j b))
                 (if (vector-ref bs j) (z (+ t j j 1))))))))))

(define (prime? n) ; spsp(2), spsp(3), and lucas pseudoprime
  (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 (digits n . args)
    (let ((b (if (null? args) 10 (car args))))
      (let loop ((n n) (d '()))
        (if (zero? n) d
            (loop (quotient n b)
                  (cons (modulo n b) d))))))
  (define (isqrt n)
    (let loop ((x n) (y (quotient (+ n 1) 2)))
      (if (<= 0 (- y x) 1) x
        (loop y (quotient (+ y (quotient n y)) 2)))))
  (define (square? n)
    (let ((n2 (isqrt n)))
      (= n (* n2 n2))))
  (define (jacobi a n)
    (if (not (and (integer? a) (integer? n) (positive? n) (odd? n)))
        (error 'jacobi "modulus must be positive odd integer")
        (let jacobi ((a a) (n n))
          (cond ((= a 0) 0)
                ((= a 1) 1)
                ((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1)))
                ((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n)))
                ((< n a) (jacobi (modulo a n) n))
                ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a)))
                (else (jacobi n a))))))
  (define legendre jacobi)
  (define (inverse x n)
    (let loop ((x (modulo x n)) (a 1))
      (cond ((zero? x) (error 'inverse "division by zero"))
            ((= x 1) a)
            (else (let ((q (- (quotient n x))))
                    (loop (+ n (* q x)) (modulo (* q a) n)))))))
  (define (miller? n a)
    (let loop ((r 0) (s (- n 1)))
      (if (even? s) (loop (+ r 1) (/ s 2))
        (if (= (expm a s n) 1) #t
          (let loop ((r r) (s s))
            (cond ((zero? r) #f)
                  ((= (expm a s n) (- n 1)) #t)
                  (else (loop (- r 1) (* s 2)))))))))
  (define (chain m f g x0 x1)
    (let loop ((ms (digits m 2)) (u x0) (v x1))
      (cond ((null? ms) (values u v))
            ((zero? (car ms)) (loop (cdr ms) (f u) (g u v)))
            (else (loop (cdr ms) (g u v) (f v))))))
  (define (lucas? n)
    (let loop ((a 11) (b 7))
      (let ((d (- (* a a) (* 4 b))))
        (cond ((square? d) (loop (+ a 2) (+ b 1)))
              ((not (= (gcd n (* 2 a b d)) 1))
                (loop (+ a 2) (+ b 2)))
              (else (let* ((x1 (modulo (- (* a a (inverse b n)) 2) n))
                           (m (quotient (- n (legendre d n)) 2))
                           (f (lambda (u) (modulo (- (* u u) 2) n)))
                           (g (lambda (u v) (modulo (- (* u v) x1) n))))
                      (let-values (((xm xm1) (chain m f g 2 x1)))
                        (zero? (modulo (- (* x1 xm) (* 2 xm1)) n)))))))))
  (cond ((or (not (integer? n)) (< n 2))
          (error 'prime? "must be integer greater than one"))
        ((even? n) (= n 2)) ((zero? (modulo n 3)) (= n 3))
        (else (and (miller? n 2) (miller? n 3) (lucas? n)))))

(define verbose? #f)

(define (factors n)
  (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 (times p q) (modulo (* p q) m))
    (let loop ((b b) (e e) (x 1))
      (if (zero? e) x
        (loop (times b b) (quotient e 2)
              (if (odd? e) (times b x) x)))))
  (define (digits n . args)
    (let ((b (if (null? args) 10 (car args))))
      (let loop ((n n) (d '()))
        (if (zero? n) d
            (loop (quotient n b)
                  (cons (modulo n b) d))))))
  (define (cycle . xs)
    (define (last-pair xs)
      (if (null? (cdr xs)) xs (last-pair (cdr xs))))
    (set-cdr! (last-pair xs) xs) xs)
  (define (msg . xs)
    (when verbose?
      (for-each display xs)
      (newline)))
  (define (seed) ; (date-and-time) => "Mon Nov 07 10:56:23 2011"
    (let ((days '("Sun" "Mon" "Tue" "Wed" "Thu" "Fri" "Sat"))
          (months '("Jan" "Feb" "Mar" "Apr" "May" "Jun"
                    "Jul" "Aug" "Sep" "Oct" "Nov" "Dec"))
          (str (date-and-time))) ; NON-PORTABLE (Chez)
      (* (length (member (substring str 0 3) days))
         (length (member (substring str 4 7) months))
         (string->number (substring str 8 10)) ; date
         (+ (string->number (substring str 11 13)) 1) ; hour
         (+ (string->number (substring str 14 16)) 1) ; minute
         (+ (string->number (substring str 17 19)) 1) ; second
         (string->number (substring str 20 24))))) ; year
  (define rand
    (let* ((a 69069) (c 1234567) (m 4294967296)
           (s (modulo (+ (* a (seed)) c) m)))
      (lambda () (set! s (modulo (+ (* a s) c) m)) s)))
  (define td-factors
    (let ((saved-b 0) (ps '()))
      (lambda (n b)
        (when (not (= b saved-b))
          (set! saved-b b) (set! ps (primes b)))
        (let loop ((n n) (ps ps) (fs '()))
          (cond ((null? ps) (values (reverse fs) n))
                ((= n 1) (values (reverse fs) n))
                ((zero? (modulo n (car ps)))
                  (loop (/ n (car ps)) ps (cons (car ps) fs)))
                (else (loop n (cdr ps) fs)))))))
  (define (rho-factor n c b)
    (define (f y) (modulo (+ (* y y) c) n))
    (define (g p x y) (modulo (* p (abs (- x y))) n))
    (let stage1 ((x 2) (y (+ 4 c)) (z (+ 4 c)) (j 1) (q 2) (p 1))
      (if (= j b) #f ; timeout
        (if (= x y) (rho-factor n (+ c 1) (- b j)) ; cycle
          (if (= j q) (let ((t (f y))) (stage1 y (f y) z (+ j 1) (* q 2) (g p y t)))
            (if (positive? (modulo j 100)) (stage1 x (f y) z (+ j 1) q (g p x y))
              (let ((d (gcd p n)))
                (if (= d 1) (stage1 x (f y) y (+ j 1) q (g p x y))
                  (if (< 1 d n) d ; factor
                    (let stage2 ((k 1) (z (f z)))
                      (if (= k 100) (rho-factor n (+ c 1) (- b j))
                        (let ((d (gcd (- z x) n)))
                          (if (= d 1) (stage2 (+ k 1) (f z))
                            (if (= d n) (rho-factor n (+ c 1) (- b j))
                              d))))))))))))))
  (define pminus1-factor
    (let ((b1-saved 0) (b2-saved 0) (ps-saved '()) (as-saved '())
          (ds-saved '()) (len 0) (dv-saved '()))
      (lambda (n b1 b2)
        (define (cycle . xs)
          (define (last-pair xs)
            (if (null? (cdr xs)) xs
              (last-pair (cdr xs))))
          (set-cdr! (last-pair xs) xs) xs)
        (when (not (and (= b1 b1-saved) (= b2 b2-saved)))
          (set! b1-saved b1) (set! b2-saved b2) (set! ps (primes b2))
          (set! as (let loop ((ps ps) (as (list)))
                     (let ((a (ilog (car ps) b1)))
                       (if (= a 1) (append (reverse as) (cycle 1))
                         (loop (cdr ps) (cons a as))))))
          (set! ds (map - (cdr ps) (reverse (cdr (reverse ps)))))
          (set! len (apply max ds)) (set! dv (make-vector (+ len 1) 0)))
        (let stage1 ((ps ps) (as as) (ds ds) (k 1) (q 2) (p 1))
          (cond ((< 100 k) (let ((g (gcd (- q 1) n)))
                  (if (< 1 g n) g (stage1 ps as ds 1 q p))))
                ((< (car ps) b1) (let ((q (expm q (expt (car ps) (car as)) n)))
                  (stage1 (cdr ps) (cdr as) (cdr ds) (+ k 1) q (modulo (* p (- q 1)) n))))
                (else (let ((g (gcd (- q 1) n))) (if (< 1 g n) g (begin
                  (do ((i 1 (+ i 1))) ((< len i)) (vector-set! dv i (expm q (+ i i) n)))
                  (let stage2 ((ds ds) (k 1) (q (expm q (car ps) n)) (p 1))
                    (cond ((null? ds) (let ((g (gcd p n))) (if (< 1 g n) g #f)))
                          ((< 100 k) (let ((g (gcd p n))) (if (< 1 g n) g (stage2 ds 1 q p))))
                          (else (let ((q (modulo (* q (vector-ref dv (car ds))) n)))
                            (stage2 (cdr ds) (+ k 1) q (modulo (* p (- q 1)) n)))))))))))))))
  (define (add p1 p2 p1-p2 n)
    (define (square x) (* x x))
    (let* ((x0 (car p1-p2)) (x1 (car p1)) (x2 (car p2))
           (z0 (cdr p1-p2)) (z1 (cdr p1)) (z2 (cdr p2))
           (t1 (modulo (* (+ x1 z1) (- x2 z2)) n))
           (t2 (modulo (* (- x1 z1) (+ x2 z2)) n)))
      (cons (modulo (* (square (+ t2 t1)) z0) n)
            (modulo (* (square (- t2 t1)) x0) n))))
  (define (double p an ad n)
    (define (square x) (* x x))
    (let* ((x (car p)) (z (cdr p))
           (x+z2 (modulo (square (+ x z)) n))
           (x-z2 (modulo (square (- x z)) n))
           (t (- x+z2 x-z2)))
      (cons (modulo (* x+z2 x-z2 4 ad) n)
            (modulo (* (+ (* t an) (* x-z2 ad 4)) t) n))))
  (define (multiply k p an ad n)
    (cond ((zero? k) (cons 0 0)) ((= k 1) p) ((= k 2) (double p an ad n))
      (else (let loop ((ks (cdr (digits k 2))) (q (double p an ad n)) (r p))
              (cond ((null? ks) r)
                    ((odd? (car ks))
                      (loop (cdr ks) (double q an ad n) (add q r p n)))
                    (else (loop (cdr ks) (add r q p n) (double r an ad n))))))))
  (define (curve12 n s)
    (let* ((u (modulo (- (* s s) 5) n))
           (v (modulo (* 4 s) n)) (v-u (- v u)))
      (values (modulo (* (* v-u v-u v-u) (+ u u u v)) n)
              (modulo (* 4 u u u v) n)
              (cons (modulo (* u u u) n)
                    (modulo (* v v v) n)))))
  (define ec-factor ; Crandall/Pomerance Algorithm 7.4.4
    (let* ((b1-saved 0) (b2-saved 0) (ps '()) (as '()) (ds '()) (d 210)
           (sv (make-vector (+ d 1) #f)) (bv (make-vector (+ d 1) #f)))
      (lambda (n b1 b2 s)
        (when (not (and (= b1 b1-saved) (= b2 b2-saved)))
          (set! b1-saved b1) (set! b2-saved b2) (set! ps (primes b2))
          (set! as (let loop ((ps ps) (as (list)))
                     (let ((a (ilog (car ps) b1)))
                       (if (= a 1) (append (reverse as) (cycle 1))
                         (loop (cdr ps) (cons a as))))))
          (set! ds (cons 0 (map (lambda (x) (/ x 2))
                     (map - (cdr ps) (reverse (cdr (reverse ps))))))))
        (let-values (((an ad q) (curve12 n s)))
          (let stage1 ((ps ps) (as as) (ds ds) (q q))
            (if (< (car ps) b1)
                (stage1 (cdr ps) (cdr as) (cdr ds)
                        (multiply (expt (car ps) (car as)) q an ad n))
                (let ((g (gcd (cdr q) n))) (if (< 1 g n) (list 1 g)
                  (let* ((v (- b1 1)) (r (multiply v q an ad n))
                         (t (multiply (- b1 1 d d) q an ad n))
                         (alpha (modulo (* (car r) (cdr r)) n)))
                    (vector-set! sv 1 (double q an ad n))
                    (vector-set! sv 2 (double (vector-ref sv 1) an ad n))
                    (do ((i 3 (+ i 1))) ((< d i))
                      (vector-set! sv i
                        (add (vector-ref sv (- i 1))
                             (vector-ref sv 1)
                             (vector-ref sv (- i 2)) n)))
                    (do ((i 1 (+ i 1))) ((< d i))
                      (vector-set! bv i
                        (modulo (* (car (vector-ref sv i))
                                   (cdr (vector-ref sv i))) n)))
                    (let stage2 ((v v) (ps ps) (ds ds) (r r) (t t) (g 1))
                      (if (null? ps) (let ((g (gcd g n))) (if (< 1 g n) (list 2 g) #f))
                        (let ((v2d (+ v d d)) (alpha (modulo (* (car r) (cdr r)) n)))
                          (let loop ((ps ps) (ds ds) (g g))
                            (if (and (pair? ps) (< (car ps) v2d))
                                (loop (cdr ps) (cdr ds) (modulo (* g (+ (*
                                      (- (car r) (car (vector-ref sv (car ds))))
                                      (+ (cdr r) (cdr (vector-ref sv (car ds)))))
                                      (- alpha) (vector-ref bv (car ds)))) n))
                                (let* ((old-r r) (r (add r (vector-ref sv d) t n)))
                                  (let ((g (gcd g n))) (if (< 1 g n) (list 2 g)
                                    (stage2 v2d ps ds r old-r g))))))))))))))))))
  (define td-limit 1000)
  (define rho-limit 500000)
  (define rho-constant 1)
  (define pminus1-b1 100000)
  (define pminus1-b2 5000000)
  (define ecf-b1 '(50000 200000 1000000 5000000))
  (define ecf-trials '(100 170 200 210))
  (define ecf-b2/b1 50)
  (let ((n n) (facts '()))
    (call-with-current-continuation (lambda (exit)
      (define (factor? method f)
        (if (not f) #f
          (let ((f (if (eq? method 'ecf) (cadr f) f))
                (stage (if (eq? method 'ecf) (car f) #f)))
            (if (prime? f)
                (begin (set! n (/ n f)) (set! facts (cons f facts))
                       (if (eq? method 'ecf)
                           (msg " In stage " stage ", found factor "
                                f ", remaining co-factor " n)
                           (msg " Found factor " f ", remaining co-factor " n))
                       (when (prime? n)
                         (msg " Factorization complete")
                         (exit (sort < (cons n facts))))
                       #t)
                (begin (if (eq? method 'ecf)
                           (msg " In stage " stage ", found non-prime factor " f)
                           (msg " Found non-prime factor " f))
                       (let ((fs (factors f)))
                         (if (or (not fs) (pair? (car fs)))
                             (exit (cons (sort < facts) (* n f)))
                             (begin (set! n (/ n f))
                                    (set! facts (append fs facts))
                                    (when (prime? n)
                                      (msg " Factorization complete")
                                      (exit (sort < (cons n facts))))
                                    #t))))))))
      (when (not (integer? n)) (error 'factors "must be integer"))
      (when (member n '(-1 0 1)) (exit (list n)))
      (when (negative? n) (exit (cons -1 (factors (- n)))))
      (when (prime? n) (msg "Input number is prime") (exit (list n)))
      (msg "Trial division: bound=" td-limit)
      (let-values (((fs cofact) (td-factors n td-limit)))
        (when (pair? fs) (msg " Found factors " fs)
          (when (< 1 cofact) (msg " Remaining co-factor " cofact)))
        (when (= cofact 1) (msg " Factorization complete") (exit fs))
        (when (prime? cofact)
          (msg " Remaining cofactor is prime")
          (msg " Factorization complete")
          (exit (append fs (list cofact))))
        (set! facts (append fs facts)) (set! n cofact))
      (let loop ((c rho-constant))
        (msg "Pollard rho: bound=" rho-limit ", constant=" c)
        (when (factor? 'rho (rho-factor n c rho-limit))
          (loop (+ c 1))))
      (let loop ()
        (msg "Pollard p-1: b1=" pminus1-b1 ", b2=" pminus1-b2)
        (when (factor? 'pm1 (pminus1-factor n pminus1-b1 pminus1-b2)) (loop)))
      (let loop ((c 1) (b1s ecf-b1) (trials ecf-trials))
        (when (pair? b1s)
          (let* ((b1 (car b1s)) (b2 (* b1 ecf-b2/b1)) (s (rand)))
            (if (< (car trials) c) (loop c (cdr b1s) (cdr trials))
              (begin
                (msg "Elliptic curve " c ": b1=" b1 ", b2=" b2 ", s=" s)
                (if (factor? 'ecf (ec-factor n b1 b2 s))
                    (loop 1 ecf-b1 ecf-trials)
                    (loop (+ c 1) b1s trials)))))))
      (cons (sort < facts) n)))))

Pages: 1 2 3

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 627 other followers

%d bloggers like this: