From octave-sources-request at bevo dot che dot wisc dot edu  Tue Jun 27 03:44:58 2000
Subject: primes.m
From: pkienzle at kienzle dot powernet dot co dot uk (Paul Kienzle)
To: octave-sources at bevo dot che dot wisc dot edu
Date: Mon, 26 Jun 2000 16:13:04 +0100 (BST)


## usage: x = primes(n)
##    lists primes up to n

function x=primes(p)
  null = zeros (1, p);
  sieve = ones (1, p);          # assume everything is prime
  sieve (1) = 0;                # 1 is not prime
  sieve (4:2:p) = null (4:2:p); # multiples of two are not prime
  for i=3:2:sqrt(p)             # check only odd i
    if (sieve (i))              # if i is prime, eliminate multiples of i
      idx = i+i : i : p;
      sieve (idx) = null (idx);
    endif
  endfor
  x = find (sieve);             # primes remaining after sieve
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
-----------------------------------------------------------------------