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