This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.32
in repository libmath-prime-util-perl.

commit 67979c337e9327bfa4c71d463fd927edb0aaf5ea
Author: Dana Jacobsen <d...@acm.org>
Date:   Fri Sep 27 15:12:01 2013 -0700

    Add miller_rabin_random; add primality testing opinions to documentation
---
 Changes                |   1 +
 lib/Math/Prime/Util.pm | 194 +++++++++++++++++++++++++++++++++++++++++--------
 t/02-can.t             |  45 ++++++++----
 3 files changed, 195 insertions(+), 45 deletions(-)

diff --git a/Changes b/Changes
index 810e5c4..2eaea16 100644
--- a/Changes
+++ b/Changes
@@ -8,6 +8,7 @@ Revision history for Perl module Math::Prime::Util
       - carmichael_lambda
       - znorder
       - prime_iterator_object
+      - miller_rabin_random
 
     - Added Math::Prime::Util::PrimeIterator.  A more feature-rich iterator
       than the simple closure one from prime_iterator.  Experimental.
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index a26a375..0c5cd10 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -23,7 +23,7 @@ our @EXPORT_OK =
       is_almost_extra_strong_lucas_pseudoprime
       is_frobenius_underwood_pseudoprime
       is_aks_prime
-      miller_rabin
+      miller_rabin miller_rabin_random
       lucas_sequence
       primes
       forprimes prime_iterator prime_iterator_object
@@ -1988,24 +1988,64 @@ sub lucas_sequence {
   return Math::Prime::Util::PP::lucas_sequence($n, $P, $Q, $k);
 }
 
+sub miller_rabin_random {
+  my($n, $k, $seed) = @_;
+  _validate_num($n) || _validate_positive_integer($n);
+  _validate_num($k) || _validate_positive_integer($k);
+
+  return 1 if $k <= 0;
+  return (is_prob_prime($n) > 0) if $n < 100;
+  return 0 unless $n & 1;
+  return Math::Prime::Util::GMP::miller_rabin_random($n, $k)
+    if $_HAVE_GMP
+    && defined &Math::Prime::Util::GMP::miller_rabin_random;
+
+  # Testing this many bases is silly, but let's pretend they have some
+  # good reason.  A composite n > 9 must have at least n/4 witnesses,
+  # hence we need to check only floor(3/4)+1 at most.  We could improve
+  # this is $_Config{'assume_rh'} is true, to 1 .. 2(logn)^2.
+  if ($k >= int(3*$n/4)) {
+    return is_strong_pseudoprime($n, 2 .. int(3*$n/4)+1+2 );
+  }
+
+  my $brange = $n-3;
+  my $irandf = _get_rand_func();
+  # Do one first before doing batches
+  return 0 unless is_strong_pseudoprime($n, $irandf->($brange)+2 );
+  $k--;
+  while ($k > 0) {
+    my $nbases = ($k >= 20) ? 20 : $k;
+    my @bases = map { $irandf->($brange)+2 } 1..$nbases;
+    if ($n <= $_XS_MAXVAL) {
+      return 0 unless _XS_miller_rabin($n, @bases);
+    } else {
+      return 0 unless is_strong_pseudoprime($n, @bases);
+    }
+    $k -= $nbases;
+  }
+  1;
+}
+
 
 #############################################################################
 
-  # Simple timings, do { is_strong_pseudoprime(2*$_+1,2) } for 1..1000000;
+  # Simple timings, with 0.32 code (Montgomery Reduction for 64-bit)
+  # my $n=2**20+1;  do { is_strong_pseudoprime($n,2); $n+=2 } for 1..1000000;
+  # my $n=2**47+1;  do { is_strong_pseudoprime($n,2); $n+=2 } for 1..1000000;
   #
   #      1.0 uS  XS 32-bit input, is_strong_pseudoprime
   #      1.8 uS  XS 64-bit input, is_strong_pseudoprime
-  #      1.4 uS  XS 32-bit input, is_strong_lucas_pseudoprime
-  #      3.0 uS  XS 64-bit input, is_strong_lucas_pseudoprime
+  #      1.3 uS  XS 32-bit input, is_strong_lucas_pseudoprime
+  #      1.9 uS  XS 64-bit input, is_strong_lucas_pseudoprime
   #      1.2 uS  XS 32-bit input, is_almost_extra_strong_lucas_pseudoprime
-  #      2.3 uS  XS 64-bit input, is_almost_extra_strong_lucas_pseudoprime
+  #      1.8 uS  XS 64-bit input, is_almost_extra_strong_lucas_pseudoprime
   #
-  #     12 uS  Perl 32-bit input, is_strong_pseudoprime
-  #   1200 uS  Perl 64-bit input, is_strong_pseudoprime
-  #   2000 uS  Perl 32-bit input, is_strong_lucas_pseudoprime
-  #   7840 uS  Perl 64-bit input, is_strong_lucas_pseudoprime
-  #    940 uS  Perl 32-bit input, is_almost_extra_strong_lucas_pseudoprime
-  #   3360 uS  Perl 64-bit input, is_almost_extra_strong_lucas_pseudoprime
+  #     13 uS  Perl 32-bit input, is_strong_pseudoprime
+  #    945 uS  Perl 64-bit input, is_strong_pseudoprime
+  #   1700 uS  Perl 32-bit input, is_strong_lucas_pseudoprime
+  #   3000 uS  Perl 64-bit input, is_strong_lucas_pseudoprime
+  #   1500 uS  Perl 32-bit input, is_almost_extra_strong_lucas_pseudoprime
+  #   3400 uS  Perl 64-bit input, is_almost_extra_strong_lucas_pseudoprime
 
 sub is_aks_prime {
   my($n) = @_;
@@ -2729,6 +2769,21 @@ and earlier with bignums.
 =back
 
 
+=head1 PRIMALITY TESTING
+
+This module provides three functions for general primality testing, as
+well as numerous specialized functions.  The three main functions are:
+L</is_prob_prime> and L</is_prime> for general use, and L</is_provable_prime>
+for proofs.  For inputs below C<2^64> the functions are identical and
+fast deterministic testing is performed.  That is, the results will always
+be correct and should take at most a few microseconds for any input.  This
+is hundreds to thousands of times faster than other CPAN modules.  For
+inputs larger than C<2^64>, an extra-strong
+L<BPSW test|http://en.wikipedia.org/wiki/Baillie-PSW_primality_test>
+is used.  See the L</PRIMALITY TESTING NOTES> section for more
+discussion.
+
+
 =head1 FUNCTIONS
 
 =head2 is_prime
@@ -3035,10 +3090,6 @@ Returning composite for all non-2 even input makes the 
function match most
 other implementations including L<Math::Primality>'s C<is_strong_pseudoprime>
 function.
 
-=head2 miller_rabin
-
-An alias for C<is_strong_pseudoprime>.  This name is deprecated.
-
 =head2 is_lucas_pseudoprime
 
 Takes a positive number as input, and returns 1 if the input is a standard
@@ -3100,6 +3151,19 @@ C<(L+2)^(n-1) = 5 + 2x mod (n, L^2 - Lx + 1)>.  The 
computational cost for this
 is between the cost of 2 and 3 strong pseudoprime tests.  There are no known
 counterexamples, but this is not a well studied test.
 
+=head2 miller_rabin_random
+
+Takes a positive number (C<n>) as input and a positive number (C<k>) of bases
+to use.  Performs C<k> Miller-Rabin tests using uniform random bases
+between 2 and C<n-2>.
+
+This should not be used in place of L</is_prob_prime>, L</is_prime>,
+or L</is_provable_prime>.  Those functions will be faster and provide
+better results than running C<k> Miller-Rabin tests.  This function can
+be used if one wants more assurances for non-proven primes, such as for
+cryptographic uses where the size is large enough that proven primes are
+not desired.
+
 
 
 =head2 is_prob_prime
@@ -4257,7 +4321,7 @@ Project Euler, problem 21 (Amicable numbers):
 
 Project Euler, problem 41 (Pandigital prime), brute force command line:
 
-  perl -MMath::Prime::Util=:all -E 'my @p = grep { 
/1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} @{primes(1000000,9999999)}; say $p[-1];
+  perl -MMath::Prime::Util=:all -E 'my @p = grep { 
/1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} @{primes(1000000,9999999)}; say $p[-1];'
 
 Project Euler, problem 47 (Distinct primes factors):
 
@@ -4268,7 +4332,7 @@ Project Euler, problem 47 (Distinct primes factors):
   $n++ while (nfactors($n) != 4 || nfactors($n+1) != 4 || nfactors($n+2) != 4 
|| nfactors($n+3) != 4);
   say $n;
 
-Project Euler, problem 69, stupid brute force solution (about 5 seconds):
+Project Euler, problem 69, stupid brute force solution (about 2 seconds):
 
   use Math::Prime::Util qw/euler_phi/;
   my ($n, $max) = (0,0);
@@ -4278,14 +4342,14 @@ Project Euler, problem 69, stupid brute force solution 
(about 5 seconds):
   } for 1..1000000;
   say "$n  $max";
 
-Here's the right way to do PE problem 69 (under 0.03s):
+Here is the right way to do PE problem 69 (under 0.03s):
 
   use Math::Prime::Util qw/pn_primorial/;
   my $n = 0;
   $n++ while pn_primorial($n+1) < 1000000;
-  say pn_primorial($n);'
+  say pn_primorial($n);
 
