[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/KgVfGvEG    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Jul 28:
; the nth prime

;;;;; small-p and small-pi for base of recursion

(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))))))

;;;;; calculation of Legendre's sum phi(x,a)

(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 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)))))))

;;;;; the prime-counting function, by Lehmer

(define (iroot k n) ; => m such that m^k <= n < (m+1)^k
  (if (= n 1) 1
    (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 (prime-pi n)
  (if (< n max-pi) (pi n)
    (let ((a (pi (iroot 4 n))) (b (pi (iroot 2 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 (iroot 2 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)))))))))))

;;;;; approximate pi

(define (logint x)
  (let ((gamma 0.57721566490153286061) (log-x (log x)))
    (let loop ((k 1) (fact 1) (num log-x)
               (sum (+ gamma (log log-x) log-x)))
      (if (< 100 k) sum
        (let* ((k (+ k 1))
               (fact (* fact k))
               (num (* num log-x))
               (sum (+ sum (/ num fact k))))
          (loop k fact num sum))))))

(define (factors n) ; trial division
  (let loop ((n n) (fs (list)))
    (if (even? n) (loop (/ n 2) (cons 2 fs))
      (if (= n 1) (if (null? fs) (list 1) fs)
        (let loop ((n n) (f 3) (fs fs))
          (cond ((< n (* f f)) (reverse (cons n fs)))
                ((zero? (modulo n f))
                  (loop (/ n f) f (cons f fs)))
                (else (loop n (+ f 2) fs))))))))

(define (mobius-mu n)
  (if (= n 1) 1
    (let loop ((fs (factors n)) (prev 0) (m 1))
      (cond ((null? fs) m)
            ((= (car fs) prev) 0)
            (else (loop (cdr fs) (car fs) (- m)))))))

(define riemann-r
  (let ((ms (let loop ((n 1) (k 1000) (ms (list)))
              (if (zero? k) (reverse ms)
                (let ((m (mobius-mu n)))
                  (if (zero? m) (loop (+ n 1) k ms)
                    (loop (+ n 1) (- k 1) (cons (* m n) ms))))))))
    (lambda (x)
      (let loop ((ms ms) (sum 0))
        (if (null? ms) sum
          (let* ((m (car ms)) (m-abs (abs m)) (m-recip (/ m)))
            (loop (cdr ms) (+ sum (* m-recip (logint (expt x (/ m-abs))))))))))))

;;;;; the nth prime

(define (prime? n)
  (letrec (
    (expm (lambda (b e m)
      (let ((times (lambda (x y) (modulo (* x y) m))))
        (cond ((zero? e) 1) ((even? e) (expm (times b b) (/ e 2) m))
        (else (times b (expm (times b b) (/ (- e 1) 2) m)))))))
    (digits (lambda (n)
      (let loop ((n n) (ds '()))
        (if (zero? n) ds (loop (quotient n 2) (cons (modulo n 2) ds))))))
    (isqrt (lambda (n)
      (let loop ((x n) (y (quotient (+ n 1) 2)))
        (if (<= 0 (- y x) 1) x (loop y (quotient (+ y (quotient n y)) 2))))))
    (square? (lambda (n) (let ((n2 (isqrt n))) (= n (* n2 n2)))))
    (jacobi (lambda (a n)
      (let loop ((a a) (n n))
        (cond ((= a 0) 0) ((= a 1) 1)
              ((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1)))
              ((even? a) (* (loop 2 n) (loop (/ a 2) n)))
              ((< n a) (loop (modulo a n) n))
              ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (loop n a)))
              (else (loop n a))))))
    (inverse (lambda (x m)
      (let loop ((a 1) (b 0) (g x) (u 0) (v 1) (w m))
        (if (zero? w) (modulo a m)
          (let ((q (quotient g w)))
            (loop u v w (- a (* q u)) (- b (* q v)) (- g (* q w))))))))
    (strong-pseudo-prime? (lambda (n a)
      (let loop ((r 0) (s (- n 1)))
        (if (even? s) (loop (+ r 1) (/ s 2))
          (if (= (expm a s n) 1) #t
            (let loop ((r r) (s s))
              (cond ((zero? r) #f) ((= (expm a s n) (- n 1)) #t)
              (else (loop (- r 1) (* s 2))))))))))
    (chain (lambda (m f g u v)
      (let loop ((ms (digits m)) (u u) (v v))
        (cond ((null? ms) (values u v))
              ((zero? (car ms)) (loop (cdr ms) (f u) (g u v)))
              (else (loop (cdr ms) (g u v) (f v)))))))
    (lucas-pseudo-prime? (lambda (n)
      (let loop ((a 11) (b 7))
        (let ((d (- (* a a) (* 4 b))))
          (cond ((square? d) (loop (+ a 2) (+ b 1)))
                ((not (= (gcd n (* 2 a b d)) 1)) (loop (+ a 2) (+ b 2)))
                (else (let* ((x1 (modulo (- (* a a (inverse b n)) 2) n))
                             (m (quotient (- n (jacobi d n)) 2))
                             (f (lambda (u) (modulo (- (* u u) 2) n)))
                             (g (lambda (u v) (modulo (- (* u v) x1) n))))
                        (let-values (((xm xm1) (chain m f g 2 x1)))
                          (zero? (modulo (- (* x1 xm) (* 2 xm1)) n)))))))))))
    (if (not (integer? n)) (error 'prime? "must be integer")
      (if (< n 2) #f (if (even? n) (= n 2) (if (zero? (modulo n 3)) (= n 3)
        (and (strong-pseudo-prime? n 2)
             (strong-pseudo-prime? n 3)
             (lucas-pseudo-prime? n))))))))

(define (next-prime n)
  (cond ((< n 2) 2) ((< n 3) 3)
  (else (let loop ((n (+ (if (even? n) 1 2) n)))
          (if (prime? n) n (loop (+ n 2)))))))

(define (prev-prime n)
  (cond ((< n 3) #f) ((= n 3) 2)
  (else (let loop ((n (- n (if (even? n) 1 2))))
          (if (prime? n) n (loop (- n 2)))))))

(define (nth-prime n)
  ; (define (approx-pi x) (inexact->exact (round (riemann-r x))))
  (define (approx-pi x) (inexact->exact (round
    (- (logint x) (/ (logint (sqrt x)) 2) (/ (logint (expt x 1/3)) 3)))))
  (if (< n 5) (vector-ref (vector 2 3 5 7) (- n 1))
    (let loop1 ((lo 1) (hi 2))
      (if (< (approx-pi hi) n) (loop1 hi (* hi 2))
        (let loop2 ((lo lo) (hi hi))
          (let* ((mid (quotient (+ lo hi) 2)) (mid-pi (approx-pi mid)))
            (cond ((< n mid-pi) (loop2 lo mid))
                  ((< mid-pi n) (loop2 mid hi))
                  (else (let* ((m (if (prime? mid) mid (prev-prime mid))) (k (prime-pi m)))
                          (cond ((< k n) (let loop3 ((k k) (m m))
                                  (if (= n k) m (loop3 (+ k 1) (next-prime m)))))
                                ((< n k) (let loop3 ((k k) (m m))
                                  (if (= k n) m (loop3 (- k 1) (prev-prime m)))))
                                (else m)))))))))))

(display (prime-pi 15485863)) (newline)

(display (nth-prime 1000000)) (newline)


Output:
1
2
1000000
15485863


Create a new paste based on this one


Comments: