[libmath-prime-util-perl] 26/59: Bigint enhancements
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:56 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list