[libmath-prime-util-perl] 144/181: Performance for PP, and a few pre-5.8 64-bit workarounds
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:15 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 7406eed3312c432ef52c9d48a9c98691edb0155b
Author: Dana Jacobsen <dana at acm.org>
Date: Wed Jan 8 19:36:34 2014 -0800
Performance for PP, and a few pre-5.8 64-bit workarounds
---
Changes | 4 +
lib/Math/Prime/Util.pm | 52 ++++++----
lib/Math/Prime/Util/PP.pm | 242 ++++++++++++++++++++++++++++++----------------
3 files changed, 200 insertions(+), 98 deletions(-)
diff --git a/Changes b/Changes
index c8e0d17..d1b324f 100644
--- a/Changes
+++ b/Changes
@@ -43,6 +43,10 @@ Revision history for Perl module Math::Prime::Util
turned off but the GMP module enabled, things will run slower since
they no longer go to GMP.
+ - Test suite should run faster. Combination of small speedups to hot
+ spots as well as pushing a few slow tasks to EXTENDED_TESTING (these
+ are generally things never used, like pure Perl AKS).
+
- Some 5.6.2-is-broken workarounds.
- Some LMO edge cases: small numbers that only show up if a #define is
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index a35e36a..b263fbb 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -88,7 +88,8 @@ BEGIN {
$_Config{'xs'} = 1;
1;
} or do {
- #carp "Using Pure Perl implementation: $@";
+ carp "Using Pure Perl implementation: $@"
+ unless defined $ENV{MPU_NO_XS} && $ENV{MPU_NO_XS} == 1;
$_Config{'xs'} = 0;
$_Config{'maxbits'} = MPU_MAXBITS;
@@ -943,8 +944,9 @@ sub primes {
}
# We know we don't have GMP and are > 2^64, so skip all the middle.
#next unless is_prob_prime($p);
- next unless Math::Prime::Util::PP::is_strong_pseudoprime($p, 2);
- next unless Math::Prime::Util::PP::is_extra_strong_lucas_pseudoprime($p);
+ #next unless Math::Prime::Util::PP::is_strong_pseudoprime($p, 2);
+ #next unless Math::Prime::Util::PP::is_extra_strong_lucas_pseudoprime($p);
+ next unless Math::Prime::Util::PP::is_bpsw_prime($p);
}
return $p;
}
@@ -1050,7 +1052,7 @@ sub primes {
}
# I've seen +0, +1, and +2 here. Maurer uses +0. Menezes uses +1.
- my ($q, $qcert) = random_maurer_prime_with_cert( ($r * $k)->bfloor + 1 );
+ my ($q, $qcert) = random_maurer_prime_with_cert( ($r * $k)->bfloor->binc );
$q = Math::BigInt->new("$q") unless ref($q) eq 'Math::BigInt';
my $I = Math::BigInt->new(2)->bpow($k-2)->bdiv($q)->bfloor->as_int();
print "r = $r k = $k q = $q I = $I\n" if $verbose && $verbose != 3;
@@ -1059,13 +1061,17 @@ sub primes {
# Big GCD's are hugely fast with GMP or Pari, but super slow with Calc.
_make_big_gcds() if $_big_gcd_use < 0;
+ my $ONE = Math::BigInt->bone;
+ my $TWO = $ONE->copy->binc;
my $loop_limit = 1_000_000 + $k * 1_000;
while ($loop_limit-- > 0) {
# R is a random number between $I+1 and 2*$I
- my $R = $I + 1 + $_RANDF->( $I - 1 );
+ #my $R = $I + 1 + $_RANDF->( $I - 1 );
+ my $R = $I->copy->binc->badd( $_RANDF->($I->copy->bdec) );
#my $n = 2 * $R * $q + 1;
- my $n = Math::BigInt->new(2)->bmul($R)->bmul($q)->binc();
+ my $nm1 = $TWO->copy->bmul($R)->bmul($q);
+ my $n = $nm1->copy->binc;
# We constructed a promising looking $n. Now test it.
print "." if $verbose > 2;
if ($_HAVE_GMP) {
@@ -1073,12 +1079,12 @@ sub primes {
next unless Math::Prime::Util::GMP::is_prob_prime($n);
} else {
# No GMP, so first do trial divisions, then a SPSP test.
- next unless Math::BigInt::bgcd($n, 111546435) == 1;
+ next unless Math::BigInt::bgcd($n, 111546435)->is_one;
if ($_big_gcd_use && $n > $_big_gcd_top) {
- next unless Math::BigInt::bgcd($n, $_big_gcd[0]) == 1;
- next unless Math::BigInt::bgcd($n, $_big_gcd[1]) == 1;
- next unless Math::BigInt::bgcd($n, $_big_gcd[2]) == 1;
- next unless Math::BigInt::bgcd($n, $_big_gcd[3]) == 1;
+ next unless Math::BigInt::bgcd($n, $_big_gcd[0])->is_one;
+ next unless Math::BigInt::bgcd($n, $_big_gcd[1])->is_one;
+ next unless Math::BigInt::bgcd($n, $_big_gcd[2])->is_one;
+ next unless Math::BigInt::bgcd($n, $_big_gcd[3])->is_one;
}
print "+" if $verbose > 2;
next unless is_strong_pseudoprime($n, 3);
@@ -1093,9 +1099,9 @@ sub primes {
# BLS75 theorem 4 (Pocklington) used by Maurer's paper.
# Check conditions -- these should be redundant.
- my $m = 2 * $R;
+ my $m = $TWO * $R;
if (! ($q->is_odd && $q > 2 && $m > 0 &&
- $m * $q + 1 == $n && 2*$q+1 > $n->copy->bsqrt()) ) {
+ $m * $q + $ONE == $n && $TWO*$q+$ONE > $n->copy->bsqrt()) ) {
carp "Maurer prime failed BLS75 theorem 3 conditions. Retry.";
next;
}
@@ -1103,8 +1109,8 @@ sub primes {
foreach my $trya (2, 3, 5, 7, 11, 13) {
my $a = Math::BigInt->new($trya);
# m/2 = R (n-1)/2 = (2*R*q)/2 = R*q
- next unless $a->copy->bmodpow($R, $n) != $n-1;
- next unless $a->copy->bmodpow($R*$q, $n) == $n-1;
+ next unless $a->copy->bmodpow($R, $n) != $nm1;
+ next unless $a->copy->bmodpow($R*$q, $n) == $nm1;
print "($k)" if $verbose > 2;
croak "Maurer prime $n=2*$R*$q+1 failed BPSW" unless is_prob_prime($n);
my $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\n" .
@@ -1187,12 +1193,24 @@ sub primorial {
return Math::BigInt->new(''.Math::Prime::Util::GMP::primorial($n))
if $_HAVE_GMP && defined &Math::Prime::Util::GMP::primorial;
- my $max = (MPU_32BIT) ? 29 : 53;
+ my $max = (MPU_32BIT) ? 29 : (OLD_PERL_VERSION) ? 43 : 53;
my $pn = (ref($_[0]) eq 'Math::BigInt') ? $_[0]->copy->bone()
: ($n >= $max) ? Math::BigInt->bone()
: 1;
if (ref($pn) eq 'Math::BigInt') {
- $pn->bmul($_) for map { Math::BigInt->new($_) } @{primes(2,$n)};
+ my $start = 2;
+ if ($n >= 97) {
+ $start = 101;
+ $pn->bdec->badd(Math::BigInt->new("2305567963945518424753102147331756070"));
+ }
+ my @plist = @{primes($start,$n)};
+ while (@plist > 2 && $plist[2] < 1625) {
+ $pn->bmul( Math::BigInt->new(shift(@plist)*shift(@plist)*shift(@plist)) );
+ }
+ while (@plist > 1 && $plist[1] < 65536) {
+ $pn->bmul( Math::BigInt->new(shift(@plist)*shift(@plist)) );
+ }
+ $pn->bmul($_) for @plist;
} else {
forprimes { $pn *= $_ } $n;
}
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 9068e37..34ab778 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -24,13 +24,18 @@ BEGIN {
use constant MPU_32BIT => MPU_MAXBITS == 32;
#use constant MPU_MAXPARAM => MPU_32BIT ? 4294967295 : 18446744073709551615;
#use constant MPU_MAXDIGITS => MPU_32BIT ? 10 : 20;
- #use constant MPU_MAXPRIME => MPU_32BIT ? 4294967291 : 18446744073709551557;
+ use constant MPU_MAXPRIME => MPU_32BIT ? 4294967291 : 18446744073709551557;
use constant MPU_MAXPRIMEIDX => MPU_32BIT ? 203280221 : 425656284035217743;
use constant MPU_HALFWORD => MPU_32BIT ? 65536 : OLD_PERL_VERSION ? 33554432 : 4294967296;
use constant UVPACKLET => MPU_32BIT ? 'L' : 'Q';
use constant MPU_INFINITY => (65535 > 0+'inf') ? 20**20**20 : 0+'inf';
use constant CONST_EULER => '0.577215664901532860606512090082402431042159335939923598805767';
use constant CONST_LI2 => '1.04516378011749278484458888919461313652261557815120157583290914407501320521';
+ use constant BZERO => Math::BigInt->bzero;
+ use constant BONE => Math::BigInt->bone;
+ use constant BTWO => Math::BigInt->new(2);
+ use constant B_PRIM759 => Math::BigInt->new("64092011671807087969");
+ use constant B_PRIM235 => Math::BigInt->new("30");
}
{
@@ -133,6 +138,13 @@ my @_prevwheel30 = (29,29,1,1,1,1,1,1,7,7,7,7,11,11,13,13,13,13,17,17,19,19,19,1
sub _is_prime7 { # n must not be divisible by 2, 3, or 5
my($n) = @_;
+ if (ref($n) eq 'Math::BigInt') {
+ return 0 unless Math::BigInt::bgcd($n, B_PRIM759)->is_one;
+ return 0 unless _miller_rabin_2($n);
+ my $is_esl_prime = is_extra_strong_lucas_pseudoprime($n);
+ return ($is_esl_prime) ? (($n <= "18446744073709551615") ? 2 : 1) : 0;
+ }
+
if ($n < 61*61) {
foreach my $i (qw/7 11 13 17 19 23 29 31 37 41 43 47 53 59/) {
return 2 if $i*$i > $n;
@@ -141,14 +153,10 @@ sub _is_prime7 { # n must not be divisible by 2, 3, or 5
return 2;
}
- if (ref($n) eq 'Math::BigInt') {
- return 0 unless Math::BigInt::bgcd($n, "64092011671807087969")->is_one;
- } else {
- return 0 if !($n % 7) || !($n % 11) || !($n % 13) || !($n % 17) ||
- !($n % 19) || !($n % 23) || !($n % 29) || !($n % 31) ||
- !($n % 37) || !($n % 41) || !($n % 43) || !($n % 47) ||
- !($n % 53) || !($n % 59);
- }
+ return 0 if !($n % 7) || !($n % 11) || !($n % 13) || !($n % 17) ||
+ !($n % 19) || !($n % 23) || !($n % 29) || !($n % 31) ||
+ !($n % 37) || !($n % 41) || !($n % 43) || !($n % 47) ||
+ !($n % 53) || !($n % 59);
if ($n <= 1_000_000) {
$n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt';
@@ -194,11 +202,8 @@ sub _is_prime7 { # n must not be divisible by 2, 3, or 5
}
# Inlined BPSW
- return 0 unless is_strong_pseudoprime($n, 2);
- if ($n <= 18446744073709551615) {
- return is_almost_extra_strong_lucas_pseudoprime($n) ? 2 : 0;
- }
- return is_extra_strong_lucas_pseudoprime($n) ? 1 : 0;
+ return 0 unless _miller_rabin_2($n);
+ return is_almost_extra_strong_lucas_pseudoprime($n) ? 2 : 0;
}
sub is_prime {
@@ -206,8 +211,12 @@ sub is_prime {
return 0 if defined $n && int($n) < 0;
_validate_positive_integer($n);
- if ($n < 7) { return ($n == 2) || ($n == 3) || ($n == 5) ? 2 : 0; }
- return 0 if !($n % 2) || !($n % 3) || !($n % 5);
+ if (ref($n) eq 'Math::BigInt') {
+ return 0 unless Math::BigInt::bgcd($n, B_PRIM235)->is_one;
+ } else {
+ if ($n < 7) { return ($n == 2) || ($n == 3) || ($n == 5) ? 2 : 0; }
+ return 0 if !($n % 2) || !($n % 3) || !($n % 5);
+ }
return _is_prime7($n);
}
@@ -221,7 +230,7 @@ sub is_prime {
sub is_bpsw_prime {
my($n) = @_;
_validate_positive_integer($n);
- return 0 unless is_strong_pseudoprime($n, 2);
+ return 0 unless _miller_rabin_2($n);
if ($n <= 18446744073709551615) {
return is_almost_extra_strong_lucas_pseudoprime($n) ? 2 : 0;
}
@@ -411,10 +420,8 @@ sub next_prime {
return $_prime_next_small[$n] if $n <= $#_prime_next_small;
- if (ref($n) ne 'Math::BigInt' && $n >= 4294967291) {
- $n = Math::BigInt->new(''.$_[0])
- if MPU_32BIT || $n >= 18446744073709551557;
- }
+ $n = Math::BigInt->new(''.$_[0])
+ if ref($n) ne 'Math::BigInt' && $n >= MPU_MAXPRIME;
#$n = ($n+1) | 1;
#while ( !($n%3) || !($n%5) || !($n%7) || !($n%11) || !($n%13)
@@ -449,16 +456,16 @@ sub prev_prime {
$n -= ($n & 1) ? 2 : 1;
my $nmod6 = $n % 6;
if ($nmod6 == 5) {
- return $n if $n % 5 && $n % 7 && _is_prime7($n);
+ return $n if ($n % 5) != 0 && ($n % 7) != 0 && _is_prime7($n);
$n -= 4;
} elsif ($nmod6 == 3) {
$n -= 2;
}
while (1) {
- return $n if $n % 5 && $n % 7 && _is_prime7($n);
+ return $n if ($n % 5) != 0 && ($n % 7) != 0 && _is_prime7($n);
$n -= 2;
- return $n if $n % 5 && $n % 7 && _is_prime7($n);
+ return $n if ($n % 5) != 0 && ($n % 7) != 0 && _is_prime7($n);
$n -= 4;
}
return $n;
@@ -536,7 +543,7 @@ sub moebius {
my($n) = @_;
return ($n == 1) ? 1 : 0 if $n <= 1;
return 0 if ($n >= 49) && (!($n % 4) || !($n % 9) || !($n % 25) || !($n%49) );
- my @factors = ($n < 1_000_000) ? trial_factor($n) : factor($n);
+ my @factors = Math::Prime::Util::factor($n);
foreach my $i (1 .. $#factors) {
return 0 if $factors[$i] == $factors[$i-1];
}
@@ -879,6 +886,7 @@ sub nth_prime {
sub _mulmod {
my($x, $y, $n) = @_;
return (($x * $y) % $n) if ($x|$y) < MPU_HALFWORD;
+ #return (($x * $y) % $n) if ($x|$y) < MPU_HALFWORD || $y == 0 || $x < int(~0/$y);
my $r = 0;
$x %= $n if $x >= $n;
$y %= $n if $y >= $n;
@@ -955,10 +963,11 @@ sub lcm {
return $lcm;
}
-# unsigned, no validation
+# no validation, x is allowed to be negative, y must be >= 0
sub _gcd_ui {
my($x, $y) = @_;
if ($y < $x) { ($x, $y) = ($y, $x); }
+ elsif ($x < 0) { $x = -$x; }
while ($y > 0) {
# y1 <- x0 % y0 ; x1 <- y0
my $t = $y;
@@ -1001,6 +1010,61 @@ sub is_pseudoprime {
return ($x == 1) ? 1 : 0;
}
+sub _miller_rabin_2 {
+ my($n, $nm1, $s, $d) = @_;
+
+ if ( ref($n) eq 'Math::BigInt' ) {
+
+ if (!defined $nm1) {
+ $nm1 = $n->copy->bdec();
+ $s = 0;
+ $d = $nm1->copy;
+ do {
+ $s++;
+ $d->brsft(BONE);
+ } while $d->is_even;
+ }
+ my $x = BTWO->copy->bmodpow($d,$n);
+ return 1 if $x->is_one || $x->bcmp($nm1) == 0;
+ foreach my $r (1 .. $s-1) {
+ $x->bmul($x)->bmod($n);
+ last if $x->is_one;
+ return 1 if $x->bcmp($nm1) == 0;
+ }
+
+ } else {
+
+ if (!defined $nm1) {
+ $nm1 = $n-1;
+ $s = 0;
+ $d = $nm1;
+ while ( ($d & 1) == 0 ) {
+ $s++;
+ $d >>= 1;
+ }
+ }
+
+ if ($n < MPU_HALFWORD) {
+ my $x = _native_powmod(2, $d, $n);
+ return 1 if $x == 1 || $x == $nm1;
+ foreach my $r (1 .. $s-1) {
+ $x = ($x*$x) % $n;
+ last if $x == 1;
+ return 1 if $x == $n-1;
+ }
+ } else {
+ my $x = _powmod(2, $d, $n);
+ return 1 if $x == 1 || $x == $nm1;
+ foreach my $r (1 .. $s-1) {
+ $x = ($x < MPU_HALFWORD) ? ($x*$x) % $n : _mulmod($x, $x, $n);
+ last if $x == 1;
+ return 1 if $x == $n-1;
+ }
+ }
+ }
+ 0;
+}
+
sub is_strong_pseudoprime {
my($n, @bases) = @_;
return 0 if defined $n && int($n) < 0;
@@ -1010,6 +1074,12 @@ sub is_strong_pseudoprime {
return 0+($n >= 2) if $n < 4;
return 0 if ($n % 2) == 0;
+ if ($bases[0] == 2) {
+ return 0 unless _miller_rabin_2($n);
+ shift @bases;
+ return 1 unless @bases;
+ }
+
# Die on invalid bases
foreach my $base (@bases) { croak "Base $base is invalid" if $base < 2 }
# Make sure we handle big bases ok.
@@ -1022,7 +1092,7 @@ sub is_strong_pseudoprime {
my $d = $nminus1->copy;
do { # n is > 3 and odd, so n-1 must be even
$s++;
- $d->brsft(1);
+ $d->brsft(BONE);
} while $d->is_even;
# Different way of doing the above. Fewer function calls, slower on ave.
#my $dbin = $nminus1->as_bin;
@@ -1078,6 +1148,8 @@ sub is_strong_pseudoprime {
}
1;
}
+
+
# Calculate Kronecker symbol (a|b). Cohen Algorithm 1.4.10.
# Extension of the Jacobi symbol, itself an extension of the Legendre symbol.
sub kronecker {
@@ -1226,45 +1298,47 @@ sub lucas_sequence {
$n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
my $ZERO = $n->copy->bzero;
- 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 - $TWO*$TWO*$Q;
+ my $D = $P*$P - BTWO*BTWO*$Q;
croak "lucas_sequence: D is zero" if $D->is_zero;
- my $U = $ONE->copy;
+ my $U = BONE->copy;
my $V = $P->copy;
my $Qk = $Q->copy;
- return ($ZERO, $TWO) if $k == 0;
+ return (BZERO->copy, BTWO->copy) 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->is_one) {
my $Dinverse = $D->copy->bmodinv($n);
- if ($P > $TWO && !$Dinverse->is_nan) {
+ if ($P > BTWO && !$Dinverse->is_nan) {
# Calculate V_k with U=V_{k+1}
- $U = $P->copy->bmul($P)->bsub($TWO)->bmod($n);
+ $U = $P->copy->bmul($P)->bsub(BTWO)->bmod($n);
while (++$bpos < length($kstr)) {
if (substr($kstr,$bpos,1)) {
$V->bmul($U)->bsub($P )->bmod($n);
- $U->bmul($U)->bsub($TWO)->bmod($n);
+ $U->bmul($U)->bsub(BTWO)->bmod($n);
} else {
$U->bmul($V)->bsub($P )->bmod($n);
- $V->bmul($V)->bsub($TWO)->bmod($n);
+ $V->bmul($V)->bsub(BTWO)->bmod($n);
}
}
# Crandall and Pomerance eq 3.13: U_n = D^-1 (2V_{n+1} - PV_n)
- $U = $Dinverse * ($TWO*$U - $P*$V);
+ $U = $Dinverse * (BTWO*$U - $P*$V);
} else {
while (++$bpos < length($kstr)) {
$U->bmul($V)->bmod($n);
- $V->bmul($V)->bsub($TWO)->bmod($n);
+ $V->bmul($V)->bsub(BTWO)->bmod($n);
if (substr($kstr,$bpos,1)) {
my $T1 = $U->copy->bmul($D);
- $U->bmul($P)->badd( $V)->badd( $U->is_odd ? $n : $ZERO )->brsft(1);
- $V->bmul($P)->badd($T1)->badd( $V->is_odd ? $n : $ZERO )->brsft(1);
+ $U->bmul($P)->badd( $V);
+ $U->badd($n) if $U->is_odd;
+ $U->brsft(BONE);
+ $V->bmul($P)->badd($T1);
+ $V->badd($n) if $V->is_odd;
+ $V->brsft(BONE);
}
}
}
@@ -1272,20 +1346,18 @@ sub lucas_sequence {
my $qsign = ($Q == -1) ? -1 : 0;
while (++$bpos < length($kstr)) {
$U->bmul($V)->bmod($n);
- if ($qsign == 1) { $V->bmul($V)->bsub($TWO)->bmod($n); }
- elsif ($qsign == -1) { $V->bmul($V)->badd($TWO)->bmod($n); }
- else { $V->bmul($V)->bsub($Qk->copy->blsft($ONE))->bmod($n); }
+ if ($qsign == 1) { $V->bmul($V)->bsub(BTWO)->bmod($n); }
+ elsif ($qsign == -1) { $V->bmul($V)->badd(BTWO)->bmod($n); }
+ else { $V->bmul($V)->bsub($Qk->copy->blsft(BONE))->bmod($n); }
if (substr($kstr,$bpos,1)) {
my $T1 = $U->copy->bmul($D);
- $U->bmul($P);
- $U->badd($V);
+ $U->bmul($P)->badd( $V);
$U->badd($n) if $U->is_odd;
- $U->brsft($ONE);
+ $U->brsft(BONE);
- $V->bmul($P);
- $V->badd($T1);
+ $V->bmul($P)->badd($T1);
$V->badd($n) if $V->is_odd;
- $V->brsft($ONE);
+ $V->brsft(BONE);
if ($qsign != 0) { $qsign = -1; }
else { $Qk->bmul($Qk)->bmul($Q)->bmod($n); }
@@ -1342,7 +1414,7 @@ sub is_strong_lucas_pseudoprime {
foreach my $r (0 .. $s-1) {
return 1 if $V->is_zero;
if ($r < ($s-1)) {
- $V->bmul($V)->bsub(2*$Qk)->bmod($n);
+ $V->bmul($V)->bsub(BTWO*$Qk)->bmod($n);
$Qk->bmul($Qk)->bmod($n);
}
}
@@ -1368,15 +1440,15 @@ sub is_extra_strong_lucas_pseudoprime {
my($s, $k) = (0, $n->copy->binc);
while ($k->is_even && !$k->is_zero) {
$s++;
- $k->brsft(1);
+ $k->brsft(BONE);
}
my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k);
- return 1 if $U->is_zero && ($V == 2 || $V == ($n-2));
+ return 1 if $U->is_zero && ($V == BTWO || $V == ($n - BTWO));
foreach my $r (0 .. $s-2) {
return 1 if $V->is_zero;
- $V->bmul($V)->bsub(2)->bmod($n);
+ $V->bmul($V)->bsub(BTWO)->bmod($n);
}
return 0;
}
@@ -1456,7 +1528,8 @@ sub is_frobenius_underwood_pseudoprime {
$fb->bsub($fa)->bmul($t)->bmod($n);
$fa->bmul($na)->bmod($n);
- if ( ($np1 >> $bit) & $ONE ) {
+ #if ( ($np1 >> $bit) & 1 ) {
+ if ( $np1->copy->brsft($bit)->is_odd ) {
$t = $fb->copy;
$fb->badd($fb)->bsub($fa);
$fa->bmul($multiplier)->badd($t);
@@ -1614,8 +1687,8 @@ sub _basic_factor {
while ( !($_[0] % 3) ) { push @factors, 3; $_[0] = int($_[0] / 3); }
while ( !($_[0] % 5) ) { push @factors, 5; $_[0] = int($_[0] / 5); }
} else {
- if (Math::BigInt::bgcd($_[0], 2*3*5) != 1) {
- while ( $_[0]->is_even) { push @factors, 2; $_[0]->brsft(1); }
+ if (!Math::BigInt::bgcd($_[0], B_PRIM235)->is_one) {
+ while ( $_[0]->is_even) { push @factors, 2; $_[0]->brsft(BONE); }
foreach my $div (3, 5) {
my ($q, $r) = $_[0]->copy->bdiv($div);
while ($r->is_zero) {
@@ -1726,7 +1799,7 @@ sub factor {
my @factors;
# Use 'n=int($n/7)' instead of 'n/=7' to not "upgrade" n to a Math::BigFloat.
if (ref($n) eq 'Math::BigInt') {
- while ($n->is_even) { push @factors, 2; $n->brsft(1); }
+ while ($n->is_even) { push @factors, 2; $n->brsft(BONE); }
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);
@@ -1770,6 +1843,7 @@ sub factor {
if (scalar @ftry > 1) {
#print " split into ", join(",", at ftry), "\n";
$n = shift @ftry;
+ $n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt' && $n <= ''.~0;
push @nstack, @ftry;
} else {
#warn "trial factor $n\n";
@@ -1828,7 +1902,7 @@ 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, $n);
+ my $f = Math::BigInt::bgcd($U-$V, $n);
if ($f->bacmp($n) == 0) {
last if $inloop++; # We've been here before
} elsif (!$f->is_one) {
@@ -1842,7 +1916,7 @@ sub prho_factor {
$U = ($U * $U + $pa) % $n;
$V = ($V * $V + $pa) % $n;
$V = ($V * $V + $pa) % $n;
- my $f = _gcd_ui( ($U > $V) ? $U-$V : $V-$U, $n );
+ my $f = _gcd_ui( $U-$V, $n );
if ($f == $n) {
last if $inloop++; # We've been here before
} elsif ($f != 1) {
@@ -1853,10 +1927,16 @@ sub prho_factor {
} else {
for my $i (1 .. $rounds) {
- $U = _mulmod($U, $U, $n); $U = (($n-$U) > $pa) ? $U+$pa : $U-$n+$pa;
- $V = _mulmod($V, $V, $n); $V = (($n-$V) > $pa) ? $V+$pa : $V-$n+$pa;
- $V = _mulmod($V, $V, $n); $V = (($n-$V) > $pa) ? $V+$pa : $V-$n+$pa;
- my $f = _gcd_ui( ($U > $V) ? $U-$V : $V-$U, $n );
+ if ($n <= (~0 >> 1)) {
+ $U = _mulmod($U, $U, $n); $U += $pa; $U -= $n if $U >= $n;
+ $V = _mulmod($V, $V, $n); $V += $pa; # Let the mulmod handle it
+ $V = _mulmod($V, $V, $n); $V += $pa; $V -= $n if $V >= $n;
+ } else {
+ $U = _mulmod($U, $U, $n); $U=$n-$U; $U = ($pa>=$U) ? $pa-$U : $n-$U+$pa;
+ $V = _mulmod($V, $V, $n); $V=$n-$V; $V = ($pa>=$V) ? $pa-$V : $n-$V+$pa;
+ $V = _mulmod($V, $V, $n); $V=$n-$V; $V = ($pa>=$V) ? $pa-$V : $n-$V+$pa;
+ }
+ my $f = _gcd_ui( $U-$V, $n );
if ($f == $n) {
last if $inloop++; # We've been here before
} elsif ($f != 1) {
@@ -1918,7 +1998,7 @@ sub pbrent_factor {
$Xi = $saveXi->copy;
do {
$Xi->bmul($Xi)->badd($pa)->bmod($n);
- $f = Math::BigInt::bgcd( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n);
+ $f = Math::BigInt::bgcd($Xm-$Xi, $n);
} while ($f != 1 && $r-- != 0);
last if $f == 1 || $f == $n;
}
@@ -1929,7 +2009,7 @@ sub pbrent_factor {
for my $i (1 .. $rounds) {
$Xi = ($Xi * $Xi + $pa) % $n;
- my $f = _gcd_ui( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n );
+ my $f = _gcd_ui($Xm-$Xi, $n);
return _found_factor($f, $n, "pbrent", @factors) if $f != 1 && $f != $n;
$Xm = $Xi if ($i & ($i-1)) == 0; # i is a power of 2
}
@@ -1940,7 +2020,7 @@ sub pbrent_factor {
# Xi^2+a % n
$Xi = _mulmod($Xi, $Xi, $n);
$Xi = (($n-$Xi) > $pa) ? $Xi+$pa : $Xi+$pa-$n;
- my $f = _gcd_ui( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n );
+ my $f = _gcd_ui($Xm-$Xi, $n);
return _found_factor($f, $n, "pbrent", @factors) if $f != 1 && $f != $n;
$Xm = $Xi if ($i & ($i-1)) == 0; # i is a power of 2
}
@@ -2012,7 +2092,7 @@ sub pminus1_factor {
my $t = $one->copy;
my $pa = $one->copy->binc();
my $savea = $pa->copy;
- my $f = 1;
+ my $f = $one->copy;
my($pc_beg, $pc_end, @bprimes);
$pc_beg = 2;
@@ -2031,13 +2111,13 @@ sub pminus1_factor {
if ($pa == 0) { push @factors, $n; return @factors; }
$f = Math::BigInt::bgcd( $pa-1, $n );
last if $f == $n;
- return _found_factor($f, $n, "pminus1", @factors) if $f != 1;
+ return _found_factor($f, $n, "pminus1", @factors) unless $f->is_one;
$saveq = $q;
$savea = $pa->copy;
}
}
$q = $bprimes[-1];
- last if $f != 1 || $pc_end >= $B1;
+ last if !$f->is_one || $pc_end >= $B1;
$pc_beg = $pc_end+1;
$pc_end += 500_000;
}
@@ -2054,12 +2134,12 @@ sub pminus1_factor {
$pa->bmodpow($k, $n);
my $f = Math::BigInt::bgcd( $pa-1, $n );
if ($f == $n) { push @factors, $n; return @factors; }
- last if $f != 1;
+ last if !$f->is_one;
$q = next_prime($q);
}
}
# STAGE 2
- if ($f == 1 && $B2 > $B1) {
+ if ($f->is_one && $B2 > $B1) {
my $bm = $pa->copy;
my $b = $one->copy;
my @precomp_bm;
@@ -2087,10 +2167,10 @@ sub pminus1_factor {
if (($j++ % 128) == 0) {
$b->bmod($n);
$f = Math::BigInt::bgcd( $b, $n );
- last if $f != 1;
+ last if !$f->is_one;
}
}
- last if $f != 1 || $pc_end >= $B2;
+ last if !$f->is_one || $pc_end >= $B2;
$pc_beg = $pc_end+1;
$pc_end += 500_000;
}
@@ -2138,7 +2218,7 @@ sub holf_factor {
next unless $mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25;
my $f = int(sqrt($m));
next unless $f*$f == $m;
- $f = _gcd_ui( ($s > $f) ? $s - $f : $f - $s, $n);
+ $f = _gcd_ui($s - $f, $n);
return _found_factor($f, $n, "HOLF ($i rounds)", @factors);
}
}
@@ -2491,25 +2571,25 @@ sub LogarithmicIntegral {
# Remember MPFR eint doesn't handle negative inputs
if ($x >= 1 && _MPFR_available()) {
my $wantbf = 0;
- my $xdigits = 17;
+ my $xdigits = 18;
if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
- do { require Math::BigFloat; Math::BigFloat->import(); }
- if !defined $Math::BigFloat::VERSION;
- $x = Math::BigFloat->new("$x") if ref($x) ne 'Math::BigFloat';
$wantbf = 1;
- $xdigits = $x->accuracy || Math::BigFloat->accuracy() || Math::BigFloat->div_scale();
+ $xdigits = $x->accuracy || Math::BigInt->accuracy() || Math::BigInt->div_scale();
}
- $x = log($x); # Use BigFloat to do the log to simplify precision tracking.
+ $xdigits += length(int(log(0.0+"$x"))) + 1;
my $rnd = 0; # MPFR_RNDN;
my $bit_precision = int($xdigits * 3.322) + 4;
my $rx = Math::MPFR->new();
Math::MPFR::Rmpfr_set_prec($rx, $bit_precision);
Math::MPFR::Rmpfr_set_str($rx, "$x", 10, $rnd);
+ Math::MPFR::Rmpfr_log($rx, $rx, $rnd);
my $lix = Math::MPFR->new();
Math::MPFR::Rmpfr_set_prec($lix, $bit_precision);
Math::MPFR::Rmpfr_eint($lix, $rx, $rnd);
my $strval = Math::MPFR::Rmpfr_get_str($lix, 10, 0, $rnd);
- return ($wantbf) ? Math::BigFloat->new($strval) : 0.0 + $strval;
+ return Math::BigFloat->new($strval, ($x->accuracy || Math::BigInt->accuracy() || Math::BigInt->div_scale()))
+ if $wantbf;
+ return 0.0 + $strval;
}
if ($x == 2) {
--
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