diff -r ba4aebad10d1 libinterp/corefcn/data.cc
--- a/libinterp/corefcn/data.cc	Wed Oct 13 10:36:20 2021 -0700
+++ b/libinterp/corefcn/data.cc	Sat Oct 16 10:14:26 2021 -0400
@@ -8471,4 +8471,151 @@
 %!error <input was not valid base64> __base64_decode_bytes__ (1)
 */
 
+/*
+  The next several functions are utility functions used by __isprimelarge__
+  which will be defined below. That function __isprimelarge__ in turn is
+  called by isprime.m. The end user should call isprime not __isprimelarge__.
+*/
+
+uint64_t safeadd (uint64_t a, uint64_t b, uint64_t modulus)
+{
+  static uint64_t IMAX = (1ULL << 63) - 1 + (1ULL << 63);  // 2^64-1
+  uint64_t retval;
+  if (b < IMAX-a)
+    retval = a+b;
+  else
+    retval = (IMAX - modulus) + b - (IMAX - a);
+
+  if (retval > modulus)
+    retval -= modulus;
+
+  return retval;
+}
+
+uint64_t safemultiply (uint64_t a, uint64_t b, uint64_t modulus)
+{
+  uint64_t retval = 0;
+  if (a < modulus / b)
+    retval = a*b;
+  else
+    {
+      while (b > 0)
+        {
+          if (b & 1)
+            retval = safeadd (retval, a, modulus);
+          b >>= 1;
+          a = safeadd (a, a, modulus);
+        }
+    }
+  return retval;
+}
+
+uint64_t safepower (uint64_t a, uint64_t b, uint64_t modulus)
+{
+  uint64_t retval = 1;
+  while (b > 0)
+    {
+      if (b & 1)
+        retval = safemultiply (retval, a, modulus);
+      b >>= 1;
+      a = safemultiply (a, a, modulus);
+    }
+  return retval;
+}
+
+bool millerrabin (uint64_t div, uint64_t d, uint64_t r, uint64_t n)
+{
+  uint64_t x = safepower (div, d, n);
+  if (x == 1 || x == n-1)
+    return true;
+
+  for (uint64_t j = 1; j < r; j++)
+    {
+      x = safemultiply (x, x, n);
+      if (x == n-1)
+        return true;
+    }
+  return false;
+}
+
+DEFUN (__isprimelarge__, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{x} =} __isprimelarge__ (@var{n})
+Use the Miller-Rabin test to find out whether the scalar N is prime or
+composite. The input N is required to be a scalar 64-bit integer. You
+probably want to call isprime(N) instead of directly calling this function.
+
+@seealso{isprime, factor}
+@end deftypefn */)
+{
+  int nargin = args.length ();
+  if (nargin != 1)
+    print_usage ();
+
+  // This function is intended for internal use by isprime.m,
+  // so the following error handling should not be necessary. But it is
+  // probably good practice for any curious users calling it directly.
+  uint64_t n = args(0).xuint64_scalar_value
+  ("__isprimelarge__: unable to convert input into uint64_t scalar. Call isprime() instead.");
+
+  if (n <= 37)  // n is too small to test reliably with this code.
+    // This case happens only if the user is directly calling this function.
+    // Calling it from isprime.m will prevent this case before it gets here.
+    error ("__isprimelarge__: input argument too small. Call isprime() instead.");
+
+  // rewrite n as d * 2^r + 1, for odd d
+  uint64_t d = n-1, r = 0;
+  while ( ! (d & 1))
+  {
+    d >>= 1;
+    r++;
+  }
+
+  // Miller-Rabin test with the first 12 primes.
+  // If the number passes all 12 tests, then it is prime.
+  // If it fails any, then it is composite.
+  // The first 12 primes suffice to test all 64-bit integers.
+  if ( ! millerrabin ( 2, d, r, n))   return ovl (false);
+  if ( ! millerrabin ( 3, d, r, n))   return ovl (false);
+  if ( ! millerrabin ( 5, d, r, n))   return ovl (false);
+  if ( ! millerrabin ( 7, d, r, n))   return ovl (false);
+  if ( ! millerrabin (11, d, r, n))   return ovl (false);
+  if ( ! millerrabin (13, d, r, n))   return ovl (false);
+  if ( ! millerrabin (17, d, r, n))   return ovl (false);
+  if ( ! millerrabin (19, d, r, n))   return ovl (false);
+  if ( ! millerrabin (23, d, r, n))   return ovl (false);
+  if ( ! millerrabin (29, d, r, n))   return ovl (false);
+  if ( ! millerrabin (31, d, r, n))   return ovl (false);
+  if ( ! millerrabin (37, d, r, n))   return ovl (false);
+  return ovl (true); // if all the way here, then it is prime
+
+/*
+  Mathematical references for the curious as to why we need only
+  the 12 smallest primes for testing all 64-bit numbers:
+  (1) https:oeis.org/A014233
+      Comment: a(12) > 2^64.  Hence the primality of numbers < 2^64 can be
+      determined by asserting strong pseudoprimality to all prime bases <= 37
+      (=prime(12)). Testing to prime bases <=31 does not suffice,
+      as a(11) < 2^64 and a(11) is a strong pseudoprime
+      to all prime bases <= 31 (=prime(11)). - Joerg Arndt, Jul 04 2012
+  (2) https:arxiv.org/abs/1509.00864
+      Strong Pseudoprimes to Twelve Prime Bases
+      Jonathan P. Sorenson, Jonathan Webster
+  
+  In addition, a source listed here: https://miller-rabin.appspot.com/
+  reports that all 64-bit numbers can be covered with only 7 divisors,
+  namely 2, 325, 9375, 28178, 450775, 9780504, and 1795265022.
+  I couldn't find a peer-reviewed article to back it up though,
+  so I am using the 12 primes <= 37 in this code.
+*/
+
+}
+
+/*
+%!assert (__isprimelarge__ (uint64 (12345)), false) 
+%!assert (__isprimelarge__ (uint64 (2147483647)), true) 
+%!assert (__isprimelarge__ (uint64 (2305843009213693951)), true) 
+%!assert (__isprimelarge__ (uint64 (18446744073709551557)), true) 
+*/
+
 OCTAVE_NAMESPACE_END
diff -r ba4aebad10d1 scripts/specfun/isprime.m
--- a/scripts/specfun/isprime.m	Wed Oct 13 10:36:20 2021 -0700
+++ b/scripts/specfun/isprime.m	Sat Oct 16 10:14:26 2021 -0400
@@ -58,10 +58,6 @@
 ## @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: @sc{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}
@@ -85,54 +81,22 @@
     return;
   endif
 
-  ## Code strategy is to build a table with the list of possible primes
+  ## Code strategy is to build a table with the list of small 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.
+  ## lookup(), then call Miller-Rabin for 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;              # some suitable transition threshold, can be small
+  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 with Miller-Rabin
+  for i = find (x>maxp)(:)'  # might as well go sequentially because __isprimelarge__ takes only microseconds.
+    t(i) = __isprimelarge__ (x(i));
+  endfor
 
 endfunction
 
@@ -162,6 +126,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])
