Project:
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 ``` ```; 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))) ```
 ```1 2 3 4 ``` ```1091346396980401 cpu time: 10 real time: 50 gc time: 10 1091346396980401 cpu time: 0 real time: 48 gc time: 0 ```