[libmath-prime-util-perl] 49/181: Add kronecker, totient, carmichael_lambda, znprimroot

[Partha P. Mukherjee]

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 272b94859bd28b9a36c9d01fbfda7d73e27e10e5
Author: Dana Jacobsen <dana at acm.org>
Date:   Fri Dec 27 18:21:50 2013 -0800

    Add kronecker, totient, carmichael_lambda, znprimroot
---
 util.c | 115 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 util.h |  18 ++++++++++-
 2 files changed, 132 insertions(+), 1 deletion(-)

diff --git a/util.c b/util.c
index 510426e..dd05c95 100644
--- a/util.c
+++ b/util.c
@@ -57,11 +57,14 @@
 
 #include "ptypes.h"
 #define FUNC_isqrt 1
+#define FUNC_lcm_ui 1
 #include "util.h"
 #include "sieve.h"
 #include "primality.h"
 #include "cache.h"
 #include "lmo.h"
+#include "factor.h"
+#include "mulmod.h"
 
 static int _verbose = 0;
 void _XS_set_verbose(int v) { _verbose = v; }
@@ -885,6 +888,118 @@ IV _XS_mertens(UV n) {
   return sum;
 }
 
+
+static int padic2(UV n) {  /* How many times does 2 divide n? */
+  UV p = 0;
+  while (n && n % 2 == 0) { n >>= 1;  p++; }
+  return p;
+}
+#define IS_MOD8_3OR5(x)  (((x)&7)==3 || ((x)&7)==5)
+
+static int kronecker_uu_sign(UV a, UV b, int s) {
+  while (a) {
+    int r = padic2(a);
+    if (r) {
+      if ((r&1)  &&  IS_MOD8_3OR5(b))  s = -s;
+      a >>= r;
+    }
+    if (a & b & 2)  s = -s;
+    { UV t = b % a;  b = a;  a = t; }
+  }
+  return (b == 1) ? s : 0;
+}
+
+int kronecker_uu(UV a, UV b) {
+  int r, s;
+  if (b & 1)   return kronecker_uu_sign(a, b, 1);
+  if (!(a&1))  return 0;
+  s = 1;
+  r = padic2(b);
+  if (r) {
+    if ((r&1) && IS_MOD8_3OR5(a))  s = -s;
+    b >>= r;
+  }
+  return kronecker_uu_sign(a, b, s);
+}
+
+int kronecker_su(IV a, UV b) {
+  int r, s;
+  if (a >= 0)  return kronecker_uu(a, b);
+  if (b == 0)  return (a == 1 || a == -1) ? 1 : 0;
+  s = 1;
+  r = padic2(b);
+  if (r) {
+    if (!(a&1))  return 0;
+    if ((r&1) && IS_MOD8_3OR5(a))  s = -s;
+    b >>= r;
+  }
+  a %= (IV) b;
+  if (a < 0)  a += b;
+  return kronecker_uu_sign(a, b, s);
+}
+
+int kronecker_ss(IV a, IV b) {
+  if (a >= 0 && b >= 0)
+    return (b & 1)  ?  kronecker_uu_sign(a, b, 1)  :  kronecker_uu(a,b);
+  if (b >= 0)
+    return kronecker_su(a, b);
+  return kronecker_su(a, -b) * ((a < 0) ? -1 : 1);
+}
+
+UV totient(UV n) {
+  UV i, facs[MPU_MAX_FACTORS+1];
+  UV nfacs = factor(n, facs);
+  UV totient = 1;
+  UV lastf = 0;
+  if (n <= 1) return n;
+  for (i = 0; i < nfacs; i++) {
+    UV f = facs[i];
+    if (f == lastf) { totient *= f;               }
+    else            { totient *= f-1;  lastf = f; }
+  }
+  return totient;
+}
+
+UV carmichael_lambda(UV n) {
+  UV fac[MPU_MAX_FACTORS+1];
+  UV exp[MPU_MAX_FACTORS+1];
+  int i, j, nfactors;
+  UV lambda = 1;
+
+  if (n < 8) return totient(n);
+  if ((n & (n-1)) == 0)  return totient(n)/2;
+
+  nfactors = factor_exp(n, fac, exp);
+  if (fac[0] == 2 && exp[0] > 2)  exp[0]--;
+  for (i = 0; i < nfactors; i++) {
+    UV pk = fac[i]-1;
+    for (j = 1; j < exp[i]; j++)
+      pk *= fac[i];
+    lambda = lcm_ui(lambda, pk);
+  }
+  return lambda;
+}
+  
+UV znprimroot(UV n) {
+  UV fac[MPU_MAX_FACTORS+1];
+  UV exp[MPU_MAX_FACTORS+1];
+  UV a, phi;
+  int i, nfactors;
+  if (n <= 4) return (n == 0) ? 0 : n-1;
+  phi = totient(n);
+  nfactors = factor_exp(phi, fac, exp);
+  for (i = 0; i < nfactors; i++)
+    exp[i] = phi / fac[i];  /* exp[i] = phi(n) / i-th-factor-of-phi(n) */
+  for (a = 2; a < n; a++) {
+    if (kronecker_uu(a, n) == 0)  continue;
+    for (i = 0; i < nfactors; i++) {
+      if (powmod(a, exp[i], n) == 1) break;
+    }
+    if (i == nfactors) return a;
+  }
+}
+
+
 double _XS_chebyshev_theta(UV n)
 {
   KAHAN_INIT(sum);
diff --git a/util.h b/util.h
index 4dbcc6e..b420c20 100644
--- a/util.h
+++ b/util.h
@@ -27,6 +27,15 @@ extern long double ld_riemann_zeta(long double x);
 extern double _XS_RiemannR(double x);
 extern UV _XS_Inverse_Li(UV x);
 
+extern int kronecker_uu(UV a, UV b);
+extern int kronecker_su(IV a, UV b);
+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 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 */
 #if (BITS_PER_WORD == 64) || HAVE_STD_U64
@@ -72,7 +81,7 @@ static UV icbrt(UV n) {
 }
 #endif
 
-#ifdef FUNC_gcd_ui
+#if defined(FUNC_gcd_ui) || defined(FUNC_lcm_ui)
 static UV gcd_ui(UV x, UV y) {
   UV t;
   if (y < x) { t = x; x = y; y = t; }
@@ -83,4 +92,11 @@ static UV gcd_ui(UV x, UV y) {
 }
 #endif
 
+#ifdef FUNC_lcm_ui
+static UV lcm_ui(UV x, UV y) {
+  /* Can overflow if lcm(x,y) > 2^64 (e.g. two primes each > 2^32) */
+  return x * (y / gcd_ui(x,y));
+}
+#endif
+
 #endif

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git