# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1274787982 -7200
# Node ID c6833d31f34eac6e786239812f2a79866a9c81a0
# Parent  ca836bcdf85ef103de064dce470135ab070de156
optimize primes and isprime

diff -r ca836bcdf85e -r c6833d31f34e scripts/ChangeLog
--- a/scripts/ChangeLog	Tue May 25 11:50:24 2010 +0200
+++ b/scripts/ChangeLog	Tue May 25 13:46:22 2010 +0200
@@ -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).
diff -r ca836bcdf85e -r c6833d31f34e scripts/specfun/isprime.m
--- a/scripts/specfun/isprime.m	Tue May 25 11:50:24 2010 +0200
+++ b/scripts/specfun/isprime.m	Tue May 25 13:46:22 2010 +0200
@@ -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
diff -r ca836bcdf85e -r c6833d31f34e scripts/specfun/primes.m
--- 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