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

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

commit 3e3b780a5a66811381d7bd314b7a8a9f4840b54b
Author: Dana Jacobsen <d...@acm.org>
Date:   Tue Jul 3 15:26:23 2012 -0600

    Bigint enhancements
---
 examples/bench-isprime-bpsw.pl |  34 ++++++++----
 lib/Math/Prime/Util.pm         |  14 +++--
 lib/Math/Prime/Util/PP.pm      | 119 ++++++++++++++++++++++-------------------
 t/17-pseudoprime.t             |   4 +-
 t/80-pp.t                      |   4 +-
 5 files changed, 100 insertions(+), 75 deletions(-)

diff --git a/examples/bench-isprime-bpsw.pl b/examples/bench-isprime-bpsw.pl
index 86c00f6..f58cd70 100755
--- a/examples/bench-isprime-bpsw.pl
+++ b/examples/bench-isprime-bpsw.pl
@@ -7,20 +7,34 @@ use Math::Prime::Util;
 use Math::Primality;
 
 srand(500);
-use bigint try=>'GMP';
+#use bigint lib=>'';
 use Math::BigInt::Random::OO;
-#my $gen = Math::BigInt::Random::OO -> new(length => 50);
-my $gen = Math::BigInt::Random::OO -> new(length => 25);
+my $gen = Math::BigInt::Random::OO -> new(length => 180);
+#my $gen = Math::BigInt::Random::OO -> new(length => 8);
 
 my @rns;
-push @rns, $gen->generate() for (1 .. 100);
+push @rns, $gen->generate() for (1 .. 50);
+#my @rns;
+#push @rns, 1000000000 + int(rand(1000000000)) for (1..100);
+
+print "Verifying";
+for my $n (@rns) {
+  die "bad MR for $n" unless Math::Prime::Util::PP::miller_rabin($n,2) == 
Math::Primality::is_strong_pseudoprime("$n","2");
+  die "bad LP for $n" unless 
Math::Prime::Util::PP::is_strong_lucas_pseudoprime($n) == 
Math::Primality::is_strong_lucas_pseudoprime("$n");
+  die "bad IP for $n" unless (Math::Prime::Util::PP::is_prime($n)?1:0) == 
(Math::Primality::is_prime("$n")?1:0);
+  print ".";
+}
+print "OK\n";
 
 use Benchmark qw/:all/;
+my $sum = 0;
 cmpthese(-.5, {
-  "MP  MR" => sub { Math::Primality::is_strong_pseudoprime("$_","2") for @rns; 
},
-  "MPU MR" => sub { Math::Prime::Util::PP::miller_rabin($_,2) for @rns; },
-  "MP  LP" => sub { Math::Primality::is_strong_lucas_pseudoprime("$_") for 
@rns;},
-  "MPU LP" => sub { Math::Prime::Util::PP::is_strong_lucas_pseudoprime($_) for 
@rns;},
-  "MP  IP" => sub { Math::Primality::is_prime("$_") for @rns;},
-  "MPU IP" => sub { Math::Prime::Util::PP::is_prime($_) for @rns;},
+  "MP  MR" => sub { $sum += Math::Primality::is_strong_pseudoprime("$_","2") 
for @rns; },
+  "MPU MR" => sub { $sum += Math::Prime::Util::PP::miller_rabin($_,2) for 
@rns; },
+  #"MPUxMR" => sub { Math::Prime::Util::miller_rabin($_,2) for @rns; },
+  "MP  LP" => sub { $sum += Math::Primality::is_strong_lucas_pseudoprime("$_") 
for @rns;},
+  "MPU LP" => sub { $sum += 
Math::Prime::Util::PP::is_strong_lucas_pseudoprime($_) for @rns;},
+  "MP  IP" => sub { $sum += Math::Primality::is_prime("$_") for @rns;},
+  "MPU IP" => sub { $sum += Math::Prime::Util::PP::is_prime($_) for @rns;},
+  #"MPUxIP" => sub { Math::Prime::Util::is_prime($_) for @rns;},
 });
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index edc3cf7..0578f0d 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1492,9 +1492,12 @@ The differences are in the implementations:
      C<prime_precalc>).  Larger inputs just need too much time and memory
      for the sieve.
 
