; pollard's p-1 factorization algorithm

(define rand
  (let* ((a 3141592653) (c 2718281829)
         (m (expt 2 35)) (x 5772156649)
         (next (lambda ()
                 (let ((x-prime (modulo (+ (* a x) c) m)))
                   (set! x x-prime) x-prime)))
         (k 103)
         (v (list->vector (reverse
              (let loop ((i k) (vs (list x)))
                (if (= i 1) vs
                  (loop (- i 1) (cons (next) vs)))))))
         (y (next))
         (init (lambda (s)
                 (set! x s) (vector-set! v 0 x)
                 (do ((i 1 (+ i 1))) ((= i k))
                   (vector-set! v i (next))))))
    (lambda seed
      (cond ((null? seed)
              (let* ((j (quotient (* k y) m))
                     (q (vector-ref v j)))
                (set! y q)
                (vector-set! v j (next)) (/ y m)))
            ((eq? (car seed) 'get) (list a c m x y k v))
            ((eq? (car seed) 'set)
              (let ((state (cadr seed)))
                (set! a (list-ref state 0))
                (set! c (list-ref state 1))
                (set! m (list-ref state 2))
                (set! x (list-ref state 3))
                (set! y (list-ref state 4))
                (set! k (list-ref state 5))
                (set! v (list-ref state 6))))
            (else (init (modulo (numerator
                    (inexact->exact (car seed))) m))
                  (rand))))))

(define (randint . args)
  (cond ((null? (cdr args))
          (inexact->exact (floor (* (rand) (car args)))))
        ((< (car args) (cadr args))
          (+ (inexact->exact (floor (* (rand) (- (cadr args) (car args))))) (car args)))
        (else (+ (inexact->exact (ceiling (* (rand) (- (cadr args) (car args))))) (car args)))))

(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 (prime? n)
  (define (witness? a n)
    (let loop ((r 0) (s (- n 1)))
      (if (even? s) (loop (+ r 1) (/ s 2))
        (if (= (expm a s n) 1) #t
          (let loop ((j 0) (s s))
            (cond ((= j r) #f)
                  ((= (expm a s n) (- n 1)) #t)
                  (else (loop (+ j 1) (* s 2)))))))))
  (cond ((< n 2) #f) ((= n 2) #t) ((even? n) #f)
        (else (let loop ((k 50))
                (cond ((zero? k) #t)
                      ((not (witness? (randint 1 n) n)) #f)
                      (else (loop (- k 1))))))))

(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 (factor n)
  (let loop1 ((b 7) (k 420))
    (let loop2 ((a (randint 2 n)))
      (let ((g (gcd a n)))
        (if (< 1 g) g
          (let ((d (gcd (- (expm a k n) 1) n)))
            (cond ((= d 1)
                    (let ((b (+ b 1)))
                      (loop1 b (lcm k b))))
                  ((= d n) (loop2 (randint 2 n)))
                  (else d))))))))

(define (factors n)
  (sort < (let fact ((n n) (fs '()))
    (cond ((= n 1) fs)
          ((even? n) (fact (/ n 2) (cons 2 fs)))
          ((prime? n) (cons n fs))
          (else (let ((f (factor n)))
                  (append fs (fact f '()) (fact (/ n f) '()))))))))

(display (factors (- (expt 2 86) 1)))

(3 431 9719 2099863 2932031007403)