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

 
 

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

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

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 25 Nov 2021 05:35:39 AM UTC, comment #25: 

Thank you again for this cleanup patch, now applied here:

https://hg.savannah.gnu.org/hgweb/octave/rev/4f4fa00edada

Confirm the changes look great, the code is blazing fast now! 🚀

Closing report.

Kai Torben Ohlhus <siko1056>
Group Member
Thu 18 Nov 2021 05:17:07 PM UTC, comment #24: 

I have attached a patch with formatting changes to align with the GNU and Octave style guides. (Moved function names to a new line, indented a while-block, removed a use of iostream that had been used for debugging in development).

(file #52293)

Arun Giridhar <arungiridhar>
Group Member
Sat 30 Oct 2021 03:05:30 PM UTC, comment #23: 

Thanks again.

I added a commit message and pushed your patch here:
http://hg.savannah.gnu.org/hgweb/octave/rev/a49c635b179d

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Sat 30 Oct 2021 12:54:57 AM UTC, comment #22: 

Thanks, Markus. I have added the tag inside the error BIST. The spacing near function names and negation operators has been checked. I edited the comments to match the style, and I reversed a for-loop in the C++ code to avoid repeated calls to numel() in case it was not cached. "Make check" passes.

Patch attached.

Re the speedup, if I run this for the very top end of the 64-bit range:

n = intmax("uint64") - 58;   # n = 2^64 - 59, the largest 64-bit prime
tic; isprime(n); toc


Results:
Old code:  Elapsed time is 20.1576 seconds.
New code:  Elapsed time is 0.000178099 seconds.
Ratio = 20.1576 / 0.000178099 = 113,182.

Saying it is 50,000 times faster for some large inputs actually might be conservative, but that's OK since it's more representative of a wider range, from 1x at the threshold (29 billion) to over 100,000x at 2^64. For inputs below the threshold, it retains the use of the conventional division technique because that is faster for smaller inputs.

(file #52173)

Arun Giridhar <arungiridhar>
Group Member
Fri 29 Oct 2021 06:02:55 PM UTC, comment #21: 

A 50,000 fold speed-up for some large numbers! That is quite an improvement to the existing implementation!

Just a few more nit-picking comments. (I hope you don't mind.)
- Use "##" for comments in .m files that start at the beginning of a line.
- Use "#" preceded by two spaces for comments that trail code in the same line.
- For the negation operator "!", follow up with a space in both C++ and Octave code.
- Use a space between function name and opening parenthesis in the actual code and in tests (unless in an [array]).
- If applicable, add a snippet from the expected error message in `%!error` tests to make sure the function is failing in the way we expect. E.g. `%!error <unable to convert input> _isprimelarge_ ({'foo'; 'bar'})`

Markus Mützel <mmuetzel>
Group administrator
Thu 28 Oct 2021 10:33:57 PM UTC, comment #20: 

@Markus: Thank you for the tip about boolNDArray. I was able to get it to work with that change. Earlier I had been trying to construct an octave value vector directly instead of converting a boolNDArray with ovl().

I'm also very surprised at exactly how much overhead was present in that last for-loop. Moving that to C++ sped up the overall performance by 2.5x from that change alone. The latest patch (Patch 12, attached to this comment) is therefore 2.5x as fast as Patch 11 in Comment #16. Thank you for the guidance to look at that part.

Patch 12 attached with the following changes compared to Patch 11:

1. Moved the final for-loop in isprime.m to C++, enabling a vector to be passed instead of a sequence of scalars.

2. Renamed internal functions in the Miller-Rabin code and made some documentation edits for the end user.

3. Since Miller-Rabin was now faster than before, the THRESHOLD value needed to be retuned in isprime.m.

Please examine the last line of isprime() and the corresponding C++ use of boolNDArray in _isprimelarge_. Is there a performance penalty to boolNDArray when the Octave function will always call it with a row vector?

(file #52166)

Arun Giridhar <arungiridhar>
Group Member
Thu 28 Oct 2021 10:16:22 AM UTC, comment #19: 

@Arun: I'm not sure I can correctly follow the changes you describe.
Could you please upload the wip patch (even if it doesn't do the correct thing currently)?

I would have thought, you could use a boolNDArray to collect the results of a C++ function `isprimelarge` (which returns a bool) and return only one single output argument from the DEFUN `__isprimelargevector__`.

Markus Mützel <mmuetzel>
Group administrator
Wed 27 Oct 2021 10:15:04 PM UTC, comment #18: 

@Markus: I tried changing the C++ function to handle n-dimensional arrays instead of scalars so that the for loop could be migrated but I keep getting an error "value on right hand side of assignment is undefined". This has nothing to do with prime numbers and everything to do with conversion to/from uint64NDArray.

If it's not a performance problem having that final for-loop, Patch 11 can be applied if it looks reliable to you.

Otherwise, this is what I'm trying to do in C++ in a new DEFUN function _isprimelargevector_ to eliminate the last for-loop, while the scalar function _isprimelarge_ has been taken out of its DEFUN. It builds but doesn't return anything and the BISTs fail.


uint64NDArray vec = args(0).xuint64_array_value
("__isprimelargevector__: unable to convert input. Call isprime() instead.");

octave_value_list retval (vec.numel());

for (octave_idx_type i = 0; i < vec.numel(); i++)
  retval(i) = __isprimelarge__ (vec(i));

return retval;


Arun Giridhar <arungiridhar>
Group Member
Wed 27 Oct 2021 04:11:13 PM UTC, comment #17: 

I haven't tested the code yet. But the patch looks very good already.

Loops in Octave's scripting language tend to be slower than in C++.
Would it make sense to move the final loop to inside `__isprimelarge__`?

Markus Mützel <mmuetzel>
Group administrator
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>
Group Member
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>
Group Member
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>
Group Member
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>
Group Member
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>
Group administrator
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.

Attachments:

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>
Group Member
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>
Group 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
  18446744073709550671


Patch 7 attached with documentation edits.

(file #52114)

Arun Giridhar <arungiridhar>
Group Member
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>
Group Member
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>
Group Member
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>
Group Member
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>
Group Member
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>
Group Member
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

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.

Arun Giridhar <arungiridhar>
Group Member
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>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #52293:  patch_clean.patch added by arungiridhar (2KiB - text/x-patch)
file #52173:  patch_clean.patch added by arungiridhar (14KiB - text/x-patch)
file #52166:  patch12.patch added by arungiridhar (14KiB - text/x-patch)
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)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 20 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-11-25 siko1056 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2021-11-18 arungiridhar Attached File- Added patch_clean.patch, #52293
    2021-10-30 mmuetzel StatusPatch Submitted Ready For Test
    2021-10-30 arungiridhar Attached File- Added patch_clean.patch, #52173
    2021-10-28 arungiridhar Attached File- Added patch12.patch, #52166
    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.13-4448.
    Corresponding source code