[libmath-prime-util-perl] 08/40: Add lucas_sequence
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:49:02 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 4ebadd564814c0670b568e51509210fdebd9fc4d
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Jun 24 14:04:19 2013 -0700
Add lucas_sequence
---
Changes | 1 +
TODO | 5 +
XS.xs | 10 ++
factor.c | 2 +-
factor.h | 1 +
lib/Math/Prime/Util.pm | 36 +++++++
lib/Math/Prime/Util/PP.pm | 233 +++++++++++++++++++++++++++-------------------
t/17-pseudoprime.t | 23 ++++-
8 files changed, 212 insertions(+), 99 deletions(-)
diff --git a/Changes b/Changes
index 200f32f..d1aa936 100644
--- a/Changes
+++ b/Changes
@@ -11,6 +11,7 @@ Revision history for Perl extension Math::Prime::Util.
- Added:
is_frobenius_underwood_pseudoprime
+ lucas_sequence
0.29 30 May 2013
diff --git a/TODO b/TODO
index 33f6bc5..f9c19c0 100644
--- a/TODO
+++ b/TODO
@@ -48,3 +48,8 @@
- Write a standalone function that demonstrates the memory leak with MULTICALL,
so we can use MULTICALL.
+
+- Tests for:
+ - F-U pseudoprime
+
+- PP F-U pseudoprime
diff --git a/XS.xs b/XS.xs
index e2f5679..502a95d 100644
--- a/XS.xs
+++ b/XS.xs
@@ -450,6 +450,16 @@ _XS_miller_rabin(IN UV n, ...)
OUTPUT:
RETVAL
+void
+_XS_lucas_sequence(IN UV n, IN IV P, IN IV Q, IN UV k)
+ PREINIT:
+ UV U, V, Qk;
+ PPCODE:
+ lucas_seq(&U, &V, &Qk, n, P, Q, k);
+ XPUSHs(sv_2mortal(newSVuv( U )));
+ XPUSHs(sv_2mortal(newSVuv( V )));
+ XPUSHs(sv_2mortal(newSVuv( Qk )));
+
int
_XS_is_lucas_pseudoprime(IN UV n, int strength)
diff --git a/factor.c b/factor.c
index 0e0f61e..61e2025 100644
--- a/factor.c
+++ b/factor.c
@@ -486,7 +486,7 @@ int _XS_is_prob_prime(UV n)
}
/* Generic Lucas sequence for any appropriate P and Q */
-static void lucas_seq(UV* Uret, UV* Vret, UV* Qkret, UV n, IV P, IV Q, UV k)
+void lucas_seq(UV* Uret, UV* Vret, UV* Qkret, UV n, IV P, IV Q, UV k)
{
UV U, V, b, Dmod, Qmod, Pmod, Qk;
IV D;
diff --git a/factor.h b/factor.h
index bce382c..b913345 100644
--- a/factor.h
+++ b/factor.h
@@ -22,6 +22,7 @@ extern int _XS_is_pseudoprime(UV n, UV a);
extern int _XS_miller_rabin(UV n, const UV *bases, int nbases);
extern int _SPRP2(UV n);
extern int _XS_is_prob_prime(UV n);
+extern void lucas_seq(UV* U, UV* V, UV* Qk, UV n, IV P, IV Q, UV k);
extern int _XS_is_lucas_pseudoprime(UV n, int strength);
extern int _XS_is_frobenius_underwood_pseudoprime(UV n);
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index c78f1fa..4278d5a 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -22,6 +22,7 @@ our @EXPORT_OK =
is_frobenius_underwood_pseudoprime
is_aks_prime
miller_rabin
+ lucas_sequence
primes
forprimes prime_iterator
next_prime prev_prime
@@ -1644,6 +1645,24 @@ sub miller_rabin {
return is_strong_pseudoprime(@_);
}
+sub lucas_sequence {
+ my($n, $P, $Q, $k) = @_;
+ _validate_num($n) || _validate_positive_integer($n);
+ _validate_num($k) || _validate_positive_integer($k);
+ { my $testP = (!defined $P || $P >= 0) ? $P : -$P;
+ _validate_num($testP) || _validate_positive_integer($testP); }
+ { my $testQ = (!defined $Q || $Q >= 0) ? $Q : -$Q;
+ _validate_num($testQ) || _validate_positive_integer($testQ); }
+ return _XS_lucas_sequence($n, $P, $Q, $k)
+ if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL
+ && ref($k) ne 'Math::BigInt' && $k <= $_XS_MAXVAL;
+ return Math::Prime::Util::GMP::lucas_sequence($n, $P, $Q, $k)
+ if $_HAVE_GMP
+ && defined &Math::Prime::Util::GMP::lucas_sequence;
+ return Math::Prime::Util::PP::lucas_sequence($n, $P, $Q, $k);
+}
+
+
#############################################################################
# Oct 2012 note: these numbers are old.
@@ -3188,6 +3207,23 @@ However AKS in general is far too slow to be of practical use. R.P. Brent,
2010: "AKS is not a practical algorithm. ECPP is much faster."
+=head2 lucas_sequence
+
+ my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k)
+
+Computes C<U_k>, C<V_k>, and C<Q_k> for the Lucas seqence defined by
+C<P>,C<Q>, modulo C<n>. The modular Lucas sequence is used in a
+number of primality tests and proofs.
+
+The following conditions must hold:
+ - C<< D = P*P - 4*Q != 0 >>
+ - C<< P > 0 >>
+ - C<< P < n >>
+ - C<< Q < n >>
+ - C<< k >= 0 >>
+ - C<< n >= 2 >>
+
+
=head2 moebius
say "$n is square free" if moebius($n) != 0;
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 07ef07e..435e5df 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -903,7 +903,7 @@ sub _jacobi {
}
# Find first D in sequence (5,-7,9,-11,13,-15,...) where (D|N) == -1
-sub _find_jacobi_d_sequence {
+sub _lucas_selfridge_params {
my($n) = @_;
# D is typically quite small: 67 max for N < 10^19. However, it is
@@ -913,39 +913,46 @@ sub _find_jacobi_d_sequence {
while (1) {
my $gcd = (ref($n) eq 'Math::BigInt') ? Math::BigInt::bgcd($d, $n)
: _gcd_ui($d, $n);
- return 0 if $gcd > 1 && $gcd != $n; # Found divisor $d
+ return (0,0,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);
+ my $D = $sign * $d;
+ my $P = 1;
+ my $Q = int( (1 - $D) / 4 );
+ ($P, $Q, $D)
}
-sub is_lucas_pseudoprime {
- return is_strong_lucas_pseudoprime($_[0], 'weak');
-}
+sub _lucas_extrastrong_params {
+ my($n) = @_;
-sub is_strong_lucas_pseudoprime {
- my($n, $doweak) = @_;
- _validate_positive_integer($n);
- $doweak = (defined $doweak && $doweak eq 'weak') ? 1 : 0;
+ my ($P, $Q, $D) = (3, 1, 5);
+ while (1) {
+ my $gcd = (ref($n) eq 'Math::BigInt') ? Math::BigInt::bgcd($D, $n)
+ : _gcd_ui($D, $n);
+ return (0,0,0) if $gcd > 1 && $gcd != $n; # Found divisor $d
+ last if _jacobi($D, $n) == -1;
+ $P++;
+ $D = $P*$P - 4;
+ }
+ ($P, $Q, $D);
+}
- return 1 if $n == 2;
- return 0 if $n < 2 || ($n % 2) == 0;
- return 0 if _is_perfect_square($n);
+# returns U_k, V_k, Q_k all mod n
+sub lucas_sequence {
+ my($n, $P, $Q, $k) = @_;
- # Determine Selfridge D, P, and Q parameters
- my $D = _find_jacobi_d_sequence($n);
- return 0 if $D == 0; # We found a divisor in the sequence
- my $P = 1;
- my $Q = int( (1 - $D) / 4 );
- # Verify we've calculated this right
- die "Selfridge error: $D, $P, $Q\n" if ($D != $P*$P - 4*$Q);
- #warn "N: $n D: $D P: $P Q: $Q\n";
+ return (0, 2) if $k == 0;
+ my $D = $P*$P - 4*$Q;
- # It's now time to perform the Lucas pseudoprimality test using $D.
+ 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;
+ croak "lucas_sequence: D is zero" if $D == 0;
if (ref($n) ne 'Math::BigInt') {
if (!defined $Math::BigInt::VERSION) {
@@ -955,46 +962,97 @@ sub is_strong_lucas_pseudoprime {
$n = Math::BigInt->new("$n");
}
- my $m = $n->copy->badd(1);
- # Traditional d,s:
- # my $d=$m->copy; my $s=0; while ($d->is_even) { $s++; $d->brsft(1); }
- # die "Invalid $m, $d, $s\n" unless $m == $d * 2**$s;
- my $dstr = substr($m->as_bin, 2);
- my $s = 0;
- if (!$doweak) {
- $dstr =~ s/(0*)$//;
- $s = length($1);
- }
-
my $ZERO = $n->copy->bzero;
my $U = $ZERO + 1;
my $V = $ZERO + $P;
my $Qk = $ZERO + $Q;
- my $bpos = 0;
- while (++$bpos < length($dstr)) {
- $U = ($U * $V) % $n;
- $V = ($V * $V - 2*$Qk) % $n;
- $Qk = ($Qk * $Qk) % $n;
- if (substr($dstr,$bpos,1)) {
- my $T1 = $U->copy->bmul($D);
- $U->badd($V);
- $U->badd($n) if $U->is_odd;
- $U->brsft(1);
- $U->bmod($n);
- $V->badd($T1);
- $V->badd($n) if $V->is_odd;
- $V->brsft(1);
- $V->bmod($n);
+ $k = Math::BigInt->new("$k") unless ref($k) eq 'Math::BigInt';
+ my $kstr = substr($k->as_bin, 2);
+ my $bpos = 0;
- $Qk = ($Qk * $Q) % $n;
+ 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);
+ }
+ }
+ } else {
+ while (++$bpos < length($kstr)) {
+ $U = ($U * $V) % $n;
+ $V = ($V * $V - 2*$Qk) % $n;
+ $Qk->bmul($Qk);
+ 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);
+
+ $Qk->bmul($Q);
+ }
+ $Qk->bmod($n);
}
}
- if ($doweak) { return $U->is_zero ? 1 : 0; }
- return 1 if $U->is_zero || $V->is_zero;
+ $U->bmod($n);
+ $V->bmod($n);
+ return ($U, $V, $Qk);
+}
- # Compute powers of V
- foreach my $r (1 .. $s-1) {
+sub is_lucas_pseudoprime {
+ my($n) = @_;
+ _validate_positive_integer($n);
+
+ 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_selfridge_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($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $n+1);
+ return $U->is_zero ? 1 : 0;
+}
+
+sub is_strong_lucas_pseudoprime {
+ my($n) = @_;
+ _validate_positive_integer($n);
+
+ 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_selfridge_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;
+ }
+ my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k);
+
+ return 1 if $U->is_zero || $V->is_zero;
+ foreach my $r (1 .. $s-1) { # Compute powers of V
$V = ($V * $V - 2*$Qk) % $n;
return 1 if $V->is_zero;
$Qk = ($Qk * $Qk) % $n if $r < ($s-1);
@@ -1010,53 +1068,18 @@ sub is_extra_strong_lucas_pseudoprime {
return 0 if $n < 2 || ($n % 2) == 0;
return 0 if _is_perfect_square($n);
- my ($P, $Q, $D) = (3, 1, 5);
- while (1) {
- my $gcd = (ref($n) eq 'Math::BigInt') ? Math::BigInt::bgcd($D, $n)
- : _gcd_ui($D, $n);
- return 0 if $gcd > 1 && $gcd != $n; # Found divisor $d
- last if _jacobi($D, $n) == -1;
- $P++;
- $D = $P*$P - 4;
- }
- die "Lucas incorrect DPQ: $D, $P, $Q\n" if ($D != $P*$P - 4*$Q);
+ 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);
- 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 $m = $n+1;
+ my($s, $k) = (0, $m);
+ while ( $k > 0 && !($k % 2) ) {
+ $s++;
+ $k >>= 1;
}
+ my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k);
- my $m = $n->copy->badd(1);
- # Traditional d,s:
- # my $d=$m->copy; my $s=0; while ($d->is_even) { $s++; $d->brsft(1); }
- # die "Invalid $m, $d, $s\n" unless $m == $d * 2**$s;
- my $dstr = substr($m->as_bin, 2);
- $dstr =~ s/(0*)$//;
- my $s = length($1);
-
- my $ZERO = $n->copy->bzero;
- my $U = $ZERO + 1;
- my $V = $ZERO + $P;
- my $bpos = 0;
- while (++$bpos < length($dstr)) {
- $U->bmul($V)->bmod($n);
- $V->bmul($V)->bsub(2)->bmod($n);
- if (substr($dstr,$bpos,1)) {
- my $U_times_D = $U->copy->bmul($D);
- $U->bmul($P)->badd($V);
- $U->badd($n) if $U->is_odd;
- $U->brsft(1);
-
- $V->bmul($P)->badd($U_times_D);
- $V->badd($n) if $V->is_odd;
- $V->brsft(1);
- }
- }
- $U->bmod($n);
- $V->bmod($n);
return 1 if $U->is_zero && ($V == 2 || $V == ($n-2));
return 1 if $V->is_zero;
@@ -2794,6 +2817,22 @@ C<n-1>, then attempts a proof using theorem 5 of Brillhart, Lehmer, and
Selfridge's 1975 paper. This can take a long time to return if given a
composite, though it should not be anywhere near as long as the Lucas test.
+=head2 lucas_sequence
+
+ my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k)
+
+Computes C<U_k>, C<V_k>, and C<Q_k> for the Lucas seqence defined by
+C<P>,C<Q>, modulo C<n>. The modular Lucas sequence is used in a
+number of primality tests and proofs.
+
+The following conditions must hold:
+ - C<< D = P*P - 4*Q != 0 >>
+ - C<< P > 0 >>
+ - C<< P < n >>
+ - C<< Q < n >>
+ - C<< k >= 0 >>
+ - C<< n >= 2 >>
+
=head1 UTILITY FUNCTIONS
diff --git a/t/17-pseudoprime.t b/t/17-pseudoprime.t
index dea5d93..3a7caed 100644
--- a/t/17-pseudoprime.t
+++ b/t/17-pseudoprime.t
@@ -6,7 +6,8 @@ use Test::More;
use Math::Prime::Util qw/is_prime
is_pseudoprime is_strong_pseudoprime
is_lucas_pseudoprime is_strong_lucas_pseudoprime
- is_extra_strong_lucas_pseudoprime/;
+ is_extra_strong_lucas_pseudoprime
+ lucas_sequence/;
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
@@ -72,6 +73,17 @@ foreach my $ppref (values %pseudoprimes) {
}
my @small_lucas_trials = (2, 9, 16, 100, 102, 2047, 2048, 5781, 9000, 14381);
+my %lucas_sequences = (
+ "323 1 1 324" => [0,2,1],
+ "323 4 1 324" => [170,308,1],
+ "323 4 5 324" => [194,156,115],
+ "323 3 1 324" => [0,2,1],
+ "323 3 1 81" => [0,287,1],
+ "323 5 -1 81" => [153,195,322],
+ "49001 25 117 24501" => [20933,18744,19141],
+ "18971 10001 -1 4743" => [5866,14421,18970],
+);
+
plan tests => 0 + 3
+ 4
+ $num_pseudoprimes
@@ -79,6 +91,7 @@ plan tests => 0 + 3
+ 1 # mr base 2 2-4k
+ 9 # mr with large bases
+ scalar @small_lucas_trials
+ + scalar(keys %lucas_sequences)
+ 1*$extra;
ok(!eval { is_strong_pseudoprime(2047); }, "MR with no base fails");
@@ -165,3 +178,11 @@ if ($extra) {
}
is($mr2fail, 0, "is_strong_pseudoprime bases 2,3 matches is_prime to 1,373,652");
}
+
+# Lucas sequences, used for quite a few tests
+sub lucas_sequence_to_native {
+ map { (ref($_) eq 'Math::BigInt') ? int($_->bstr) : $_ } lucas_sequence(@_);
+}
+while (my($params, $expect) = each (%lucas_sequences)) {
+ is_deeply( [lucas_sequence_to_native(split(' ', $params))], $expect, "Lucas sequence $params" );
+}
--
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