; 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 (inverse x n)
(let loop ((x (modulo x n)) (a 1))
(cond ((zero? x) (error 'inverse "division by zero"))
((= x 1) a)
(else (let ((q (- (quotient n x))))
(loop (+ n (* q x)) (modulo (* q a) n)))))))
(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))