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

[Michael Leitner]

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

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.

    _______________________________________________________

Reply to this item at:

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

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