; more prime-counting functions

(define (make-hash hash eql? oops size)
  (let ((table (make-vector size '())))
    (lambda (message . args)
      (if (eq? message 'enlist)
          (let loop ((k 0) (result '()))
            (if (= size k)
                result
                (loop (+ k 1) (append (vector-ref table k) result))))
          (let* ((key (car args))
                 (index (modulo (hash key) size))
                 (bucket (vector-ref table index)))
            (case message
              ((lookup fetch get ref recall)
                (let loop ((bucket bucket))
                  (cond ((null? bucket) oops)
                        ((eql? (caar bucket) key) (cdar bucket))
                        (else (loop (cdr bucket))))))
              ((insert insert! ins ins! set set! store store! install install!)
                (vector-set! table index
                  (let loop ((bucket bucket))
                    (cond ((null? bucket)
                            (list (cons key (cadr args))))
                          ((eql? (caar bucket) key)
                            (cons (cons key (cadr args)) (cdr bucket)))
                          (else (cons (car bucket) (loop (cdr bucket))))))))
              ((delete delete! del del! remove remove!)
                (vector-set! table index
                  (let loop ((bucket bucket))
                    (cond ((null? bucket) '())
                          ((eql? (caar bucket) key)
                            (cdr bucket))
                          (else (cons (car bucket) (loop (cdr bucket))))))))
              ((update update!)
                (vector-set! table index
                  (let loop ((bucket bucket))
                    (cond ((null? bucket)
                            (list (cons key (caddr args))))
                          ((eql? (caar bucket) key)
                            (cons (cons key ((cadr args) key (cdar bucket))) (cdr bucket)))
                          (else (cons (car bucket) (loop (cdr bucket))))))))
              (else (error 'hash-table "unrecognized message")) ))))))

(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 (iroot k n) ; => m such that m^k <= n < (m+1)^k
  (let loop ((hi 1))
    (if (< (expt hi k) n) (loop (* hi 2))
      (let loop ((lo (/ hi 2)) (hi hi))
        (if (= (- hi lo) 1)
            (if (= (expt hi k) n) hi lo)
            (let* ((mid (quotient (+ lo hi) 2))
                   (mid^k (expt mid k)))
              (cond ((< mid^k n) (loop mid hi))
                    ((< n mid^k) (loop lo mid))
                    (else mid))))))))

(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 max-pi #e32e5)

(define ps (list->vector (cons #f (primes max-pi))))

(define max-p (- (vector-length ps) 1))

(define (p n)
 (if (< max-p n) (nth-prime n)
   (vector-ref ps n)))

(define (pi n)
  (if (< max-pi n) (prime-pi n)
    (let loop ((lo 1) (hi (- (vector-length ps) 1)))
      (let ((mid (quotient (+ lo hi) 2)))
        (cond ((< hi lo) mid)
              ((< n (p mid)) (loop lo (- mid 1)))
              ((< (p mid) n) (loop (+ mid 1) hi))
              (else mid))))))

(define phi
  (let ((memo (make-hash (lambda (x) (+ (* 100000 (car x)) (cadr x))) equal? #f 999983)))
    (lambda args
      (let ((x (car args)) (a (cadr args)) (t (memo 'lookup args)))
        (cond (t t) ; return memoized value
              ((= a 1) (let ((t (quotient (+ x 1) 2))) (memo 'insert args t) t))
              (else (let ((t (- (phi x (- a 1)) (phi (quotient x (p a)) (- a 1)))))
                      (memo 'insert args t) t)))))))

(define (legendre-pi n)
  (if (< n max-pi) (pi n)
    (let ((a (pi (isqrt n))))
      (+ (phi n a) a -1))))

(define (meissel-pi n)
  (if (< n max-pi) (pi n)
    (let ((b (pi (isqrt n))) (c (pi (iroot 3 n))))
      (let loop ((i (+ c 1)) (sum (+ (phi n c) (quotient (* (+ b c -2) (- b c -1)) 2))))
        (if (< b i) sum
          (loop (+ i 1) (- sum (pi (quotient n (p i))))))))))

(define (lehmer-pi n)
  (if (< n max-pi) (pi n)
    (let ((a (pi (iroot 4 n))) (b (pi (isqrt n))) (c (pi (iroot 3 n))))
      (let i-loop ((i (+ a 1)) (sum (+ (phi n a) (quotient (* (+ b a -2) (- b a -1)) 2))))
        (if (< b i) sum
          (let* ((w (quotient n (p i))) (lim (pi (isqrt w))) (sum (- sum (pi w))))
            (if (< c i) (i-loop (+ i 1) sum)
              (let j-loop ((j i) (sum sum))
                (if (< lim j) (i-loop (+ i 1) sum)
                  (j-loop (+ j 1) (- sum (pi (quotient w (p j))) (- j) 1)))))))))))

(define prime-pi legendre-pi)
(time (display (prime-pi #e1e7)) (newline))

(define prime-pi meissel-pi)
(time (display (prime-pi #e1e8)) (newline))

(define prime-pi lehmer-pi)
(time (display (prime-pi #e1e9)) (newline))

664579
cpu time: 220 real time: 1147 gc time: 70
5761455
cpu time: 160 real time: 1058 gc time: 80
50847534
cpu time: 250 real time: 1934 gc time: 70