; the factorization of f7

(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 f7 (+ (expt 2 (expt 2 7) 1))
(define fb (make-factor-base f7 257 60000 2700))
(define residues (residue f7 257 fb 1330000))

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

(display (answer f7 fb residues))