codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
; prime-count and prime-sum ; http://www.icecreambreakfast.com/primecount/PrimeCountingSurvey.pdf (define (factors n) ; trial division (let loop ((n n) (f 2) (fs (list))) (cond ((< n (* f f)) (reverse (if (= n 1) fs (cons n fs)))) ((zero? (modulo n f)) (loop (/ n f) f (cons f fs))) (else (loop n (+ f 1) fs))))) (define (mobius n) (if (= n 1) 1 (let loop ((fs (factors n)) (mob 1) (prev 0)) (cond ((null? fs) mob) ((= (car fs) prev) 0) (else (loop (cdr fs) (- mob) (car fs))))))) (define (choose n k) ; binomial theorem (let loop ((n n) (k k) (c 1)) (if (zero? k) c (loop (- n 1) (- k 1) (* c n (/ k)))))) (define (iroot k n) ; largest integer x such that x^k <= n (let ((k-1 (- k 1))) (let loop ((u n) (s (+ n 1))) (if (<= s u) s (loop (quotient (+ (* k-1 u) (quotient n (expt u k-1))) k) u))))) (define (ilog b n) ; largest integer x such that b^x <= n (let loop1 ((lo 0) (b^lo 1) (hi 1) (b^hi b)) (if (< b^hi n) (loop1 hi b^hi (* hi 2) (* b^hi b^hi)) (let loop2 ((lo lo) (b^lo b^lo) (hi hi) (b^hi b^hi)) (if (<= (- hi lo) 1) (if (= b^hi n) hi lo) (let* ((mid (quotient (+ lo hi) 2)) (b^mid (* b^lo (expt b (- mid lo))))) (cond ((< n b^mid) (loop2 lo b^lo mid b^mid)) ((< b^mid n) (loop2 mid b^mid hi b^hi)) (else mid)))))))) (define-syntax sigma (syntax-rules () ((sigma var lo hi func ...) (let loop ((var lo) (sum 0)) (if (< hi var) sum (loop (+ var 1) (+ sum ((lambda (var) func ...) var)))))))) (define (prime-count n) (define (d k a n) (if (zero? k) 1 (if (= k 1) (- n a -1) (sigma j 1 k (* (choose k j) (sigma m a (iroot k n) (d (- k j) (+ m 1) (quotient n (expt m j))))))))) (define (p n) (sigma k 1 (ilog 2 n) (quotient (* (expt -1 (+ k 1)) (d k 2 n)) k))) (sigma j 1 (ilog 2 n) (quotient (* (mobius j) (p (iroot j n))) j))) (define (prime-sum n) (define (d s k a n) (if (zero? k) 1 (if (= k 1) (- (* n (+ n 1) 1/2) (* a (- a 1) 1/2)) (sigma j 1 k (* (choose k j) (sigma m a (iroot k n) (* (expt m (* s j)) (d s (- k j) (+ m 1) (quotient n (expt m j)))))))))) (define (p s n) (sigma k 1 (ilog 2 n) (quotient (* (expt -1 (+ k 1)) (d s k 2 n)) k))) (sigma j 1 (ilog 2 n) (quotient (* (mobius j) (p j (iroot j n))) j))) (display (prime-count 100)) (newline) ; 25 (display (prime-count 1000)) (newline) ; 168 (display (prime-count 1000000)) (newline) ; 78498 (display (prime-sum 100)) (newline) ; 1060 (display (prime-sum 1000)) (newline) ; 76127
Private
[
?
]
Run code
Submit