From octave-sources-request at bevo dot che dot wisc dot edu  Wed Jun 28 18:12:13 2000
Subject: Re: primes.m
From: pkienzle at kienzle dot powernet dot co dot uk (Paul Kienzle)
To: F dot Potorti at cnuce dot cnr dot it, pkienzle@kienzle.powernet.co.uk
Cc: octave-sources at bevo dot che dot wisc dot edu
Date: Wed, 28 Jun 2000 01:48:12 +0100 (BST)

... and fixed again :(

	- Paul

function x=primes(p)
  if (p > 2)
    len = floor((p-1)/2);         # length of the sieve
    sieve = ones (1, len);        # assume every odd number is prime
    for i=1:(sqrt(p)-1)/2         # check up to sqrt(p)
      if (sieve(i))               # if i is prime, eliminate multiples of i
      	sieve(3*i+1:2*i+1:len) = 0; # do it
      endif
    endfor
    x = [2, 1+2*find(sieve)];     # primes remaining after sieve
  elseif (p==2)
    x = 2;
  else
    x = [];
  endif
endfunction


Note that if you prefer list_primes(n), which gives the first n primes, you
can use the fact that the density of prime numbers around x is about
log (x) to make a much faster version of list_primes [as in a factor of
700 faster for n=1000] using the following:

function x=list_primes(n)
  ## density of primes at x is approx. log(x), so the nth prime
  ## is located around S_1^n[log(x)]dx = n log(n/e) + 1.  By trial
  ## and error, n log(5n) is big enough for n less than 300,000.
  x = primes ( 2 + n * log ( 5*n+1 ) ); 

  k = length(x);
  if ( k >= n )                 # if too many, keep the first n
    x = x (1:n);
  else                          # if not enough, generate the last few
    x = [ x, inf*ones (1, n-k) ];
    p = x (k);                  # start checking after the last prime
    while ( k < n )
      p = p+2;                  # skip to the next odd number
      idx = find(x < sqrt(p));
      q = p./x(idx);
      if !any(q==fix(q))
   	x(++k) = p;
      endif
    endwhile
  endif
endfunction



-----------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
Subscription information:  http://www.che.wisc.edu/octave/archive.html
-----------------------------------------------------------------------