[ create a new paste ] login | about

Link: http://codepad.org/49axmL16    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Aug 24:
; two sieving problems

(define (primes n)
  (let ((ps (list)) (sieve (make-vector (+ n 1) #t)))
    (do ((p 2 (+ p 1))) ((< n p) (reverse ps))
      (when (vector-ref sieve p)
        (set! ps (cons p ps))
        (do ((i (* p p) (+ i p))) ((< n i))
          (vector-set! sieve i #f))))))

(define prime?
  (let ((seed 3141592654))
    (lambda (n)
      (define (rand)
        (set! seed (modulo (+ (* 69069 seed) 1234567) 4294967296))
        (+ (quotient (* seed (- n 2)) 4294967296) 2))
      (define (expm b e m)
        (define (times x y) (modulo (* x y) m))
        (let loop ((b b) (e e) (r 1))
          (if (zero? e) r
            (loop (times b b) (quotient e 2)
                  (if (odd? e) (times b r) r)))))
      (define (spsp? n a)
        (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
            ((odd? d)
              (let ((t (expm a d n)))
                (if (or (= t 1) (= t (- n 1))) #t
                  (do ((s (- s 1) (- s 1))
                       (t (expm t 2 n) (expm t 2 n)))
                      ((or (zero? s) (= t (- n 1)))
                        (positive? s))))))))
      (if (not (integer? n)) (error 'prime? "must be integer")
        (if (< n 2) #f (if (= n 2) #t
            (do ((a (rand) (rand)) (k 10 (- k 1)))
                ((or (zero? k) (not (spsp? n a)))
                  (zero? k)))))))))

(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 (inverse x n)
  (let loop ((x (modulo x n)) (a 1))
    (cond ((zero? x) (error 'inverse "division by zero"))
          ((= x 1) a)
          (else (let ((q (- (quotient n x))))
                  (loop (+ n (* q x)) (modulo (* q a) n)))))))

(define (big-primes lo hi delta limit)
  (let* ((output (list))
         (sieve (make-vector delta #t))
         (ps (cdr (primes limit)))
         (qs (map (lambda (p) (modulo (* -1/2 (+ lo p 1)) p)) ps)))
    (let loop ((lo lo) (qs qs))
      (if (not (< lo hi)) (reverse output)
        (begin
          (do ((i 0 (+ i 1))) ((= i delta)) (vector-set! sieve i #t))
          (do ((ps ps (cdr ps)) (qs qs (cdr qs))) ((null? ps))
            (do ((j (car qs) (+ j (car ps)))) ((<= delta j))
              (vector-set! sieve j #f)))
          (do ((i 0 (+ i 1)) (t (+ lo 1) (+ t 2)))
              ((or (<= delta i) (<= hi t)))
            (if (and (vector-ref sieve i) (prime? t))
              (set! output (cons t output))))
          (loop (+ lo (* 2 delta))
                (map (lambda (p q) (modulo (- q delta) p)) ps qs)))))))

(define (primes1mod4 lo hi delta)
  (let* ((output (list))
         (sieve (make-vector delta #t))
         (ps (cdr (primes (isqrt hi))))
         (qs (map (lambda (p) (modulo (* -1 (inverse 4 p) (+ lo p 1)) p)) ps)))
    (let loop ((lo lo) (qs qs))
      (if (not (< lo hi)) (reverse output)
        (begin
          (do ((i 0 (+ i 1))) ((= i delta)) (vector-set! sieve i #t))
          (do ((ps ps (cdr ps)) (qs qs (cdr qs))) ((null? ps))
            (do ((j (car qs) (+ j (car ps)))) ((<= delta j))
              (vector-set! sieve j #f)))
          (do ((i 0 (+ i 1)) (t (+ lo 1) (+ t 4)))
              ((or (<= delta i) (<= hi t)))
            (if (vector-ref sieve i) (set! output (cons t output))))
          (loop (+ lo (* 4 delta))
                (map (lambda (p q) (modulo (- q delta) p)) ps qs)))))))

(display (length (big-primes #e1e25 (+ #e1e25 20000) 5000 250000))) (newline)
(display (length (primes1mod4 1000000 2000000 25000))) (newline)


Output:
1
2
327
35241


Create a new paste based on this one


Comments: