Mercurial > hg > octave-terminal
changeset 10657:c6833d31f34e
optimize primes and isprime
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Tue, 25 May 2010 13:46:22 +0200 |
parents | ca836bcdf85e |
children | c66a4657d764 |
files | scripts/ChangeLog scripts/specfun/isprime.m scripts/specfun/primes.m |
diffstat | 3 files changed, 25 insertions(+), 22 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2010-05-25 Jaroslav Hajek <highegg@gmail.com> + + * specfun/primes.m: Use logical masks rather than numeric. + * specfun/isprime.m: Rewrite using isprime. + 2010-05-25 Jaroslav Hajek <highegg@gmail.com> * miscellaneous/unimplemented.m: Don't mention onCleanup (supported).
--- a/scripts/specfun/isprime.m +++ b/scripts/specfun/isprime.m @@ -36,20 +36,18 @@ function t = isprime (n) if (nargin == 1) - if (! isscalar (n)) - nel = numel (n); - t = zeros (size (n), "logical"); - for i = 1:nel - t(i) = isprime (n(i)); - endfor - elseif (n != fix (n) || n < 2) - t = logical (0); - elseif (n < 9) - t = all (n != [4, 6, 8]); - else - q = n./[2, 3:2:sqrt(n)]; - t = all (q != fix (q)); - endif + n = n(:); + idx = 1:numel (n); + for p = primes (sqrt (max (n(:)))) + if (isempty (idx)) + break; + endif + mask = rem (n, p) != 0; + n = n(mask); + idx = idx(mask); + endfor + t = false (size (n)); + t(idx) = true; else print_usage (); endif
--- a/scripts/specfun/primes.m +++ b/scripts/specfun/primes.m @@ -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