--- a/scripts/specfun/isprime.m	Wed Aug 13 21:29:53 2014 +0200
+++ b/scripts/specfun/isprime.m	Mon Aug 25 21:30:23 2014 -0700
@@ -22,9 +22,21 @@
 ## Return a logical array which is true where the elements of @var{x} are
 ## prime numbers and false where they are not.
 ##
-## @code{isprime} is appropriate if the maximum value in @var{x} is not too
-## large (< 1e15).  For larger values special purpose factorization code
-## should be used.
+## A prime number is conventionally defined as a positive integer (1, 2, 3,
+## @dots{}) which is divisible only by itself and 1.  Octave extends this
+## definition to include both negative integers and complex values.  A
+## negative integer is prime if its positive counterpart is prime.  This is
+## equivalent to @code{isprime (abs (x))}.
+## 
+## If @code{class (@var{x})} is complex, then primality is tested in the domain
+## of Gaussian integers (@url{http://en.wikipedia.org/wiki/Gaussian_integer}).
+## Some non-complex integers are prime in the ordinary sense, but not in the
+## domain of Gaussian integers.  For example, @math{5 = (1+2i)*(1-2i)} shows
+## that 5 is not prime because it has a factor other than itself and 1.
+## Exercise caution when testing complex and real values together in the same
+## matrix.
+##
+## Examples:
 ##
 ## @example
 ## @group
@@ -33,17 +45,19 @@
 ## @end group
 ## @end example
 ##
-## If @code{class (@var{x})} is complex, then primality will be tested
-## in the domain of Gaussian integers. This means that some real prime
-## numbers will fail the primality test in this case. For example, 5 =
-## (1+2i)*(1-2i) and 2 = i*(1-i)^2.
-##
 ## @example
 ## @group
 ## isprime ([i, 2, 3, 5])
 ##     @result{} [0, 0, 1, 0]
 ## @end group
 ## @end example
+##
+## Programming Note: @code{isprime} is appropriate if the maximum value in
+## @var{x} is not too large (< 1e15).  For larger values special purpose
+## factorization code should be used.
+##
+## Compatibility Note: @var{matlab} does not extend the definition of prime
+## numbers and will produce an error if given negative or complex inputs.
 ## @seealso{primes, factor, gcd, lcm}
 ## @end deftypefn
 
@@ -65,19 +79,23 @@
     return;
   endif
 
-  ## Handle negative entries
-  x = abs (x);
+  ## Code strategy is to build a table with the list of possible primes
+  ## and then quickly compare entries in x with the table of primes using
+  ## lookup().  The table size is limited to save memory and computation
+  ## time during its creation.  All entries larger than the maximum in the
+  ## table are checked by straightforward division.
 
+  x = abs (x);              # handle negative entries
   maxn = max (x(:));
   ## generate prime table of suitable length.
-  maxp = min (maxn, max (sqrt (maxn), 1e7)); # FIXME: threshold not optimized.
+  ## 1e7 threshold requires ~0.15 seconds of computation, 1e8 requires 1.8.
+  maxp = min (maxn, max (sqrt (maxn), 1e7));
   pr = primes (maxp);
-  ## quick search for table matches.
-  t = lookup (pr, x, "b");
-  ## take the rest.
+  t = lookup (pr, x, "b");  # quick search for table matches.
+
+  ## process any remaining large entries
   m = x(x > maxp);
   if (! isempty (m))
-    ## there are still possible primes. filter them out by division.
     if (maxn <= intmax ("uint32"))
       m = uint32 (m);
     elseif (maxn <= intmax ("uint64"))
@@ -85,16 +103,25 @@
     else
       warning ("isprime: X contains integers too large to be tested");
     endif
+
+    ## Start by dividing through by the small primes until the remaining
+    ## list of entries is small (and most likely prime themselves).
     pr = cast (pr(pr <= sqrt (maxn)), class (m));
     for p = pr
       m = m(rem (m, p) != 0);
-      if (length (m) < length (pr) / 10)
+      if (numel (m) < numel (pr) / 10)
         break;
       endif
     endfor
+
+    ## Check the remaining list of possible primes against the
+    ## remaining prime factors which were not tested in the for loop.
+    ## This is just an optimization to use arrayfun over for loo
     pr = pr(pr > p);
     mm = arrayfun (@(x) all (rem (x, pr)), m);
     m = m(mm);
+
+    ## Add any remaining entries, which are truly prime, to the results.
     if (! isempty (m))
       m = cast (sort (m), class (x));
       t |= lookup (m, x, "b");
@@ -104,7 +131,7 @@
 endfunction
 
 function t = isgaussianprime (z)
-  ## Assumed prime unless proven otherwise
+  ## Assume prime unless proven otherwise
   t = true (size (z));
 
   x = real (z);
@@ -118,15 +145,13 @@
   t(xidx) &= isprime (x(xidx));
   t(yidx) &= isprime (y(yidx));
 
-  ## Otherwise, prime if x^2 + y^2 is prime and NOT 3 mod 4.
-  zabs = x.^2 + y.^2;
-  zidx = mod (zabs, 4) != 3 & !xidx & !yidx;
-
-  t(zidx) &= isprime (zabs (zidx));
+  ## Otherwise, prime if x^2 + y^2 is prime
+  zidx = ! (xidx | yidx);          # Skip entries that were already evaluated
+  zabs = x(zidx).^2 + y(zidx).^2;
+  t(zidx) &= isprime (zabs);
 endfunction
 
 
-
 %!assert (isprime (3), true)
 %!assert (isprime (4), false)
 %!assert (isprime (5i), false)
@@ -144,3 +169,4 @@
 %!error isprime (1, 2)
 %!error <X contains non-integer entries> isprime (0.5i)
 %!error <X contains non-integer entries> isprime (0.5)
+