```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 ``` ```; 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)) (newline) (time (display (length (primes1mod4 1000000 2000000 25000))) (newline)) ```
 ```1 2 3 ``` ```(101 109 113 137 149 157 173 181 193 197 229 233 241 257 269 277 281 289 293) 35241 cpu time: 50 real time: 954 gc time: 0 ```