Revamp znprimroot, can remove all powers.

Changes

@@ -248,6 +248,8 @@ Revision history for Perl module Math::Prime::Util
 
     - is_primitive_root faster.  Much faster for p^k and 2p^k with k >= 2.
 
+    - znprimroot faster.  2-20x faster for p^k and 2p^k with k >= 2.
+
     [OTHER]
 
     - Legendre, Meissel, and Lehmer prime counting got rewritten to use our

t/19-primroots.t

@@ -37,6 +37,10 @@ if ($use64) {
   $primroots{14123555781055773271} = 6;     # bmodpow hits RT 71548
   $primroots{89637484042681} = 335;         # smallest root is large
   $primroots{9223372036854775837} = 5;      # Pari #905
+  $primroots{36002292036481} = 13;
+  $primroots{72004584072962} = 13;
+  $primroots{2067900233973681742} = 17;
+  $primroots{8000468009126059319} = 13;
 }
 
 plan tests => 0

util.c

@@ -1363,48 +1363,58 @@ UV znorder(UV a, UV n) {
 UV znprimroot(UV n) {
   UV fac[MPU_MAX_FACTORS+1];
   UV phi_div_fac[MPU_MAX_FACTORS+1];
-  UV a, phi, on, r;
-  int i, nfactors;
+  UV k, p, phi, a, psquared;
+  int i, nfactors, isneven, ispow;
 
   if (n <= 4) return (n == 0) ? 0 : n-1;
   if (n % 4 == 0)  return 0;
 
-  on = (n&1) ? n : (n>>1);
-  a = powerof(on);
-  r = rootint(on, a);
-  if (!is_prob_prime(r)) return 0;        /* c^a or 2c^a */
-  phi = (r-1) * (on/r);                   /* p^a or 2p^a */
+  isneven = !(n & 1);
+  if (isneven)  n >>= 1;
+
+  k = powerof(n);
+  p = rootint(n, k);
+  if (p == 3 && isneven) return 5;
+  if (!is_prob_prime(p)) return 0;
+  ispow = (k > 1);
+
+  phi = p-1;  /* p an odd prime */
+  psquared = ispow ? p*p : 0;
 
   nfactors = factor_exp(phi, fac, 0);
-  for (i = 0; i < nfactors; i++)
+  for (i = 1; i < nfactors; i++)
     phi_div_fac[i] = phi / fac[i];
 
 #if USE_MONTMATH
-  if (n & 1) {
-    const uint64_t npi = mont_inverse(n),  mont1 = mont_get1(n);
-    for (a = 2; a < n; a++) {
-      if (a == 4 || a == 8 || a == 9) continue;  /* Skip some perfect powers */
-      /* Skip values we know can't be right: (a|n) = 0 (or 1 for odd primes) */
-      if (phi == n-1) { if (kronecker_uu(a, n) != -1) continue; }
-      else            { if (gcd_ui(a,n) != 1) continue; }
-      r = mont_geta(a, n);
-      for (i = 0; i < nfactors; i++)
-        if (mont_powmod(r, phi_div_fac[i], n) == mont1)
-          break;
-      if (i == nfactors) return a;
-    }
-  } else
-#endif
-  for (a = 2; a < n; a++) {
+ {
+  UV r;
+  const uint64_t npi = mont_inverse(p),  mont1 = mont_get1(p);
+  for (a = 2; a < p; a++) {
+    if (isneven && !(a&1)) continue;
     if (a == 4 || a == 8 || a == 9) continue;  /* Skip some perfect powers */
-    /* Skip values we know can't be right: (a|n) = 0 (or 1 for odd primes) */
-    if (phi == n-1) { if (kronecker_uu(a, n) != -1) continue; }
-    else            { if (gcd_ui(a,n) != 1) continue; }
-    for (i = 0; i < nfactors; i++)
-      if (powmod(a, phi_div_fac[i], n) == 1)
+    if (kronecker_uu(a, p) != -1) continue;
+    r = mont_geta(a, p);
+    for (i = 1; i < nfactors; i++)
+      if (mont_powmod(r, phi_div_fac[i], p) == mont1)
         break;
-    if (i == nfactors) return a;
+    if (i == nfactors)
+      if (!ispow || powmod(a, phi, psquared) != 1)
+        return a;
   }
+ }
+#else
+  for (a = 2; a < p; a++) {
+    if (isneven && !(a&1)) continue;
+    if (a == 4 || a == 8 || a == 9) continue;  /* Skip some perfect powers */
+    if (kronecker_uu(a, p) != -1) continue;
+    for (i = 1; i < nfactors; i++)
+      if (powmod(a, phi_div_fac[i], p) == 1)
+        break;
+    if (i == nfactors)
+      if (!ispow || powmod(a, phi, psquared) != 1)
+        return a;
+  }
+#endif
   return 0;
 }