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

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

commit 8dffb63afdee53cc803f6aa7faed6a03ee7ddd31
Author: Dana Jacobsen <d...@acm.org>
Date:   Fri Jul 19 14:53:17 2013 -0700

    Clean up Lucas, update PP is_prime, remove an test func that got put in 
XS.xs, add documentation
---
 Changes                   |   3 +-
 TODO                      |  11 ----
 XS.xs                     |  27 ---------
 lib/Math/Prime/Util.pm    |  54 ++++++++++++------
 lib/Math/Prime/Util/PP.pm | 140 ++++++++++++++++++++++++++++++++++++----------
 5 files changed, 150 insertions(+), 85 deletions(-)

diff --git a/Changes b/Changes
index 7adfc64..b809319 100644
--- a/Changes
+++ b/Changes
@@ -7,7 +7,8 @@ Revision history for Perl extension Math::Prime::Util.
 
     - Fixed a rare refcount / bignum / callback issue.
 
-    - Small mulmod speedup for non-gcc/x86_64 platforms.
+    - Small mulmod speedup for non-gcc/x86_64 platforms, and for any platform
+      with gcc 4.4 or newer.
 
     - Add more conditions to ECPP block verification.
 
diff --git a/TODO b/TODO
index 598ea96..c89456e 100644
--- a/TODO
+++ b/TODO
@@ -30,13 +30,6 @@
 
 - Implement S2 calculation for LMO prime count.
 
