From octave-sources-request at bevo dot che dot wisc dot edu  Thu Jun 29 06:31:13 2000
Subject: Re: primes.m
From: Dirk Laurie <dirk at calvyn dot puk dot ac dot za>
To: Paul Kienzle <pkienzle at kienzle dot powernet dot co dot uk>
Cc: F dot Potorti at cnuce dot cnr dot it, octave-sources@bevo.che.wisc.edu
Date: Thu, 29 Jun 2000 13:30:59 +0200

Paul Kienzle skryf:
> ... 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
> 
A further optimization is possible at the expense of complicating
the code.  The main point is a 33% memory reduction, but there is
also a system-dependent speed improvement.

On my system the speed difference is slight until we come close to
the memory limit, when swapping takes over.  E.g. primes(2e7)
takes 95 seconds real time on the above code whereas the code below 
takes 24 seconds real time.  At p=1e7 the speed improvement is only
about 20%.

The trick can be iterated (separate sieves for 30k+-1, 30k+-7, 
30k+-11, 30k+-13, etc) but there is a law of diminishing returns.

-- Dirk Laurie

---------------------------------------------------------------------

function x=primes(p)
  if (p > 3)
    lenm = floor((p+1)/6);        # length of the 6n-1 sieve
    lenp = floor((p-1)/6);        # length of the 6n+1 sieve
    sievem = ones (1, lenm);      # assume every number of form 6n-1 is prime
    sievep = ones (1, lenp);      # assume every number of form 6n+1 is prime
    for i=1:(sqrt(p)+1)/6         # check up to sqrt(p)
      if (sievem(i))              # if i is prime, eliminate multiples of i
        sievem(7*i-1:6*i-1:lenm) = 0; 
        sievep(5*i-1:6*i-1:lenp) = 0;  
      endif                       # if i is prime, eliminate multiples of i
      if (sievep(i))
        sievep(7*i+1:6*i+1:lenp) = 0; 
        sievem(5*i+1:6*i+1:lenm) = 0;  
      endif
    endfor
    x = sort([2, 3, 6*find(sievem)-1 6*find(sievep)+1]);     
  else
    x=[2 3]; x=x(find(x<=p));
  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
-----------------------------------------------------------------------