; the factorization of f7, part 1

(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) ; new version
  (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)))
          (if fs
              (loop a-1 a q r (- (* 2 g) r) q+1 (+ i 1) (cons (cons* i q a-1 fs) qs) (- lim 1))
              (loop a-1 a q r (- (* 2 g) r) q+1 (+ i 1) qs (- lim 1))))))))

(for-each
  (lambda (x) (display x) (newline))
  (residue 13290059 1 '(2 5 31 41 43 53 113) 44))

(5 2050 171341 2 41)
(10 1333 6700527 31 43)
(22 4633 5235158 41 113)
(23 226 1914221 2 113)
(26 3286 11455708 2 31 53)
(31 5650 1895246 2 113)
(40 4558 3213960 2 43 53)