-Project Euler, problem 187, stupid brute force solution, ~3 minutes:
+Project Euler, problem 187, stupid brute force solution, ~4 minutes:
 
   use Math::Prime::Util qw/factor -nobigint/;
   my $nsemis = 0;
@@ -4293,7 +4357,7 @@ Project Euler, problem 187, stupid brute force solution, 
~3 minutes:
      for 1 .. int(10**8)-1;
   say $nsemis;
 
-Here's the best way for PE187.  Under 30 milliseconds from the command line:
+Here is the best way for PE187.  Under 30 milliseconds from the command line:
 
   use Math::Prime::Util qw/primes prime_count -nobigint/;
   use List::Util qw/sum/;
@@ -4303,6 +4367,78 @@ Here's the best way for PE187.  Under 30 milliseconds 
from the command line:
                1 .. scalar @primes );
 
 
+=head1 PRIMALITY TESTING NOTES
+
+Above C<2^64>, L</is_prob_prime> performs an extra-strong
+L<BPSW test|http://en.wikipedia.org/wiki/Baillie-PSW_primality_test>
+which is fast (a little less than the time to perform 3 Miller-Rabin
+tests) and has no known counterexamples.  If you believe the primality
+testing done by Pari, Maple, SAGE, FLINT, etc., then this function
+should be appropriate for you.  L</is_prime> will do the same BPSW
+test as well as some additional testing, making it slightly more time
+consuming but less likely to produce a false result.  This is a little
+more stringent than Mathematica.  L</is_provable_prime> constructs a
+primality proof.  If a certificate is requested, then either BLS75
+theorem 5 or ECPP is performed.  Without a certificate, the method
+is implementation specific (currently it is identical, but later
+releases may use APRCL).  With L<Math::Prime::Util::GMP> installed,
+this is quite fast through 300 or so digits.
+
+Math systems 30 years ago used to use Miller-Rabin tests with C<k>
+bases (typically fixed bases, sometimes random) for primality
+testing, but these have generally been replaced by some form of BPSW
+as used in this module.  See Pinch's 1993 paper for examples of why
+using C<k> M-R tests leads to poor results.  The three exceptions in
+common contemporary use I am aware of are:
+
+=over 4
+
+=item libtommath
+
+Uses the first C<k> prime bases.  This is unacceptable for
+cryptographic use, as there are known methods (e.g. Arnault 1994) for
+constructing counterexamples.  The number of bases required to avoid
+false results is unreasonably high, hence performance is slow even
+if one ignores counterexamples.  Unfortunately this is the
+multi-precision math library used for Perl 6 and at least one CPAN
+Crypto module.
+
+=item GMP/MPIR
+
+Uses a set of C<k> static-random bases.  The bases are randomly chosen
+using a PRNG that is seeded identically each call (the seed changes
+with each release).  This offers a very slight advantage over using
+the first C<k> prime bases, but not much.  See, for example, Nicely's
+L<mpz_probab_prime_p pseudoprimes|http://www.trnicely.net/misc/mpzspsp.html>
+page.
+
+=item L<Math::Pari>
+
+Pari 2.1.7 is the default version installed with the L<Math::Pari>
+module.  It uses 10 random M-R bases (the PRNG uses a fixed seed
+set at compile time).  Pari 2.3.0 was released in May 2006 and it,
+like all later releases through at least 2.6.1, use BPSW / APRCL,
+after complaints of false results from using M-R tests.
+
+=back
+
+Basically the problem is that it is just too easy to get counterexamples
+from running C<k> M-R tests, forcing one to use a very large number of
+tests (at least 20) to avoid frequent false results.  Using the BPSW test
+results in no known counterexamples after 30+ years and runs much faster.
+It can be enhanced with one or more random bases if one desires, and
+will I<still> be much faster.
+
+Using C<k> fixed bases has another problem, which is that in any
+adversarial situation we can assume the inputs will be selected such
+that they are one of our counterexamples.  Now we need absurdly large
+numbers of tests.  This is like playing "pick my number" but the
+number is fixed forever at the start, the guesser gets to know
+everyone else's guesses and results who has every played, and can
+keep playing as long as they like.  It's only valid if the players
+are completely oblivious to what is happening.
+
+
 =head1 LIMITATIONS
 
 Perl versions earlier than 5.8.0 have a rather broken 64-bit implementation,
@@ -4458,7 +4594,7 @@ Some of the highlights:
 
 =item C<isprime>
 
-Similar to MPU's L<is_prob_prime> or L<is_prime> functions.
+Similar to MPU's L</is_prob_prime> or L</is_prime> functions.
 MPU is deterministic for native integers, and uses a strong
 BPSW test for bigints (with a quick primality proof tried as well).  The
 default version of Pari used by L<Math::Pari> (2.1.7) uses 10 random M-R
@@ -4479,7 +4615,7 @@ function will perform a BPSW test similar to 
L</is_prob_prime>.
 
 =item C<primepi>
 
-Only available with version 2.3 of Pari.  Similar to MPU's L<prime_count>
+Only available with version 2.3 of Pari.  Similar to MPU's L</prime_count>
 function in API, but uses a naive counting algorithm with its precalculated
 primes, so is not of practical use.  Incidently, Pari 2.6 (not usable from
 Perl) has fixed the pre-calculation requirement so it is more useful, but is
@@ -4495,7 +4631,7 @@ doesn't support segmenting.
 
 =item C<factorint>
 
-Similar to MPU's L<factor> though with a different return (I
+Similar to MPU's L</factor> though with a different return (I
 find the result value quite inconvenient to work with, but others may like
 its vector of factor/exponent format).  Slower than MPU for all 64-bit inputs
 on an x86_64 platform, it may be faster for large values on other platforms.
@@ -4504,27 +4640,27 @@ faster in MPU.
 
 =item C<eulerphi>
 
-Similar to MPU's L<euler_phi>.  MPU is 2-20x faster for native integers.
+Similar to MPU's L</euler_phi>.  MPU is 2-20x faster for native integers.
 There is also support for a range, which can be much more efficient.
 Without L<Math::Prime::Util::GMP> installed, MPU is very slow with bigints.
 With it installed, it is about 2x slower than Math::Pari.
 
 =item C<moebius>
 
-Similar to MPU's L<moebius>.  Comparisons are similar to C<eulerphi>.
+Similar to MPU's L</moebius>.  Comparisons are similar to C<eulerphi>.
 
 =item C<sigma>
 
-Similar to MPU's L<divisor_sum>.  MPU is ~10x faster for native integers
+Similar to MPU's L</divisor_sum>.  MPU is ~10x faster for native integers
 and about 2x slower for bigints.
 
 =item C<eint1>
 
-Similar to MPU's L<ExponentialIntegral>.
+Similar to MPU's L</ExponentialIntegral>.
 
 =item C<zeta>
 
-A more feature-rich version MPU's L<RiemannZeta> function (supports negative
+A more feature-rich version MPU's L</RiemannZeta> function (supports negative
 and complex inputs).
 
 =back
diff --git a/t/02-can.t b/t/02-can.t
index 89498c0..160daa9 100644
--- a/t/02-can.t
+++ b/t/02-can.t
@@ -6,21 +6,34 @@ use Math::Prime::Util;
 use Test::More  tests => 1;
 
 my @functions =  qw(
-                     prime_get_config prime_set_config
-                     prime_precalc prime_memfree
-                     is_prime is_prob_prime is_provable_prime
-                     is_strong_pseudoprime is_strong_lucas_pseudoprime
-                     is_aks_prime
-                     miller_rabin
-                     primes
-                     next_prime  prev_prime
-                     prime_count prime_count_lower prime_count_upper 
prime_count_approx
-                     nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
-                     random_prime random_ndigit_prime random_nbit_prime 
random_maurer_prime
-                     primorial pn_primorial
-                     factor all_factors
-                     moebius euler_phi jordan_totient
-                     divisor_sum
-                     ExponentialIntegral LogarithmicIntegral RiemannZeta 
RiemannR
+      prime_get_config prime_set_config
+      prime_precalc prime_memfree
+      is_prime is_prob_prime is_provable_prime is_provable_prime_with_cert
+      prime_certificate verify_prime
+      is_pseudoprime is_strong_pseudoprime
+      is_lucas_pseudoprime
+      is_strong_lucas_pseudoprime
+      is_extra_strong_lucas_pseudoprime
+      is_almost_extra_strong_lucas_pseudoprime
+      is_frobenius_underwood_pseudoprime
+      is_aks_prime
+      miller_rabin_random
+      lucas_sequence
+      primes
+      forprimes prime_iterator prime_iterator_object
+      next_prime  prev_prime
+      prime_count
+      prime_count_lower prime_count_upper prime_count_approx
+      nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
+      random_prime random_ndigit_prime random_nbit_prime random_strong_prime
+      random_proven_prime random_proven_prime_with_cert
+      random_maurer_prime random_maurer_prime_with_cert
+      primorial pn_primorial consecutive_integer_lcm
+      factor all_factors
+      moebius mertens euler_phi jordan_totient exp_mangoldt
+      chebyshev_theta chebyshev_psi
+      divisor_sum
+      carmichael_lambda znorder
+      ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR
                    );
 can_ok( 'Math::Prime::Util', @functions);

-- 
Alioth's /usr/local/bin/git-commit-notice on 
/srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git

_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits

Reply via email to