[Octave-bug-tracker] [bug #61312] Extending isprime() with Miller-Rabin test

[Arun Giridhar]

Follow-up Comment #2, bug #61312 (project octave):

Much better patch attached (Patch 2).

Due to problems with idivide() for large inputs (see Bug #61319), which
affected the reliability of Patch 1, this patch moves all the Miller-Rabin
code to C++, which is then called by scripts/specfun/isprime.m for any inputs
that are larger than a modest threshold.

Note to devs: The new C++ code has been added to the end of
libinterp/corefcn/data.cc. If you need to move it to a separate new file
(analogous to gcd.cc) please do so.

Benchmark test:

# create input vector over several orders of magnitude
n = uint64(2) .^ (1:63) - 1;        # some are Mersenne primes, some are not
n(end+1) = intmax("uint64");        # == 2^64-1
n(end+1) = intmax("uint64") - 58;   # == 2^64 - 59 == largest 64-bit prime

# actual test
tic, x = isprime(n); toc

assert (all (find(x) == [2 3 5 7 13 17 19 31 61 65]))     # check correctness

# prime n == 3  7  31 127 8191 131071 524287 2147483647 2305843009213693951
18446744073709551557
# all are Mersenne primes except the last, which is the largest 64-bit prime

whos


Benchmark results:
Unpatched Octave:  65.5 seconds
Patched Octave:    21.8 milliseconds

Patch is on average about 3000X faster than unpatched over the whole range of
64-bit integers, increasing to about 30,000X for inputs around 10^19. It
typically reduces 20-second runtimes for the upper end to about 0.7
milliseconds.


    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?61312>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/