; sieving a polynomial

(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 (primes n) ; list of primes not exceeding n
  (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
    (let loop ((i 0) (p 3) (ps (list 2)))
      (cond ((< n (* p p))
              (do ((i i (+ i 1)) (p p (+ p 2))
                   (ps ps (if (vector-ref bits i) (cons p ps) ps)))
                  ((= i len) (reverse ps))))
            ((vector-ref bits i)
              (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
                  ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
                (vector-set! bits j #f)))
            (else (loop (+ i 1) (+ p 2) ps))))))

(define (factors-1mod4 n)

    (let* ((len (quotient (- n 1) 4))
           (factors (make-vector (+ len 1) (list))))

      (define (sieve p pp) ; p = prime, pp = prime-power
        (let* ((k (if (= (modulo pp 4) 1) pp (* 3 pp)))
               (k (quotient (- k 1) 4)))
          (do ((i k (+ i pp))) ((< len i))
            (vector-set! factors i
              (cons p (vector-ref factors i))))))

      (do ((ps (cdr (primes (isqrt n))) (cdr ps))) ((null? ps))
        (do ((pp (car ps) (* pp (car ps)))) ((< n pp))
          (sieve (car ps) pp)))

      (do ((k 0 (+ k 1))) ((< len k) factors)
        (let* ((n (+ (* 4 k) 1))
               (fs (vector-ref factors k))
               (f (/ n (apply * fs))))
          (vector-set! factors k
            (cond ((null? fs) (list n))
                  ((= f 1) (reverse fs))
                  (else (reverse (cons f fs)))))))))

(define fs1mod4 (factors-1mod4 10001))

(do ((k 1 (+ k 1))) ((= (vector-length fs1mod4) k))
  (display (+ (* 4 k) 1)) (display #\tab)
  (display (vector-ref fs1mod4 k)) (newline))