diff -r 33d895260fa4 NEWS
--- a/NEWS Sun Oct 17 12:28:51 2021 +0200
+++ b/NEWS Mon Oct 18 12:26:50 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 libinterp/corefcn/__isprimelarge__.cc
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/__isprimelarge__.cc Mon Oct 18 12:26:50 2021 -0400
@@ -0,0 +1,161 @@
+////////////////////////////////////////////////////////////////////////
+//
+// Copyright (C) 1994-2021 The Octave Project Developers
+//
+// See the file COPYRIGHT.md in the top-level directory of this
+// distribution or .
+//
+// This file is part of Octave.
+//
+// Octave is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// Octave is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with Octave; see the file COPYING. If not, see
+// .
+//
+////////////////////////////////////////////////////////////////////////
+
+#include "defun.h"
+#include "error.h"
+#include "ovl.h"
+
+OCTAVE_NAMESPACE_BEGIN
+
+// This function implements the Schrage technique for modular multiplication.
+// The returned value is equivalent to "mod (a*b, modulus)"
+// but calculated without overflow.
+uint64_t safemultiply (uint64_t a, uint64_t b, uint64_t modulus)
+{
+ uint64_t q = modulus / a;
+ uint64_t r = modulus - q*a;
+ uint64_t term1 = a * (b % q);
+ uint64_t term2 = (r < q) ? r * (b/q) : safemultiply (r, b / q, modulus);
+ return (term1 > term2) ? (term1 - term2) : (term1 + modulus - term2);
+}
+
+// This function returns "mod (a^b, modulus)"
+// but calculated without overflow.
+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;
+}
+
+// This function implements a single round of Miller-Rabin primality testing.
+// Returns false if composite, true if pseudoprime for this divisor.
+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;
+ uint64_t 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)
+
+%!error (__isprimelarge__ ([1 2 3]))
+%!error (__isprimelarge__ (36))
+%!error (__isprimelarge__ (-1))
+*/
+
+OCTAVE_NAMESPACE_END
diff -r 33d895260fa4 libinterp/corefcn/module.mk
--- a/libinterp/corefcn/module.mk Sun Oct 17 12:28:51 2021 +0200
+++ b/libinterp/corefcn/module.mk Mon Oct 18 12:26:50 2021 -0400
@@ -119,6 +119,7 @@
%reldir%/__gammainc__.cc \
%reldir%/__ichol__.cc \
%reldir%/__ilu__.cc \
+ %reldir%/__isprimelarge__.cc \
%reldir%/__lin_interpn__.cc \
%reldir%/__magick_read__.cc \
%reldir%/__pchip_deriv__.cc \
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 12:26:50 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 12:26:50 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])