codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
; 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)))
Private
[
?
]
Run code
Submit