; pollard rho, revisited
(define (cons* first . rest)
(let loop ((curr first) (rest rest))
(if (null? rest) curr
(cons curr (loop (car rest) (cdr rest))))))
(define (last-pair xs)
(if (pair? (cdr xs))
(last-pair (cdr xs))
xs))
(define sort #f)
(define merge #f)
(let ()
(define dosort
(lambda (pred? ls n)
(if (= n 1)
(list (car ls))
(let ((i (quotient n 2)))
(domerge pred?
(dosort pred? ls i)
(dosort pred? (list-tail ls i) (- n i)))))))
(define domerge
(lambda (pred? l1 l2)
(cond
((null? l1) l2)
((null? l2) l1)
((pred? (car l2) (car l1))
(cons (car l2) (domerge pred? l1 (cdr l2))))
(else (cons (car l1) (domerge pred? (cdr l1) l2))))))
(set! sort
(lambda (pred? l)
(if (null? l) l (dosort pred? l (length l)))))
(set! merge
(lambda (pred? l1 l2)
(domerge pred? l1 l2))))
(define (prime? n)
(letrec ((expm (lambda (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))))))
(digits (lambda (n)
(let loop ((n n) (ds '()))
(if (zero? n) ds
(loop (quotient n 2) (cons (modulo n 2) ds))))))
(isqrt (lambda (n)
(let loop ((x n) (y (quotient (+ n 1) 2)))
(if (<= 0 (- y x) 1) x
(loop y (quotient (+ y (quotient n y)) 2))))))
(square? (lambda (n)
(let ((n2 (isqrt n))) (= n (* n2 n2)))))
(jacobi (lambda (a n)
(let loop ((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) (* (loop 2 n) (loop (/ a 2) n)))
((< n a) (loop (modulo a n) n))
((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (loop n a)))
(else (loop n a))))))
(inverse (lambda (x m)
(let loop ((a 1) (b 0) (g x) (u 0) (v 1) (w m))
(if (zero? w) (modulo a m)
(let ((q (quotient g w)))
(loop u v w (- a (* q u)) (- b (* q v)) (- g (* q w))))))))
(strong-pseudo-prime? (lambda (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))))))))))
(chain (lambda (m f g u v)
(let loop ((ms (digits m)) (u u) (v v))
(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)))))))
(lucas-pseudo-prime? (lambda (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 (jacobi 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)))))))))))
(if (not (integer? n)) (error 'prime? "non-integer")
(if (< n 2) #f (if (even? n) (= n 2)
(if (zero? (modulo n 3)) (= n 3)
(and (strong-pseudo-prime? n 2)
(strong-pseudo-prime? n 3)
(lucas-pseudo-prime? n))))))))
(define (wheel-factors n limit)
(define (wheel . xs) (set-cdr! (last-pair xs) xs) xs)
(let loop ((n n) (f 2) (fs '())
(ws (cons* 1 2 2 (wheel 4 2 4 2 4 6 2 6))))
(cond ((< limit f) (values n (reverse fs)))
((or (= n 1) (< n (* f f))) (values 1 (reverse (cons n fs))))
((zero? (modulo n f)) (loop (/ n f) f (cons f fs) ws))
(else (loop n (+ f (car ws)) fs (cdr ws))))))
(define (pollard factor n)
(let-values (((n fs) (wheel-factors n 1000)))
(sort < (let fact ((n n) (fs fs))
(cond ((= n 1) fs)
((prime? n) (cons n fs))
(else (let ((f (factor n 1 100000000)))
(append fs (fact f '()) (fact (/ n f) '())))))))))
(define (floyd n c limit)
(define (f x) (modulo (+ (* x x) c) n))
(define (g p t h) (modulo (* p (abs (- t h))) n))
(let loop1 ((j 1) (t 2) (h (f 2)) (x 2) (y (f 2)) (p 1))
(if (= j limit) (error 'floyd "timeout")
(if (= t h) (floyd n (+ c 1) (- limit j))
(if (positive? (modulo j 100)) (loop1 (+ j 1) (f t) (f (f h)) x y (g p t h))
(let ((d (gcd p n)))
(if (= d 1) (loop1 (+ j 1) (f t) (f (f h)) t h 1)
(if (< 1 d n) d
(let loop2 ((k 1) (x x) (y y) (d (gcd (- x y) n)))
(if (= k 100) (floyd n (+ c 1) (- limit j))
(if (= d 1) (loop2 (+ k 1) (f x) (f (f y)) (gcd (- x y) n))
(if (= d n) (floyd n (+ c 1) (- limit j))
d))))))))))))
(define (brent n c limit)
(define (f y) (modulo (+ (* y y) c) n))
(define (g p x y) (modulo (* p (abs (- x y))) n))
(let loop1 ((x 2) (y (+ 4 c)) (z (+ 4 c)) (j 1) (q 2) (p 1))
(if (= j limit) (error 'brent "timeout")
(if (= x y) (brent n (+ c 1) (- limit j)) ; cycle
(if (= j q) (let ((t (f y))) (loop1 y (f y) z (+ j 1) (* q 2) (g p y t)))
(if (positive? (modulo j 100)) (loop1 x (f y) z (+ j 1) q (g p x y))
(let ((d (gcd p n)))
(if (= d 1) (loop1 x (f y) y (+ j 1) q (g p x y))
(if (< 1 d n) d ; factor
(let loop2 ((k 1) (z (f z)))
(if (= k 100) (brent n (+ c 1) (- limit j))
(let ((d (gcd (- z x) n)))
(if (= d 1) (loop2 (+ k 1) (f z))
(if (= d n) (brent n (+ c 1) (- limit j))
d))))))))))))))
(define (fib n)
(define (square x) (* x x))
(cond ((= n 0) 0) ((= n 1) 1) ((= n 2) 1)
((odd? n) (let ((n2 (+ (quotient n 2) 1)))
(+ (square (fib n2)) (square (fib (- n2 1))))))
(else (let ((n2 (quotient n 2)))
(* (fib n2) (+ (* (fib (- n2 1)) 2) (fib n2)))))))
(time (do ((i 3 (+ i 1))) ((< 100 i))
(display i) (display " ")
(display (pollard floyd (fib i)))
(newline)))
(time (do ((i 3 (+ i 1))) ((< 100 i))
(display i) (display " ")
(display (pollard brent (fib i)))
(newline)))