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 <d...@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

_______________________________________________
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

Reply via email to