isl_val_gcdext: revert to open-coded version of isl_int_gcdext
authorSven Verdoolaege <skimo@kotnet.org>
Wed, 28 May 2014 20:51:49 +0000 (28 22:51 +0200)
committerSven Verdoolaege <skimo@kotnet.org>
Thu, 29 May 2014 08:17:19 +0000 (29 10:17 +0200)
In 10a7f8a (isl_ast_build.c: use isl_int_gcdext instead of open-coded version,
Wed Apr 17 17:51:55 2013 +0200), the open-coded version was replaced
by a call to the isl_int_gcdext macro introduced in e191579
(add isl_int_gcdext, Wed Apr 17 17:51:28 2013 +0200).
In 0a83197 (try and avoid depending on mpz_gcdext internals,
Wed Jun 19 11:22:18 2013 +0200), we introduced a isl_gmp_gcdext
to try and mask differences in the outcomes of the underlying
mpz_gcdext over different versions of GMP.

However, there is no guarantee that we caught all the differences.
Furthermore, we will be introducing an alternative to GMP soon and
the differences across libraries are likely to be even greater than
those between different versions of the same library.

We therefore revert back to our (slightly modified) original implementation.

Signed-off-by: Sven Verdoolaege <skimo@kotnet.org>
configure.ac
isl_gmp.c
isl_int.h
isl_val.c

index 9b44f07..48e53ee 100644 (file)
@@ -80,30 +80,10 @@ need_get_memory_functions=false
 AC_CHECK_DECLS(mp_get_memory_functions,[],[
        need_get_memory_functions=true
 ],[#include <gmp.h>])
-AC_RUN_IFELSE([AC_LANG_PROGRAM([[#include <gmp.h>]], [[
-       mpz_t x,y,g,a,b;
-       mpz_init(x);
-       mpz_init(y);
-       mpz_init(g);
-       mpz_init(a);
-       mpz_init(b);
-       mpz_set_si(x, -1);
-       mpz_set_si(y, 9);
-       mpz_gcdext(g, a, b, x, y);
-       if (mpz_get_si(a) == -1 && mpz_get_si(b) == 0)
-               return 0;
-       else
-               return 1;
-]])], [need_normalized_gcdext=false], [need_normalized_gcdext=true],
-[need_normalized_gcdext=true])
 CPPFLAGS="$SAVE_CPPFLAGS"
 LDFLAGS="$SAVE_LDFLAGS"
 LIBS="$SAVE_LIBS"
 AM_CONDITIONAL(NEED_GET_MEMORY_FUNCTIONS, test x$need_get_memory_functions = xtrue)
-if test $need_normalized_gcdext = true; then
-AC_DEFINE([GMP_NORMALIZE_GCDEXT], [],
-       [result of mpz_gcdext needs to be normalized])
-fi
 AC_CHECK_DECLS(ffs,[],[],[#include <strings.h>])
 AC_CHECK_DECLS(__builtin_ffs,[],[],[])
 
index cf2ddee..d488fd2 100644 (file)
--- a/isl_gmp.c
+++ b/isl_gmp.c
@@ -22,36 +22,3 @@ uint32_t isl_gmp_hash(mpz_t v, uint32_t hash)
                isl_hash_byte(hash, *data);
        return hash;
 }
-
-/* This function tries to produce outputs that do not depend on
- * the version of GMP that is being used.
- *
- * In particular, when computing the extended gcd of -1 and 9,
- * some versions will produce
- *
- *     1 = -1 * -1 + 0 * 9
- *
- * while other versions will produce
- *
- *     1 = 8 * -1 + 1 * 9
- *
- * If configure detects that we are in the former case, then
- * mpz_gcdext will be called directly.  Otherwise, this function
- * is called and then we try to mimic the behavior of the other versions.
- */
-void isl_gmp_gcdext(mpz_t G, mpz_t S, mpz_t T, mpz_t A, mpz_t B)
-{
-       if (mpz_divisible_p(B, A)) {
-               mpz_set_si(S, mpz_sgn(A));
-               mpz_set_si(T, 0);
-               mpz_abs(G, A);
-               return;
-       }
-       if (mpz_divisible_p(A, B)) {
-               mpz_set_si(S, 0);
-               mpz_set_si(T, mpz_sgn(B));
-               mpz_abs(G, B);
-               return;
-       }
-       mpz_gcdext(G, S, T, A, B);
-}
index 520b486..71198f2 100644 (file)
--- a/isl_int.h
+++ b/isl_int.h
@@ -71,12 +71,6 @@ typedef void (*isl_int_print_gmp_free_t)(void *, size_t);
 #define isl_int_submul_ui(r,i,j)       mpz_submul_ui(r,i,j)
 
 #define isl_int_gcd(r,i,j)     mpz_gcd(r,i,j)
-#ifdef GMP_NORMALIZE_GCDEXT
-void isl_gmp_gcdext(mpz_t G, mpz_t S, mpz_t T, mpz_t A, mpz_t B);
-#define isl_int_gcdext(g,x,y,i,j)      isl_gmp_gcdext(g,x,y,i,j)
-#else
-#define isl_int_gcdext(g,x,y,i,j)      mpz_gcdext(g,x,y,i,j)
-#endif
 #define isl_int_lcm(r,i,j)     mpz_lcm(r,i,j)
 #define isl_int_divexact(r,i,j)        mpz_divexact(r,i,j)
 #define isl_int_divexact_ui(r,i,j)     mpz_divexact_ui(r,i,j)
index 8a9c0cb..38b4c97 100644 (file)
--- a/isl_val.c
+++ b/isl_val.c
@@ -942,6 +942,48 @@ error:
        return NULL;
 }
 
+/* Compute x, y and g such that g = gcd(a,b) and a*x+b*y = g.
+ */
+static void isl_int_gcdext(isl_int g, isl_int x, isl_int y,
+       isl_int a, isl_int b)
+{
+       isl_int d, tmp;
+       isl_int a_copy, b_copy;
+
+       isl_int_init(a_copy);
+       isl_int_init(b_copy);
+       isl_int_init(d);
+       isl_int_init(tmp);
+       isl_int_set(a_copy, a);
+       isl_int_set(b_copy, b);
+       isl_int_abs(g, a_copy);
+       isl_int_abs(d, b_copy);
+       isl_int_set_si(x, 1);
+       isl_int_set_si(y, 0);
+       while (isl_int_is_pos(d)) {
+               isl_int_fdiv_q(tmp, g, d);
+               isl_int_submul(x, tmp, y);
+               isl_int_submul(g, tmp, d);
+               isl_int_swap(g, d);
+               isl_int_swap(x, y);
+       }
+       if (isl_int_is_zero(a_copy))
+               isl_int_set_si(x, 0);
+       else if (isl_int_is_neg(a_copy))
+               isl_int_neg(x, x);
+       if (isl_int_is_zero(b_copy))
+               isl_int_set_si(y, 0);
+       else {
+               isl_int_mul(tmp, a_copy, x);
+               isl_int_sub(tmp, g, tmp);
+               isl_int_divexact(y, tmp, b_copy);
+       }
+       isl_int_clear(d);
+       isl_int_clear(tmp);
+       isl_int_clear(a_copy);
+       isl_int_clear(b_copy);
+}
+
 /* Given two integer values v1 and v2, return their greatest common divisor g,
  * as well as two integers x and y such that x * v1 + y * v2 = g.
  */