summaryrefslogtreecommitdiffstats
path: root/libtommath/bn_mp_gcd.c
diff options
context:
space:
mode:
authorKevin B Kenny <kennykb@acm.org>2005-01-19 22:41:26 (GMT)
committerKevin B Kenny <kennykb@acm.org>2005-01-19 22:41:26 (GMT)
commitef78ca64ce6ba6a8786f083318fe536f2bd52925 (patch)
tree47f8ad0d7291237c7f9af988c5e05275ed9286ee /libtommath/bn_mp_gcd.c
parentb23d942a1e86ddee18c2309afd7fa7e9afa79ef8 (diff)
downloadtcl-ef78ca64ce6ba6a8786f083318fe536f2bd52925.zip
tcl-ef78ca64ce6ba6a8786f083318fe536f2bd52925.tar.gz
tcl-ef78ca64ce6ba6a8786f083318fe536f2bd52925.tar.bz2
Import of libtommath 0.33
Diffstat (limited to 'libtommath/bn_mp_gcd.c')
-rw-r--r--libtommath/bn_mp_gcd.c109
1 files changed, 109 insertions, 0 deletions
diff --git a/libtommath/bn_mp_gcd.c b/libtommath/bn_mp_gcd.c
new file mode 100644
index 0000000..6265df1
--- /dev/null
+++ b/libtommath/bn_mp_gcd.c
@@ -0,0 +1,109 @@
+#include <tommath.h>
+#ifdef BN_MP_GCD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is a library that provides multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library was designed directly after the MPI library by
+ * Michael Fromberger but has been written from scratch with
+ * additional optimizations in place.
+ *
+ * The library is free for all purposes without any express
+ * guarantee it works.
+ *
+ * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
+ */
+
+/* Greatest Common Divisor using the binary method */
+int mp_gcd (mp_int * a, mp_int * b, mp_int * c)
+{
+ mp_int u, v;
+ int k, u_lsb, v_lsb, res;
+
+ /* either zero than gcd is the largest */
+ if (mp_iszero (a) == 1 && mp_iszero (b) == 0) {
+ return mp_abs (b, c);
+ }
+ if (mp_iszero (a) == 0 && mp_iszero (b) == 1) {
+ return mp_abs (a, c);
+ }
+
+ /* optimized. At this point if a == 0 then
+ * b must equal zero too
+ */
+ if (mp_iszero (a) == 1) {
+ mp_zero(c);
+ return MP_OKAY;
+ }
+
+ /* get copies of a and b we can modify */
+ if ((res = mp_init_copy (&u, a)) != MP_OKAY) {
+ return res;
+ }
+
+ if ((res = mp_init_copy (&v, b)) != MP_OKAY) {
+ goto LBL_U;
+ }
+
+ /* must be positive for the remainder of the algorithm */
+ u.sign = v.sign = MP_ZPOS;
+
+ /* B1. Find the common power of two for u and v */
+ u_lsb = mp_cnt_lsb(&u);
+ v_lsb = mp_cnt_lsb(&v);
+ k = MIN(u_lsb, v_lsb);
+
+ if (k > 0) {
+ /* divide the power of two out */
+ if ((res = mp_div_2d(&u, k, &u, NULL)) != MP_OKAY) {
+ goto LBL_V;
+ }
+
+ if ((res = mp_div_2d(&v, k, &v, NULL)) != MP_OKAY) {
+ goto LBL_V;
+ }
+ }
+
+ /* divide any remaining factors of two out */
+ if (u_lsb != k) {
+ if ((res = mp_div_2d(&u, u_lsb - k, &u, NULL)) != MP_OKAY) {
+ goto LBL_V;
+ }
+ }
+
+ if (v_lsb != k) {
+ if ((res = mp_div_2d(&v, v_lsb - k, &v, NULL)) != MP_OKAY) {
+ goto LBL_V;
+ }
+ }
+
+ while (mp_iszero(&v) == 0) {
+ /* make sure v is the largest */
+ if (mp_cmp_mag(&u, &v) == MP_GT) {
+ /* swap u and v to make sure v is >= u */
+ mp_exch(&u, &v);
+ }
+
+ /* subtract smallest from largest */
+ if ((res = s_mp_sub(&v, &u, &v)) != MP_OKAY) {
+ goto LBL_V;
+ }
+
+ /* Divide out all factors of two */
+ if ((res = mp_div_2d(&v, mp_cnt_lsb(&v), &v, NULL)) != MP_OKAY) {
+ goto LBL_V;
+ }
+ }
+
+ /* multiply by 2**k which we divided out at the beginning */
+ if ((res = mp_mul_2d (&u, k, c)) != MP_OKAY) {
+ goto LBL_V;
+ }
+ c->sign = MP_ZPOS;
+ res = MP_OKAY;
+LBL_V:mp_clear (&u);
+LBL_U:mp_clear (&v);
+ return res;
+}
+#endif