; 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)))