-   - L<Math::Primality> uses a GMP BPSW test which is overkill for our 64-bit
-     range.  It's generally an order of magnitude or two slower than any
-     of the others.
+   - L<Math::Primality> uses GMP for all work.  Under ~32-bits it uses 2 or 3
+     MR tests, while above 4759123141 it performs a BPSW test.  This is is
+     fantastic for bigints over 2^64, but it is significantly slower than
+     native precision tests.  With 64-bit numbers it is generally an order of
+     magnitude or more slower than any of the others.  This reverses when
+     numbers get larger.
 
    - L<Math::Pari> has some very effective code, but it has some overhead to 
get
      to it from Perl.  That means for small numbers it is relatively slow: an
@@ -1509,8 +1512,9 @@ The differences are in the implementations:
    - L<Math::Prime::Util> looks in the sieve for a fast bit lookup if that
      exists (default up to 30,000 but it can be expanded, e.g.
      C<prime_precalc>), uses trial division for numbers higher than this but
-     not too large (0.1M on 64-bit machines, 100M on 32-bit machines), and a
-     deterministic set of Miller-Rabin tests for large numbers.
+     not too large (0.1M on 64-bit machines, 100M on 32-bit machines), a
+     deterministic set of Miller-Rabin tests for 64-bit and smaller numbers,
+     and a BPSW test for bigints.
 
 
 
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index bb1b98a..ff39e71 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -609,15 +609,17 @@ sub _find_jacobi_d_sequence {
   my($n) = @_;
 
   # D is typically quite small: 67 max for N < 10^19.  However, it is
-  # theoretically possible D could grow unreasonably.  We are ignoring this.
+  # theoretically possible D could grow unreasonably.  I'm giving up at 4000M.
   my $d = 5;
   my $sign = 1;
   while (1) {
     my $gcd = _gcd_ui($d, $n);
+    #my $gcd = Math::BigInt::bgcd($d, $n);
     return 0 if $gcd > 1 && $gcd != $n;  # Found divisor $d
     my $j = _jacobi($d * $sign, $n);
     last if $j == -1;
     $d += 2;
+    croak "Could not find Jacobi sequence for $n" if $d > 4_000_000_000;
     $sign = -$sign;
   }
   return ($sign * $d);
@@ -628,6 +630,10 @@ sub is_strong_lucas_pseudoprime {
   my($n) = @_;
   _validate_positive_integer($n);
 
+  # We're trying to limit the bignum calculations as much as possible.
+  # It's also important to try to match whatever they passed in.  For instance
+  # if they use a GMP or Pari object, we must do the same.
+
   return 1 if $n == 2;
   return 0 if $n < 2 || ($n % 2) == 0;
 
@@ -636,27 +642,25 @@ sub is_strong_lucas_pseudoprime {
   #     Math::Primality
 
   # Check for perfect square
-  {
+  if (ref($n) eq 'Math::BigInt') {
+    my $mcheck = ($n & 127)->numify;
+    if (($mcheck*0x8bc40d7d) & ($mcheck*0xa1e2f5d1) & 0x14020a) {
+      # ~82% of non-squares were rejected by the bloom filter
+      my $sq = int($n->copy()->bsqrt());
+      return 0 if ($sq*$sq) == $n;
+    }
+  } else {
     my $mcheck = $n & 127;
     if (($mcheck*0x8bc40d7d) & ($mcheck*0xa1e2f5d1) & 0x14020a) {
-      # ... 82% of non-squares were rejected by the bloom filter
-      # For bigints we should put in the rest of this filter, which prunes
-      # 99.92% of nonsquares, for the cost of one bigint mod and some
-      # non-bigint operations.
       my $sq = int(sqrt($n));
       return 0 if ($sq*$sq) == $n;
     }
   }
 
-  # We have to make sure we use bigints
-  my $result = 0;
- {
-  use bigint;
-  $n = Math::BigInt->new($n) unless ref($n) eq 'Math::BigInt';
-
   # Determine Selfridge D, P, and Q parameters
   my $D = _find_jacobi_d_sequence($n);
-  $D = $D->numify if $D <= ~0 && ref($D) eq 'Math::BigInt';
+  #die if ref($D) eq 'Math::BigInt';
+  #$D = $D->numify if $D <= ~0 && ref($D) eq 'Math::BigInt';
   return 0 if $D == 0;  # We found a divisor in the sequence
   my $P = 1;
   my $Q = int( (1 - $D) / 4 );
@@ -664,11 +668,16 @@ sub is_strong_lucas_pseudoprime {
   die "Selfridge error: $D, $P, $Q\n" if ($D != $P*$P - 4*$Q);
   #warn "N: $n  D: $D  P: $P  Q: $Q\n";
 
-  # It's now time to perform the Lucase pseudoprimality test using $D.
+  # It's now time to perform the Lucas pseudoprimality test using $D.
 
-  my $m = $n + 1;
+  if (ref($n) ne 'Math::BigInt') {
+    require Math::BigInt;
+    $n = Math::BigInt->new($n);
+  }
+
+  my $m = $n->copy() + 1;
   my $s = 0;
-  my $d = $m;
+  my $d = $m->copy();
   while ( ($d & 1) == 0 ) {
     $s++;
     $d >>= 1;
@@ -676,61 +685,59 @@ sub is_strong_lucas_pseudoprime {
   # $m == $d * 2 ** $s
   #die "Invalid $m, $d, $s\n" unless $m == $d * 2**$s;
 
-  my $U = 1;
-  my $V = $P;
-  my $U_2m = 1;
-  my $V_2m = $P;
-  my $Q_m = $Q;
-  my $Q_m2 = 2 * $Q;
-  my $Qkd = $Q;
+  my $ZERO = $n->copy->bzero();
+  my $U = $ZERO + 1;
+  my $V = $ZERO + $P;
+  my $U_2m = $ZERO + 1;
+  my $V_2m = $ZERO + $P;
+  my $Q_m = $ZERO + $Q;
+  my $Q_m2 = ($ZERO + $Q) * 2;
+  my $Qkd = $ZERO + $Q;
   $d >>= 1;
+  $d = $d->numify if $d <= ~0 && ref($d) eq 'Math::BigInt';
   #my $i = 1;
   while ($d != 0) {
     #warn "U=$U  V=$V  Qm=$Q_m  Qm2=$Q_m2\n";
-    $U_2m = ($U_2m * $V_2m) % $n;
-    $V_2m = ($V_2m * $V_2m - $Q_m2) % $n;
+    $U_2m->bmul($V_2m);             $U_2m->bmod($n);
+    $V_2m->bmuladd($V_2m, -$Q_m2);  $V_2m->bmod($n);
     #warn "  $i  U2m=$U_2m  V2m=$V_2m\n";  $i++;
-    $Q_m  = ($Q_m * $Q_m) % $n;
-    $Q_m2 = 2 * $Q_m;   # no mod
+    $Q_m->bmul($Q_m);             $Q_m->bmod($n);
+    $Q_m2 = 2 * $Q_m; # don't bother to mod n
     if ($d & 1) {
-      my $Uold = $U;
-      $U = $U_2m * $V  +  $V_2m * $U;
-      $U += $n if $U & 1;
-      $U = int($U / 2) % $n;
-
-      $V = $V_2m * $V  +  $U_2m * $Uold * $D;
-      $V += $n if $V & 1;
-      $V = int($V / 2) % $n;
-
-      $Qkd = ($Qkd * $Q_m) % $n;
+      my $Uold = $U->copy();
+      $U->bmuladd($V_2m, $U_2m * $V);  # U = U*V_2m + V*U_2m
+      $U->badd($n) if $U->is_odd;      # U += n if U & 1
+      $U->brsft(1,2);                  # U = floor(U / 2)
+      $U->bmod($n);                    # U = U % n
+
+      $V->bmuladd($V_2m, $U_2m * $Uold * $D);
+      $V->badd($n) if $V->is_odd;
+      $V->brsft(1,2);
+      $V->bmod($n);
+
+      $Qkd->bmul($Q_m);
+      $Qkd->bmod($n);
     }
     $d >>= 1;
   }
   #warn "l0 U=$U  V=$V\n";
-  if ($U == 0 || $V == 0) {
-    $result = 1;
-  } else {
-    # Compute powers of V
-    my $Qkd2 = 2 * $Qkd;
-    foreach my $r (1 .. $s-1) {
-      #warn "l$r U=$U  V=$V  Qkd2=$Qkd2\n";
-      $V = ($V * $V - $Qkd2) % $n;
-      if ($V == 0) {
-        $result = 1;
-        last;
-      }
-      if ($r < ($s-1)) {
-        $Qkd = ($Qkd * $Qkd) % $n;
-        $Qkd2 = 2 * $Qkd;
-      }
+  return 1 if $U == 0 || $V == 0;
+
+  # Compute powers of V
+  my $Qkd2 = 2 * $Qkd;
+  foreach my $r (1 .. $s-1) {
+    #warn "l$r U=$U  V=$V  Qkd2=$Qkd2\n";
+    $V->bmuladd($V, -$Qkd2);  $V->bmod($n);
+    return 1 if $V == 0;
+    if ($r < ($s-1)) {
+      $Qkd->bmul($Qkd);  $Qkd->bmod($n);
+      $Qkd2 = 2 * $Qkd;
     }
   }
-  no bigint;
- }
   #warn "Math::BigInt is loaded\n" if defined $Math::BigInt::VERSION;
   #warn "bigint is loaded\n" if defined $bigint::VERSION;
   #warn "bigint is in effect\n" if bigint::in_effect();
-  return $result;
+  return 0;
 }
 
 
diff --git a/t/17-pseudoprime.t b/t/17-pseudoprime.t
index f2f04ed..b842a94 100644
--- a/t/17-pseudoprime.t
+++ b/t/17-pseudoprime.t
@@ -77,8 +77,8 @@ ok(!eval { miller_rabin(2047,1); }, "MR base 1 fails");
 
 is( miller_rabin(0, 2), 0, "MR with 0 shortcut composite");
 is( miller_rabin(1, 2), 0, "MR with 0 shortcut composite");
-is( miller_rabin(2, 2), 2, "MR with 2 shortcut prime");
-is( miller_rabin(3, 2), 2, "MR with 3 shortcut prime");
+is( miller_rabin(2, 2), 1, "MR with 2 shortcut prime");
+is( miller_rabin(3, 2), 1, "MR with 3 shortcut prime");
 
 
 # Check that each strong pseudoprime base b makes it through MR with that base
diff --git a/t/80-pp.t b/t/80-pp.t
index 99c5888..79d4f52 100644
--- a/t/80-pp.t
+++ b/t/80-pp.t
@@ -332,8 +332,8 @@ while (my($n, $nth) = each (%nthprimes_small)) {
 
 is( miller_rabin(0, 2), 0, "MR with 0 shortcut composite");
 is( miller_rabin(1, 2), 0, "MR with 0 shortcut composite");
-is( miller_rabin(2, 2), 2, "MR with 2 shortcut prime");
-is( miller_rabin(3, 2), 2, "MR with 3 shortcut prime");
+is( miller_rabin(2, 2), 1, "MR with 2 shortcut prime");
+is( miller_rabin(3, 2), 1, "MR with 3 shortcut prime");
 
 while (my($base, $ppref) = each (%pseudoprimes)) {
   foreach my $p (@$ppref) {

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