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 <d...@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 _______________________________________________ Pkg-perl-cvs-commits mailing list Pkg-perl-cvs-commits@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits