; 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