diff -r 33d895260fa4 NEWS
--- a/NEWS	Sun Oct 17 12:28:51 2021 +0200
+++ b/NEWS	Mon Oct 18 10:26:23 2021 -0400
@@ -82,7 +82,10 @@
 (A*x - b)`.  Previously, Octave computed a minimum-norm solution.
 
 - The `factor` function has been overhauled for speed.  For large
-composite numbers > 1e14, it can be up to 10,000 times faster.
+inputs > 1e14, it can be up to 10,000 times faster.
+
+- The `isprime` function uses a new primality testing algorithm
+that is up to 50,000 times faster for all inputs < 2^64.
 
 - The `betainc` function now calculates an exact output for the
 important special cases where a or b are 1.
diff -r 33d895260fa4 scripts/specfun/factor.m
--- a/scripts/specfun/factor.m	Sun Oct 17 12:28:51 2021 +0200
+++ b/scripts/specfun/factor.m	Mon Oct 18 10:26:23 2021 -0400
@@ -51,13 +51,16 @@
     error ("factor: Q must be a real non-negative integer");
   endif
 
-  ## Special case of no primes less than sqrt (q).
-  if (q < 4)
+  ## Special case if q is prime, because isprime() is now much faster than factor().
+  ## This also absorbs the case of q < 4, where there are no primes less than sqrt(q).
+  if (q < 4 || isprime (q))
     pf = q;
     n = 1;
     return;
   endif
 
+  ## If we are here, then q is composite.
+
   cls = class (q);  # store class
   if (isfloat (q) && q > flintmax (q))
     error ("factor: Q too large to factor (> flintmax)");
diff -r 33d895260fa4 scripts/specfun/isprime.m
--- a/scripts/specfun/isprime.m	Sun Oct 17 12:28:51 2021 +0200
+++ b/scripts/specfun/isprime.m	Mon Oct 18 10:26:23 2021 -0400
@@ -58,9 +58,8 @@
 ## @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.
+## Programming Note: @code{isprime} is suitable for all @var{x}
+## in the range abs(@var{x}) < 2^64.
 ##
 ## Compatibility Note: @sc{matlab} does not extend the definition of prime
 ## numbers and will produce an error if given negative or complex inputs.
@@ -85,54 +84,21 @@
     return;
   endif
 
-  ## 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.
+  ## Code strategy is to quickly compare entries in x with small primes
+  ## using lookup(), then call Miller-Rabin for all bigger numbers.
 
   x = abs (x);              # handle negative entries
-  maxn = max (x(:));
+
   ## generate prime table of suitable length.
-  ## 1e7 threshold requires ~0.15 seconds of computation, 1e8 requires 1.8.
-  maxp = min (maxn, max (sqrt (maxn), 1e7));
-  pr = primes (maxp);
+  maxp = 37;  # this needs to be at least 37,
+              # because of the method used by __isprimelarge__ below.
+  pr = [2 3 5 7 11 13 17 19 23 29 31 37]; # list of primes up to maxp
   t = lookup (pr, x, "b");  # quick search for table matches.
 
-  ## process any remaining large entries
-  m = x(x > maxp);
-  if (! isempty (m))
-    if (maxn <= intmax ("uint32"))
-      m = uint32 (m);
-    elseif (maxn <= intmax ("uint64"))
-      m = uint64 (m);
-    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 (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");
-    endif
-  endif
+  ## process remaining entries sequentially with Miller-Rabin
+  for i = find (x>maxp)(:)'
+    t(i) = __isprimelarge__ (x(i));
+  endfor
 
 endfunction
 
@@ -162,6 +128,7 @@
 
 %!assert (isprime (3), true)
 %!assert (isprime (4), false)
+%!assert (isprime (uint64 (18446744073709551557)), true)
 %!assert (isprime (5i), false)
 %!assert (isprime (7i), true)
 %!assert (isprime ([1+2i, (2+3i)*(-1+2i)]), [true, false])
