[libmath-prime-util-perl] 52/181: Move znorder to XS; add tests for znprimroot
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:06 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.36
in repository libmath-prime-util-perl.
commit ba82ffc64622c4df8f378cc845b4f4486fc785ea
Author: Dana Jacobsen <dana at acm.org>
Date: Sat Dec 28 01:24:22 2013 -0800
Move znorder to XS; add tests for znprimroot
---
XS.xs | 36 +++++++++++++++++++++++++++++++++
factor.c | 31 ++++++++++++-----------------
lib/Math/Prime/Util.pm | 54 +++++++++++++++++++++++++-------------------------
t/19-moebius.t | 41 +++++++++++++++++++++++++++++++++++++-
util.c | 29 +++++++++++++++++++++++++++
util.h | 1 +
6 files changed, 146 insertions(+), 46 deletions(-)
diff --git a/XS.xs b/XS.xs
index 21903a5..1576e72 100644
--- a/XS.xs
+++ b/XS.xs
@@ -739,6 +739,42 @@ carmichael_lambda(IN SV* svn)
_vcallsub("Math::Prime::Util::_generic_carmichael_lambda");
XSRETURN(1);
+int
+znorder(IN SV* sva, IN SV* svn)
+ PREINIT:
+ int astatus, nstatus;
+ PPCODE:
+ astatus = _validate_int(sva, 0);
+ nstatus = _validate_int(svn, 0);
+ if (astatus == 1 && nstatus == 1) {
+ UV a, n, order;
+ set_val_from_sv(a, sva);
+ set_val_from_sv(n, svn);
+ order = znorder(a, n);
+ if (order == 0) XSRETURN_UNDEF;
+ XSRETURN_UV(order);
+ } else {
+ dTHX;
+ dSP;
+ int count;
+ ENTER;
+ SAVETMPS;
+ (void) POPs; (void) POPs;
+ PUSHMARK(SP);
+ XPUSHs(sva);
+ XPUSHs(svn);
+ PUTBACK;
+ count = call_pv("Math::Prime::Util::_generic_znorder", G_SCALAR);
+ SPAGAIN;
+ if (count != 1) croak("callback sub should return one value");
+ TOPs = SvREFCNT_inc(TOPs);
+ PUTBACK;
+ FREETMPS;
+ LEAVE;
+ TOPs = sv_2mortal(TOPs);
+ XSRETURN(1);
+ }
+
void
_XS_moebius(IN UV lo, IN UV hi = 0)
PREINIT:
diff --git a/factor.c b/factor.c
index 4e0c337..9315961 100644
--- a/factor.c
+++ b/factor.c
@@ -396,31 +396,26 @@ UV _XS_divisor_sum(UV n, UV k)
}
} else if (k == 1) {
for (i = 0; i < nfac; i++) {
- UV e = 1, f = factors[i];
- UV fmult = 1 + f;
- while (i+1 < nfac && f == factors[i+1]) { e++; i++; }
- if (e > 1) {
- UV pke = f;
- for (j = 1; j < (int)e; j++) {
- pke *= f;
- fmult += pke;
- }
+ UV f = factors[i];
+ UV pke = f, fmult = 1 + f;
+ while (i+1 < nfac && f == factors[i+1]) {
+ pke *= f;
+ fmult += pke;
+ i++;
}
product *= fmult;
}
} else {
for (i = 0; i < nfac; i++) {
- UV e = 1, f = factors[i];
- UV fmult, pk = f;
+ UV f = factors[i];
+ UV fmult, pke, pk = f;
for (j = 1; j < (int)k; j++) pk *= f;
- while (i+1 < nfac && f == factors[i+1]) { e++; i++; }
fmult = 1 + pk;
- if (e > 1) {
- UV pke = pk;
- for (j = 1; j < (int)e; j++) {
- pke *= pk;
- fmult += pke;
- }
+ pke = pk;
+ while (i+1 < nfac && f == factors[i+1]) {
+ pke *= pk;
+ fmult += pke;
+ i++;
}
product *= fmult;
}
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 5f43c73..4fb7262 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -107,6 +107,7 @@ BEGIN {
*exp_mangoldt = \&Math::Prime::Util::_generic_exp_mangoldt;
*carmichael_lambda = \&Math::Prime::Util::_generic_carmichael_lambda;
*kronecker = \&Math::Prime::Util::_generic_kronecker;
+ *znorder = \&Math::Prime::Util::_generic_znorder;
*znprimroot = \&Math::Prime::Util::_generic_znprimroot;
*forprimes = sub (&$;$) { _generic_forprimes(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
*fordivisors = sub (&$) { _generic_fordivisors(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
@@ -1678,7 +1679,7 @@ sub _generic_carmichael_lambda {
return $lcm;
}
-sub znorder {
+sub _generic_znorder {
my($a, $n) = @_;
_validate_num($a) || _validate_positive_integer($a);
_validate_num($n) || _validate_positive_integer($n);
@@ -1688,33 +1689,25 @@ sub znorder {
# Sadly, Calc/FastCalc are horrendously slow for this function.
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 (divisors(euler_phi($n), $n).
- # Much better, but slow with many divisors.
- # Method 3: check all k in (divisors(carmichael_lambda($n), $n).
- # Stronger result, but still can have too many divisors.
- # Method 4: Abhijit Das, alg 1.7, just enough multiples of p.
- # Combine this with method 3.
- #
- # 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.
+ # The answer is one of the divisors of phi(n) and lambda(n).
+ my $lambda = carmichael_lambda($n);
$a = Math::BigInt->new("$a") unless ref($a) eq 'Math::BigInt';
- # Method 3:
- # foreach my $k (divisors(carmichael_lambda($n))) {
- # return $k if $k != 1 && $a->copy->bmodpow("$k", $n)->is_one;
- # }
- # return;
+ # This is easy and usually fast, but can bog down with too many divisors.
+ if ($lambda <= $_XS_MAXVAL) {
+ foreach my $k (_XS_divisors($lambda)) {
+ return $k if $a->copy->bmodpow("$k", $n)->is_one;
+ }
+ return;
+ }
- # Method 4 (Das algorithm 1.7 applied to Carmichael Lambda)
- #my $phi = Math::BigInt->new('' . euler_phi($n));
- my $phi = Math::BigInt->new('' . carmichael_lambda($n));
- my @pe = ($phi <= $_XS_MAXVAL) ? _XS_factor_exp($phi) : factor_exp($phi);
+ # Algorithm 1.7 from A. Das applied to Carmichael Lambda.
+ $lambda = Math::BigInt->new("$lambda") unless ref($lambda) eq 'Math::BigInt';
my $k = Math::BigInt->bone;
- foreach my $f (@pe) {
+ foreach my $f (factor_exp($lambda)) {
my($pi, $ei, $enum) = (Math::BigInt->new("$f->[0]"), $f->[1], 0);
- my $phidiv = $phi / ($pi**$ei);
+ my $phidiv = $lambda / ($pi**$ei);
my $b = $a->copy->bmodpow($phidiv, $n);
while ($b != 1) {
return if $enum++ >= $ei;
@@ -3856,11 +3849,18 @@ C<MultiplicativeOrder[n]> function.
=head2 znprimroot
-Given a positive integer C<n>, returns a primitive root of C<(Z/nZ)^*>.
-A root exists when C<euler_phi($n) == carmichael_lambda($n)>, which will
-be true for all prime C<n> and some composites.
-(L<OEIS A033948|http://oeis.org/A033948>) is the list of such values.
-If a primitive root does not exist for C<n>, this function returns undef.
+Given a positive integer C<n>, returns the smallest primitive root
+of C<(Z/nZ)^*>, or C<undef> if no root exists. A root exists when
+C<euler_phi($n) == carmichael_lambda($n)>, which will be true for
+all prime C<n> and some composites.
+
+L<OEIS A033948|http://oeis.org/A033948> is a sequence of integers where
+the primitive root exists, while L<OEIS A046145|http://oeis.org/A046145>
+is a list of the smallest primitive roots, which is what this function
+produces.
+
+This will always produce the smallest primitive root. Pari does not
+necessarily produce the smallest result for composites.
=head2 random_prime
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 451ecc1..1423ccb 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -5,7 +5,9 @@ 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 znorder liouville/;
+ chebyshev_theta chebyshev_psi carmichael_lambda znorder liouville
+ znprimroot
+ /;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
@@ -192,6 +194,38 @@ my @mult_orders = (
[31407,2147475467,266],
);
+my %primroots = (
+ -11 => 2,
+ 0 => undef,
+ 1 => 0,
+ 2 => 1,
+ 3 => 2,
+ 4 => 3,
+ 5 => 2,
+ 6 => 5,
+ 7 => 3,
+ 8 => undef,
+ 9 => 2,
+ 10 => 3, # 3 is the smallest root. Pari gives the other root 7.
+ 1729 => undef,
+ 5109721 => 94,
+ 17551561 => 97,
+ 90441961 => 113,
+1407827621 => 2,
+1520874431 => 17,
+1685283601 => 164,
+);
+if ($use64) {
+ $primroots{2232881419280027} = 6;
+ $primroots{14123555781055773271} = 6;
+ $primroots{89637484042681} = 335;
+}
+
+# znprimroot:
+# 1729
+# 2232881419280027 factor divide goes to FP
+# 14123555781055773271 bmodpow hits RT 71548
+
# These are slow with XS, and *really* slow with PP.
if (!$usexs) {
%big_mertens = map { $_ => $big_mertens{$_} }
@@ -219,6 +253,7 @@ plan tests => 0 + 1
+ 7 + scalar(keys %totients)
+ 1 # Small Carmichael Lambda
+ scalar(@mult_orders)
+ + scalar(keys %primroots)
+ scalar(keys %jordan_totients)
+ 2 # Dedekind psi calculated two ways
+ 2 # Calculate J5 two different ways
@@ -365,6 +400,10 @@ foreach my $moarg (@mult_orders) {
my $zn = znorder($a, $n);
is( $zn, $exp, "znorder($a, $n) = " . ((defined $exp) ? $exp : "<undef>") );
}
+###### znprimroot
+while (my($n, $root) = each (%primroots)) {
+ is( znprimroot($n), $root, "znprimroot($n) == " . ((defined $root) ? $root : "<undef>") );
+}
###### liouville
foreach my $i (@liouville_pos) {
is( liouville($i), 1, "liouville($i) = 1" );
diff --git a/util.c b/util.c
index dd05c95..9032028 100644
--- a/util.c
+++ b/util.c
@@ -979,6 +979,35 @@ UV carmichael_lambda(UV n) {
}
return lambda;
}
+
+UV znorder(UV a, UV n) {
+ UV fac[MPU_MAX_FACTORS+1];
+ UV exp[MPU_MAX_FACTORS+1];
+ int i, j, nfactors;
+ UV phi, k = 1;
+
+ if (n == 0 || a == 0) return 0;
+ if (n == 1 || a == 1) return 1;
+ if (gcd_ui(a,n) > 1) return 0;
+
+ /* Abhijit Das, algorithm 1.7, applied to Carmichael Lambda */
+ phi = carmichael_lambda(n);
+ nfactors = factor_exp(phi, fac, exp);
+ for (i = 0; i < nfactors; i++) {
+ UV b, ek;
+ UV pi = fac[i], ei = exp[i];
+ UV phidiv = phi / pi;
+ for (j = 1; j < ei; j++) phidiv /= pi;
+ b = powmod(a, phidiv, n);
+ ek = 0;
+ while (b != 1) {
+ if (ek++ >= ei) return 0;
+ b = powmod(b, pi, n);
+ k *= pi;
+ }
+ }
+ return k;
+}
UV znprimroot(UV n) {
UV fac[MPU_MAX_FACTORS+1];
diff --git a/util.h b/util.h
index 520ea76..1cb51d5 100644
--- a/util.h
+++ b/util.h
@@ -34,6 +34,7 @@ extern int kronecker_ss(IV a, IV b);
extern UV totient(UV n);
extern UV carmichael_lambda(UV n);
extern UV znprimroot(UV n);
+extern UV znorder(UV a, UV n);
/* Above this value, is_prime will do deterministic Miller-Rabin */
/* With 64-bit math, we can do much faster mulmods from 2^16-2^32 */
--
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