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 4d3882cb7ef54ef73370fdc6c40af3338cf9a773
Author: Dana Jacobsen <d...@acm.org>
Date: Fri Jan 10 14:11:08 2014 -0800
Add simple (and likely buggy) Pollard Rho DLP
---
TODO | 9 ++++++++
factor.c | 53 ++++++++++++++++++++++++++++++++++++++++++++
factor.h | 3 +++
t/19-moebius.t | 1 +
util.c | 69 +++++++++++++++++++++++++++++++++++++++++++++++++---------
util.h | 3 +++
6 files changed, 128 insertions(+), 10 deletions(-)
diff --git a/TODO b/TODO
index 6c8682a..2150c53 100644
--- a/TODO
+++ b/TODO
@@ -61,3 +61,12 @@
- look at sieve.c style prime walking
- Add Inverse Li and Legendre Phi to API?
+
+- Iterators speedup:
+ 1) config option for sieved next_prime. Very general, would make
+ next_prime run fast when called many times sequentially. Nasty
+ mixing with threads.
+ 2) iterator, PrimeIterator, or PrimeArray in XS using segment sieve.
+
+- Perhaps have main segment know the filled in range. That would allow
+ a sieved next_prime, and might speed up some counts and the like.
diff --git a/factor.c b/factor.c
index bbb8931..636aebc 100644
--- a/factor.c
+++ b/factor.c
@@ -920,3 +920,56 @@ int squfof_factor(UV n, UV *factors, UV rounds)
factors[0] = n;
return 1;
}
+
+UV dlp_trial(UV a, UV g, UV p, UV maxrounds) {
+ UV t, n = 1;
+ if (maxrounds > p) maxrounds = p;
+ for (n = 1; n < maxrounds; n++) {
+ t = powmod(g, n, p);
+ if (t == a)
+ return n;
+ }
+ return 0;
+}
+
+#define pollard_rho_cycle(u,v,w,p,n,a,g) \
+ switch (u % 3) { \
+ case 0: u = mulmod(u,u,p); v = mulmod(v,2,n); w = mulmod(w,2,n);
break;\
+ case 1: u = mulmod(u,a,p); v = addmod(v,1,n);
break;\
+ case 2: u = mulmod(u,g,p); w = addmod(w,1,n);
break;\
+ }
+
+UV dlp_prho(UV a, UV g, UV p, UV maxrounds) {
+ UV i;
+ UV n = znorder(g, p);
+ UV u=1, v=0, w=0;
+ UV U=u, V=v, W=w;
+ int const verbose = _XS_get_verbose();
+
+ if (verbose > 1 && n != p-1) printf("for g=%lu p=%lu, order is %lu\n", g, p,
n);
+ if (maxrounds > n) maxrounds = n;
+ for (i = 1; i < maxrounds; i++) {
+ pollard_rho_cycle(u,v,w,p,n,a,g); /* xi, ai, bi */
+ pollard_rho_cycle(U,V,W,p,n,a,g);
+ pollard_rho_cycle(U,V,W,p,n,a,g); /* x2i, a2i, b2i */
+ if (verbose > 3) printf( "%3lu %4lu %3lu %3lu %4lu %3lu %3lu\n", i, u,
v, w, U, V, W );
+ if (u == U) {
+ UV r1, r2, k;
+ r1 = submod(v, V, n);
+ if (r1 == 0) {
+ if (verbose) printf("DLP Rho failure, r=0\n");
+ return 0;
+ }
+ r2 = submod(W, w, n);
+ k = divmod(r2, r1, n);
+ if (powmod(g,k,p) != a) {
+ if (verbose > 2) printf("r1 = %lu r2 = %lu k = %lu\n", r1, r2, k);
+ if (verbose) printf("Incorrect DLP Rho solution: %lu\n", k);
+ return 0;
+ }
+ if (verbose) printf("DLP Rho solution found after %lu steps\n", i);
+ return k;
+ }
+ }
+ return 0;
+}
diff --git a/factor.h b/factor.h
index 7c4682c..03d324c 100644
--- a/factor.h
+++ b/factor.h
@@ -21,4 +21,7 @@ extern int squfof_factor(UV n, UV *factors, UV rounds);
extern UV* _divisor_list(UV n, UV *num_divisors);
+extern UV dlp_trial(UV a, UV g, UV p, UV maxrounds);
+extern UV dlp_prho(UV a, UV g, UV p, UV maxrounds);
+
#endif
diff --git a/t/19-moebius.t b/t/19-moebius.t
index dd56fe6..9ee1c56 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -332,6 +332,7 @@ my @znlogs = (
[ [3,3,8], 1],
[ [10,2,101], 25],
[ [2,55,101], 73], # 2 = 55^73 mod 101
+ [ [228,2,383], 110],
[ [3061666278, 499998, 3332205179], 22],
);
if ($usexs) {
diff --git a/util.c b/util.c
index 611fb94..9a3515c 100644
--- a/util.c
+++ b/util.c
@@ -1062,17 +1062,66 @@ UV znprimroot(UV n) {
return 0;
}
-/* Find smallest n where a = g^n mod p */
-/* This implementation is just a stupid placeholder. */
+/* Calculate 1/a mod p. From William Hart. */
+UV modinverse(UV a, UV p) {
+ IV u1 = 1, u3 = a;
+ IV v1 = 0, v3 = p;
+ IV t1 = 0, t3 = 0;
+ IV quot;
+ while (v3) {
+ quot = u3 - v3;
+ if (u3 < (v3<<2)) {
+ if (quot < v3) {
+ if (quot < 0) {
+ t1 = u1; u1 = v1; v1 = t1;
+ t3 = u3; u3 = v3; v3 = t3;
+ } else {
+ t1 = u1 - v1; u1 = v1; v1 = t1;
+ t3 = u3 - v3; u3 = v3; v3 = t3;
+ }
+ } else if (quot < (v3<<1)) {
+ t1 = u1 - (v1<<1); u1 = v1; v1 = t1;
+ t3 = u3 - (v3<<1); u3 = v3; v3 = t3;
+ } else {
+ t1 = u1 - v1*3; u1 = v1; v1 = t1;
+ t3 = u3 - v3*3; u3 = v3; v3 = t3;
+ }
+ } else {
+ quot = u3 / v3;
+ t1 = u1 - v1*quot; u1 = v1; v1 = t1;
+ t3 = u3 - v3*quot; u3 = v3; v3 = t3;
+ }
+ }
+ if (u1 < 0) u1 += p;
+ return u1;
+}
+
+UV divmod(UV a, UV b, UV n) { /* a / b mod n */
+ UV binv = modinverse(b, n);
+ if (binv == 0) return 0;
+ return mulmod(a, binv, n);
+}
+
+/* Find smallest n where a = g^n mod p
+ * This implementation is just a stupid placeholder.
+ * When prho or bsgs gets working well, lower the trial limit
+ */
+#define DLP_TRIAL_NUM 1000000
UV znlog(UV a, UV g, UV p) {
- UV t, n = 1;
- if (a == 0 || g == 0 || p < 2) return 0;
- for (n = 1; n < p; n++) {
- t = powmod(g, n, p);
- if (t == a)
- return n;
- }
- return 0;
+ UV k;
+ const int verbose = _XS_get_verbose();
+ if (a == 0 || g == 0 || p < 2)
+ return 0;
+ k = dlp_trial(a, g, p, DLP_TRIAL_NUM);
+ if (k != 0 || p <= DLP_TRIAL_NUM)
+ return k;
+ if (verbose) printf(" dlp trial failed. Trying prho\n");
+ k = dlp_prho(a, g, p, 1000000);
+ if (k != 0)
+ return k;
+ if (verbose) printf(" dlp prho failed. Back to trial\n");
+ k = dlp_trial(a, g, p, p);
+ return k;
}
long double _XS_chebyshev_theta(UV n)
diff --git a/util.h b/util.h
index 0afdff7..e214598 100644
--- a/util.h
+++ b/util.h
@@ -31,6 +31,9 @@ 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 modinverse(UV a, UV p); /* Returns 1/a mod p */
+extern UV divmod(UV a, UV b, UV n); /* Returns a/b mod n */
+
extern UV totient(UV n);
extern int moebius(UV n);
extern UV exp_mangoldt(UV n);
--
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
