From octave-sources-request at bevo dot che dot wisc dot edu  Thu Jun 29 08:49:51 2000
Subject: Re: primes.m
From: pkienzle at kienzle dot powernet dot co dot uk (Paul Kienzle)
To: dirk at calvyn dot puk dot ac dot za
Cc: F dot Potorti at cnuce dot cnr dot it, octave-sources@bevo.che.wisc.edu
Date: Thu, 29 Jun 2000 03:40:18 +0100 (BST)

Care to debug the following?

  octave> primes(32)
  error: A([]) = X: X must also be an empty matrix
  error: evaluating assignment expression near line 9, column 34
  error: evaluating if command near line 8, column 7
  error: evaluating for command near line 7, column 5
  error: evaluating if command near line 2, column 3
  error: called from `primes' in file `/home/pkienzle/octave/primes.m'

From: Dirk Laurie <dirk at calvyn dot puk dot ac dot za>
> 
>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
-----------------------------------------------------------------------