[libmath-prime-util-perl] 126/181: Performance updates for no-GMP, focusing on making test suite run faster.
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:13 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.36
in repository libmath-prime-util-perl.
commit dbf2bd304023c5fe8f3308fd39bd32a49cdecc9d
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Jan 6 19:26:10 2014 -0800
Performance updates for no-GMP, focusing on making test suite run faster.
---
lib/Math/Prime/Util.pm | 12 +-
lib/Math/Prime/Util/PP.pm | 282 ++++++++++++++++++++++++++++------------------
t/23-primality-proofs.t | 2 +
t/81-bignum.t | 45 +++++---
4 files changed, 209 insertions(+), 132 deletions(-)
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index deb0b58..3655e0b 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1193,7 +1193,11 @@ sub primorial {
my $pn = (ref($_[0]) eq 'Math::BigInt') ? $_[0]->copy->bone()
: ($n >= $max) ? Math::BigInt->bone()
: 1;
- forprimes { $pn *= $_ } $n;
+ if (ref($pn) eq 'Math::BigInt') {
+ $pn->bmul($_) for map { Math::BigInt->new($_) } @{primes(2,$n)};
+ } else {
+ forprimes { $pn *= $_ } $n;
+ }
return $pn;
}
@@ -1758,9 +1762,9 @@ sub is_provable_prime_with_cert {
my $header = "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\n";
if ($n <= $_XS_MAXVAL) {
- my $isp = is_prime("$n");
+ my $isp = is_prime($n);
return ($isp, '') unless $isp == 2;
- return (2, "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\nType Small\nN $n\n");
+ return (2, $header . "Type Small\nN $n\n");
}
if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::is_provable_prime_with_cert) {
@@ -1783,7 +1787,7 @@ sub is_provable_prime_with_cert {
{
my $isp = is_prob_prime($n);
return ($isp, '') if $isp == 0;
- return (2, "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\nType Small\nN $n\n") if $isp == 2;
+ return (2, $header . "Type Small\nN $n\n") if $isp == 2;
}
# Choice of methods for proof:
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index b855b17..c75817e 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -156,6 +156,7 @@ sub _is_prime7 { # n must not be divisible by 2, 3, or 5
}
if ($n <= 1_000_000) {
+ $n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt';
my $limit = int(sqrt($n));
my $i = 61;
while (($i+30) <= $limit) {
@@ -320,12 +321,10 @@ sub _sieve_segment {
foreach my $s (split("0", substr($$primesieveref, 3), -1)) {
$p += 2 + 2 * length($s);
my $p2 = $p*$p;
- last if $p2 > $end;
if ($p2 < $beg) {
- $p2 = int($beg / $p) * $p;
- $p2 += $p if $p2 < $beg;
- $p2 += $p if ($p2 % 2) == 0; # Make sure p2 is odd
- }
+ my $f = 1+int(($beg-1)/$p);
+ $p2 = $p * ($f + !($f & 1));
+ } elsif ($p2 > $end) { last; }
# With large bases and small segments, it's common to find we don't hit
# the segment at all. Skip all the setup if we find this now.
if ($p2 <= $end) {
@@ -442,15 +441,32 @@ sub next_prime {
sub prev_prime {
my($n) = @_;
_validate_positive_integer($n);
- if ($n <= 7) {
- return ($n <= 2) ? 0 : ($n <= 3) ? 2 : ($n <= 5) ? 3 : 5;
+ if ($n <= 11) {
+ return ($n <= 2) ? 0 : ($n <= 3) ? 2 : ($n <= 5) ? 3 : ($n <= 7) ? 5 : 7;
}
- $n++ if ($n % 2) == 0;
- do {
+ #$n++ if ($n % 2) == 0;
+ #do {
+ # $n -= 2;
+ #} while ( (($n % 3) == 0) || (($n % 5) == 0) || (!_is_prime7($n)) );
+ #return $n;
+
+ $n -= ($n & 1) ? 2 : 1;
+ my $nmod6 = $n % 6;
+ if ($nmod6 == 5) {
+ return $n if $n % 5 && $n % 7 && _is_prime7($n);
+ $n -= 4;
+ } elsif ($nmod6 == 3) {
$n -= 2;
- } while ( (($n % 3) == 0) || (($n % 5) == 0) || (!_is_prime7($n)) );
- $n;
+ }
+
+ while (1) {
+ return $n if $n % 5 && $n % 7 && _is_prime7($n);
+ $n -= 2;
+ return $n if $n % 5 && $n % 7 && _is_prime7($n);
+ $n -= 4;
+ }
+ return $n;
# This is faster for larger intervals, slower for short ones.
#my $base = 30 * int($n/30);
@@ -873,17 +889,19 @@ sub _mulmod {
$y %= $n if $y >= $n;
($x,$y) = ($y,$x) if $x < $y;
if ($n <= (~0 >> 1)) {
- while ($y > 0) {
+ while ($y > 1) {
if ($y & 1) { $r += $x; $r -= $n if $r >= $n; }
$y >>= 1;
- if ($y) { $x += $x; $x -= $n if $x >= $n; }
+ $x += $x; $x -= $n if $x >= $n;
}
+ if ($y & 1) { $r += $x; $r -= $n if $r >= $n; }
} else {
- while ($y > 0) {
+ while ($y > 1) {
if ($y & 1) { $r = $n-$r; $r = ($x >= $r) ? $x-$r : $n-$r+$x; }
$y >>= 1;
- if ($y) { $x = ($x > ($n - $x)) ? ($x - $n) + $x : $x + $x; }
+ $x = ($x > ($n - $x)) ? ($x - $n) + $x : $x + $x;
}
+ if ($y & 1) { $r = $n-$r; $r = ($x >= $r) ? $x-$r : $n-$r+$x; }
}
$r;
}
@@ -994,7 +1012,8 @@ sub is_strong_pseudoprime {
_validate_positive_integer($n);
croak "No bases given to miller_rabin" unless @bases;
- return $n >> 1 if $n <= 3;
+ return 0+($n >= 2) if $n < 4;
+ return 0 if ($n % 2) == 0;
# Die on invalid bases
foreach my $base (@bases) { croak "Base $base is invalid" if $base < 2 }
@@ -1003,14 +1022,13 @@ sub is_strong_pseudoprime {
if ( ref($n) eq 'Math::BigInt' ) {
- return 0 if $n->is_even;
my $nminus1 = $n->copy->bdec();
my $s = 0;
my $d = $nminus1->copy;
- while ($d->is_even) {
+ do { # n is > 3 and odd, so n-1 must be even
$s++;
$d->brsft(1);
- }
+ } while $d->is_even;
# Different way of doing the above. Fewer function calls, slower on ave.
#my $dbin = $nminus1->as_bin;
#my $last1 = rindex($dbin, '1');
@@ -1019,7 +1037,7 @@ sub is_strong_pseudoprime {
foreach my $ma (@bases) {
my $x = $n->copy->bzero->badd($ma)->bmodpow($d,$n);
- next if ($x->is_one) || ($x->bcmp($nminus1) == 0);
+ next if $x->is_one || $x->bcmp($nminus1) == 0;
foreach my $r (1 .. $s-1) {
$x->bmul($x); $x->bmod($n);
return 0 if $x->is_one;
@@ -1030,7 +1048,6 @@ sub is_strong_pseudoprime {
} else {
- return 0 unless $n & 1;
my $s = 0;
my $d = $n - 1;
while ( ($d & 1) == 0 ) {
@@ -1214,22 +1231,22 @@ sub lucas_sequence {
$n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
my $ZERO = $n->copy->bzero;
- my $ONE = $ZERO + 1;
- my $TWO = $ZERO + 2;
+ my $ONE = $ZERO->copy->binc;
+ my $TWO = $ONE->copy->binc;
$P = $ZERO+$P unless ref($P) eq 'Math::BigInt';
$Q = $ZERO+$Q unless ref($Q) eq 'Math::BigInt';
- my $D = $P*$P - 4*$Q;
- croak "lucas_sequence: D is zero" if $D == 0;
+ my $D = $P*$P - $TWO*$TWO*$Q;
+ croak "lucas_sequence: D is zero" if $D->is_zero;
my $U = $ONE->copy;
my $V = $P->copy;
my $Qk = $Q->copy;
- return ($ZERO, $ZERO+2) if $k == 0;
+ return ($ZERO, $TWO) 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) {
+ if ($Q->is_one) {
my $Dinverse = $D->copy->bmodinv($n);
if ($P > $TWO && !$Dinverse->is_nan) {
# Calculate V_k with U=V_{k+1}
@@ -1295,9 +1312,8 @@ sub is_lucas_pseudoprime {
return 0 if defined $n && int($n) < 0;
_validate_positive_integer($n);
- return 1 if $n == 2;
- return 0 if $n < 2 || ($n % 2) == 0;
- return 0 if _is_perfect_square($n);
+ return 0+($n >= 2) if $n < 4;
+ return 0 if ($n % 2) == 0 || _is_perfect_square($n);
my ($P, $Q, $D) = _lucas_selfridge_params($n);
return 0 if $D == 0; # We found a divisor in the sequence
@@ -1312,9 +1328,8 @@ sub is_strong_lucas_pseudoprime {
return 0 if defined $n && int($n) < 0;
_validate_positive_integer($n);
- return 1 if $n == 2;
- return 0 if $n < 2 || ($n % 2) == 0;
- return 0 if _is_perfect_square($n);
+ return 0+($n >= 2) if $n < 4;
+ return 0 if ($n % 2) == 0 || _is_perfect_square($n);
my ($P, $Q, $D) = _lucas_selfridge_params($n);
return 0 if $D == 0; # We found a divisor in the sequence
@@ -1344,23 +1359,23 @@ sub is_extra_strong_lucas_pseudoprime {
return 0 if defined $n && int($n) < 0;
_validate_positive_integer($n);
- return 1 if $n == 2;
- return 0 if $n < 2 || ($n % 2) == 0;
- return 0 if _is_perfect_square($n);
+ return 0+($n >= 2) if $n < 4;
+ return 0 if ($n % 2) == 0 || _is_perfect_square($n);
my ($P, $Q, $D) = _lucas_extrastrong_params($n);
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;
- }
# We have to convert n to a bigint or Math::BigInt::GMP's stupid set_si bug
# (RT 71548) will hit us and make the test $V == $n-2 always return false.
$n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
+
+ my($s, $k) = (0, $n->copy->binc);
+ while ($k->is_even && !$k->is_zero) {
+ $s++;
+ $k->brsft(1);
+ }
+
my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k);
return 1 if $U->is_zero && ($V == 2 || $V == ($n-2));
@@ -1381,9 +1396,8 @@ sub is_almost_extra_strong_lucas_pseudoprime {
$increment = 1;
}
- return 1 if $n == 2;
- return 0 if $n < 2 || ($n % 2) == 0;
- return 0 if _is_perfect_square($n);
+ return 0+($n >= 2) if $n < 4;
+ return 0 if ($n % 2) == 0 || _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
@@ -1392,26 +1406,27 @@ sub is_almost_extra_strong_lucas_pseudoprime {
$n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
my $ZERO = $n->copy->bzero;
- my $V = $ZERO + $P; # V_{k}
- my $W = $ZERO + $P*$P-2; # V_{k+1}
+ my $TWO = $ZERO->copy->binc->binc;
+ my $V = $ZERO + $P; # V_{k}
+ my $W = $ZERO + $P*$P-$TWO; # V_{k+1}
my $kstr = substr($n->copy->binc()->as_bin, 2);
$kstr =~ s/(0*)$//;
my $s = length($1);
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);
+ $V->bmul($W)->bsub($P )->bmod($n);
+ $W->bmul($W)->bsub($TWO)->bmod($n);
} else {
- $W->bmul($V)->bsub($P)->bmod($n);
- $V->bmul($V)->bsub( 2)->bmod($n);
+ $W->bmul($V)->bsub($P )->bmod($n);
+ $V->bmul($V)->bsub($TWO)->bmod($n);
}
}
- return 1 if $V == 2 || $V == ($n-2);
+ return 1 if $V == 2 || $V == ($n-$TWO);
foreach my $r (0 .. $s-2) {
return 1 if $V->is_zero;
- $V->bmul($V)->bsub(2)->bmod($n);
+ $V->bmul($V)->bsub($TWO)->bmod($n);
}
return 0;
}
@@ -1420,14 +1435,13 @@ sub is_frobenius_underwood_pseudoprime {
my($n) = @_;
return 0 if defined $n && int($n) < 0;
_validate_positive_integer($n);
- return 0 if $n < 2;
- return 1 if $n < 4;
- return 0 if ($n % 2) == 0;
- return 0 if _is_perfect_square($n);
+ return 0+($n >= 2) if $n < 4;
+ return 0 if ($n % 2) == 0 || _is_perfect_square($n);
$n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
my $ZERO = $n->copy->bzero;
+ my $ONE = $ZERO->copy->binc;
my $fa = $ZERO + 1;
my $fb = $ZERO + 2;
@@ -1442,13 +1456,15 @@ sub is_frobenius_underwood_pseudoprime {
$result %= $n if $result > $n;
$multiplier %= $n if $multiplier > $n;
foreach my $bit (reverse 0 .. $len-2) {
- $na = $fa * (($fa*$x) + ($fb+$fb));
- $fb = ( ($fb + $fa) * ($fb - $fa) ) % $n;
- $fa = $na % $n;
- if ( ($np1 >> $bit) & 1 ) {
- $na = $fb + ($fa * $multiplier);
- $fb += ($fb - $fa);
- $fa = $na;
+ $na = ($fa * $x) + ($fb + $fb);
+ $t = $fb + $fa;
+ $fb->bsub($fa)->bmul($t)->bmod($n);
+ $fa->bmul($na)->bmod($n);
+
+ if ( ($np1 >> $bit) & $ONE ) {
+ $t = $fb->copy;
+ $fb->badd($fb)->bsub($fa);
+ $fa->bmul($multiplier)->badd($t);
}
}
$fa->bmod($n);
@@ -1486,25 +1502,22 @@ sub _poly_mod_mul {
my $px_degree = $#$px;
my $py_degree = $#$py;
- my @res;
+ my @res = map { $_poly_bignum ? Math::BigInt->bzero : 0 } 0 .. $r-1;
# convolve(px, py) mod (X^r-1,n)
my @indices_y = grep { $py->[$_] } (0 .. $py_degree);
- for (my $ix = 0; $ix <= $px_degree; $ix++) {
+ foreach my $ix (0 .. $px_degree) {
my $px_at_ix = $px->[$ix];
next unless $px_at_ix;
if ($_poly_bignum) {
foreach my $iy (@indices_y) {
my $rindex = ($ix + $iy) % $r; # reduce mod X^r-1
- $res[$rindex] = Math::BigInt->bzero unless defined $res[$rindex];
$res[$rindex]->badd($px_at_ix->copy->bmul($py->[$iy]))->bmod($n);
}
} else {
foreach my $iy (@indices_y) {
my $rindex = ($ix + $iy) % $r; # reduce mod X^r-1
- $res[$rindex] = 0 unless defined $res[$rindex];
- my $py_px = $px_at_ix * $py->[$iy];
- $res[$rindex] = ($res[$rindex] + $py_px) % $n;
+ $res[$rindex] = ($res[$rindex] + $px_at_ix * $py->[$iy]) % $n;
}
}
}
@@ -1540,21 +1553,24 @@ sub is_aks_prime {
return 0 if defined $n && int($n) < 0;
_validate_positive_integer($n);
- $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
-
return 0 if $n < 2;
return 0 if _is_perfect_power($n);
- do { require Math::BigFloat; Math::BigFloat->import(); }
- if !defined $Math::BigFloat::VERSION;
- # limit = floor( log2(n) * log2(n) ). o_r(n) must be larger than this
- my $floatn = Math::BigFloat->new($n);
- my $sqrtn = _bigint_to_int($floatn->copy->bsqrt->bfloor);
- # The following line seems to trigger a memory leak in Math::BigFloat::blog
- # (the part where $MBI is copied to $int) if $n is a Math::BigInt::GMP.
- my $log2n = $floatn->copy->blog(2);
- my $log2_squared_n = $log2n * $log2n;
- my $limit = _bigint_to_int($log2_squared_n->bfloor);
+ my($log2n, $limit);
+ if ($n > 2**48) {
+ do { require Math::BigFloat; Math::BigFloat->import(); }
+ if !defined $Math::BigFloat::VERSION;
+ # limit = floor( log2(n) * log2(n) ). o_r(n) must be larger than this
+ my $floatn = Math::BigFloat->new("$n");
+ #my $sqrtn = _bigint_to_int($floatn->copy->bsqrt->bfloor);
+ # The following line seems to trigger a memory leak in Math::BigFloat::blog
+ # (the part where $MBI is copied to $int) if $n is a Math::BigInt::GMP.
+ $log2n = $floatn->copy->blog(2);
+ $limit = _bigint_to_int( ($log2n * $log2n)->bfloor );
+ } else {
+ $log2n = log($n)/log(2) + 0.0001; # Error on large side.
+ $limit = int( $log2n*$log2n + 0.0001 );
+ }
my $r = next_prime($limit);
foreach my $f (@{primes(0,$r-1)}) {
@@ -1572,13 +1588,17 @@ sub is_aks_prime {
return 1 if $r >= $n;
# Since r is a prime, phi(r) = r-1
- my $rlimit = _bigint_to_int( Math::BigFloat->new("$r")->bdec()
- ->bsqrt->bmul($log2n)->bfloor);
+ my $rlimit = (ref($log2n) eq 'Math::BigFloat')
+ ? _bigint_to_int( Math::BigFloat->new("$r")->bdec()
+ ->bsqrt->bmul($log2n)->bfloor)
+ : int( (sqrt(($r-1)) * $log2n) + 0.001 );
$_poly_bignum = 1;
if ( $n < (MPU_HALFWORD-1) ) {
$_poly_bignum = 0;
- $n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt';
+ #$n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt';
+ } else {
+ $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
}
for (my $a = 1; $a <= $rlimit; $a++) {
@@ -1630,26 +1650,49 @@ sub trial_factor {
@factors = ($n == 1) ? () : ($n);
return @factors;
}
- while ( !($n % 2) ) { push @factors, 2; $n = int($n / 2); }
- while ( !($n % 3) ) { push @factors, 3; $n = int($n / 3); }
- while ( !($n % 5) ) { push @factors, 5; $n = int($n / 5); }
+ if (ref($n) ne 'Math::BigInt' || !Math::BigInt::bgcd($n, 30030)->is_one) {
+ while ( !($n % 2) ) { push @factors, 2; $n = int($n / 2); }
+ while ( !($n % 3) ) { push @factors, 3; $n = int($n / 3); }
+ while ( !($n % 5) ) { push @factors, 5; $n = int($n / 5); }
+ while ( !($n % 7) ) { push @factors, 7; $n = int($n / 7); }
+ while ( !($n %11) ) { push @factors,11; $n = int($n /11); }
+ while ( !($n %13) ) { push @factors,13; $n = int($n /13); }
+ }
$n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt' && $n <= ''.~0;
return @factors if $n < 4;
my $limit = int(sqrt($n) + 0.001);
$limit = $maxlim if $limit > $maxlim;
- my $f = 7;
- SEARCH: while ($f <= $limit) {
- foreach my $finc (4, 2, 4, 2, 4, 6, 2, 6) {
- if ( (($n % $f) == 0) && ($f <= $limit) ) {
- do { push @factors, $f; $n = int($n/$f); } while (($n % $f) == 0);
- $n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt' && $n <= ''.~0;
- #last SEARCH if $n == 1 || Math::Prime::Util::is_prob_prime($n);
- last SEARCH if $n == 1;
- $limit = int( sqrt($n) + 0.001);
- $limit = $maxlim if $limit > $maxlim;
+
+ if (ref($n) eq 'Math::BigInt') {
+ my $f = Math::BigInt->new(17);
+ $limit = Math::BigInt->new("$limit");
+ my @incs = map { Math::BigInt->new($_) } (2, 4, 6, 2, 6, 4, 2, 4);
+ SEARCH: while ($f <= $limit) {
+ foreach my $finc (@incs) {
+ if ($n->copy->bmod($f)->is_zero && $f->bacmp($limit) <= 0) {
+ my $sf = ($f <= ''.~0) ? _bigint_to_int($f) : $f;
+ do { push @factors, $sf; $n = int($n/$f); } while (($n % $f) == 0);
+ last SEARCH if $n->is_one;
+ $limit = int( sqrt($n) + 0.001);
+ $limit = $maxlim if $limit > $maxlim;
+ $limit = Math::BigInt->new("$limit");
+ }
+ $f->badd($finc);
+ }
+ }
+ } else {
+ my $f = 17;
+ SEARCH: while ($f <= $limit) {
+ foreach my $finc (2, 4, 6, 2, 6, 4, 2, 4) {
+ if ( (($n % $f) == 0) && ($f <= $limit) ) {
+ do { push @factors, $f; $n = int($n/$f); } while (($n % $f) == 0);
+ last SEARCH if $n == 1;
+ $limit = int( sqrt($n) + 0.001);
+ $limit = $maxlim if $limit > $maxlim;
+ }
+ $f += $finc;
}
- $f += $finc;
}
}
push @factors, $n if $n > 1;
@@ -1687,16 +1730,30 @@ sub factor {
my @factors;
# Use 'n=int($n/7)' instead of 'n/=7' to not "upgrade" n to a Math::BigFloat.
- while (($n % 2) == 0) { push @factors, 2; $n = int($n / 2); }
- while (($n % 3) == 0) { push @factors, 3; $n = int($n / 3); }
- while (($n % 5) == 0) { push @factors, 5; $n = int($n / 5); }
- while (($n % 7) == 0) { push @factors, 7; $n = int($n / 7); }
- while (($n % 11) == 0) { push @factors, 11; $n = int($n / 11); }
- while (($n % 13) == 0) { push @factors, 13; $n = int($n / 13); }
- while (($n % 17) == 0) { push @factors, 17; $n = int($n / 17); }
- while (($n % 19) == 0) { push @factors, 19; $n = int($n / 19); }
- while (($n % 23) == 0) { push @factors, 23; $n = int($n / 23); }
- while (($n % 29) == 0) { push @factors, 29; $n = int($n / 29); }
+ if (ref($n) eq 'Math::BigInt') {
+ while ($n->is_even) { push @factors, 2; $n->brsft(1); }
+ if (!Math::BigInt::bgcd($n, "3234846615")->is_one) {
+ foreach my $div (3, 5, 7, 11, 13, 17, 19, 23, 29) {
+ my ($q, $r) = $n->copy->bdiv($div);
+ while ($r->is_zero) {
+ push @factors, $div;
+ $n = $q;
+ ($q, $r) = $n->copy->bdiv($div);
+ }
+ }
+ }
+ } else {
+ while (($n % 2) == 0) { push @factors, 2; $n = int($n / 2); }
+ while (($n % 3) == 0) { push @factors, 3; $n = int($n / 3); }
+ while (($n % 5) == 0) { push @factors, 5; $n = int($n / 5); }
+ while (($n % 7) == 0) { push @factors, 7; $n = int($n / 7); }
+ while (($n % 11) == 0) { push @factors, 11; $n = int($n / 11); }
+ while (($n % 13) == 0) { push @factors, 13; $n = int($n / 13); }
+ while (($n % 17) == 0) { push @factors, 17; $n = int($n / 17); }
+ while (($n % 19) == 0) { push @factors, 19; $n = int($n / 19); }
+ while (($n % 23) == 0) { push @factors, 23; $n = int($n / 23); }
+ while (($n % 29) == 0) { push @factors, 29; $n = int($n / 29); }
+ }
if ($n < (31*31)) {
push @factors, $n if $n != 1;
return @factors;
@@ -1768,6 +1825,7 @@ sub prho_factor {
if ( ref($n) eq 'Math::BigInt' ) {
+ $pa = Math::BigInt->new("$pa");
$U = $n->copy->bzero->badd($U);
$V = $n->copy->bzero->badd($V);
for my $i (1 .. $rounds) {
@@ -1775,10 +1833,10 @@ sub prho_factor {
$U->bmul($U)->badd($pa)->bmod($n);
$V->bmul($V)->badd($pa);
$V->bmul($V)->badd($pa)->bmod($n);
- my $f = Math::BigInt::bgcd( ($U > $V) ? $U-$V : $V-$U, $n);
- if ($f == $n) {
+ my $f = Math::BigInt::bgcd( $U-$V, $n);
+ if ($f->bacmp($n) == 0) {
last if $inloop++; # We've been here before
- } elsif ($f != 1) {
+ } elsif (!$f->is_one) {
return _found_factor($f, $n, "prho", @factors);
}
}
diff --git a/t/23-primality-proofs.t b/t/23-primality-proofs.t
index 9eb5f49..25edfde 100644
--- a/t/23-primality-proofs.t
+++ b/t/23-primality-proofs.t
@@ -60,6 +60,8 @@ foreach my $p (@plist) {
ok( is_prime($p), "$p is prime" );
skip "These take a long time on non-64-bit. Skipping", 5
if !$use64 && !$extra && $p =~ /^(6778|9800)/;
+ skip "Skipping a certificate without GMP", 5
+ if !prime_get_config->{'gmp'} && !$extra && $p =~ /^9800/;
my($isp, $cert) = is_provable_prime_with_cert($p);
is( $isp, 2, " is_provable_prime_with_cert returns 2" );
ok( defined($cert) && $cert =~ /^Type/m,
diff --git a/t/81-bignum.t b/t/81-bignum.t
index ece2356..5e21c6c 100644
--- a/t/81-bignum.t
+++ b/t/81-bignum.t
@@ -76,7 +76,7 @@ plan tests => 0
+ 6*2*$extra # more PC tests
+ 2*scalar(keys %factors)
+ scalar(keys %allfactors)
- + 13 # moebius, euler_phi, jordan totient, divsum, znorder, etc.
+ + 13+4*$extra # moebius, euler_phi, jordan totient, divsum, etc.
+ 2 # liouville
+ 3 # gcd
+ 3 # lcm
@@ -136,12 +136,12 @@ use Math::Prime::Util qw/
# LogarithmicIntegral
# RiemannR
+my $usegmp = Math::Prime::Util::prime_get_config->{gmp};
my $bignumver = $bigint::VERSION;
my $bigintver = $Math::BigInt::VERSION;
my $bigintlib = Math::BigInt->config()->{lib};
-$bigintlib =~ s/^Math::BigInt:://;
-my $mpugmpver = Math::Prime::Util::prime_get_config->{gmp}
- ? $Math::Prime::Util::GMP::VERSION : "<none>";
+ $bigintlib =~ s/^Math::BigInt:://;
+my $mpugmpver = $usegmp ? $Math::Prime::Util::GMP::VERSION : "<none>";
diag "BigInt $bignumver/$bigintver, lib: $bigintlib. MPU::GMP $mpugmpver\n";
@@ -160,8 +160,10 @@ foreach my $n (@composites) {
foreach my $n (@proveprimes) {
ok( is_prime($n), "$n is prime" );
SKIP: {
- skip "Large proof on 32-bit machine.", 1
+ skip "Large proof on 32-bit machine without EXTENDED_TESTING.", 1
if !$use64 && !$extra && $n > 2**66;
+ skip "Large proof without GMP or EXTENDED_TESTING.", 1
+ if !$usegmp && !$extra && $n > 2**66;
skip "Skipping provable primes on broken 64-bit", 1 if $broken64;
ok( is_provable_prime($n), "$n is provably prime" );
}
@@ -223,23 +225,34 @@ SKIP: {
###############################################################################
SKIP: {
- skip "Your 64-bit Perl is broken, skipping moebius, totient, etc.", 13 if $broken64;
+ skip "Your 64-bit Perl is broken, skipping moebius, totient, etc.", 13+4*$extra if $broken64;
my $n;
$n = 618970019642690137449562110;
is( moebius($n), -1, "moebius($n)" );
is( euler_phi($n), 145857122964987051805507584, "euler_phi($n)" );
is( carmichael_lambda($n), 3271601336256, "carmichael_lambda($n)" );
- $n = 48981631802481400359696467;
- is( jordan_totient(5,$n), 281946200770875813001683560563488308767928594805846855593191749929654015729263525162226378019837608857421063724603387506651820000, "jordan_totient(5,$n)" );
- is( divisor_sum( $n, sub { my $d=shift; $d**5 * moebius($n/$d); }), 281946200770875813001683560563488308767928594805846855593191749929654015729263525162226378019837608857421063724603387506651820000, "jordan totient using divisor_sum and moebius" );
+
+ $n = 2188536338969724335807;
+ is( jordan_totient(5,$n), 50207524710890617788554288878260755791080217791665431423557510096680804997771551711694188532723268222129800, "jordan_totient(5,$n)" );
+ is( divisor_sum( $n, sub { my $d=shift; $d**5 * moebius($n/$d); }), 50207524710890617788554288878260755791080217791665431423557510096680804997771551711694188532723268222129800, "jordan totient using divisor_sum and moebius" );
+
+ if ($extra) {
+ $n = 48981631802481400359696467;
+ is( jordan_totient(5,$n), "281946200770875813001683560563488308767928594805846855593191749929654015729263525162226378019837608857421063724603387506651820000", "jordan_totient(5,$n)" );
+ is( divisor_sum( $n, sub { my $d=shift; $d**5 * moebius($n/$d); }), "281946200770875813001683560563488308767928594805846855593191749929654015729263525162226378019837608857421063724603387506651820000", "jordan totient using divisor_sum and moebius" );
+ }
+
# Done wrong, the following will have a bunch of extra zeros.
my $hundredfac = Math::BigInt->new(100)->bfac;
is( divisor_sum($hundredfac), 774026292208877355243820142464115597282472420387824628823543695735957009720184359087194959566149232506852422409529601312686157396490982598473425595924480000000, "Divisor sum of 100!" );
# These should yield bigint results.
# Quoted 0 to prevent error in perl 5.8.2 + bigint 0.23 (0 turns into NaN)
- is( divisor_sum(pn_primorial(71),"0"), 2361183241434822606848, "Divisor count(353#)" );
- is( divisor_sum(pn_primorial(71),1), 592169807666179080336898884075191344863843751107274613826065194910163387683715846870630955555390054490059876013007363004327526400000000000000000, "Divisor sum(353#)" );
- is( divisor_sum(pn_primorial(71),2), "12949784465615028275107011121945805610528825503288465119226912396970037707579655747291137846306343809131200618880146749230653882973421307691846381555612687582146340434261447200658536708625570145324567757917046739100833453606420350207262720000000000000000000000000000000000000000000000000", "sigma_2(353#)" );
+ is( divisor_sum(pn_primorial(27),"0"), 134217728, "Divisor count(103#)" );
+ is( divisor_sum(pn_primorial(27),1), "123801167235014219383860918985791897600000", "Divisor sum(103#)" );
+ is( divisor_sum(pn_primorial(27),2), "872887488619258559049272439859735080160421720974947767918289356800000000000000000", "sigma_2(103#)" );
+ if ($extra) {
+ is( divisor_sum(pn_primorial(71),"0"), 2361183241434822606848, "Divisor count(353#)" );
+ }
# Calc/FastCalc are slugs with this function, so tone things down.
#is( znorder(82734587234,927208363107752634625923555185111613055040823736157),
# 4360156780036190093445833597286118936800,
@@ -287,10 +300,10 @@ cmp_ok( $randprime, '>', 2**79, "random 80-bit prime is not too small");
cmp_ok( $randprime, '<', 2**80, "random 80-bit prime is not too big");
ok( is_prime($randprime), "random 80-bit prime is just right");
-$randprime = random_strong_prime(256);
-cmp_ok( $randprime, '>', 2**255, "random 256-bit strong prime is not too small");
-cmp_ok( $randprime, '<', 2**256, "random 256-bit strong prime is not too big");
-ok( is_prime($randprime), "random 256-bit strong prime is just right");
+$randprime = random_strong_prime(190);
+cmp_ok( $randprime, '>', 2**189, "random 190-bit strong prime is not too small");
+cmp_ok( $randprime, '<', 2**190, "random 190-bit strong prime is not too big");
+ok( is_prime($randprime), "random 190-bit strong prime is just right");
$randprime = random_maurer_prime(80);
cmp_ok( $randprime, '>', 2**79, "random 80-bit Maurer prime is not too small");
--
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