```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 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 ``` ```; two sieving problems (define (primes n) (let ((ps (list)) (sieve (make-vector (+ n 1) #t))) (do ((p 2 (+ p 1))) ((< n p) (reverse ps)) (when (vector-ref sieve p) (set! ps (cons p ps)) (do ((i (* p p) (+ i p))) ((< n i)) (vector-set! sieve i #f)))))) (define prime? (let ((seed 3141592654)) (lambda (n) (define (rand) (set! seed (modulo (+ (* 69069 seed) 1234567) 4294967296)) (+ (quotient (* seed (- n 2)) 4294967296) 2)) (define (expm b e m) (define (times x y) (modulo (* x y) m)) (let loop ((b b) (e e) (r 1)) (if (zero? e) r (loop (times b b) (quotient e 2) (if (odd? e) (times b r) r))))) (define (spsp? n a) (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1))) ((odd? d) (let ((t (expm a d n))) (if (or (= t 1) (= t (- n 1))) #t (do ((s (- s 1) (- s 1)) (t (expm t 2 n) (expm t 2 n))) ((or (zero? s) (= t (- n 1))) (positive? s)))))))) (if (not (integer? n)) (error 'prime? "must be integer") (if (< n 2) #f (if (= n 2) #t (do ((a (rand) (rand)) (k 10 (- k 1))) ((or (zero? k) (not (spsp? n a))) (zero? k))))))))) (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 (big-primes lo hi delta limit) (let* ((output (list)) (sieve (make-vector delta #t)) (ps (cdr (primes limit))) (qs (map (lambda (p) (modulo (* -1/2 (+ 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 2))) ((or (<= delta i) (<= hi t))) (if (and (vector-ref sieve i) (prime? t)) (set! output (cons t output)))) (loop (+ lo (* 2 delta)) (map (lambda (p q) (modulo (- q delta) p)) ps qs))))))) (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 (length (big-primes #e1e25 (+ #e1e25 20000) 5000 250000))) (newline) (display (length (primes1mod4 1000000 2000000 25000))) (newline) (define (omega n k) (let ((sieve (make-vector (+ n 1) 0)) (zs (list))) (do ((p 2 (+ p 1))) ((< n p)) (when (zero? (vector-ref sieve p)) (do ((i p (+ i p))) ((< n i)) (vector-set! sieve i (+ (vector-ref sieve i) 1))))) (do ((p 2 (+ p 1))) ((< n p) (reverse zs)) (when (= (vector-ref sieve p) k) (set! zs (cons p zs)))))) (display (omega 25 1)) ```
 ```1 2 3 ``` ```327 35241 (2 3 4 5 7 8 9 11 13 16 17 19 23 25)```