; primes congruent to 1 (mod 4)

(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 (euclid x y)
  (let loop ((a 1) (b 0) (g x) (u 0) (v 1) (w y))
    (if (zero? w) (values a b g)
      (let ((q (quotient g w)))
        (loop u v w (- a (* q u)) (- b (* q v)) (- g (* q w)))))))

(define (inverse x m)
  (if (not (= (gcd x m) 1))
      (error 'inverse "divisor must be coprime to modulus")
      (call-with-values
        (lambda () (euclid x m))
        (lambda (a b g) (modulo a m)))))

(define (primes n)
  (let ((sieve (make-vector (+ n 1) #t)))
    (let loop ((p 2) (ps (list)))
      (cond ((= p n) (reverse ps))
            ((vector-ref sieve p)
              (do ((i (* p p) (+ i p))) ((< n i))
                (vector-set! sieve i #f))
              (loop (+ p 1) (cons p ps)))
            (else (loop (+ p 1) ps))))))

(define (primes1mod4 lo hi delta)
  (let* ((output (list))
         (sieve (make-vector delta #t))
         (ps (cdr (primes (isqrt hi))))
         (qs (map (lambda (p) (modulo (* -1 (inverse 4 p) (+ lo p 1)) p)) ps)))
    (let loop ((lo lo) (qs qs))
      (if (not (< lo hi)) (reverse output)
        (begin
          (do ((i 0 (+ i 1))) ((= i delta)) (vector-set! sieve i #t))
          (do ((ps ps (cdr ps)) (qs qs (cdr qs))) ((null? ps))
            (do ((j (car qs) (+ j (car ps)))) ((<= delta j))
              (vector-set! sieve j #f)))
          (do ((i 0 (+ i 1)) (t (+ lo 1) (+ t 4)))
              ((or (<= delta i) (<= hi t)))
            (if (vector-ref sieve i) (set! output (cons t output))))
          (loop (+ lo (* 4 delta))
                (map (lambda (p q) (modulo (- q delta) p)) ps qs)))))))

(display (primes1mod4 100 300 25))

(101 109 113 137 149 157 173 181 193 197 229 233 241 257 269 277 281 289 293)