Pseudo prime tests

[Torbjorn Granlund]

bodrato at mail.dm.unipi.it writes:

  I'd like to keep the "Return 2 if n is definitely prime" feature of the
  current _probab_prime_ function.
  
Yes, that would be nice.

  > 2. mpz_notdiv_pprime_p, like mpz_pprime_p but without trial dividing.
  
  Do we really need a new function for this? An additional parameter should
  be enough.
  
I am not sure that's easier for users.

  > 3. mpz_millerrabin, replacing the current internal function with the
  >    same name (but with different parameters).
  
  This can be useful if we decide to export it, but it (slightly) slows down
  the repeated use, because nm1, one, nredc, q and k must be recomputed at
  each call.

Indeed.  But I see no really good solution to that.
Fortunately, I think the slowdown is almost unmeasurable, except
when tested numbers are tiny.

  Moreover, the comment to this function says: "FIXME: Stay in REDC,
  requires reimplementation of mpz_powm". In this case we should require an
  already redcified base argument.

I agree.  (Perhaps the current powm code could then invoke that code.)

  > The new code makes use of redc functions at the mpz level, which I wrote
  > several years ago.  (Such functions could also at some point be made
  > part of the public GMP interface.)
  
  But I fear mpz_millerrabin doesn't use them the best way: nm1 is redcfied
  too early, then it's compared to the non redcfied result of _powm,
  possibly giving a wrong "composite" answer on prime numbers! Moreover we
  don't need to use the redcfy function on both 1 and (n-1). I attach a
  working copy.
  
My function could falsely claim primes are composite, not just that
composites are prime.  I thought that was more symmetric.  :-)

  About the redcfy interface: _modredc can be used both after a
  multiplication, or after a (sequence of) subtraction/addition?
  
I don't think that would work, since it multiplies by B^(-n) mod m.
Only after multiplying two residues in the special B^n format, it will
work as expected to divide out B^n.

Thanks for reviewing my code, and thanks for the working version of it!

-- 
Torbjörn