[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Sep 27:
; sieve of sundaram

(define (primes n)
  (let* ((m (quotient n 2)) (pv (make-vector (+ m 1) #t)))
    (do ((i 1 (+ i 1))) ((< (quotient m 4) i))
      (do ((j i (+ j 1))) ((< (quotient (- m i) (+ i i 1)) j))
        (vector-set! pv (+ i j (* 2 i j)) #f)))
    (let loop ((i 1) (ps (list 2)))
      (cond ((= i m) (reverse ps))
            ((vector-ref pv i) (loop (+ i 1) (cons (+ i i 1) ps)))
            (else (loop (+ i 1) ps))))))

(display "Sundaram:     ") (time (display (length (primes 2000000))) (display " "))

(define (primes limit)
  (define (exact x) (inexact->exact (floor x)))
  (let ((sieve (make-vector (+ (quotient limit 2) (modulo limit 2)) #f))
        (primes (list 3 2)))
    (define (flip! m) (vector-set! sieve m (not (vector-ref sieve m))))

    (let ((x-max (exact (sqrt (/ (- limit 1) 4)))) (x2 0))
      (do ((xd 4 (+ xd 8))) ((<= (+ (* x-max 8) 2) xd))
        (set! x2 (+ x2 xd))
        (let* ((y-max (exact (sqrt (- limit x2))))
               (n (+ x2 (* y-max y-max)))
               (n-diff (+ y-max y-max -1)))
          (when (even? n) (set! n (- n n-diff)) (set! n-diff (- n-diff 2)))
          (do ((d (* (- n-diff 1) 2) (- d 8))) ((<= d -1))
            (when (member (modulo n 12) (list 1 5)) (flip! (quotient n 2)))
            (set! n (- n d))))))

    (let ((x-max (exact (sqrt (/ (- limit 1) 3)))) (x2 0))
      (do ((xd 3 (+ xd 6))) ((<= (+ (* x-max 6) 2) xd))
        (set! x2 (+ x2 xd))
        (let* ((y-max (exact (sqrt (- limit x2))))
               (n (+ x2 (* y-max y-max)))
               (n-diff (+ y-max y-max -1)))
          (when (even? n) (set! n (- n n-diff)) (set! n-diff (- n-diff 2)))
          (do ((d (* (- n-diff 1) 2) (- d 8))) ((<= d -1))
            (when (= (modulo n 12) 7) (flip! (quotient n 2)))
            (set! n (- n d))))))

    (let ((x-max (exact (/ (+ (sqrt (- 4 (* (- 1 limit) 8))) 2) 4)))
          (y-min -1) (x2 0) (xd 3))
      (do ((x 1 (+ x 1))) ((<= (+ x-max 1) x))
        (set! x2 (+ x2 xd)) (set! xd (+ xd 6))
        (when (<= limit x2)
          (set! y-min (* (- (* (- (inexact->exact (ceiling (sqrt (- x2 limit)))) 1) 2) 2) 2)))
        (let ((n (- (* (+ (* x x) x) 2) 1))
              (n-diff (* (- (* (- x 1) 2) 2) 2)))
          (do ((d n-diff (- d 8))) ((<= d y-min))
            (when (= (modulo n 12) 11) (flip! (quotient n 2)))
            (set! n (+ n d))))))

    (do ((n 2 (+ n 1))) ((<= (quotient (+ (exact (sqrt limit)) 1) 2) n))
      (when (vector-ref sieve n)
        (let* ((p (+ n n 1)) (p2 (* p p)))
          (set! primes (cons p primes))
          (do ((k p2 (+ k (+ p2 p2)))) ((<= limit k))
            (vector-set! sieve (quotient k 2) #f)))))

    (do ((p (+ (exact (sqrt limit)) 1) (+ p 2))) ((<= limit p))
      (when (vector-ref sieve (quotient p 2))
        (set! primes (cons p primes))))

    (reverse primes)))

(display "Atkin:        ") (time (display (length (primes 2000000))) (display " "))

(define (primes . args) ; (primes [lo] hi) inclusive at both ends
  (let* ((lo (if (null? (cdr args)) 0 (car args)))
         (hi (if (null? (cdr args)) (car args) (cadr args))))
    (cond ((and (<= lo 100000) (<= hi 1000000)) ; simple sieve
           (let* ((max-index (quotient (- hi 3) 2))
                  (v (make-vector (+ max-index 1) #t)))
             (let loop ((i 0) (ps (list 2)))
               (let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
                 (cond ((< hi (* p p))
                        (let loop ((j i) (ps ps))
                          (cond ((< max-index j)
                                 (let loop ((ps (reverse ps)))
                                   (if (<= lo (car ps)) ps
                                     (loop (cdr 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))
                          (when (<= j max-index)
                            (vector-set! v j #f) (loop (+ j p))))
                        (loop (+ i 1) (cons p ps)))
                       (else (loop (+ i 1) ps)))))))
          ((< lo (sqrt hi))
           (let* ((r (inexact->exact (ceiling (sqrt hi))))
                  (r (if (even? r) r (+ r 1))))
             (append (primes lo r) (primes r hi))))
          (else ; segmented sieve
           (let* ((lo (if (even? lo) lo (- lo 1)))
                  (b 50000) (bs (make-vector b #t))
                  (r (inexact->exact (ceiling (sqrt hi))))
                  (ps (cdr (primes r)))
                  (qs (map (lambda (p)
                             (modulo (* -1/2 (+ lo 1 p)) p)) ps))
                  (zs (list)) (z (lambda (p) (set! zs (cons p zs)))))
             (do ((t lo (+ t b b))
                  (qs qs (map (lambda (p q) (modulo (- q b) p))
                              ps qs)))
                 ((<= hi t)
                   (let loop ((zs zs))
                     (if (<= (car zs) hi) (reverse zs)
                       (loop (cdr zs)))))
               (do ((i 0 (+ i 1))) ((= i b)) (vector-set! bs i #t))
               (do ((ps ps (cdr ps)) (qs qs (cdr qs))) ((null? qs))
                 (do ((j (car qs) (+ j (car ps)))) ((<= b j))
                   (vector-set! bs j #f)))
               (do ((j 0 (+ j 1))) ((= j b))
                 (if (vector-ref bs j) (z (+ t j j 1))))))))))

(display "Segmented:    ") (time (display (length (primes 2000000))) (display " "))

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

(display "Eratosthenes: ") (time (display (length (primes 2000000))) (display " "))


Output:
1
2
3
4
Sundaram:     148933 cpu time: 720 real time: 3564 gc time: 0
Atkin:        148933 cpu time: 300 real time: 1506 gc time: 10
Segmented:    148933 cpu time: 210 real time: 959 gc time: 0
Eratosthenes: 148933 cpu time: 170 real time: 855 gc time: 0


Create a new paste based on this one


Comments: