Mercurial > octave-antonio
diff scripts/specfun/primes.m @ 10657:c6833d31f34e
optimize primes and isprime
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Tue, 25 May 2010 13:46:22 +0200 |
parents | 95c3e38098bf |
children | 3140cb7a05a1 |
line wrap: on
line diff
--- a/scripts/specfun/primes.m Tue May 25 11:50:24 2010 +0200 +++ b/scripts/specfun/primes.m Tue May 25 13:46:22 2010 +0200 @@ -58,26 +58,26 @@ ## p = 3e6. Hardly worthwhile, but Dirk reports better numbers. 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 + sievem = true (1, lenm); # assume every number of form 6n-1 is prime + sievep = true (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; + sievem(7*i-1:6*i-1:lenm) = false; + sievep(5*i-1:6*i-1:lenp) = false; 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; + sievep(7*i+1:6*i+1:lenp) = false; + sievem(5*i+1:6*i+1:lenm) = false; endif endfor x = sort([2, 3, 6*find(sievem)-1, 6*find(sievep)+1]); elseif (p > 352) # nothing magical about 352; must be >2 len = floor ((p-1)/2); # length of the sieve - sieve = ones (1, len); # assume every odd number is prime + sieve = true (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 + sieve(3*i+1:2*i+1:len) = false; # do it endif endfor x = [2, 1+2*find(sieve)]; # primes remaining after sieve