; pollard's p-1 factorization algorithm, revisited
(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 (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 (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 pminus1 ; no backtracking
(letrec ((diffs (lambda (xs)
(let loop ((x (car xs)) (xs (cdr xs)) (zs (list)))
(if (null? (cdr xs))
(reverse (cons (- (car xs) x) zs))
(loop (car xs) (cdr xs) (cons (- (car xs) x) zs))))))
(last-pair (lambda (xs)
(if (null? (cdr xs)) xs (last-pair (cdr xs)))))
(cycle (lambda xs (set-cdr! (last-pair xs) xs) xs)))
(let* ((b1 #e16e4) (b2 #e32e5) (ps (primes b2))
(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))))))
(ds (map (lambda (x) (/ x 2)) (diffs ps))))
(lambda (n)
(let loop1 ((ps ps) (as as) (ds ds) (q 2))
; (display "loop1: ") (display q) (newline)
(let ((g (gcd (- q 1) n)))
(if (< 1 g n) g
(if (< (car ps) b1)
(loop1 (cdr ps) (cdr as) (cdr ds)
(expm q (expt (car ps) (car as)) n))
(let* ((len (apply max ds))
(dv (make-vector (+ len 1) 0)))
(do ((i 1 (+ i 1))) ((< len i))
(vector-set! dv i (expm q (+ i i) n)))
(let loop2 ((ds ds) (q (expm q (car ps) n)) (p 1))
; (display "loop2: ") (display q)
; (display " ") (display p) (newline)
(let ((g (gcd p n)))
(if (< 1 g n) g
(if (null? ds) #f
(let ((q (modulo (* q (vector-ref dv (car ds))) n)))
(loop2 (cdr ds) q (modulo (* p (- q 1)) n))))))))))))))))
(time (begin (display (pminus1 1355145973557738676476835291149239036254674001)) (newline)))
(define pminus1 ; with backtracking
(letrec ((diffs (lambda (xs)
(let loop ((x (car xs)) (xs (cdr xs)) (zs (list)))
(if (null? (cdr xs))
(reverse (cons (- (car xs) x) zs))
(loop (car xs) (cdr xs) (cons (- (car xs) x) zs))))))
(last-pair (lambda (xs)
(if (null? (cdr xs)) xs (last-pair (cdr xs)))))
(cycle (lambda xs (set-cdr! (last-pair xs) xs) xs)))
(let* ((b1 #e16e4) (b2 #e32e5) (ps (primes b2))
(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))))))
(ds (map (lambda (x) (/ x 2)) (diffs ps))))
(lambda (n)
(let loop1 ((ps ps) (as as) (ds ds) (j 1) (q 2) (p 1) (ps-saved ps) (as-saved as) (q-saved 2))
; (display "loop1: ") (display q) (newline)
(if (and (positive? (modulo j 100)) (< (car ps) b1))
(let ((q (expm q (expt (car ps) (car as)) n)))
(loop1 (cdr ps) (cdr as) (cdr ds) (+ j 1) q (modulo (* p (- q 1)) n) ps-saved as-saved q-saved))
(let ((g (gcd p n)))
(cond ((< 1 g n) g)
((= g n) (let loop ((ps ps-saved) (as as-saved) (q q-saved))
(let ((g (gcd (- q 1) n)))
(if (= g n) #f (if (< 1 g) g
(loop (cdr ps) (cdr as) (expm q (expt (car ps) (car as)) n)))))))
((< (car ps) b1) (loop1 ps as ds 1 q 1 ps as q))
(else (let* ((len (apply max ds)) (dv (make-vector (+ len 1) 0)))
(do ((i 1 (+ i 1))) ((< len i)) (vector-set! dv i (expm q (+ i i) n)))
(let loop2 ((ds ds) (j 1) (q (expm q (car ps) n)) (p 1) (ds-saved ds) (q-saved (expm q (car ps) n)))
; (display "loop2: ") (display q) (display " ") (display p) (newline)
(if (and (positive? (modulo j 100)) (pair? ds))
(let ((q (modulo (* q (vector-ref dv (car ds))) n)))
(loop2 (cdr ds) (+ j 1) q (modulo (* p (- q 1)) n) ds-saved q-saved))
(let ((g (gcd p n)))
(cond ((< 1 g n) g)
((= g n) (let loop ((ds ds-saved) (q q-saved))
(let ((g (gcd (- q 1) n)))
(if (= g n) #f (if (< 1 g) g
(loop (cdr ds) (modulo (* q (vector-ref dv (car ds))) n)))))))
((pair? ds) (loop2 ds 1 q 1 ds q))
(else #f)))))))))))))))
(time (begin (display (pminus1 1355145973557738676476835291149239036254674001)) (newline)))
programmingpraxis
-
Scheme,
pasted
on Sep 14: