; dixon's factorization algorithm

(define (cons* first . rest)
  (let loop ((curr first) (rest rest))
    (if (null? rest) curr
        (cons curr (loop (car rest) (cdr rest))))))

(define (isqrt n)
  (if (not (and (positive? n) (integer? n)))
      (error 'isqrt "must be positive integer")
      (let loop ((x n))
        (let ((y (quotient (+ x (quotient n x)) 2)))
          (if (< y x) (loop y) x)))))

(define (jacobi a n)
  (if (not (and (integer? a) (integer? n) (positive? n) (odd? n)))
      (error 'jacobi "modulus must be positive odd integer")
      (let jacobi ((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) (* (jacobi 2 n) (jacobi (quotient a 2) n)))
              ((< n a) (jacobi (modulo a n) n))
              ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a)))
              (else (jacobi n a))))))

(define (primes n)
  (let* ((max-index (quotient (- n 3) 2))
         (v (make-vector (+ 1 max-index) #t)))
    (let loop ((i 0) (ps '(2)))
      (let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
        (cond ((>= (* p p) n)
               (let loop ((j i) (ps ps))
                  (cond ((> j max-index) (reverse 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))
                  (if (<= j max-index)
                      (begin (vector-set! v j #f)
                             (loop (+ j p)))))
                      (loop (+ 1 i) (cons p ps)))
              (else (loop (+ 1 i) ps)))))))

(define (make-factor-base n k bound lim)
  (let loop ((i (- lim 1)) (ps (cdr (primes bound))) (fb '(2)))
    (cond ((or (zero? i) (null? ps)) (reverse fb))
          ((not (negative? (jacobi (* k n) (car ps))))
            (loop (- i 1) (cdr ps) (cons (car ps) fb)))
          (else (loop i (cdr ps) fb)))))

(define (smooth n q fb)
  (let loop ((q q) (fb fb) (fs '()) (i 0))
    (cond ((= q 1) (reverse fs))
          ((null? fb) #f)
          ((zero? (modulo q (car fb)))
            (if (and (pair? fs) (= (car fb) (car fs)))
                (loop (/ q (car fb)) fb (cdr fs) (+ i 1))
                (loop (/ q (car fb)) fb (cons (car fb) fs) (+ i 1))))
          (else (loop q (cdr fb) fs (+ i 1))))))

;(define (residue n k fb lim)
;  (let* ((kn (* k n)) (g (isqrt kn)) (a-3 0) (a-2 1)
;         (q-2 kn) (r-2 g) (g+p-1 g) (q-1 1)
;         (t-1 (quotient g+p-1 q-1)) (r-1 (- g+p-1 (* q-1 t-1)))
;         (q0 (+ q-2 (* t-1 (- r-1 r-2)))))
;    (let loop ((a-2 1) (a-1 g) (q-1 1) (r-1 0) (g+p (+ g g)) (q q0) (i 1) (qs '()) (lim lim))
;      ; (for-each display `(,i " " ,g+p " " ,q " " ,a-1 #\newline))
;      (if (or (= q 1) (zero? lim)) (reverse qs)
;        (let* ((t (quotient g+p q))
;               (r (- g+p (* q t)))
;               (a (modulo (+ (* t a-1) a-2) n))
;               (q+1 (+ q-1 (* t (- r r-1))))
;               (fs (smooth n q fb)))
;          (cond ((null? fs)
;                  (let ((d (gcd (- a-1 (isqrt q)) n)))
;                    (if (< 1 d n) d
;                      (loop a-1 a q r (- (* 2 g) r) q+1 (+ i 1) qs (- lim 1)))))
;                (fs (loop a-1 a q r (- (* 2 g) r) q+1 (+ i 1) (cons (cons* i q a-1 fs) qs) (- lim 1)))
;                (else (loop a-1 a q r (- (* 2 g) r) q+1 (+ i 1) qs (- lim 1)))))))))

(define (residue n k fb lim)
  (let loop ((r (isqrt n)) (qs '()))
    (if (zero? lim) (reverse qs)
      (let* ((r2n (modulo (* r r) n))
             (fs (smooth n r2n fb)))
        (cond ((null? fs)
                (let ((d (gcd (- r (isqrt r2n)) n)))
                  (if (< 1 d n) d (loop (+ r 1) qs))))
              (fs (loop (+ r 1) (cons (cons* 0 r2n r fs) qs)))
              (else (loop (+ r 1) qs)))))))

; (fibonacci 73) => 9375829 * 86020717
(define fb (make-factor-base 806515533049393 1 60000 270))
(define residues (residue 806515533049393 1 fb 2000000))

(define (make-row rs fb len)
  (let ((vec (make-vector (+ len 1) 0)))
    (vector-set! vec 0 (if (odd? (car rs)) 1 0))
    (let loop ((i 1) (rs (cdddr rs)) (fb fb))
      ; (for-each display `(,i " " ,rs " " ,fb " " #\newline))
      (cond ((or (null? rs) (null? fb)) vec)
            ((= (car fb) (car rs))
              (vector-set! vec i 1)
              (loop (+ i 1) (cdr rs) (cdr fb)))
            (else (loop (+ i 1) rs (cdr fb)))))))

(define (make-ev fb len residues)
  (let loop ((rss residues) (es '()) (as '()) (qs '()))
    (if (null? rss)
        (values (list->vector (reverse es))
                (list->vector (reverse as))
                (list->vector (reverse qs)))
        (loop (cdr rss)
              (cons (make-row (car rss) fb len) es)
              (cons (caddar rss) as)
              (cons (cadar rss) qs)))))

(define (make-hv ev-len)
  (let ((hv (make-vector ev-len)))
    (do ((i 0 (+ i 1))) ((= i ev-len) hv)
      (let ((h (make-vector ev-len 0)))
        (vector-set! h i 1)
        (vector-set! hv i h)))))

(define (right-most-one rv start)
  (let loop ((i start))
    (cond ((negative? i) i)
          ((= (vector-ref rv i) 1) i)
          (else (loop (- i 1))))))

(define (vector-add v1 v2)
  (let* ((len (vector-length v1))
         (v3 (make-vector len)))
    (do ((i 0 (+ i 1))) ((= i len) v3)
      (vector-set! v3 i
        (if (= (vector-ref v1 i)
               (vector-ref v2 i)) 0 1)))))

(define (vector-sum vec)
  (do ((i 0 (+ i 1)) (sum 0 (+ sum (vector-ref vec i))))
      ((= i (vector-length vec)) sum)))

(define (a-minus-q n vec ev-len av qv)
  (let loop ((i 0) (a-prod 1) (q-prod 1))
    (cond ((= i ev-len) (modulo (- a-prod (isqrt q-prod)) n))
          ((= (vector-ref vec i) 1)
            (loop (+ i 1)
                  (modulo (* a-prod (vector-ref av i)) n)
                  (* q-prod (vector-ref qv i))))
          (else (loop (+ i 1) a-prod q-prod)))))

(define (answer n fb residues)
  (let-values (((ev av qv) (make-ev fb (length fb) residues)))
    (let* ((fb-len (length fb)) (ev-len (vector-length ev))
           (hv (make-hv ev-len)) (ptr (make-vector ev-len)))
      (do ((i 0 (+ i 1))) ((= i ev-len))
        (vector-set! ptr i (right-most-one (vector-ref ev i) fb-len)))
      ; reduce ev and hv
      (do ((j fb-len (- j 1))) ((negative? j))
        (let loop ((i 0))
          (when (< i ev-len)
            (if (= (vector-ref ptr i) j)
                (do ((m (+ i 1) (+ m 1))) ((= m ev-len))
                  (when (= (vector-ref ptr m) j)
                    (vector-set! ev m (vector-add (vector-ref ev i) (vector-ref ev m)))
                    (vector-set! hv m (vector-add (vector-ref hv i) (vector-ref hv m)))
                    (vector-set! ptr m (right-most-one (vector-ref ev m) j))))
                (loop (+ i 1))))))
      ; examine s-sets
      (let loop ((i 0))
        (cond ((= i ev-len) #f)
              ((zero? (vector-sum (vector-ref ev i)))
                (let ((d (gcd (a-minus-q n (vector-ref hv i) ev-len av qv) n)))
                  (if (< 1 d n) d (loop (+ i 1)))))
              (else (loop (+ i 1))))))))

(answer 806515533049393 fb residues)