; 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

26
170
78501
1174
78517