bugGNU Octave - Bugs: bug #61312, Extending isprime() with...


bug #61312: Extending isprime() with Miller-Rabin test

Submitted by:  Arun Giridhar <arungiridhar>
Submitted on:  Fri 08 Oct 2021 10:43:58 PM UTC  
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Performance
Status:  Patch Submitted Assigned to:  None
Originator Name:  Open/Closed:  Open
Release:  dev Operating System:  Any

Add a New Comment (Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


( Jump to the original submission)

Mon 18 Oct 2021 10:30:22 PM UTC, comment #16: 

Patch 11 attached with the following changes:

1. Tuned value of the transition in isprime.m from division to Miller-Rabin. This had to be empirical, so I've added a FIXME to that part of isprime.m explaining how. The tuned value is 380e9 == 3.8e11. I found the casting from double to uint32 etc was slower than simply using double, so I've removed the casts for m and x. Please check for performance regressions.

2. Edited NEWS in view of the performance differences between old and new for small and big.

3. Optimizations to Schrage based on empirical testing. Alternatives tried before converging: the code in comment #14, various combinations and sequences for the opening code (if-else vs ternary operators), and whether to swap a and b or whether to call the function with a and b already ordered. The differences between them were minor though, compared to earlier patches.

Patch 11 attached.

(file #52122)

Arun Giridhar <arungiridhar>
Mon 18 Oct 2021 08:51:28 PM UTC, comment #15: 

I found a bigger problem with the patch. Let me work on that first before changing the Schrage part.

Try this stress test for 10 million numbers:

n = 4294967295 - (0:9999999); tic; isprime(n); toc

With unpatched Octave (straight division), it takes 26.4 seconds.
With the patched Octave (Miller-Rabin), it takes 69.9 seconds.

Clearly there needs to be some tuning of the threshold so that numbers smaller than say 2^32 or 10^10 are tested by straight division. But raising the prime threshold above 37 causes the Miller-Rabin performance for large numbers to suffer, so it may be prudent to split the inputs into two parts, those below say 2^32 and those above, like in Patch 1 all the way at the bottom. I don't know the better answer yet.

For inputs that include both numbers less than 2^32 and numbers larger than say 2^53, the Miller Rabin test is very much faster. For numbers near 10^19, it's much faster. But the threshold cannot be as low as 37 evidently.

Arun Giridhar <arungiridhar>
Mon 18 Oct 2021 05:59:58 PM UTC, comment #14: 

I just tested a bit, and it seems that the Schrage recursion often (always?) ends up with b equal to 0, but r>=q, thus going further in the recursion, where obviously it could already return. So I seem to get about a factor of two improvement if in the body of savemultiply I define

  if (a < 2) return a * b;
  uint64_t q = modulus / a;
  uint64_t r = modulus - q*a;
  uint64_t term1 = a * (b % q);
  uint64_t d = b / q;
  uint64_t term2 = (r < q) ? r * d : r < d ? safemultiply (r, d, modulus) : safemultiply (d, r, modulus);
  return (term1 > term2) ? (term1 - term2) : (term1 + modulus - term2);

You see that additionally I indeed call safemultiply (at least in the recursion) so that the first factor is the smaller as I proposed below -- this seems to give also about 10%, but perhaps it would be worth to test this further.

Michael Leitner <mleitner>
Mon 18 Oct 2021 04:59:42 PM UTC, comment #13: 

Apologies for missing a really obvious check from the code. ("if input is even, return false immediately") which makes a 40% difference to the speed. Please use this latest attached patch that rectifies the oversight.

Patch 10 attached.

(file #52120)

Arun Giridhar <arungiridhar>
Mon 18 Oct 2021 04:31:10 PM UTC, comment #12: 

Thank you for the guidance in adding it to the build system. After pruning the headers, it turned out it needed only three headers (defun, error and ovl). Added a few more BISTs. Make check passes.

Patch 9 attached.

(file #52119)

Arun Giridhar <arungiridhar>
Mon 18 Oct 2021 03:25:18 PM UTC, comment #11: 

I agree with Kai. Always good work!

1. Sound reasonable.

2. Looks good to me.

3. Adding the file to `COREFCN_SRC` in `libinterp/corefcn/module.mk` should be enough. (I believe we try to sort that list in lexical order.)

4. It would be nice if you could weed out the headers and only include those that are actually needed.

Markus Mützel <mmuetzel>
Project Member
Mon 18 Oct 2021 02:29:57 PM UTC, comment #10: 

@Kai: Thank you. Here are my thoughts:

1. I think isprime.m should continue to exist instead of being replaced completely with isprime.cc, because it handles the cases where the input is a vector or matrix much more easily than reimplementing that part in C++. It also acts as a good front-end for negative and complex inputs which won't be sped up much in C++. Only the core primality calculation has been moved to C++.

2. I have moved all the newly added C++ functions from libinterp/corefcn/data.cc to a new file called libinterp/corefcn/__isprimelarge__.cc but I have not added it to the Mercurial repository. Is it only "hg add [filename]" or are there other steps as well?

3. I also need guidance on how to add the new file to Octave's build system. I expect it should be added to Makefile.in but I see that gcd.cc for example is mentioned in Makefile.in in multiple places. If you can make the changes to Makefile.in this time I can use it as a template for my future patches. Does it also need to be added to libinterp/corefcn/module.mk?

4. I copied over all the header files from data.cc into the new _isprimelarge_.cc. I expect many of them can be removed, Once I learn from you how to add it to the build system, I will prune the list of included files.


1. New file _isprimelarge_.cc to be added to libinterp/corefcn/ and to the build system with Kai's guidance.

2. Updated patch (Patch 8) that undoes the additions to data.cc from earlier patches, but retains the changes to isprime.m, factor.m and NEWS. This patch won't be functional until the new _isprimelarge_.cc is added to the build system.

(file #52117, file #52118)

Arun Giridhar <arungiridhar>
Mon 18 Oct 2021 02:16:50 AM UTC, comment #9: 

@arungiridhar, recently I feel excited, when I see you submitted a new patch 😉  Same this time:

>> tic; n = n(isprime(n)); toc      # patch 7
Elapsed time is 0.0182581 seconds.
>> tic; n = n(isprime2(n)); toc     # Octave 7
Elapsed time is 128.464 seconds.

Octave's current version has indeed an enormous memory footprint for large inputs.

Your changes look very good to me.

As lots of code was moved from Octave to C++, I am almost tempted to suggest a new clean and complete C++ re-implementation in a new file `libinterp/corefun/isprime.cc` (like `gcd.cc` in the same directory)?

Are you interested in further work into this direction?

Otherwise, I can suggest a patch (in a few days) [<-- dangerous words in this bug tracker, that might turn out to months sometimes 😇].

Kai Torben Ohlhus <siko1056>
Project Member
Sun 17 Oct 2021 06:03:24 PM UTC, comment #8: 

Patch 7 attached. Only edits to documentation this time, not code.

New stress test / benchmark and results:

n = intmax("uint64") - (0:999);  ## 1000 largest 64-bit integers
tic; n = n(isprime(n)); toc

Unpatched Octave:    Elapsed time is 923.134 seconds.   (and also takes 3.2 GB of memory.)

Patched Octave:      Elapsed time is 0.038898 seconds.

In both cases the result is the same. There are 21 primes among the 1000 largest 64-bit integers:

n =
  18446744073709551557  18446744073709551533  18446744073709551521  18446744073709551437
  18446744073709551427  18446744073709551359  18446744073709551337  18446744073709551293
  18446744073709551263  18446744073709551253  18446744073709551191  18446744073709551163
  18446744073709551113  18446744073709550873  18446744073709550791  18446744073709550773
  18446744073709550771  18446744073709550719  18446744073709550717  18446744073709550681

Patch 7 attached with documentation edits.

(file #52114)

Arun Giridhar <arungiridhar>
Sun 17 Oct 2021 12:49:14 PM UTC, comment #7: 

Patch 6 attached, same as Patch 5 with a minor bugfix to factor() that was causing it to fail for inputs 0 and 1.

(file #52113)

Arun Giridhar <arungiridhar>
Sun 17 Oct 2021 12:28:07 PM UTC, comment #6: 

Patch 5 attached with two more improvements:

1. This version of isprime uses Schrage multiplication and eliminates the safeadd utility function. The average execution time is decreased from 30 microseconds to 18 microsecongs, and the code is shorter. An entry has been made in NEWS if it is appropriate.

2. The related function factor(), which was recently sped up for composite numbers but was still slow for large prime numbers, now uses the patched version of isprime() to return immediately if its input is prime, meaning that factor() now becomes several orders of magnitude faster for ALL inputs, not just composite ones. The corresponding NEWS item has been edited.

Patch 5 attached. Please test and benchmark using the benchmark code in comment #2.

(file #52112)

Arun Giridhar <arungiridhar>
Sat 16 Oct 2021 10:19:37 PM UTC, comment #5: 

Patch 4 attached with a comment added about performance experiments and removal of trailing spaces. In short, there was no benefit to swapping a and b in multiplication based on value.

Patch 4 attached.

(file #52111)

Arun Giridhar <arungiridhar>
Sat 16 Oct 2021 02:23:19 PM UTC, comment #4: 

Minor tweak to Patch 2. Lowered the threshold inside isprime.m and replaced the call to primes() with a manual list of small primes.

Benchmark performance improves from 21.8 milliseconds for Patch 2 to 2 milliseconds fpr Patch 3! Overall speedup compared to current primes is 30,000X for the whole range

(file #52108)

Arun Giridhar <arungiridhar>
Sat 16 Oct 2021 03:42:12 AM UTC, comment #3: 

Patch 2 attached to comment #2 but disappeared on previewing the comment!

Patch 2 attached again!

(file #52103)

Arun Giridhar <arungiridhar>
Sat 16 Oct 2021 03:40:56 AM UTC, comment #2: 

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


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.

Arun Giridhar <arungiridhar>
Sat 09 Oct 2021 08:49:34 AM UTC, comment #1: 

I did not test, I just looked at the code.

A first quick improvement: savemultiply takes time proportional to log2(b), so you could swap a and b so that a is larger than b. However, this will not help much, because with the repeated multiplications you will quickly be where each bit (including the leading one if the modulus is large) of both a and b will be set with probability 1/2.

So savemultiply will take on average 63 rounds of saveadd. You could try Schrage's algorithm of modular multiplication: Schrage, L. ``A More Portable Fortran Random Number Generator.'' ACM Trans. Math. Software 5, 132-138, 1979.
See also https://stackoverflow.com/questions/20971888/modular-multiplication-of-large-numbers-in-c -- see the comment to the first answer, that this potentially has to be done iteratively. Note that it should read z2=[z/q]. And as you want that optimally r<q, it will help if a is small. Thus, also here I would say that one should swap a and z so that a is the smaller one of the two, making q large and r small, so that r*[z/q] is as small as possible. And note that this is actually the criterion: r*[z/q] must not overflow for the iteration to terminate.

I do not see immediately how fast this will terminate, but also in the worst case I would guess rather in log2(64) rounds than in 63 rounds. In the average case perhaps even in two rounds.

As to the FIXME: no, I don't think that any step can be vectorized. Thus, writing this as compiled code would probably give a performance increase on the order of 1000.

Michael Leitner <mleitner>
Fri 08 Oct 2021 10:43:58 PM UTC, original submission:  

I have attempted to extend isprime() to go to the full range of 64-bit integers. Currently isprime() tends to slow down for any input over about 1e15 with the prime generation and direct division technique, so I've implemented a Miller-Rabin primality test for larger n, up to 2^64-1.

Please test the attached function and test/benchmark script. For development it's temporarily named isprime2().

Feedback solicited, please!

Arun Giridhar <arungiridhar>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files
file #52122:  patch11.patch added by arungiridhar (13KiB - text/x-patch)
file #52120:  patch10.patch added by arungiridhar (11KiB - text/x-patch)
file #52119:  patch9.patch added by arungiridhar (11KiB - text/x-patch)
file #52117:  patch8.patch added by arungiridhar (5KiB - text/x-patch)
file #52118:  __isprimelarge__.cc added by arungiridhar (6KiB - text/x-c++src)
file #52114:  patch7.patch added by arungiridhar (10KiB - text/x-patch)
file #52113:  patch6.patch added by arungiridhar (9KiB - text/x-patch)
file #52112:  patch5.patch added by arungiridhar (9KiB - text/x-patch)
file #52111:  patch4.patch added by arungiridhar (10KiB - text/x-patch)
file #52108:  patch3.patch added by arungiridhar (8KiB - text/x-patch)
file #52103:  patch2.patch added by arungiridhar (8KiB - text/x-patch)
file #52077:  primality_benchmark.m added by arungiridhar (472B - text/x-objcsrc)
file #52078:  isprime2.m added by arungiridhar (11KiB - text/x-objcsrc)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by siko1056 (Updated the item)
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by arungiridhar (Submitted the item)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.




    Follow 14 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-10-18 arungiridhar Attached File- => Added patch11.patch, #52122
    2021-10-18 arungiridhar Attached File- => Added patch10.patch, #52120
    2021-10-18 arungiridhar Attached File- => Added patch9.patch, #52119
    2021-10-18 arungiridhar Attached File- => Added patch8.patch, #52117
        Attached File- => Added __isprimelarge__.cc, #52118
    2021-10-17 arungiridhar Attached File- => Added patch7.patch, #52114
    2021-10-17 arungiridhar Attached File- => Added patch6.patch, #52113
    2021-10-17 arungiridhar Attached File- => Added patch5.patch, #52112
    2021-10-17 siko1056 StatusNone => Patch Submitted
    2021-10-16 arungiridhar Attached File- => Added patch4.patch, #52111
    2021-10-16 arungiridhar Attached File- => Added patch3.patch, #52108
    2021-10-16 arungiridhar Attached File- => Added patch2.patch, #52103
    2021-10-08 arungiridhar Attached File- => Added primality_benchmark.m, #52077
        Attached File- => Added isprime2.m, #52078

    Back to the top

    Powered by Savane 3.7