codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
; counting primes (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) (let* ((max-index (quotient (- n 3) 2)) (v (make-vector (+ 1 max-index) #t))) (let loop ((i 0) (ps '(2))) (let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3))) (cond ((>= (* p p) n) (let loop ((j i) (ps ps)) (cond ((> j max-index) (reverse ps)) ((vector-ref v j) (loop (+ j 1) (cons (+ j j 3) ps))) (else (loop (+ j 1) ps))))) ((vector-ref v i) (let loop ((j startj)) (if (<= j max-index) (begin (vector-set! v j #f) (loop (+ j p))))) (loop (+ 1 i) (cons p ps))) (else (loop (+ 1 i) ps))))))) (define (make-pi-vector) (with-output-to-file "pi-vector.ss" (lambda () (display "#(") (newline) (display 0) (newline) (let* ((ps (cdr (primes (+ (isqrt #e1e12) 1)))) (qs (map (lambda (p) (modulo (* -1/2 (+ #e1e8 1 p)) p)) ps)) (z (length (primes #e1e8)))) (do ((t #e1e8 (+ t #e1e5 #e1e5)) (qs qs (map (lambda (p q) (modulo (- q #e1e5) p)) ps qs))) ((> t #e1e12) (display ")") (newline)) (when (zero? (modulo t #e1e8)) (display z) (newline)) (let ((bv (make-vector #e1e5 #t))) (do ((qs qs (cdr qs)) (ps ps (cdr ps))) ((null? qs)) (do ((j (car qs) (+ j (car ps)))) ((<= #e1e5 j)) (vector-set! bv j #f))) (do ((j 0 (+ j 1))) ((= j #e1e5)) (if (vector-ref bv j) (set! z (+ z 1)))))))))) (define (prime-pi n) (if (< #e1e12 n) (error 'prime-pi "out of range") (if (< n #e1e8) (length (primes n)) (let* ((n (if (even? n) (- n 1) n)) (b (quotient n #e1e8)) (l (* b #e1e8)) (ps (cdr (primes (+ (isqrt n) 1)))) (qs (map (lambda (p) (modulo (* -1/2 (+ l 1 p)) p)) ps)) (z (vector-ref pi-vector b)) (pi #f)) (do ((t l (+ t #e2e5)) (qs qs (map (lambda (p q) (modulo (- q #e1e5) p)) ps qs))) (pi pi) (let ((bv (make-vector #e1e5 #t))) (do ((qs qs (cdr qs)) (ps ps (cdr ps))) ((null? qs)) (do ((j (car qs) (+ j (car ps)))) ((<= #e1e5 j)) (vector-set! bv j #f))) (do ((j 0 (+ j 1))) ((= j #e1e5)) (when (vector-ref bv j) (set! z (+ z 1))) (when (= (+ t j j 1) n) (set! pi z))))))))) (define (nth-prime n) (define (bsearch n) (let loop ((lo 0) (hi 10000)) (if (< hi lo) hi (let ((mid (quotient (+ lo hi) 2))) (cond ((< (vector-ref pi-vector mid) n) (loop (+ mid 1) hi)) ((< n (vector-ref pi-vector mid)) (loop lo (- mid 1))) (else (- mid 1))))))) (if (< 37607912018 n) (error 'nth-prime "out of range") (if (< n #e1e8) (list-ref (primes #e1e8) (- n 1)) (let* ((b (bsearch n)) (l (* b #e1e8)) (ps (cdr (primes (+ (isqrt (+ l #e1e8)) 1)))) (qs (map (lambda (p) (modulo (* -1/2 (+ l 1 p)) p)) ps)) (z (vector-ref pi-vector b)) (nth #f)) (do ((t l (+ t #e2e5)) (qs qs (map (lambda (p q) (modulo (- q #e1e5) p)) ps qs))) (nth nth) (let ((bv (make-vector #e1e5 #t))) (do ((qs qs (cdr qs)) (ps ps (cdr ps))) ((null? qs)) (do ((j (car qs) (+ j (car ps)))) ((<= #e1e5 j)) (vector-set! bv j #f))) (do ((j 0 (+ j 1))) ((= j #e1e5)) (when (vector-ref bv j) (set! z (+ z 1))) (when (and (not nth) (= z n)) (set! nth (+ t j j 1))))))))))
Private
[
?
]
Run code
Submit