This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.32 in repository libmath-prime-util-perl.
commit 3c14ce7d1e5f574be4aa0c5fc756b1300d2d4084 Author: Dana Jacobsen <d...@acm.org> Date: Thu Sep 12 17:55:57 2013 -0700 Add znorder function for multiplicative order --- Changes | 1 + TODO | 5 --- lib/Math/Prime/Util.pm | 107 ++++++++++++++++++++++++++++++++++------------ lib/Math/Prime/Util/PP.pm | 27 +++++++----- t/19-moebius.t | 34 ++++++++++++++- t/81-bignum.t | 8 +++- 6 files changed, 134 insertions(+), 48 deletions(-) diff --git a/Changes b/Changes index ef8973f..66d1318 100644 --- a/Changes +++ b/Changes @@ -6,6 +6,7 @@ Revision history for Perl module Math::Prime::Util - is_proven_prime - is_proven_prime_with_cert - carmichael_lambda + - znorder - random_nbit_prime now uses Fouque and Tibouchi A1. Slightly better uniformity and typically a bit faster. diff --git a/TODO b/TODO index 97bb00f..0ffcdba 100644 --- a/TODO +++ b/TODO @@ -42,8 +42,3 @@ - Rewrite 23-primality-proofs.t for new format (keep some of the old tests?). - Use Montgomery routines in more places: Factoring and Lucas tests. - -- Add a group order function. We have a naieve method used in AKS (where the - values are so small it doesn't matter). See: - http://cs-haven.blogspot.com/2012/02/efficiently-computing-order-of-element.html - for some discussion of different methods. diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index a2d7743..9f971da 100644 --- a/lib/Math/Prime/Util.pm +++ b/lib/Math/Prime/Util.pm @@ -40,7 +40,7 @@ our @EXPORT_OK = moebius mertens euler_phi jordan_totient exp_mangoldt chebyshev_theta chebyshev_psi divisor_sum - carmichael_lambda + carmichael_lambda znorder ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR ); our %EXPORT_TAGS = (all => [ @EXPORT_OK ]); @@ -1414,8 +1414,6 @@ sub jordan_totient { # Mathematica and Pari both have functions like this. sub divisor_sum { my($n, $sub) = @_; - # I really need to get cracking on an XS validator. - #return _XS_divisor_sum($n) if !defined $sub && defined $n && $n <= $_XS_MAXVAL && $_Config{'nobigint'}; return (0,1)[$n] if defined $n && $n <= 1; _validate_num($n) || _validate_positive_integer($n); @@ -1437,8 +1435,9 @@ sub divisor_sum { return $product; } + # TODO: Alternately the argument could be the k for sigma, so this + # should be 0, the above should be 1, etc. if (ref($sub) ne 'CODE' && int($sub) == 1) { - return 1 if $n == 1; my ($product, $exponent) = (1, 1); my @factors = ($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n); while (@factors) { @@ -1558,16 +1557,64 @@ sub carmichael_lambda { eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; } or do { croak "Cannot load Math::BigInt"; }; } + @factors = map { Math::BigInt->new("$_") } @factors; my $lcm = Math::BigInt::blcm( - map { Math::BigInt->new("$_") - ->bpow($factor_mult{$_}-1) - ->bmul($_-1) - } @factors + map { $_->copy->bpow($factor_mult{$_}-1)->bmul($_-1) } @factors ); $lcm = int($lcm->bstr) if $lcm->bacmp(''.~0) <= 0; return $lcm; } +sub znorder { + my($a, $n) = @_; + _validate_num($a) || _validate_positive_integer($a); + _validate_num($n) || _validate_positive_integer($n); + return if $n <= 0; + return (undef,1)[$a] if $a <= 1; + return 1 if $n == 1; + if (!defined $Math::BigInt::VERSION) { + eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; } + or do { croak "Cannot load Math::BigInt"; }; + } + return if Math::BigInt::bgcd($a, $n) > 1; + # Method 1: check all a^k 1 .. $n-1. + # Naive and terrible slow. + # Method 2: check all k in (all_factors(euler_phi($n), $n). + # Method 3: check all k in (all_factors(carmichael_lambda($n), $n). + # Good for most cases, but slow when factor quantity is large. + # Method 4: Das algorithm 1.7, just enough multiples of p + # Fastest. + # + # Most of the time is spent factoring $n and $phi. We could do the phi + # construction here and build its factors to save a little more time. + + $a = Math::BigInt->new("$a") unless ref($a) eq 'Math::BigInt'; + my $phi = Math::BigInt->new('' . euler_phi($n)); + my %e; + my @p = grep { !$e{$_}++ } + ($phi <= $_XS_MAXVAL) ? _XS_factor($phi) : factor($phi); + my $k = Math::BigInt->bone; + foreach my $i (0 .. $#p) { + my($pi, $ei, $enum) = (Math::BigInt->new("$p[$i]"), $e{$p[$i]}, 0); + my $phidiv = $phi / ($pi**$ei); + my $b = $a->copy->bmodpow($phidiv, $n); + while ($b != 1) { + return if $enum++ >= $ei; + $b = $b->copy->bmodpow($pi, $n); + $k *= $pi; + } + } + return $k; + + # Method 3: + # my $cl = carmichael_lambda($n); + # foreach my $k (all_factors($cl), $cl) { + # my $t = $a->copy->bmodpow("$k", $n); + # return $k if $t->is_one; + # } + # return; +} + @@ -1850,27 +1897,21 @@ sub lucas_sequence { ############################################################################# - # Oct 2012 note: these numbers are old. + # Simple timings, do { is_strong_pseudoprime(2*$_+1,2) } for 1..1000000; # - # Timings for various combinations, given the current possibilities of: - # 1) XS MR optimized (either x86-64, 32-bit on 64-bit mach, or half-word) - # 2) XS MR non-optimized (big input not on 64-bit machine) - # 3) PP MR with small input (non-bigint Perl) - # 4) PP MR with large input (using functions for mulmod) - # 5) PP MR with full bigints - # 6) PP Lucas with small input - # 7) PP Lucas with large input - # 8) PP Lucas with full bigints + # 1.0 uS XS 32-bit input, is_strong_pseudoprime + # 1.8 uS XS 64-bit input, is_strong_pseudoprime + # 1.4 uS XS 32-bit input, is_strong_lucas_pseudoprime + # 3.0 uS XS 64-bit input, is_strong_lucas_pseudoprime + # 1.2 uS XS 32-bit input, is_almost_extra_strong_lucas_pseudoprime + # 2.3 uS XS 64-bit input, is_almost_extra_strong_lucas_pseudoprime # - # Time for one test: - # 0.5uS XS MR with small input - # 0.8uS XS MR with large input - # 7uS PP MR with small input - # 400uS PP MR with large input - # 5000uS PP MR with bigint - # 2700uS PP LP with small input - # 6100uS PP LP with large input - # 7400uS PP LP with bigint + # 12 uS Perl 32-bit input, is_strong_pseudoprime + # 1200 uS Perl 64-bit input, is_strong_pseudoprime + # 2000 uS Perl 32-bit input, is_strong_lucas_pseudoprime + # 7840 uS Perl 64-bit input, is_strong_lucas_pseudoprime + # 940 uS Perl 32-bit input, is_almost_extra_strong_lucas_pseudoprime + # 3360 uS Perl 64-bit input, is_almost_extra_strong_lucas_pseudoprime sub is_aks_prime { my($n) = @_; @@ -2340,7 +2381,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 superset sqrt(N) gcd(A^M +=for stopwords forprimes Möbius Deléglise totient moebius mertens znorder irand primesieve uniqued k-tuples von SoE pari yafu fonction qui compte le nombre nombres voor PhD superset sqrt(N) gcd(A^M =head1 NAME @@ -3501,6 +3542,16 @@ positive integer m such that a^m = 1 mod n for every integer a coprime to n. This is L<OEIS series A002322|http://oeis.org/A002322>. +=head2 znorder + + $order = znorder(2, next_prime(10**19)-6); + +Given two positive integers C<a> and C<n>, returns the multiplicative order +of C<a> modulo <n>. This is the smallest positive integer C<k> such that +C<a^k ≡ 1 mod n>. Returns 1 if C<a = 1>. Return undef if C<a = 0> or if +C<a> and C<n> are not coprime, since no value will result in 1 mod n. + + =head2 random_prime my $small_prime = random_prime(1000); # random prime <= limit diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm index e2036ce..551d5a7 100644 --- a/lib/Math/Prime/Util/PP.pm +++ b/lib/Math/Prime/Util/PP.pm @@ -198,8 +198,10 @@ sub _is_prime7 { # n must not be divisible by 2, 3, or 5 # Slow since it's all in PP and uses bigints. return 0 unless miller_rabin($n, 2); - return 0 unless is_extra_strong_lucas_pseudoprime($n); - return ($n <= 18446744073709551615) ? 2 : 1; + if ($n <= 18446744073709551615) { + return is_almost_extra_strong_lucas_pseudoprime($n) ? 2 : 0; + } + return is_extra_strong_lucas_pseudoprime($n) ? 1 : 0; } sub is_prime { @@ -1097,6 +1099,15 @@ sub is_extra_strong_lucas_pseudoprime { $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. + 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($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k); return 1 if $U->is_zero && ($V == 2 || $V == ($n-2)); @@ -1121,13 +1132,6 @@ sub is_almost_extra_strong_lucas_pseudoprime { 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; } @@ -1139,8 +1143,9 @@ sub is_almost_extra_strong_lucas_pseudoprime { 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 $kstr = substr($n->copy->badd(1)->as_bin, 2); + $kstr =~ s/(0*)$//; + my $s = length($1); my $bpos = 0; while (++$bpos < length($kstr)) { if (substr($kstr,$bpos,1)) { diff --git a/t/19-moebius.t b/t/19-moebius.t index 2e3796e..b44b888 100644 --- a/t/19-moebius.t +++ b/t/19-moebius.t @@ -5,7 +5,7 @@ use warnings; use Test::More; use Math::Prime::Util qw/moebius mertens euler_phi jordan_totient divisor_sum exp_mangoldt - chebyshev_theta chebyshev_psi carmichael_lambda/; + chebyshev_theta chebyshev_psi carmichael_lambda znorder/; my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING}; my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32; @@ -157,6 +157,30 @@ my %chebyshev2 = ( my @A002322 = (0,1,1,2,2,4,2,6,2,6,4,10,2,12,6,4,4,16,6,18,4,6,10,22,2,20,12,18,6,28,4,30,8,10,16,12,6,36,18,12,4,40,6,42,10,12,22,46,4,42,20,16,12,52,18,20,6,18,28,58,4,60,30,6,16,12,10,66,16,22,12,70,6,72,36,20,18,30,12,78,4,54,40,82,6,16,42,28,10,88,12,12,22,30,46,36,8,96,42,30,20,100,16,102,12,12,52,106,18,108,20,36,12,112,18,44,28,12,58,48,4,110,60,40,30,100,6,126,32,42,12,130,10,18,66,36,16,136,22,138,12,46,70,60,12,28,72,42,36,148,20,150,18,48,30,60,12,156,78,52,8,66,54,162,40,20, [...] +my @mult_orders = ( + [1, 35, 1], + [2, 35, 12], + [4, 35, 6], + [6, 35, 2], + [7, 35], + #[2,1000000000000031,81788975100], + [1, 1, 1], + [0, 0], + [1, 0], + [25, 0], + [1, 1, 1], + [19, 1, 1], + [1, 19, 1], + [2, 19, 18], + [3, 20, 4], + [57,1000000003,189618], + [67,999999749,30612237], + [22,999991815,69844], + [10,2147475467,31448382], + [141,2147475467,1655178], + [7410,2147475467,39409], + [31407,2147475467,266], +); plan tests => 0 + 1 + 1 # Small Moebius @@ -165,6 +189,7 @@ plan tests => 0 + 1 + 2 # Small Phi + 7 + scalar(keys %totients) + 1 # Small Carmichael Lambda + + scalar(@mult_orders) + scalar(keys %jordan_totients) + 2 # Dedekind psi calculated two ways + 1 # Calculate J5 two different ways @@ -291,7 +316,12 @@ while (my($n, $c2) = each (%chebyshev2)) { my @lambda = map { carmichael_lambda($_) } (0 .. $#A002322); is_deeply( \@lambda, \@A002322, "carmichael_lambda with range: 0, $#A000010" ); } - +###### znorder +foreach my $moarg (@mult_orders) { + my ($a, $n, $exp) = @$moarg; + my $zn = znorder($a, $n); + is( $zn, $exp, "znorder($a, $n) = " . ((defined $exp) ? $exp : "<undef>") ); +} sub cmp_closeto { my $got = shift; diff --git a/t/81-bignum.t b/t/81-bignum.t index a2f2c86..a422b3d 100644 --- a/t/81-bignum.t +++ b/t/81-bignum.t @@ -75,7 +75,7 @@ plan tests => 0 + 6*2*$extra # more PC tests + scalar(keys %factors) + scalar(keys %allfactors) - + 5 # moebius, euler_phi, jordan totient + + 6 # moebius, euler_phi, jordan totient + 15 # random primes + 0; @@ -98,6 +98,7 @@ use Math::Prime::Util qw/ euler_phi jordan_totient divisor_sum + znorder ExponentialIntegral LogarithmicIntegral RiemannR @@ -208,7 +209,7 @@ SKIP: { ############################################################################### SKIP: { - skip "Your 64-bit Perl is broken, skipping moebius and euler_phi tests", 5 if $broken64; + skip "Your 64-bit Perl is broken, skipping moebius and euler_phi tests", 6 if $broken64; my $n; $n = 618970019642690137449562110; is( moebius($n), -1, "moebius($n)" ); @@ -219,6 +220,9 @@ SKIP: { # 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!" ); + is( znorder(82734587234,927208363107752634625923555185111613055040823736157), + 4360156780036190093445833597286118936800, + "znorder" ); } ############################################################################### -- 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