Ciao,
Il Mar, 5 Giugno 2012 11:34 pm, Torbjorn Granlund ha scritto:
> 1. mpz_pprime_p, for doing some trial dividing and unspecified
> probabilistic primality tests to specified probability.
I'd like to keep the "Return 2 if n is definitely prime" feature of the
current _probab_prime_ function.
> 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.
> 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.
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. A loop like the one in _pprime_p will
then use the random numbers generated by _urandomm as redcified...
Or we can write a mpn_redc_millerrabin, and use it for all the mpz
functions above.
> 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.
About the redcfy interface: _modredc can be used both after a
multiplication, or after a (sequence of) subtraction/addition?
Regards,
Marco
--
http://bodrato.it/software/strassen.html
/* mpz_pprime_p and mpz_notdiv_pprime_p.
Copyright 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2005,
2008, 2012 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
#include "gmp.h"
#include "gmp-impl.h"
#include "redcify.h"
#undef mpz_millerrabin
int mpz_millerrabin (mpz_srcptr, mpz_srcptr);
/* Test n for probabilistic primality tests to probability 2^(-logprob),
using rs for generating random numbers. Performs trial dividing.
To be used when the likelihood for small prime factors are high. */
int
mpz_pprime_p (mpz_srcptr n, int logprob, gmp_randstate_t rs)
{
mp_size_t nprimes;
int where;
mp_ptr np = PTR(n);
mp_size_t nn = ABSIZ(n);
mp_limb_t divisor;
mpz_t r, nm3;
int isprime;
TMP_DECL;
/* Handle even numbers. */
if ((np[0] & 1) == 0)
return nn == 1 && np[0] == 2;
/* Do some trial dividing. FIXME: We should probably trial divide to sqrt(n)
for small enough n. */
where = 0;
nprimes = GMP_NUMB_BITS * nn; /* FIXME: Compute more carefully */
divisor = mpn_trialdiv (np, nn, nprimes, &where);
if (divisor != 0)
return 0;
TMP_MARK;
MPZ_TMP_INIT (nm3, nn + 1);
mpz_sub_ui (nm3, n, 3);
MPZ_TMP_INIT (r, nn + 1);
/* Perform a number of Miller-Rabin tests. */
isprime = 1;
while (logprob >= 0 && isprime)
{
mpz_urandomm (r, rs, nm3);
mpz_add_ui (r, r, 2);
isprime = mpz_millerrabin (n, r);
logprob -= 2;
}
TMP_FREE;
return isprime;
}
/* Test n for probabilistic primality tests to probability 2^(-logprob),
using rs for generating random numbers. Performs trial dividing.
To be used when the likelihood for small prime factors are low. */
int
mpz_notdiv_pprime_p (mpz_srcptr n, int logprob, gmp_randstate_t rs)
{
mp_ptr np = PTR(n);
mp_size_t nn = ABSIZ(n);
mpz_t r, nm3;
int isprime;
TMP_DECL;
/* Handle even numbers. */
if ((np[0] & 1) == 0)
return nn == 1 && np[0] == 2;
TMP_MARK;
MPZ_TMP_INIT (nm3, nn + 1);
mpz_sub_ui (nm3, n, 3);
MPZ_TMP_INIT (r, nn + 1);
/* Perform a number of Miller-Rabin tests. */
isprime = 1;
while (logprob >= 0 && isprime)
{
mpz_urandomm (r, rs, nm3);
mpz_add_ui (r, r, 2);
isprime = mpz_millerrabin (n, r);
logprob -= 2;
}
TMP_FREE;
return isprime;
}
/* Perform a Miller-Rabin test on n for the base r.
We use REDC representation, except that we call mpz_powm using plain
residues. FIXME: Stay in REDC, requires reimplementation of mpz_powm. */
int
mpz_millerrabin (mpz_srcptr n, mpz_srcptr r)
{
mpz_t nm1, one, y, q;
mp_bitcnt_t k, i;
mpz_redc_t nredc;
int retval;
TMP_DECL;
TMP_MARK;
ASSERT (mpz_cmp_ui (n, 2) > 0 && mpz_odd_p (n));
MPZ_TMP_INIT (nm1, SIZ(n) + 1);
mpz_sub_ui (nm1, n, 1);
MPZ_TMP_INIT (q, 2 * SIZ(n));
MPZ_TMP_INIT (y, SIZ(n));
/* Find q and k, where q is odd and n = 1 + 2**k * q. */
k = mpz_scan1 (nm1, 0); /* There is a shortcut for scan1(., 0), use it. */
mpz_tdiv_q_2exp (q, n, k);
mpz_powm (y, r, q, n);
if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0)
{
retval = 1;
}
else {
retval = 0;
if (k > 1) { /* if n = 3 mod 4, i.e. k==1 Miller-Rabin is done. */
mpz_redc_init_set (nredc, n);
MPZ_TMP_INIT (one, SIZ(n) + 1);
mpz_set_ui (one, 1);
mpz_redcify (one, one, nredc);
mpz_sub (nm1, n, one); /* mpz_redcify (nm1, nm1, nredc); */
mpz_redcify (y, y, nredc);
for (i = 1; i < k; k--) /* FIXME: use a do {} while, get rid of i. */
{
mpz_mul (q, y, y);
mpz_modredc (y, q, nredc);
if (mpz_cmp (y, nm1) == 0)
{
retval = 1;
break;
}
if (mpz_cmp (y, one) == 0)
break;
}
mpz_redc_clear (nredc);
}}
TMP_FREE;
return retval;
}
int
main ()
{
mpz_t t;
gmp_randstate_t rands;
gmp_randinit_default (rands);
mpz_init (t);
for (;;)
{
mpz_urandomb (t, rands, 1000);
gmp_printf ("%40Zd %d\n", t, mpz_pprime_p (t, 1, rands));
mpz_nextprime (t, t);
gmp_printf ("%40Zd (%d)\n", t, mpz_pprime_p (t, 1, rands));
}
}_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
http://gmplib.org/mailman/listinfo/gmp-devel
