[libmath-prime-util-perl] 26/72: Add znorder function for multiplicative order
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:49:38 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list