-- Adding a Primo certificate parser (see WraithX's code) would
-  be awesome, but may be a lot more work.  It would still be nice to have yet
-  another independent codebase for this.
-
-- Similar to the previous, having a method to convert our certificates into
-  Primo-style ones would be useful for easy independent verification.
-
 - Big features:
    - LMO prime count
    - QS factoring
@@ -51,7 +44,3 @@
   broken on 5.16 and previous, or if this was just a 5.17 issue.
 
 - Switch to new text proofs.
-
-- Sync all the Lucas pseudoprimes with GMP.
-
-- Update Perl code for Lucas, use fastest available (probably AES).
diff --git a/XS.xs b/XS.xs
index e12c7c0..5953f31 100644
--- a/XS.xs
+++ b/XS.xs
@@ -707,33 +707,6 @@ _validate_num(SV* n, ...)
   OUTPUT:
     RETVAL
 
-UV
-testfac (IN UV bits, IN UV num)
-  PREINIT:
-    UV n, i, nf, base, numfacs;
-    UV factors[64];
-  CODE:
-    numfacs = 0;
-    base = 1 + (1UL << bits);
-    for (i = 0; i < num; i++) {
-      n = base+2*i;
-      nf = trial_factor(n, factors, 1001);
-      n = factors[--nf];
-      if (_XS_is_prob_prime(n)) continue;
-      //numfacs += holf_factor(n, factors, 20000) - 1;
-      //numfacs += pminus1_factor(n, factors, 5000, 80000) - 1;
-      //numfacs += pminus1_factor(n, factors, 300, 60000) - 1;
-      numfacs += pplus1_factor(n, factors, 5000) - 1;
-      //numfacs += squfof_factor(n, factors, 20000) - 1;
-      //numfacs += racing_squfof_factor(n, factors, 10000) - 1;
-      //numfacs += pbrent_factor(n, factors, 1500, 3) - 1;
-      //numfacs += pplus1_factor(n, factors, 200) - 1;
-      //numfacs += lehman_factor(n, factors) - 1;
-    }
-    RETVAL = numfacs;
-  OUTPUT:
-    RETVAL
-
 void
 forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
   PROTOTYPE: &$;$
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index ba2412e..f2e552b 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -2483,7 +2483,7 @@ __END__
 
 =encoding utf8
 
-=for stopwords forprimes Möbius Deléglise totient moebius mertens irand 
primesieve uniqued k-tuples von SoE pari yafu fonction qui compte le nombre 
nombres voor PhD
+=for stopwords forprimes Möbius Deléglise totient moebius mertens irand 
primesieve uniqued k-tuples von SoE pari yafu fonction qui compte le nombre 
nombres voor PhD superset
 
 
 =head1 NAME
@@ -2740,10 +2740,11 @@ a lot.  There are some brittle behaviors on 5.12.4 and 
earlier with bignums.
   print "$n is prime" if is_prime($n);
 
 Returns 2 if the number is prime, 0 if not.  For numbers larger than C<2^64>
-it will return 0 for composite and 1 for probably prime, using a strong BPSW
-test.  If L<Math::Prime::Util::GMP> is installed, some additional
-primality tests are also performed on large inputs, and a quick attempt is
-made to perform a primality proof, so it will return 2 for many other inputs.
+it will return 0 for composite and 1 for probably prime, using an
+extra-strong BPSW test.  If L<Math::Prime::Util::GMP> is installed, some
+additional primality tests are also performed on large inputs, and a
+quick attempt is made to perform a primality proof, so it will
+return 2 for many other inputs.
 
 Also see the L</is_prob_prime> function, which will never do additional
 tests, and the L</is_provable_prime> function which will construct a proof
@@ -2753,8 +2754,9 @@ expense of speed).
 For native precision numbers (anything smaller than C<2^64>, all three
 functions are identical and use a deterministic set of tests (selected
 Miller-Rabin bases or BPSW).  For larger inputs both L</is_prob_prime> and
-L</is_prime> return probable prime results using the strong Baillie-PSW test,
-which has had no counterexample found since it was published in 1980.
+L</is_prime> return probable prime results using the extra-strong
+Baillie-PSW test, which has had no counterexample found since it was
+published in 1980.
 
 
 
@@ -3029,17 +3031,37 @@ L<OEIS A217255|http://oeis.org/A217255>.
 
 Takes a positive number as input, and returns 1 if the input passes the extra
 strong Lucas test (as defined in
-L<Grantham 2000|http://www.ams.org/mathscinet-getitem?mr=1680879>).  Parameter
-selection is done by incrementing C<P> from C<3> until C<jacobi(D,n) = -1>.
-This has slightly more restrictive conditions than the strong Lucas test,
-but uses different starting parameters so is not directly comparable.
+L<Grantham 2000|http://www.ams.org/mathscinet-getitem?mr=1680879>).  This
+test has more stringent conditions than the strong Lucas test, and produces
+about 60% fewer pseudoprimes.  Performance is typically 20-30% I<faster>
+than the strong Lucas test.
+
+The parameters are selected using the
+L<Baillie-OEIS method|http://oeis.org/A217719>
+method: increment C<P> from C<3> until C<jacobi(D,n) = -1>.
 Removing primes, this produces the sequence
 L<OEIS A217719|http://oeis.org/A217719>.
 
-The extra strong Lucas test typically performs 20 to 30% faster than the
-strong Lucas test, and produces fewer pseudoprimes.  There are no
-counterexamples below C<2^64> with BPSW using any of the Lucas tests, and
-no published counterexamples of any size.
+=head2 is_almost_extra_strong_lucas_pseudoprime
+
+This is similar to the L</is_extra_strong_lucas_pseudoprime> function, but
+does not calculate C<U>, so is a little faster, but also weaker.
+With the current implementations, there is little reason to prefer this unless
+trying to reproduce specific results.  The extra-strong implementation has been
+optimized to use similar features, removing most of the performance advantage.
+
+An optional second argument (an integer between 1 and 256) indicates the
+increment amount for C<P> parameter selection.  The default value of 1 yields
+the parameter selection described in L</is_extra_strong_lucas_pseudoprime>,
+creating a pseudoprime sequence which is a superset of the latter's
+pseudoprime sequence L<OEIS A217719|http://oeis.org/A217719>.
+A value of 2 yields the method used by
+L<Pari|http://pari.math.u-bordeaux.fr/faq.html#primetest>.
+
+Because the C<U = 0> condition is ignored, this produces about 5% more
+pseudoprimes than the extra-strong Lucas test.  However this is still only
+66% of the number produced by the strong Lucas-Selfridge test.  No BPSW
+counterexamples have been found with any of the Lucas tests described.
 
 =head2 is_frobenius_underwood_pseudoprime
 
@@ -3066,7 +3088,7 @@ This has been verified with Jan Feitsma's 2-PSP database 
to produce no false
 results for 64-bit inputs.  Hence the result will always be 0 (composite) or
 2 (prime).
 
-For inputs larger than C<2^64>, a strong Baillie-PSW primality test is
+For inputs larger than C<2^64>, an extra-strong Baillie-PSW primality test is
 performed (also called BPSW or BSW).  This is a probabilistic test, so only
 0 (composite) and 1 (probably prime) are returned.  There is a possibility that
 composites may be returned marked prime, but since the test was published in
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index e868d73..00e2f2a 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -150,7 +150,7 @@ sub _is_prime7 {  # n must not be divisible by 2, 3, or 5
               !($n % 37) || !($n % 41) || !($n % 43) || !($n % 47) ||
               !($n % 53) || !($n % 59);
 
-  if ($n <= 10_000_000) {
+  if ($n <= 1_000_000) {
     my $limit = int(sqrt($n));
     my $i = 61;
     while (($i+30) <= $limit) {
@@ -176,18 +176,18 @@ sub _is_prime7 {  # n must not be divisible by 2, 3, or 5
     return 2;
   }
 
-  if ($n < 105936894253) {   # BPSW seems to be faster after this
+  if ($n < 47636622961201) {   # BPSW seems to be faster after this
     # Deterministic set of Miller-Rabin tests.  If the MR routines can handle
     # bases greater than n, then this can be simplified.
     my @bases;
-    if    ($n <          9080191) { @bases = (31,       73); }
-    elsif ($n <         19471033) { @bases = ( 2,   299417); }
+    # n > 1_000_000 because of the previous block.
+    if    ($n <         19471033) { @bases = ( 2,  299417); }
     elsif ($n <         38010307) { @bases = ( 2,  9332593); }
     elsif ($n <        316349281) { @bases = ( 11000544, 31481107); }
     elsif ($n <       4759123141) { @bases = ( 2, 7, 61); }
-    elsif ($n <     105936894253) { @bases = ( 2, 1005905886, 1340600841); }
-    elsif ($n <   31858317218647) { @bases = ( 2, 642735, 553174392, 
3046413974); }
-    elsif ($n < 3071837692357849) { @bases = ( 2, 75088, 642735, 203659041, 
3613982119); }
+    elsif ($n <     154639673381) { @bases = ( 15, 176006322, 4221622697); }
+    elsif ($n <   47636622961201) { @bases = ( 2, 2570940, 211991001, 
3749873356); }
+    elsif ($n < 3770579582154547) { @bases = ( 2, 2570940, 880937, 610386380, 
4130785767); }
     else                          { @bases = ( 2, 325, 9375, 28178, 450775, 
9780504, 1795265022); }
     return miller_rabin($n, @bases)  ?  2  :  0;
   }
@@ -195,10 +195,10 @@ sub _is_prime7 {  # n must not be divisible by 2, 3, or 5
   # BPSW probable prime.  No composites are known to have passed this test
   # since it was published in 1980, though we know infinitely many exist.
   # It has also been verified that no 64-bit composite will return true.
-  # Slow since it's all in PP, but it's the Right Thing To Do.
+  # Slow since it's all in PP and uses bigints.
 
   return 0 unless miller_rabin($n, 2);
-  return 0 unless is_strong_lucas_pseudoprime($n);
+  return 0 unless is_extra_strong_lucas_pseudoprime($n);
   return ($n <= 18446744073709551615)  ?  2  :  1;
 }
 
@@ -927,7 +927,8 @@ sub _lucas_selfridge_params {
 }
 
 sub _lucas_extrastrong_params {
-  my($n) = @_;
+  my($n, $increment) = @_;
+  $increment = 1 unless defined $increment;
 
   my ($P, $Q, $D) = (3, 1, 5);
   while (1) {
@@ -935,7 +936,8 @@ sub _lucas_extrastrong_params {
                                           : _gcd_ui($D, $n);
     return (0,0,0) if $gcd > 1 && $gcd != $n;  # Found divisor $d
     last if _jacobi($D, $n) == -1;
-    $P++;
+    $P += $increment;
+    croak "Could not find Jacobi sequence for $n" if $P > 65535;
     $D = $P*$P - 4;
   }
   ($P, $Q, $D);
@@ -945,13 +947,12 @@ sub _lucas_extrastrong_params {
 sub lucas_sequence {
   my($n, $P, $Q, $k) = @_;
 
-  return (0, 2) if $k == 0;
-  my $D = $P*$P - 4*$Q;
-
   croak "lucas_sequence: n must be >= 2" if $n < 2;
   croak "lucas_sequence: k must be >= 0" if $k < 0;
   croak "lucas_sequence: P out of range" if $P < 0 || $P >= $n;
   croak "lucas_sequence: Q out of range" if $Q >= $n;
+
+  my $D = $P*$P - 4*$Q;
   croak "lucas_sequence: D is zero" if $D == 0;
 
   if (ref($n) ne 'Math::BigInt') {
@@ -967,25 +968,36 @@ sub lucas_sequence {
   my $V = $ZERO + $P;
   my $Qk = $ZERO + $Q;
 
+  return ($ZERO, $ZERO+2) if $k == 0;
   $k = Math::BigInt->new("$k") unless ref($k) eq 'Math::BigInt';
   my $kstr = substr($k->as_bin, 2);
   my $bpos = 0;
 
   if ($Q == 1) {
-    while (++$bpos < length($kstr)) {
-      $U = ($U * $V) % $n;
-      $V = ($V * $V - 2) % $n;
-      if (substr($kstr,$bpos,1)) {
-        my $T1 = $U->copy->bmul($D);
-        $U->bmul($P);
-        $U->badd($V);
-        $U->badd($n) if $U->is_odd;
-        $U->brsft(1);
-
-        $V->bmul($P);
-        $V->badd($T1);
-        $V->badd($n) if $V->is_odd;
-        $V->brsft(1);
+    my $Dinverse = ($ZERO+$D)->bmodinv($n);
+    if ($P > 2 && !$Dinverse->is_nan) {
+      # Calculate V_k with U=V_{k+1}
+      $U = $ZERO + ($P*$P - 2);
+      while (++$bpos < length($kstr)) {
+        if (substr($kstr,$bpos,1)) {
+          $V->bmul($U)->bsub($P)->bmod($n);
+          $U->bmul($U)->bsub( 2)->bmod($n);
+        } else {
+          $U->bmul($V)->bsub($P)->bmod($n);
+          $V->bmul($V)->bsub( 2)->bmod($n);
+        }
+      }
+      # Crandall and Pomerance eq 3.13: U_n = D^-1 (2V_{n+1} - PV_n)
+      $U = $Dinverse * (2*$U - $P*$V);
+    } else {
+      while (++$bpos < length($kstr)) {
+        $U = ($U * $V) % $n;
+        $V = ($V * $V - 2) % $n;
+        if (substr($kstr,$bpos,1)) {
+          my $T1 = $U->copy->bmul($D);
+          $U->bmul($P)->badd( $V)->badd( $U->is_odd ? $n : 0 )->brsft(1);
+          $V->bmul($P)->badd($T1)->badd( $V->is_odd ? $n : 0 )->brsft(1);
+        }
       }
     }
   } else {
@@ -1090,6 +1102,59 @@ sub is_extra_strong_lucas_pseudoprime {
   return 0;
 }
 
+sub is_almost_extra_strong_lucas_pseudoprime {
+  my($n, $increment) = @_;
+  _validate_positive_integer($n);
+  $increment = 1 unless defined $increment;
+  _validate_positive_integer($increment, 1, 256);
+
+  return 1 if $n == 2;
+  return 0 if $n < 2 || ($n % 2) == 0;
+  return 0 if _is_perfect_square($n);
+
+  my ($P, $Q, $D) = _lucas_extrastrong_params($n, $increment);
+  return 0 if $D == 0;  # We found a divisor in the sequence
+  die "Lucas parameter error: $D, $P, $Q\n" if ($D != $P*$P - 4*$Q);
+
+  my $m = $n+1;
+  my($s, $k) = (0, $m);
+  while ( $k > 0 && !($k % 2) ) {
+    $s++;
+    $k >>= 1;
+  }
+
+  if (ref($n) ne 'Math::BigInt') {
+    if (!defined $Math::BigInt::VERSION) {
+      eval { require Math::BigInt;  Math::BigInt->import(try=>'GMP,Pari'); 1; }
+      or do { croak "Cannot load Math::BigInt "; }
+    }
+    $n = Math::BigInt->new("$n");
+  }
+
+  my $ZERO = $n->copy->bzero;
+  my $V = $ZERO + $P;        # V_{k}
+  my $W = $ZERO + $P*$P-2;   # V_{k+1}
+  $k = Math::BigInt->new("$k") unless ref($k) eq 'Math::BigInt';
+  my $kstr = substr($k->as_bin, 2);
+  my $bpos = 0;
+  while (++$bpos < length($kstr)) {
+    if (substr($kstr,$bpos,1)) {
+      $V->bmul($W)->bsub($P)->bmod($n);
+      $W->bmul($W)->bsub( 2)->bmod($n);
+    } else {
+      $W->bmul($V)->bsub($P)->bmod($n);
+      $V->bmul($V)->bsub( 2)->bmod($n);
+    }
+  }
+
+  return 1 if $V == 2 || $V == ($n-2);
+  foreach my $r (0 .. $s-2) {
+    return 1 if $V->is_zero;
+    $V->bmul($V)->bsub(2)->bmod($n);
+  }
+  return 0;
+}
+
 sub is_frobenius_underwood_pseudoprime {
   my($n) = @_;
   _validate_positive_integer($n);
@@ -2724,8 +2789,9 @@ functions.  Alternately, L<Math::Pari> has a lot of these 
types of functions.
   print "$n is prime" if is_prime($n);
 
 Returns 2 if the number is prime, 0 if not.  For numbers larger than C<2^64>
-it will return 0 for composite and 1 for probably prime, using a strong BPSW
-test.  Also note there are probabilistic prime testing functions available.
+it will return 0 for composite and 1 for probably prime, using an
+extra-strong BPSW test.  Also note there are probabilistic prime testing
+functions available.
 
 
 =head2 primes
@@ -2842,6 +2908,20 @@ L<Grantham 
2000|http://www.ams.org/mathscinet-getitem?mr=1680879>).
 This has slightly more restrictive conditions than the strong Lucas test,
 but uses different starting parameters so is not directly comparable.
 
+=head2 is_almost_extra_strong_lucas_pseudoprime
+
+Takes a positive number as input, and returns 1 if the input passes the extra
+strong Lucas test, ignoring the C<U_n = 0> condition.  This produces about 5%
+more pseudoprimes than the extra strong test, but 66% fewer than the strong
+test.
+
+An optional second argument (an integer between 1 and 256) indicates the
+increment used when selecting the C<P> parameter.  A default value of 1
+creates the same parameters as the extra strong test, so the pseudoprime
+sequence from this function will contain all the extra strong Lucas
+pseudoprimes.  A value of 2 yields the method used by
+L<Pari|http://pari.math.u-bordeaux.fr/faq.html#primetest>.
+
 =head2 is_frobenius_underwood_pseudoprime
 
 Takes a positive number as input, and returns 1 if the input passes the minimal

-- 
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