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