From posting-system@google.com Tue Nov 13 18:07:39 2001
Date: Tue, 13 Nov 2001 15:07:34 -0800
From: oleg-at-pobox.com
Newsgroups: comp.lang.scheme
X-comment: added the second algorithm.
Subject: Re: arbritary precision rationals
References: <87668rm2wf.fsf@stuff.gen.nz>
Message-ID: <7eb8ac3e.0111131507.650021e5@posting.google.com>
Status: OR

In penance for the thread about the probability of two random integers
being mutually prime, I humbly submit a message only about Scheme.
This article shows a not-so-naive realization of the Eratosthenes
Sieve. The code computes a sequence of prime numbers not exceeding a
given number N, rather fast. The sophistication of the algorithm is in
avoiding wasting time and space on even numbers. This message also
shows the example of computing with exact rationals.
	
The Eratosthenes sieve algorithm is described in
        http://mathworld.wolfram.com/SieveofEratosthenes.html
	
In the naive implementation, we would use a vector v of size N+1. v[i]
is a boolean indicating the status of integer i, i=0..N. v[i] is #f if
integer i has been eliminated (dropped through the sieve). Initially v
is set to #t, with the exception v[0]=v[1]=#f (obviously 0 and 1 are
not primes). We start by locating the first i (greater than the last
found prime) so that v[i] is #t. We add v[i] to the list of found
primes, and eliminate all multiples of i. That is, we set v[2i],
v[3i], etc. to #f.  We repeat the process until there is nothing to
eliminate.
	
Since all primes except 2 are odd, it is wasteful to even consider even
numbers. Now, we will interpret v[i] as an indication of the status of
a number n(i) = (2*i)+3. Initially, all elements of the vector are set
to #t. We locate the first i in the yet-to-be-examined portion of v 
so that v[i] is #t. Such i corresponds to a number p(i)= (2*i)+3. We
add this number to the list of found primes and eliminate all its
multiples. Because we do not consider even numbers, we need to
eliminate only odd multiples, (2k+1)*(2*i+3). k=1.... That is, we need to set
	v[ ( (2k+1)*(2*i+3) - 3 ) / 2 ] to #f
or
	v[ k*p(i) + i] = #f
	
The algorithm becomes
	
(define (prime-sieve N)
  (let* ((max-index (quotient (- N 3) 2))
	 (v (make-vector (+ 1 max-index) #t)))
    ; i is the current index on the tape
    ; primes is the list of found primes, in reverse order
    (let loop ((i 0) (primes '(2)))
      (cond
       ((> i max-index) (reverse primes))
       ((vector-ref v i)
	(let ((prime (+ i i 3)))	; newly found prime
	  (do ((j (+ i prime) (+ j prime)))
	      ((> j max-index))
	    (vector-set! v j #f))
	  (loop (+ 1 i) (cons prime primes))))
       (else
	(loop (+ 1 i) primes))))))

This algorithm has no multiplications (let alone divisions) in the
main loop. 

George Kangas from the Univ. of Kansas has pointed out that if we allow
general multiplications, we can eliminate redundant assignments of
the elements of v to #f. Once we detected the prime
p(i)= (2*i)+3, we can start the elimination pass, the setting of
      v[ k*p(i) + i] = #f
starting with k=i+1 rather than k=1. That is, we start the elimination
from the element corresponding to prime^2. All composite elements
of v before that index have already been eliminated in the earlier
passes. He furthermore pointed out that once the element corresponding
to prime^2 is outside the array `v', no productive eliminations are
possible. Therefore, we can stop and scan the vector `v' for elements
still containing #t, collecting the corresponding primes.

(define (prime-sieve N)
  (let* ((max-index (max 0 (quotient (- N 3) 2)))
	 (v (make-vector (+ 1 max-index) #t)))
    (let loop ((i 0))                   ; i is the current index on the tape
      (cond
       ((vector-ref v i)
	(let* ((prime (+ i i 3))	; newly found prime
	       (p2i (+ i (* (+ i 1) prime)))) ; the index of prime^2
	  (do ((j p2i (+ j prime)))
	      ((> j max-index))
	    (vector-set! v j #f))
	  (if (> p2i max-index)  ; no point to continue elimination, so
	                         ; we collect the primes, scanning backwards
	     (let collect ((i max-index) (primes '()))
	       (cond
	        ((negative? i) (cons 2 primes))
		((vector-ref v i) (collect (- i 1) (cons (+ i i 3) primes)))
		(else (collect (- i 1) primes))))
	     (loop (+ 1 i)))))
       (else (loop (+ 1 i)))))))

	
For example, (prime-sieve 10000) gives the list of all prime numbers
not exceeding 10000. The first 20 are
> (take 20 (prime-sieve 10000))
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71)
and the last 10 are
> (reverse (take 10 (reverse (prime-sieve 10000))))
(9887 9901 9907 9923 9929 9931 9941 9949 9967 9973)
	
It took only 0.579 CPU seconds for a Petite Chez _interpreter_ to compute all
primes not exceeding 1 million (Pentium IV 2GHz). The second algorithm 
takes only 0.378 CPU seconds.
BTW, there are 78498 such primes, well in agreement with the approximation
of the prime count function, pi(1000000) = (/ 1000000 (log 1000000)) =
72382
	
As an illustration of exact rational arithmetics (implemented in
Gambit), let's compute P(K) = Prod{ (1 - 1/p[i]^2), i=1..K} where
p[i] is the i-th prime number (p[1]=2).
	
(define (approx N)
  (exact->inexact
   (apply *
	  (map 
	   (lambda (p) (- 1 (/ 1 (* p p))))
	   (prime-sieve N)))))
	
All computation except the final exact->inexact is exact.
	
> (approx 10)
.626938775510204
> (approx 100)
.6090337253995166
> (approx 200)
.6083818935932919
> (approx 1000)
.6080043073061253
> (approx 10000)
.6079330691140551
	
	
We have argued that P(K) converges to 1/Z[2] = 6/Pi^2:
> (let ((pi (* 4 (atan 1))))
    (/ 6 (* pi pi)))
.6079271018540267
	
That does seem to be the case indeed.