From octave-sources-request at bevo dot che dot wisc dot edu  Tue Jun 27 08:45:05 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: Tue, 27 Jun 2000 03:37:55 +0100 (BST)

Your optimizations, fixed.

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)/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


Paul Kienzle
pkienzle at kienzle dot powernet dot co dot uk



-----------------------------------------------------------------------
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
-----------------------------------------------------------------------