[libmath-prime-util-perl] 06/23: Move mulmod, powmod, etc. to separate file. A few changes.
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:45:54 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.14
in repository libmath-prime-util-perl.
commit 73b3a6e6141e8fbd9a881d9f0bb9a7f11532862f
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Nov 22 17:56:40 2012 -0800
Move mulmod, powmod, etc. to separate file. A few changes.
---
MANIFEST | 1 +
factor.c | 101 +--------------------------------------------------------
mulmod.h | 110 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 files changed, 112 insertions(+), 100 deletions(-)
diff --git a/MANIFEST b/MANIFEST
index 66357cb..4cf35c5 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -10,6 +10,7 @@ README
TODO
XS.xs
ptypes.h
+mulmod.h
aks.h
aks.c
cache.h
diff --git a/factor.c b/factor.c
index b89ca99..deedbae 100644
--- a/factor.c
+++ b/factor.c
@@ -7,6 +7,7 @@
#include "factor.h"
#include "util.h"
#include "sieve.h"
+#include "mulmod.h"
/*
* You need to remember to use UV for unsigned and IV for signed types that
@@ -69,106 +70,6 @@ int trial_factor(UV n, UV *factors, UV maxtrial)
}
-#if (BITS_PER_WORD == 32) && HAVE_STD_U64
-
- /* We have 64-bit available, but UV is 32-bit. Do the math in 64-bit.
- * Even if it is emulated, it should be as fast or faster than us doing it.
- */
- #define addmod(n,a,m) (UV)(((uint64_t)(n)+(uint64_t)(a)) % ((uint64_t)(m)))
- #define mulmod(a,b,m) (UV)(((uint64_t)(a)*(uint64_t)(b)) % ((uint64_t)(m)))
- #define sqrmod(n,m) (UV)(((uint64_t)(n)*(uint64_t)(n)) % ((uint64_t)(m)))
-
-#elif defined(__GNUC__) && defined(__x86_64__)
-
- /* GCC on a 64-bit Intel x86 */
- static UV _mulmod(UV a, UV b, UV c) {
- UV d; /* to hold the result of a*b mod c */
- /* calculates a*b mod c, stores result in d */
- asm ("mov %1, %%rax;" /* put a into rax */
- "mul %2;" /* mul a*b -> rdx:rax */
- "div %3;" /* (a*b)/c -> quot in rax remainder in rdx */
- "mov %%rdx, %0;" /* store result in d */
- :"=r"(d) /* output */
- :"r"(a), "r"(b), "r"(c) /* input */
- :"%rax", "%rdx" /* clobbered registers */
- );
- return d;
- }
- #define mulmod(a,b,m) _mulmod(a,b,m)
- #define sqrmod(n,m) _mulmod(n,n,m)
-
-#else
-
- /* UV is the largest integral type available (that we know of). */
-
- /* Do it by hand */
- static UV _mulmod(UV a, UV b, UV m) {
- UV r = 0;
- while (b > 0) {
- if (b & 1) {
- if (r == 0) {
- r = a;
- } else {
- r = m - r;
- r = (a >= r) ? a-r : m-r+a;
- }
- }
- a = (a > (m-a)) ? (a-m)+a : a+a;
- b >>= 1;
- }
- return r;
- }
-
- /* if n is smaller than this, you can multiply without overflow */
- #define HALF_WORD (UVCONST(1) << (BITS_PER_WORD/2))
- #define mulmod(a,b,m) (((a)|(b)) < HALF_WORD) ? ((a)*(b))%(m):_mulmod(a,b,m)
- #define sqrmod(n,m) ((n) < HALF_WORD) ? ((n)*(n))%(m):_mulmod(n,n,m)
-
-#endif
-
-#ifndef addmod
- #define addmod(n,a,m) ((((m)-(n)) > (a)) ? ((n)+(a)) : ((n)+(a)-(m)))
-#endif
-
-/* n^power mod m */
-#ifndef HALF_WORD
- static UV powmod(UV n, UV power, UV m) {
- UV t = 1;
- n %= m;
- while (power) {
- if (power & 1) t = mulmod(t, n, m);
- power >>= 1;
- if (power) n = sqrmod(n, m);
- }
- return t;
- }
-#else
- static UV powmod(UV n, UV power, UV m) {
- UV t = 1;
- n %= m;
- if (m < HALF_WORD) {
- while (power) {
- if (power & 1) t = (t*n)%m;
- power >>= 1;
- if (power) n = (n*n)%m;
- }
- } else {
- while (power) {
- if (power & 1) t = mulmod(t, n, m);
- power >>= 1;
- if (power) n = sqrmod(n,m);
- }
- }
- return t;
- }
-#endif
-
-/* n^power + a mod m */
-#define powaddmod(n, p, a, m) addmod(powmod(n,p,m),a,m)
-
-/* n^2 + a mod m */
-#define sqraddmod(n, a, m) addmod(sqrmod(n,m), a,m)
-
/* Return 0 if n is not a perfect square. Set sqrtn to int(sqrt(n)) if so.
*
* Some simple solutions:
diff --git a/mulmod.h b/mulmod.h
new file mode 100644
index 0000000..bce8eb7
--- /dev/null
+++ b/mulmod.h
@@ -0,0 +1,110 @@
+#ifndef MPU_MULMOD_H
+#define MPU_MULMOD_H
+
+#include "ptypes.h"
+
+#if defined(__GNUC__) || defined(_MSC_VER)
+ #define INLINE inline
+#else
+ #define INLINE
+#endif
+
+/* if n is smaller than this, you can multiply without overflow */
+#define HALF_WORD (UVCONST(1) << (BITS_PER_WORD/2))
+
+#if (BITS_PER_WORD == 32) && HAVE_STD_U64
+
+ /* We have 64-bit available, but UV is 32-bit. Do the math in 64-bit.
+ * Even if it is emulated, it should be as fast or faster than us doing it.
+ */
+ #define addmod(n,a,m) (UV)(((uint64_t)(n)+(uint64_t)(a)) % ((uint64_t)(m)))
+ #define mulmod(a,b,m) (UV)(((uint64_t)(a)*(uint64_t)(b)) % ((uint64_t)(m)))
+ #define sqrmod(n,m) (UV)(((uint64_t)(n)*(uint64_t)(n)) % ((uint64_t)(m)))
+
+#elif defined(__GNUC__) && defined(__x86_64__)
+
+ /* GCC on a 64-bit Intel x86 */
+ static INLINE UV _mulmod(UV a, UV b, UV c) {
+ UV d; /* to hold the result of a*b mod c */
+ /* calculates a*b mod c, stores result in d */
+ asm ("mov %1, %%rax;" /* put a into rax */
+ "mul %2;" /* mul a*b -> rdx:rax */
+ "div %3;" /* (a*b)/c -> quot in rax remainder in rdx */
+ "mov %%rdx, %0;" /* store result in d */
+ :"=r"(d) /* output */
+ :"r"(a), "r"(b), "r"(c) /* input */
+ :"%rax", "%rdx" /* clobbered registers */
+ );
+ return d;
+ }
+ #define mulmod(a,b,m) _mulmod(a,b,m)
+ #define sqrmod(n,m) _mulmod(n,n,m)
+
+#else
+
+ /* UV is the largest integral type available (that we know of). */
+
+ /* Do it by hand */
+ static INLINE UV _mulmod(UV a, UV b, UV m) {
+ UV r = 0;
+ a %= m; /* These are wasteful given that careful attention from the */
+ b %= m; /* caller should make them unnecessary. */
+ while (b > 0) {
+ if (b & 1) r = ((m-r) > a) ? r+a : r+a-m; /* r = (r + a) % m */
+ b >>= 1;
+ if (b) a = ((m-a) > a) ? a+a : a+a-m; /* a = (a + a) % m */
+ }
+ return r;
+ }
+
+ #define mulmod(a,b,m) (((a)|(b)) < HALF_WORD) ? ((a)*(b))%(m):_mulmod(a,b,m)
+ #define sqrmod(n,m) ((n) < HALF_WORD) ? ((n)*(n))%(m):_mulmod(n,n,m)
+
+#endif
+
+#ifndef addmod
+ #define addmod(n,a,m) ((((m)-(n)) > (a)) ? ((n)+(a)) : ((n)+(a)-(m)))
+#endif
+
+/* n^2 + a mod m */
+#define sqraddmod(n, a, m) addmod(sqrmod(n,m), a,m)
+/* i*j + a mod m */
+#define muladdmod(i, j, a, m) addmod(mulmod(i,j,m), a, m)
+
+/* n^power mod m */
+#ifndef HALF_WORD
+ static INLINE UV powmod(UV n, UV power, UV m) {
+ UV t = 1;
+ n %= m;
+ while (power) {
+ if (power & 1) t = mulmod(t, n, m);
+ power >>= 1;
+ if (power) n = sqrmod(n, m);
+ }
+ return t;
+ }
+#else
+ static INLINE UV powmod(UV n, UV power, UV m) {
+ UV t = 1;
+ n %= m;
+ if (m < HALF_WORD) {
+ while (power) {
+ if (power & 1) t = (t*n)%m;
+ power >>= 1;
+ if (power) n = (n*n)%m;
+ }
+ } else {
+ while (power) {
+ if (power & 1) t = mulmod(t, n, m);
+ power >>= 1;
+ if (power) n = sqrmod(n,m);
+ }
+ }
+ return t;
+ }
+#endif
+
+/* n^power + a mod m */
+#define powaddmod(n, p, a, m) addmod(powmod(n,p,m),a,m)
+
+#endif
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list