[libmath-prime-util-perl] 27/40: Clean up Lucas, update PP is_prime, remove an test func that got put in XS.xs, add documentation
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:49:05 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list