Initial version of the integer set library
authorSven Verdoolaege <skimo@kotnet.org>
Thu, 7 Aug 2008 18:45:58 +0000 (20:45 +0200)
committerSven Verdoolaege <skimo@kotnet.org>
Thu, 7 Aug 2008 23:18:28 +0000 (01:18 +0200)
This version is very incomplete and unoptimized.

31 files changed:
.gitignore [new file with mode: 0644]
.gitmodules [new file with mode: 0644]
Makefile.am [new file with mode: 0644]
autogen.sh [new file with mode: 0755]
configure.ac [new file with mode: 0644]
include/isl_blk.h [new file with mode: 0644]
include/isl_ctx.h.in [new file with mode: 0644]
include/isl_int.h [new file with mode: 0644]
include/isl_lp.h [new file with mode: 0644]
include/isl_lp_piplib.h [new file with mode: 0644]
include/isl_map.h [new file with mode: 0644]
include/isl_map_piplib.h [new file with mode: 0644]
include/isl_map_polylib.h [new file with mode: 0644]
include/isl_piplib.h [new file with mode: 0644]
include/isl_polylib.h [new file with mode: 0644]
include/isl_seq.h [new file with mode: 0644]
include/isl_set.h [new file with mode: 0644]
include/isl_set_polylib.h [new file with mode: 0644]
isl_blk.c [new file with mode: 0644]
isl_ctx.c [new file with mode: 0644]
isl_gmp.c [new file with mode: 0644]
isl_lp.c [new file with mode: 0644]
isl_lp_piplib.c [new file with mode: 0644]
isl_map.c [new file with mode: 0644]
isl_map_affine_hull.c [new file with mode: 0644]
isl_map_piplib.c [new file with mode: 0644]
isl_map_polylib.c [new file with mode: 0644]
isl_map_private.h [new file with mode: 0644]
isl_piplib.c [new file with mode: 0644]
isl_seq.c [new file with mode: 0644]
piplib [new submodule]

diff --git a/.gitignore b/.gitignore
new file mode 100644 (file)
index 0000000..1151f52
--- /dev/null
@@ -0,0 +1,11 @@
+Makefile.in
+aclocal.m4
+autom4te.cache
+config.guess
+config.h.in
+config.sub
+configure
+depcomp
+install-sh
+ltmain.sh
+missing
diff --git a/.gitmodules b/.gitmodules
new file mode 100644 (file)
index 0000000..d7912b0
--- /dev/null
@@ -0,0 +1,3 @@
+[submodule "piplib"]
+       path = piplib
+       url = git://repo.or.cz/piplib.git
diff --git a/Makefile.am b/Makefile.am
new file mode 100644 (file)
index 0000000..f46706a
--- /dev/null
@@ -0,0 +1,77 @@
+if BUNDLED_PIPLIB
+MAYBE_PIPLIB = piplib
+endif
+
+SUBDIRS = $(MAYBE_PIPLIB) .
+
+lib_LTLIBRARIES = libisl.la
+
+if HAVE_POLYLIB
+ISL_POLYLIB = \
+       isl_map_polylib.c \
+       isl_map_polylib.h \
+       isl_polylib.h
+endif
+
+if HAVE_PIPLIB
+ISL_PIPLIB = \
+       isl_lp_piplib.h \
+       isl_lp_piplib.c \
+       isl_map_piplib.c \
+       isl_map_piplib.h \
+       isl_piplib.h \
+       isl_piplib.c
+endif
+if BUNDLED_PIPLIB
+PIPLIB_LA = $(top_builddir)/piplib/libpiplibMP.la
+endif
+
+libisl_la_SOURCES = \
+       $(ISL_PIPLIB) \
+       $(ISL_POLYLIB) \
+       isl_blk.c \
+       isl_blk.h \
+       isl_ctx.c \
+       isl_ctx.h \
+       isl_gmp.c \
+       isl_int.h \
+       isl_lp.h \
+       isl_lp.c \
+       isl_map.c \
+       isl_map.h \
+       isl_map_affine_hull.c \
+       isl_map_private.h \
+       isl_set.h \
+       isl_seq.c \
+       isl_seq.h
+EXTRA_libisl_la_SOURCES = \
+       isl_lp_piplib.h \
+       isl_lp_piplib.c \
+       isl_map_piplib.c \
+       isl_map_piplib.h \
+       isl_map_polylib.c \
+       isl_map_polylib.h \
+       isl_piplib.h \
+       isl_piplib.c \
+       isl_polylib.h
+libisl_la_LIBADD = $(PIPLIB_LA) @PIPLIB_LIBS@ @POLYLIB_LIBS@ -lgmp
+libisl_la_LDFLAGS = -release @VERSION@ @PIPLIB_LDFLAGS@ @POLYLIB_LDFLAGS@ \
+       @GMP_LDFLAGS@
+libisl_la_CPPFLAGS = -I$(srcdir)/include -Iinclude/ \
+       @PIPLIB_CPPFLAGS@ @POLYLIB_CPPFLAGS@ \
+       @GMP_CPPFLAGS@
+
+pkginclude_HEADERS = \
+       include/isl_blk.h \
+       include/isl_ctx.h \
+       include/isl_int.h \
+       include/isl_lp.h \
+       include/isl_lp_piplib.h \
+       include/isl_map.h \
+       include/isl_map_piplib.h \
+       include/isl_map_polylib.h \
+       include/isl_piplib.h \
+       include/isl_polylib.h \
+       include/isl_seq.h \
+       include/isl_set.h \
+       include/isl_set_polylib.h
diff --git a/autogen.sh b/autogen.sh
new file mode 100755 (executable)
index 0000000..505354d
--- /dev/null
@@ -0,0 +1,9 @@
+#!/bin/sh
+libtoolize -c
+aclocal
+autoheader
+automake -a -c --foreign
+autoconf
+if test -f piplib/autogen.sh; then
+       (cd piplib; ./autogen.sh)
+fi
diff --git a/configure.ac b/configure.ac
new file mode 100644 (file)
index 0000000..9a246fe
--- /dev/null
@@ -0,0 +1,186 @@
+AC_INIT
+AM_INIT_AUTOMAKE(isl, 0.00)
+
+AC_PROG_CC
+
+AC_PROG_LIBTOOL
+
+AC_ARG_WITH(gmp_prefix,
+       [AS_HELP_STRING([--with-gmp-prefix=DIR],
+                       [Location of GMP installation])])
+
+AC_SUBST(GMP_CPPFLAGS)
+AC_SUBST(GMP_LDFLAGS)
+if test "x$with_gmp_prefix" != "x"; then
+       isl_configure_args="$isl_configure_args --with-gmp=$with_gmp_prefix"
+       GMP_CPPFLAGS="-I$with_gmp_prefix/include"
+       GMP_LDFLAGS="-L$with_gmp_prefix/lib"
+fi
+
+AC_DEFUN([ISL_SUBMODULE],[
+       AC_ARG_WITH($1_prefix,
+               [AS_HELP_STRING([--with-$1-prefix=DIR],
+                               [Location of system $1])])
+       AC_ARG_WITH($1_builddir,
+                   [AS_HELP_STRING([--with-$1-builddir=DIR],
+                                   [Location of $1 builddir])])
+       if test "x$with_$1_prefix" != "x"; then
+               if test "x$with_$1" != "x" -a "x$with_$1" != "xsystem"; then
+                       AC_MSG_ERROR([Setting $with_$1_prefix implies use of system $1])
+               fi
+               with_$1="system"
+       fi
+       if test "x$with_$1_builddir" != "x"; then
+               if test "x$with_$1" != "x" -a "x$with_$1" != "xbuild"; then
+                       AC_MSG_ERROR([Setting $with_$1_builddir implies use of build $1])
+               fi
+               with_$1="build"
+               $1_srcdir=`echo @abs_srcdir@ | $with_$1_builddir/config.status --file=-`
+               AC_MSG_NOTICE($1 sources in $$1_srcdir)
+       fi
+       case "$with_$1" in
+       no|bundled|system|build)
+               ;;
+       *)
+               if test -d $srcdir/.git -a \
+                       -d $srcdir/$1 -a \
+                       ! -d $srcdir/$1/.git; then
+                       AC_MSG_WARN(
+[git repo detected, but submodule $1 not initialized])
+                       AC_MSG_WARN([You may want to run])
+                       AC_MSG_WARN([   git submodule init])
+                       AC_MSG_WARN([   git submodule update])
+                       AC_MSG_WARN([   sh autogen.sh])
+               fi
+               if test -f $srcdir/$1/configure; then
+                       with_$1="bundled"
+               fi
+               ;;
+       esac
+])
+
+AC_ARG_WITH(polylib,
+       [AS_HELP_STRING([--with-polylib=build|system|no],
+                       [Which PolyLib to use])])
+ISL_SUBMODULE(polylib)
+if test "x$with_polylib" = "x"; then
+       with_polylib="no"
+fi
+AC_MSG_CHECKING([which polylib to use])
+AC_MSG_RESULT($with_polylib)
+
+have_polylib=false
+AC_SUBST(POLYLIB_CPPFLAGS)
+AC_SUBST(POLYLIB_LDFLAGS)
+AC_SUBST(POLYLIB_LIBS)
+case "$with_polylib" in
+       build)
+               polylibs=`echo @polylibs@ | $with_polylib_builddir/config.status --file=-`
+               AC_MSG_NOTICE(Configured polylibs: $polylibs)
+               isl_cv_polylib=missing
+               for bits in $polylibs; do
+                       if test "$bits" = "libpolylibgmp.la"; then
+                               isl_cv_polylib=ok
+                       fi
+               done
+               if test "$isl_cv_polylib" = "missing"; then
+                       AC_MSG_ERROR(no gmp polylib configured)
+               fi
+               POLYLIB_CPPFLAGS="-I$with_polylib_builddir/include -I$polylib_srcdir/include"
+       ;;
+       system)
+               POLYLIB_LIBS="-lpolylibgmp"
+               if test "x$with_polylib_prefix" != "x"; then
+                       POLYLIB_CPPFLAGS="-I$with_polylib_prefix/include"
+                       POLYLIB_LDFLAGS="-L$with_polylib_prefix/lib"
+               fi
+               SAVE_CPPFLAGS="$CPPFLAGS"
+               SAVE_LDFLAGS="$LDFLAGS"
+               CPPFLAGS="$POLYLIB_CPPFLAGS $CPPFLAGS"
+               LDFLAGS="$POLYLIB_LDFLAGS $LDFLAGS"
+               AC_CHECK_LIB(polylibgmp, PolyhedronTSort,[ true ],[
+                       AC_MSG_ERROR(Need polylib)
+               ])
+               CPPFLAGS="$SAVE_CPPFLAGS"
+               LDFLAGS="$SAVE_LDFLAGS"
+       ;;
+       no)
+       ;;
+       *)
+               AC_MSG_ERROR(unsupported)
+       ;;
+esac
+if test "$with_polylib" != "no"; then
+       AC_DEFINE(ISL_POLYLIB,,polylib is available)
+       have_polylib=true
+fi
+AM_CONDITIONAL(HAVE_POLYLIB, test x$have_polylib = xtrue)
+
+AC_ARG_WITH(piplib,
+       [AS_HELP_STRING([--with-piplib=build|bundled|system|no],
+                       [Which piplib to use])])
+ISL_SUBMODULE(piplib)
+if test "x$with_piplib" = "x"; then
+       with_piplib="system"
+fi
+AC_MSG_CHECKING([which piplib to use])
+AC_MSG_RESULT($with_piplib)
+
+have_piplib=false
+AC_SUBST(PIPLIB_CPPFLAGS)
+AC_SUBST(PIPLIB_LDFLAGS)
+AC_SUBST(PIPLIB_LIBS)
+case "$with_piplib" in
+       bundled)
+               PIPLIB_CPPFLAGS="-I$srcdir/piplib/include"
+               isl_configure_args="$isl_configure_args --enable-mp-version"
+       ;;
+       build)
+               PIPLIB_CPPFLAGS="-I$piplib_srcdir/include"
+               PIPLIB_LIBS="$with_piplib_builddir/libpiplibMP.la"
+       ;;
+       system)
+               PIPLIB_LIBS="-lpiplibMP"
+               if test "x$with_piplib_prefix" != "x"; then
+                       PIPLIB_CPPFLAGS="-I$with_piplib_prefix/include"
+                       PIPLIB_LDFLAGS="-L$with_piplib_prefix/lib"
+               fi
+               SAVE_CPPFLAGS="$CPPFLAGS"
+               SAVE_LDFLAGS="$LDFLAGS"
+               CPPFLAGS="$PIPLIB_CPPFLAGS $CPPFLAGS"
+               LDFLAGS="$PIPLIB_LDFLAGS $LDFLAGS"
+               AC_CHECK_LIB(piplibMP, pip_solve,[
+                       AC_CHECK_MEMBER(PipOptions.Urs_parms, [], [
+                               AC_MSG_ERROR([Piplib too old; please install version 1.3.6 or newer])
+                       ],[#include <piplib/piplibMP.h>])
+               ],[
+                       AC_MSG_ERROR([Piplib not found])
+               ])
+               CPPFLAGS="$SAVE_CPPFLAGS"
+               LDFLAGS="$SAVE_LDFLAGS"
+       ;;
+       no)
+       ;;
+       *)
+               AC_MSG_ERROR(unsupported)
+       ;;
+esac
+if test "$with_piplib" != "no"; then
+       AC_DEFINE(ISL_PIPLIB,,piplib is available)
+       have_piplib=true
+fi
+AM_CONDITIONAL(HAVE_PIPLIB, test x$have_piplib = xtrue)
+AM_CONDITIONAL(BUNDLED_PIPLIB, test $with_piplib = bundled)
+
+AC_CONFIG_HEADERS(config.h)
+AC_CONFIG_HEADERS(include/isl_ctx.h)
+AC_CONFIG_FILES(Makefile)
+if test $with_piplib = bundled; then
+       AC_CONFIG_SUBDIRS(piplib)
+fi
+AC_CONFIG_COMMANDS_POST([
+       dnl pass on arguments to subdir configures, but don't
+       ac_configure_args="$ac_configure_args $isl_configure_args"
+       ac_configure_args="$ac_configure_args $isl_configure_args"
+])
+AC_OUTPUT
diff --git a/include/isl_blk.h b/include/isl_blk.h
new file mode 100644 (file)
index 0000000..f46a0f1
--- /dev/null
@@ -0,0 +1,25 @@
+#ifndef ISL_BLK_H
+#define ISL_BLK_H
+
+#include <isl_int.h>
+#include <isl_ctx.h>
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+struct isl_blk {
+       size_t size;
+       isl_int *data;
+};
+
+struct isl_blk isl_blk_alloc(struct isl_ctx *ctx, size_t n);
+struct isl_blk isl_blk_extend(struct isl_ctx *ctx, struct isl_blk block,
+                               size_t new_n);
+void isl_blk_free(struct isl_ctx *ctx, struct isl_blk block);
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
diff --git a/include/isl_ctx.h.in b/include/isl_ctx.h.in
new file mode 100644 (file)
index 0000000..4b58276
--- /dev/null
@@ -0,0 +1,80 @@
+#ifndef ISL_CTX_H
+#define ISL_CTX_H
+
+#include <assert.h>
+#include <stdlib.h>
+
+#include "isl_int.h"
+
+#undef ISL_POLYLIB
+#undef ISL_PIPLIB
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+/* Nearly all isa functions require a struct isl_ctx allocated using
+ * isl_ctx_alloc.  This ctx contains (or will contain) options that
+ * control the behavior of the library and some caches.
+ *
+ * An object allocated within a given ctx should never be used inside
+ * another ctx.  Functions for moving objects from one ctx to another
+ * will be added as the need arises.
+ *
+ * A given context should only be used inside a single thread.
+ * A global context for synchronization between different threads
+ * as well as functions for moving a context to a different thread
+ * will be added as the need arises.
+ *
+ * If anything goes wrong (out of memory, failed assertion), then
+ * the library will currently simply abort.  This will be made
+ * configurable in the future.
+ * Users of the library should expect functions that return
+ * a pointer to a structure, to return NULL, indicating failure.
+ * Any function accepting a pointer to a structure will treat
+ * a NULL argument as a failure, resulting in the function freeing
+ * the remaining structures (if any) and returning NULL itself
+ * (in case of pointer return type).
+ * The only exception is the isl_ctx argument, which shoud never be NULL.
+ */
+struct isl_ctx {
+       isl_int         one;
+#ifdef ISL_POLYLIB
+       unsigned        MaxRays;
+#endif
+};
+
+/* Some helper macros */
+
+#define FL_INIT(l, f)   (l) = (f)               /* Specific flags location. */
+#define FL_SET(l, f)    ((l) |= (f))
+#define FL_CLR(l, f)    ((l) &= ~(f))
+#define FL_ISSET(l, f)  ((l) & (f))
+
+#define F_INIT(p, f)    FL_INIT((p)->flags, f)  /* Structure element flags. */
+#define F_SET(p, f)     FL_SET((p)->flags, f)
+#define F_CLR(p, f)     FL_CLR((p)->flags, f)
+#define F_ISSET(p, f)   FL_ISSET((p)->flags, f)
+
+#define isl_alloc(ctx,type,size)       (type *)malloc(size)
+#define isl_realloc(ctx,ptr,type,size) (type *)realloc(ptr,size)
+#define isl_alloc_type(ctx,type)       isl_alloc(ctx,type,sizeof(type))
+#define isl_realloc_type(ctx,ptr,type) isl_realloc(ctx,ptr,type,sizeof(type))
+#define isl_alloc_array(ctx,type,n)    isl_alloc(ctx,type,(n)*sizeof(type))
+#define isl_realloc_array(ctx,ptr,type,n) \
+                                   isl_realloc(ctx,ptr,type,(n)*sizeof(type))
+
+#define isl_assert(ctx,test,code)      assert(test);
+
+#define isl_min(a,b)                   ((a < b) ? (a) : (b))
+
+/* struct isl_ctx functions */
+
+struct isl_ctx *isl_ctx_alloc();
+void isl_ctx_free(struct isl_ctx *ctx);
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
diff --git a/include/isl_int.h b/include/isl_int.h
new file mode 100644 (file)
index 0000000..6e39a5d
--- /dev/null
@@ -0,0 +1,71 @@
+#ifndef ISL_INT_H
+#define ISL_INT_H
+
+#include <sys/types.h>
+#include <string.h>
+#include <gmp.h>
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+/* isl_int is the basic integer type.  It currently always corresponds
+ * to a gmp mpz_t, but in the future, different types such as long long
+ * or cln::cl_I will be supported.
+ */
+typedef mpz_t  isl_int;
+
+#define isl_int_init(i)                mpz_init(i)
+#define isl_int_clear(i)       mpz_clear(i)
+
+#define isl_int_set(r,i)       mpz_set(r,i)
+#define isl_int_set_si(r,i)    mpz_set_si(r,i)
+#define isl_int_abs(r,i)       mpz_abs(r,i)
+#define isl_int_neg(r,i)       mpz_neg(r,i)
+#define isl_int_swap(i,j)      mpz_swap(i,j)
+#define isl_int_sub_ui(r,i,j)  mpz_sub_ui(r,i,j)
+
+#define isl_int_add(r,i,j)     mpz_add(r,i,j)
+#define isl_int_sub(r,i,j)     mpz_sub(r,i,j)
+#define isl_int_mul(r,i,j)     mpz_mul(r,i,j)
+#define isl_int_addmul(r,i,j)  mpz_addmul(r,i,j)
+
+#define isl_int_gcd(r,i,j)     mpz_gcd(r,i,j)
+#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_cdiv_q(r,i,j)  mpz_cdiv_q(r,i,j)
+
+#define isl_int_print(out,i)                                           \
+       do {                                                            \
+               char *s;                                                \
+               void (*gmp_free) (void *, size_t);                      \
+               s = mpz_get_str(0, 10, i);                              \
+               fprintf(out, "%s", s);                                  \
+               mp_get_memory_functions(NULL, NULL, &gmp_free);         \
+               (*gmp_free)(s, strlen(s)+1);                            \
+       } while (0)
+
+#define isl_int_sgn(i)         mpz_sgn(i)
+#define isl_int_cmp_si(i,si)   mpz_cmp_si(i,si)
+#define isl_int_eq(i,j)                (mpz_cmp(i,j) == 0)
+#define isl_int_ne(i,j)                (mpz_cmp(i,j) != 0)
+#define isl_int_lt(i,j)                (mpz_cmp(i,j) < 0)
+#define isl_int_abs_lt(i,j)    (mpz_cmpabs(i,j) < 0)
+
+
+#define isl_int_is_zero(i)     (isl_int_sgn(i) == 0)
+#define isl_int_is_one(i)      (isl_int_cmp_si(i,1) == 0)
+#define isl_int_is_negone(i)   (isl_int_cmp_si(i,-1) == 0)
+#define isl_int_is_pos(i)      (isl_int_sgn(i) > 0)
+#define isl_int_is_neg(i)      (isl_int_sgn(i) < 0)
+#define isl_int_is_nonpos(i)   (isl_int_sgn(i) <= 0)
+#define isl_int_is_nonneg(i)   (isl_int_sgn(i) >= 0)
+#define isl_int_is_divisible_by(i,j)   mpz_divisible_p(i,j)
+
+#define isl_int_hash(v,h)      isl_gmp_hash(v,h)
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
diff --git a/include/isl_lp.h b/include/isl_lp.h
new file mode 100644 (file)
index 0000000..5179cac
--- /dev/null
@@ -0,0 +1,24 @@
+#ifndef ISL_LP_H
+#define ISL_LP_H
+
+#include <isl_map.h>
+
+enum isl_lp_result {
+       isl_lp_error = -1,
+       isl_lp_ok = 0,
+       isl_lp_unbounded,
+       isl_lp_empty
+};
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+enum isl_lp_result isl_solve_lp(struct isl_basic_map *bmap, int maximize,
+                                     isl_int *f, isl_int denom, isl_int *opt);
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
diff --git a/include/isl_lp_piplib.h b/include/isl_lp_piplib.h
new file mode 100644 (file)
index 0000000..4374d3b
--- /dev/null
@@ -0,0 +1,18 @@
+#ifndef ISL_LP_PIPLIB_H
+#define ISL_LP_PIPLIB_H
+
+#include <isl_lp.h>
+#include <isl_piplib.h>
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+enum isl_lp_result isl_pip_solve_lp(struct isl_basic_map *bmap, int maximize,
+                                     isl_int *f, isl_int denom, isl_int *opt);
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
diff --git a/include/isl_map.h b/include/isl_map.h
new file mode 100644 (file)
index 0000000..f707cd7
--- /dev/null
@@ -0,0 +1,219 @@
+#ifndef ISL_MAP_H
+#define ISL_MAP_H
+
+#include <stdio.h>
+
+#include <isl_int.h>
+#include <isl_ctx.h>
+#include <isl_blk.h>
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+/* General notes:
+ *
+ * All structures are reference counted to allow reuse without duplication.
+ * A *_copy operation will increase the reference count, while a *_free
+ * operation will decrease the reference count and only actually release
+ * the structures when the reference count drops to zero.
+ *
+ * Functions that return an isa structure will in general _destroy_
+ * all argument isa structures (the obvious execption begin the _copy
+ * functions).  A pointer passed to such a function may therefore
+ * never be used after the function call.  If you want to keep a
+ * reference to the old structure(s), use the appropriate _copy function.
+ */
+
+/* A "basic map" is a relation between two sets of variables,
+ * called the "in" and "out" variables.
+ *
+ * It is implemented as a set with two extra fields:
+ * n_in is the number of in variables
+ * n_out is the number of out variables
+ * n_in + n_out should be equal to set.dim
+ */
+struct isl_basic_map {
+       int ref;
+#define ISL_PRIMITIVE_MAP_FINAL                (1 << 0)
+       unsigned flags;
+
+       unsigned nparam;
+       unsigned n_in;
+       unsigned n_out;
+       unsigned extra;
+
+       unsigned n_eq;
+       unsigned n_ineq;
+
+       size_t c_size;
+       isl_int **eq;
+       isl_int **ineq;
+
+       unsigned n_div;
+
+       isl_int **div;
+
+       struct isl_blk block;
+       struct isl_blk block2;
+};
+struct isl_basic_set;
+
+/* A "map" is a (disjoint) union of basic maps.
+ *
+ * Currently, the isl_set structure is identical to the isl_map structure
+ * and the library depends on this correspondence internally.
+ * However, users should not depend on this correspondence.
+ */
+struct isl_map {
+       int ref;
+#define ISL_MAP_DISJOINT               (1 << 0)
+       unsigned flags;
+
+       unsigned nparam;
+       unsigned n_in;
+       unsigned n_out;
+
+       int n;
+
+       size_t size;
+       struct isl_basic_map *p[0];
+};
+struct isl_set;
+
+struct isl_basic_map *isl_basic_map_alloc(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out, unsigned extra,
+               unsigned n_eq, unsigned n_ineq);
+struct isl_basic_map *isl_basic_map_identity(struct isl_ctx *ctx,
+               unsigned nparam, unsigned dim);
+struct isl_basic_map *isl_basic_map_finalize(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+void isl_basic_map_free(struct isl_ctx *ctx, struct isl_basic_map *bmap);
+struct isl_basic_map *isl_basic_map_copy(struct isl_ctx *ctx,
+                                       struct isl_basic_map *bmap);
+struct isl_basic_map *isl_basic_map_extend(struct isl_ctx *ctx,
+               struct isl_basic_map *base,
+               unsigned nparam, unsigned n_in, unsigned n_out, unsigned extra,
+               unsigned n_eq, unsigned n_ineq);
+struct isl_basic_map *isl_basic_map_equal(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out, unsigned n_equal);
+struct isl_basic_map *isl_basic_map_less_at(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out, unsigned pos);
+struct isl_basic_map *isl_basic_map_more_at(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out, unsigned pos);
+struct isl_basic_map *isl_basic_map_empty(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out);
+
+struct isl_basic_map *isl_basic_map_intersect_domain(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap,
+               struct isl_basic_set *bset);
+struct isl_basic_map *isl_basic_map_intersect(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2);
+struct isl_map *isl_basic_map_union(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2);
+struct isl_basic_map *isl_basic_map_apply_domain(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2);
+struct isl_basic_map *isl_basic_map_apply_range(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2);
+struct isl_basic_map *isl_basic_map_reverse(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+struct isl_basic_set *isl_basic_map_domain(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+struct isl_basic_set *isl_basic_map_range(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+struct isl_basic_map *isl_basic_map_from_basic_set(
+               struct isl_ctx *ctx, struct isl_basic_set *bset,
+               unsigned n_in, unsigned n_out);
+struct isl_basic_set *isl_basic_set_from_basic_map(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap);
+struct isl_basic_map *isl_basic_map_simplify(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap);
+
+struct isl_map *isl_basic_map_lexmax(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, struct isl_basic_set *dom,
+               struct isl_set **empty);
+struct isl_map *isl_basic_map_lexmin(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, struct isl_basic_set *dom,
+               struct isl_set **empty);
+
+void isl_basic_map_dump(struct isl_ctx *ctx, struct isl_basic_map *bmap,
+                               FILE *out, int indent);
+
+int isl_basic_map_is_empty(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+int isl_basic_map_is_subset(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2);
+int isl_basic_map_is_strict_subset(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2);
+
+struct isl_map *isl_map_alloc(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out, int n,
+               unsigned flags);
+struct isl_map *isl_map_empty(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out);
+struct isl_map *isl_map_dup(struct isl_ctx *ctx, struct isl_map *map);
+struct isl_map *isl_map_add(struct isl_ctx *ctx, struct isl_map *map,
+                               struct isl_basic_map *bmap);
+struct isl_map *isl_map_identity(struct isl_ctx *ctx,
+               unsigned nparam, unsigned dim);
+struct isl_map *isl_map_finalize(struct isl_ctx *ctx, struct isl_map *map);
+void isl_map_free(struct isl_ctx *ctx, struct isl_map *map);
+struct isl_map *isl_map_copy(struct isl_ctx *ctx, struct isl_map *map);
+struct isl_map *isl_map_extend(struct isl_ctx *ctx, struct isl_map *base,
+               unsigned nparam, unsigned n_in, unsigned n_out);
+struct isl_map *isl_map_reverse(struct isl_ctx *ctx, struct isl_map *map);
+struct isl_map *isl_map_union(struct isl_ctx *ctx,
+                       struct isl_map *map1, struct isl_map *map2);
+struct isl_map *isl_map_union_disjoint(struct isl_ctx *ctx,
+                       struct isl_map *map1, struct isl_map *map2);
+struct isl_map *isl_map_intersect_domain(
+               struct isl_ctx *ctx, struct isl_map *map,
+               struct isl_set *set);
+struct isl_map *isl_map_intersect_range(
+               struct isl_ctx *ctx, struct isl_map *map,
+               struct isl_set *set);
+struct isl_map *isl_map_apply_domain(
+               struct isl_ctx *ctx, struct isl_map *map1,
+               struct isl_map *map2);
+struct isl_map *isl_map_apply_range(
+               struct isl_ctx *ctx, struct isl_map *map1,
+               struct isl_map *map2);
+struct isl_map *isl_map_intersect(struct isl_ctx *ctx, struct isl_map *map1,
+               struct isl_map *map2);
+struct isl_map *isl_map_subtract(struct isl_ctx *ctx, struct isl_map *map1,
+               struct isl_map *map2);
+struct isl_map *isl_map_fix_input_si(struct isl_ctx *ctx, struct isl_map *map,
+               unsigned input, int value);
+struct isl_basic_set *isl_basic_map_deltas(struct isl_ctx *ctx,
+                               struct isl_basic_map *bmap);
+struct isl_set *isl_map_range(struct isl_ctx *ctx, struct isl_map *map);
+struct isl_basic_map *isl_map_affine_hull(struct isl_ctx *ctx,
+               struct isl_map *map);
+
+struct isl_set *isl_map_domain(struct isl_ctx *ctx, struct isl_map *bmap);
+struct isl_map *isl_map_from_basic_map(struct isl_ctx *ctx,
+                               struct isl_basic_map *bmap);
+struct isl_map *isl_map_from_set(
+               struct isl_ctx *ctx, struct isl_set *set,
+               unsigned n_in, unsigned n_out);
+
+int isl_map_is_empty(struct isl_ctx *ctx, struct isl_map *map);
+int isl_map_is_subset(struct isl_ctx *ctx, struct isl_map *map1,
+               struct isl_map *map2);
+int isl_map_is_equal(struct isl_ctx *ctx,
+               struct isl_map *map1, struct isl_map *map2);
+
+void isl_map_dump(struct isl_ctx *ctx, struct isl_map *map, FILE *out,
+                 int indent);
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
diff --git a/include/isl_map_piplib.h b/include/isl_map_piplib.h
new file mode 100644 (file)
index 0000000..fb12c01
--- /dev/null
@@ -0,0 +1,31 @@
+#ifndef ISL_MAP_PIPLIB_H
+#define ISL_MAP_PIPLIB_H
+
+#include <isl_map.h>
+#include <isl_piplib.h>
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+PipMatrix *isl_basic_map_to_pip(struct isl_basic_map *bmap, unsigned pip_param,
+                        unsigned extra_front, unsigned extra_back);
+
+struct isl_map *isl_map_from_quast(struct isl_ctx *ctx, PipQuast *q,
+               unsigned keep,
+               struct isl_basic_set *context,
+               struct isl_set **rest);
+struct isl_map *pip_isl_basic_map_lexmax(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, struct isl_basic_set *dom,
+               struct isl_set **empty);
+struct isl_map *pip_isl_basic_map_lexmin(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, struct isl_basic_set *dom,
+               struct isl_set **empty);
+struct isl_map *pip_isl_basic_map_compute_divs(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
diff --git a/include/isl_map_polylib.h b/include/isl_map_polylib.h
new file mode 100644 (file)
index 0000000..411d52e
--- /dev/null
@@ -0,0 +1,25 @@
+#ifndef ISL_MAP_POLYLIB_H
+#define ISL_MAP_POLYLIB_H
+
+#include <isl_map.h>
+#include <isl_polylib.h>
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+struct isl_basic_map *isl_basic_map_new_from_polylib(
+                       struct isl_ctx *ctx, Polyhedron *P,
+                       unsigned nparam, unsigned in, unsigned out);
+struct isl_map *isl_map_new_from_polylib(struct isl_ctx *ctx,
+                       Polyhedron *D,
+                       unsigned nparam, unsigned in, unsigned out);
+Polyhedron *isl_basic_map_to_polylib(struct isl_ctx *ctx,
+                       struct isl_basic_map *bmap);
+Polyhedron *isl_map_to_polylib(struct isl_ctx *ctx, struct isl_map *map);
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
diff --git a/include/isl_piplib.h b/include/isl_piplib.h
new file mode 100644 (file)
index 0000000..f075b57
--- /dev/null
@@ -0,0 +1,14 @@
+#ifndef ISL_PIPLIB_H
+#define ISL_PIPLIB_H
+
+#include <isl_ctx.h>
+#include <isl_int.h>
+#ifndef ISL_PIPLIB
+#error "no piplib"
+#endif
+
+#include <piplib/piplibMP.h>
+
+void isl_seq_cpy_to_pip(Entier *dst, isl_int *src, unsigned len);
+
+#endif
diff --git a/include/isl_polylib.h b/include/isl_polylib.h
new file mode 100644 (file)
index 0000000..b2c500d
--- /dev/null
@@ -0,0 +1,11 @@
+#ifndef ISL_POLYLIB_H
+#define ISL_POLYLIB_H
+
+#include <isl_ctx.h>
+#ifndef ISL_POLYLIB
+#error "no polylib"
+#endif
+
+#include <polylib/polylibgmp.h>
+
+#endif
diff --git a/include/isl_seq.h b/include/isl_seq.h
new file mode 100644 (file)
index 0000000..98864b7
--- /dev/null
@@ -0,0 +1,23 @@
+#ifndef ISL_SEQ_H
+#define ISL_SEQ_H
+
+#include <sys/types.h>
+#include <isl_int.h>
+
+/* Some common operations on sequences of isl_int's */
+
+void isl_seq_clr(isl_int *p, unsigned len);
+void isl_seq_neg(isl_int *dat, isl_int *src, unsigned len);
+void isl_seq_cpy(isl_int *dst, isl_int *src, unsigned len);
+void isl_seq_scale(isl_int *dst, isl_int *src, isl_int f, unsigned len);
+void isl_seq_scale_down(isl_int *dst, isl_int *src, isl_int f, unsigned len);
+void isl_seq_elim(isl_int *dst, isl_int *src, unsigned pos, unsigned len,
+                 isl_int *m);
+void isl_seq_gcd(isl_int *p, unsigned len, isl_int *gcd);
+int isl_seq_first_non_zero(isl_int *p, unsigned len);
+int isl_seq_abs_min_non_zero(isl_int *p, unsigned len);
+int isl_seq_eq(isl_int *p1, isl_int *p2, unsigned len);
+
+u_int32_t isl_seq_hash(isl_int *p, unsigned len, unsigned bits);
+
+#endif
diff --git a/include/isl_set.h b/include/isl_set.h
new file mode 100644 (file)
index 0000000..eb49fd9
--- /dev/null
@@ -0,0 +1,130 @@
+#ifndef ISL_SET_H
+#define ISL_SET_H
+
+#include "isl_map.h"
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+/* A "basic set" is a basic map with a zero-dimensional
+ * domain.
+ */
+struct isl_basic_set {
+       int ref;
+#define ISL_PRIMITIVE_SET_FINAL                (1 << 0)
+       unsigned flags;
+
+       unsigned nparam;
+       unsigned zero;
+       unsigned dim;
+       unsigned extra;
+
+       unsigned n_eq;
+       unsigned n_ineq;
+
+       size_t c_size;
+       isl_int **eq;
+       isl_int **ineq;
+
+       unsigned n_div;
+
+       isl_int **div;
+
+       struct isl_blk block;
+       struct isl_blk block2;
+};
+
+/* A "set" is a (possibly disjoint) union of basic sets.
+ *
+ * See the documentation of isl_map.
+ */
+struct isl_set {
+       int ref;
+#define ISL_SET_DISJOINT               (1 << 0)
+       unsigned flags;
+
+       unsigned nparam;
+       unsigned zero;
+       unsigned dim;
+
+       int n;
+
+       size_t size;
+       struct isl_basic_set *p[0];
+};
+
+struct isl_basic_set *isl_basic_set_alloc(struct isl_ctx *ctx,
+               unsigned nparam, unsigned dim, unsigned extra,
+               unsigned n_eq, unsigned n_ineq);
+struct isl_basic_set *isl_basic_set_extend(struct isl_ctx *ctx,
+               struct isl_basic_set *base,
+               unsigned nparam, unsigned dim, unsigned extra,
+               unsigned n_eq, unsigned n_ineq);
+struct isl_basic_set *isl_basic_set_finalize(struct isl_ctx *ctx,
+               struct isl_basic_set *bset);
+void isl_basic_set_free(struct isl_ctx *ctx, struct isl_basic_set *bset);
+struct isl_basic_set *isl_basic_set_copy(struct isl_ctx *ctx,
+                                       struct isl_basic_set *bset);
+void isl_basic_set_dump(struct isl_ctx *ctx, struct isl_basic_set *bset,
+                               FILE *out, int indent);
+struct isl_basic_set *isl_basic_set_swap_vars(struct isl_ctx *ctx,
+               struct isl_basic_set *bset, unsigned n);
+struct isl_basic_set *isl_basic_set_drop_vars(struct isl_ctx *ctx,
+               struct isl_basic_set *bset, unsigned first, unsigned n);
+struct isl_basic_set *isl_basic_set_intersect(
+               struct isl_ctx *ctx, struct isl_basic_set *bset1,
+               struct isl_basic_set *bset2);
+struct isl_basic_set *isl_basic_set_affine_hull(struct isl_ctx *ctx,
+                                               struct isl_basic_set *bset);
+struct isl_basic_set *isl_basic_set_simplify(
+               struct isl_ctx *ctx, struct isl_basic_set *bset);
+
+struct isl_set *isl_basic_set_lexmin(struct isl_ctx *ctx,
+               struct isl_basic_set *bset);
+struct isl_set *isl_basic_set_union(
+               struct isl_ctx *ctx, struct isl_basic_set *bset1,
+               struct isl_basic_set *bset2);
+
+struct isl_set *isl_set_alloc(struct isl_ctx *ctx,
+               unsigned nparam, unsigned dim, int n, unsigned flags);
+struct isl_set *isl_set_empty(struct isl_ctx *ctx,
+               unsigned nparam, unsigned dim);
+struct isl_set *isl_set_add(struct isl_ctx *ctx, struct isl_set *set,
+                               struct isl_basic_set *bset);
+struct isl_set *isl_set_finalize(struct isl_ctx *ctx, struct isl_set *set);
+struct isl_set *isl_set_copy(struct isl_ctx *ctx, struct isl_set *set);
+void isl_set_free(struct isl_ctx *ctx, struct isl_set *set);
+struct isl_set *isl_set_dup(struct isl_ctx *ctx, struct isl_set *set);
+struct isl_set *isl_set_from_basic_set(struct isl_ctx *ctx,
+                               struct isl_basic_set *bset);
+struct isl_basic_set *isl_set_affine_hull(struct isl_ctx *ctx,
+               struct isl_set *set);
+
+struct isl_set *isl_set_union_disjoint(struct isl_ctx *ctx,
+                       struct isl_set *set1, struct isl_set *set2);
+struct isl_set *isl_set_union(struct isl_ctx *ctx,
+                       struct isl_set *set1, struct isl_set *set2);
+struct isl_set *isl_set_subtract(struct isl_ctx *ctx, struct isl_set *set1,
+               struct isl_set *set2);
+struct isl_set *isl_set_apply(struct isl_ctx *ctx, struct isl_set *set,
+               struct isl_map *map);
+
+void isl_set_dump(struct isl_ctx *ctx, struct isl_set *set, FILE *out,
+                 int indent);
+struct isl_set *isl_set_swap_vars(struct isl_ctx *ctx,
+               struct isl_set *set, unsigned n);
+int isl_set_is_empty(struct isl_ctx *ctx, struct isl_set *set);
+int isl_set_is_subset(struct isl_ctx *ctx,
+               struct isl_set *set1, struct isl_set *set2);
+int isl_set_is_equal(struct isl_ctx *ctx,
+               struct isl_set *set1, struct isl_set *set2);
+
+struct isl_set *isl_basic_set_compute_divs(struct isl_ctx *ctx,
+               struct isl_basic_set *bset);
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
diff --git a/include/isl_set_polylib.h b/include/isl_set_polylib.h
new file mode 100644 (file)
index 0000000..7e01ae4
--- /dev/null
@@ -0,0 +1,25 @@
+#ifndef ISL_SET_POLYLIB_H
+#define ISL_SET_POLYLIB_H
+
+#include <isl_set.h>
+#include <isl_polylib.h>
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+struct isl_basic_set *isl_basic_set_new_from_polylib(
+                       struct isl_ctx *ctx,
+                       Polyhedron *P, unsigned nparam, unsigned dim);
+Polyhedron *isl_basic_set_to_polylib(struct isl_ctx *ctx,
+                                               struct isl_basic_set *bset);
+
+struct isl_set *isl_set_new_from_polylib(struct isl_ctx *ctx,
+                       Polyhedron *D, unsigned nparam, unsigned dim);
+Polyhedron *isl_set_to_polylib(struct isl_ctx *ctx, struct isl_set *set);
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
diff --git a/isl_blk.c b/isl_blk.c
new file mode 100644 (file)
index 0000000..a7a4496
--- /dev/null
+++ b/isl_blk.c
@@ -0,0 +1,45 @@
+#include "isl_blk.h"
+
+struct isl_blk isl_blk_alloc(struct isl_ctx *ctx, size_t n)
+{
+       struct isl_blk block;
+
+       block.data = isl_alloc_array(ctx, isl_int, n);
+       if (!block.data)
+               block.size = -1;
+       else {
+               int i;
+               block.size = n;
+               for (i = 0; i < n; ++i)
+                       isl_int_init(block.data[i]);
+       }
+
+       return block;
+}
+
+struct isl_blk isl_blk_extend(struct isl_ctx *ctx, struct isl_blk block,
+                               size_t new_n)
+{
+       if (block.size >= new_n)
+               return block;
+       block.data = isl_realloc_array(ctx, block.data, isl_int, new_n);
+       if (!block.data)
+               block.size = -1;
+       else {
+               int i;
+               for (i = block.size; i < new_n; ++i)
+                       isl_int_init(block.data[i]);
+               block.size = new_n;
+       }
+
+       return block;
+}
+
+void isl_blk_free(struct isl_ctx *ctx, struct isl_blk block)
+{
+       int i;
+
+       for (i = 0; i < block.size; ++i)
+               isl_int_clear(block.data[i]);
+       free(block.data);
+}
diff --git a/isl_ctx.c b/isl_ctx.c
new file mode 100644 (file)
index 0000000..686a9be
--- /dev/null
+++ b/isl_ctx.c
@@ -0,0 +1,30 @@
+#include "isl_ctx.h"
+#ifdef ISL_POLYLIB
+#include <polylib/polylibgmp.h>
+#endif
+
+struct isl_ctx *isl_ctx_alloc()
+{
+       struct isl_ctx *ctx;
+
+       ctx = isl_alloc_type(NULL, struct isl_ctx);
+       if (!ctx)
+               return NULL;
+
+       isl_int_init(ctx->one);
+       isl_int_set_si(ctx->one, 1);
+
+#ifdef ISL_POLYLIB
+       ctx->MaxRays = POL_NO_DUAL | POL_INTEGER;
+#endif
+
+       return ctx;
+}
+
+void isl_ctx_free(struct isl_ctx *ctx)
+{
+       if (!ctx)
+               return;
+       isl_int_clear(ctx->one);
+       free(ctx);
+}
diff --git a/isl_gmp.c b/isl_gmp.c
new file mode 100644 (file)
index 0000000..0622763
--- /dev/null
+++ b/isl_gmp.c
@@ -0,0 +1,19 @@
+#include "isl_int.h"
+
+u_int32_t isl_gmp_hash(mpz_t v, u_int32_t hash)
+{
+       int sa = v[0]._mp_size;
+       int abs_sa = sa < 0 ? -sa : sa;
+       unsigned char *data = (unsigned char *)v[0]._mp_d;
+       unsigned char *end = data + abs_sa * sizeof(v[0]._mp_d[0]);
+
+       if (sa < 0) {
+               hash *= 16777619;
+               hash ^= 0xFF;
+       }
+       for (; data < end; ++data) {
+               hash *= 16777619;
+               hash ^= *data;
+       }
+       return hash;
+}
diff --git a/isl_lp.c b/isl_lp.c
new file mode 100644 (file)
index 0000000..2e3d5fc
--- /dev/null
+++ b/isl_lp.c
@@ -0,0 +1,15 @@
+#include "isl_ctx.h"
+#include "isl_lp.h"
+#ifdef ISL_PIPLIB
+#include "isl_lp_piplib.h"
+#endif
+
+enum isl_lp_result isl_solve_lp(struct isl_basic_map *bmap, int maximize,
+                                     isl_int *f, isl_int denom, isl_int *opt)
+{
+#ifdef ISL_PIPLIB
+       return isl_pip_solve_lp(bmap, maximize, f, denom, opt);
+#else
+       return isl_lp_error;
+#endif
+}
diff --git a/isl_lp_piplib.c b/isl_lp_piplib.c
new file mode 100644 (file)
index 0000000..0072295
--- /dev/null
@@ -0,0 +1,52 @@
+#include "isl_map.h"
+#include "isl_lp.h"
+#include "isl_piplib.h"
+#include "isl_map_piplib.h"
+
+enum isl_lp_result isl_pip_solve_lp(struct isl_basic_map *bmap, int maximize,
+                                     isl_int *f, isl_int denom, isl_int *opt)
+{
+       enum isl_lp_result res = isl_lp_ok;
+       PipMatrix       *domain = NULL;
+       PipOptions      *options;
+       PipQuast        *sol;
+       unsigned         total;
+
+       total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+       domain = isl_basic_map_to_pip(bmap, 0, 1, 0);
+       if (!domain)
+               goto error;
+       entier_set_si(domain->p[0][1], -1);
+       isl_seq_cpy_to_pip(domain->p[0]+2, f, total);
+
+       options = pip_options_init();
+       if (!options)
+               goto error;
+       options->Urs_unknowns = -1;
+       options->Maximize = maximize;
+       options->Nq = 0;
+       sol = pip_solve(domain, NULL, -1, options);
+       pip_options_free(options);
+       if (!sol)
+               goto error;
+
+       if (!sol->list)
+               res = isl_lp_empty;
+       else if (entier_zero_p(sol->list->vector->the_deno[0]))
+               res = isl_lp_unbounded;
+       else {
+               if (maximize)
+                       mpz_fdiv_q(*opt, sol->list->vector->the_vector[0],
+                                        sol->list->vector->the_deno[0]);
+               else
+                       mpz_cdiv_q(*opt, sol->list->vector->the_vector[0],
+                                        sol->list->vector->the_deno[0]);
+       }
+       pip_matrix_free(domain);
+       pip_quast_free(sol);
+       return res;
+error:
+       if (domain)
+               pip_matrix_free(domain);
+       return isl_lp_error;
+}
diff --git a/isl_map.c b/isl_map.c
new file mode 100644 (file)
index 0000000..12d6780
--- /dev/null
+++ b/isl_map.c
@@ -0,0 +1,2769 @@
+#include <string.h>
+#include <strings.h>
+#include "isl_ctx.h"
+#include "isl_blk.h"
+#include "isl_seq.h"
+#include "isl_set.h"
+#include "isl_map.h"
+#include "isl_map_private.h"
+#ifdef ISL_PIPLIB
+#include "isl_map_piplib.h"
+#endif
+
+static struct isl_basic_map *basic_map_init(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap,
+               unsigned nparam, unsigned n_in, unsigned n_out, unsigned extra,
+               unsigned n_eq, unsigned n_ineq)
+{
+       int i;
+       size_t row_size = 1 + nparam + n_in + n_out + extra;
+
+       bmap->block = isl_blk_alloc(ctx, (n_eq + n_ineq) * row_size);
+       if (!bmap->block.data) {
+               free(bmap);
+               return NULL;
+       }
+
+       bmap->eq = isl_alloc_array(ctx, isl_int *, n_eq + n_ineq);
+       if (!bmap->eq) {
+               isl_blk_free(ctx, bmap->block);
+               free(bmap);
+               return NULL;
+       }
+
+       if (extra == 0) {
+               bmap->block2.size = 0;
+               bmap->block2.data = NULL;
+               bmap->div = NULL;
+       } else {
+               bmap->block2 = isl_blk_alloc(ctx, extra * (1 + row_size));
+               if (!bmap->block2.data) {
+                       free(bmap->eq);
+                       isl_blk_free(ctx, bmap->block);
+                       free(bmap);
+                       return NULL;
+               }
+
+               bmap->div = isl_alloc_array(ctx, isl_int *, extra);
+               if (!bmap->div) {
+                       isl_blk_free(ctx, bmap->block2);
+                       free(bmap->eq);
+                       isl_blk_free(ctx, bmap->block);
+                       free(bmap);
+                       return NULL;
+               }
+       }
+
+       for (i = 0; i < n_eq + n_ineq; ++i)
+               bmap->eq[i] = bmap->block.data + i * row_size;
+
+       for (i = 0; i < extra; ++i)
+               bmap->div[i] = bmap->block2.data + i * (1 + row_size);
+
+       bmap->ref = 1;
+       bmap->flags = 0;
+       bmap->c_size = n_eq + n_ineq;
+       bmap->ineq = bmap->eq + n_eq;
+       bmap->nparam = nparam;
+       bmap->n_in = n_in;
+       bmap->n_out = n_out;
+       bmap->extra = extra;
+       bmap->n_eq = 0;
+       bmap->n_ineq = 0;
+       bmap->n_div = 0;
+
+       return bmap;
+}
+
+struct isl_basic_set *isl_basic_set_alloc(struct isl_ctx *ctx,
+               unsigned nparam, unsigned dim, unsigned extra,
+               unsigned n_eq, unsigned n_ineq)
+{
+       struct isl_basic_map *bmap;
+       bmap = isl_basic_map_alloc(ctx, nparam, 0, dim, extra, n_eq, n_ineq);
+       return (struct isl_basic_set *)bmap;
+}
+
+struct isl_basic_map *isl_basic_map_alloc(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out, unsigned extra,
+               unsigned n_eq, unsigned n_ineq)
+{
+       struct isl_basic_map *bmap;
+
+       bmap = isl_alloc_type(ctx, struct isl_basic_map);
+       if (!bmap)
+               return NULL;
+
+       return basic_map_init(ctx, bmap,
+                              nparam, in, out, extra, n_eq, n_ineq);
+}
+
+static void dup_constraints(struct isl_ctx *ctx,
+               struct isl_basic_map *dst, struct isl_basic_map *src)
+{
+       int i;
+       unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
+
+       for (i = 0; i < src->n_eq; ++i) {
+               int j = isl_basic_map_alloc_equality(ctx, dst);
+               isl_seq_cpy(dst->eq[j], src->eq[i], 1+total);
+       }
+
+       for (i = 0; i < src->n_ineq; ++i) {
+               int j = isl_basic_map_alloc_inequality(ctx, dst);
+               isl_seq_cpy(dst->ineq[j], src->ineq[i], 1+total);
+       }
+
+       for (i = 0; i < src->n_div; ++i) {
+               int j = isl_basic_map_alloc_div(ctx, dst);
+               isl_seq_cpy(dst->div[j], src->div[i], 1+1+total);
+       }
+       F_SET(dst, ISL_PRIMITIVE_SET_FINAL);
+}
+
+struct isl_basic_map *isl_basic_map_dup(struct isl_ctx *ctx,
+                                       struct isl_basic_map *bmap)
+{
+       struct isl_basic_map *dup;
+
+       dup = isl_basic_map_alloc(ctx, bmap->nparam,
+                       bmap->n_in, bmap->n_out,
+                       bmap->n_div, bmap->n_eq, bmap->n_ineq);
+       if (!dup)
+               return NULL;
+       dup_constraints(ctx, dup, bmap);
+       return dup;
+}
+
+struct isl_basic_set *isl_basic_set_dup(struct isl_ctx *ctx,
+                                       struct isl_basic_set *bset)
+{
+       struct isl_basic_map *dup;
+
+       dup = isl_basic_map_dup(ctx, (struct isl_basic_map *)bset);
+       return (struct isl_basic_set *)dup;
+}
+
+struct isl_basic_set *isl_basic_set_copy(struct isl_ctx *ctx,
+                                       struct isl_basic_set *bset)
+{
+       if (!bset)
+               return NULL;
+
+       if (F_ISSET(bset, ISL_PRIMITIVE_SET_FINAL)) {
+               bset->ref++;
+               return bset;
+       }
+       return isl_basic_set_dup(ctx, bset);
+}
+
+struct isl_set *isl_set_copy(struct isl_ctx *ctx, struct isl_set *set)
+{
+       if (!set)
+               return NULL;
+
+       set->ref++;
+       return set;
+}
+
+struct isl_basic_map *isl_basic_map_copy(struct isl_ctx *ctx,
+                                       struct isl_basic_map *bmap)
+{
+       if (!bmap)
+               return NULL;
+
+       if (F_ISSET(bmap, ISL_PRIMITIVE_SET_FINAL)) {
+               bmap->ref++;
+               return bmap;
+       }
+       return isl_basic_map_dup(ctx, bmap);
+}
+
+struct isl_map *isl_map_copy(struct isl_ctx *ctx, struct isl_map *map)
+{
+       if (!map)
+               return NULL;
+
+       map->ref++;
+       return map;
+}
+
+void isl_basic_map_free(struct isl_ctx *ctx, struct isl_basic_map *bmap)
+{
+       if (!bmap)
+               return;
+
+       if (--bmap->ref > 0)
+               return;
+
+       free(bmap->div);
+       isl_blk_free(ctx, bmap->block2);
+       free(bmap->eq);
+       isl_blk_free(ctx, bmap->block);
+       free(bmap);
+}
+
+void isl_basic_set_free(struct isl_ctx *ctx, struct isl_basic_set *bset)
+{
+       isl_basic_map_free(ctx, (struct isl_basic_map *)bset);
+}
+
+int isl_basic_map_alloc_equality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       if (!bmap)
+               return -1;
+       isl_assert(ctx, bmap->n_eq + bmap->n_ineq < bmap->c_size, return -1);
+       isl_assert(ctx, bmap->eq + bmap->n_eq <= bmap->ineq, return -1);
+       if (bmap->eq + bmap->n_eq == bmap->ineq) {
+               isl_int *t;
+               int j = isl_basic_map_alloc_inequality(ctx, bmap);
+               if (j < 0)
+                       return -1;
+               t = bmap->ineq[0];
+               bmap->ineq[0] = bmap->ineq[j];
+               bmap->ineq[j] = t;
+               bmap->n_ineq--;
+               bmap->ineq++;
+       }
+       return bmap->n_eq++;
+}
+
+int isl_basic_map_free_equality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned n)
+{
+       isl_assert(ctx, n <= bmap->n_eq, return -1);
+       bmap->n_eq -= n;
+       return 0;
+}
+
+int isl_basic_map_drop_equality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned pos)
+{
+       isl_int *t;
+       isl_assert(ctx, pos < bmap->n_eq, return -1);
+
+       if (pos != bmap->n_eq - 1) {
+               t = bmap->eq[pos];
+               bmap->eq[pos] = bmap->eq[bmap->n_eq - 1];
+               bmap->eq[bmap->n_eq - 1] = t;
+       }
+       bmap->n_eq--;
+       return 0;
+}
+
+int isl_basic_map_alloc_inequality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       isl_assert(ctx, (bmap->ineq - bmap->eq) + bmap->n_ineq < bmap->c_size,
+                       return -1);
+       return bmap->n_ineq++;
+}
+
+int isl_basic_map_free_inequality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned n)
+{
+       isl_assert(ctx, n <= bmap->n_ineq, return -1);
+       bmap->n_ineq -= n;
+       return 0;
+}
+
+int isl_basic_map_drop_inequality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned pos)
+{
+       isl_int *t;
+       isl_assert(ctx, pos < bmap->n_ineq, return -1);
+
+       if (pos != bmap->n_ineq - 1) {
+               t = bmap->ineq[pos];
+               bmap->ineq[pos] = bmap->ineq[bmap->n_ineq - 1];
+               bmap->ineq[bmap->n_ineq - 1] = t;
+       }
+       bmap->n_ineq--;
+       return 0;
+}
+
+int isl_basic_map_alloc_div(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       isl_assert(ctx, bmap->n_div < bmap->extra, return -1);
+       return bmap->n_div++;
+}
+
+int isl_basic_map_free_div(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned n)
+{
+       isl_assert(ctx, n <= bmap->n_div, return -1);
+       bmap->n_div -= n;
+       return 0;
+}
+
+/* Copy constraint from src to dst, putting the vars of src at offset
+ * dim_off in dst and the divs of src at offset div_off in dst.
+ * If both sets are actually map, then dim_off applies to the input
+ * variables.
+ */
+static void copy_constraint(struct isl_basic_map *dst_map, isl_int *dst,
+                           struct isl_basic_map *src_map, isl_int *src,
+                           unsigned in_off, unsigned out_off, unsigned div_off)
+{
+       isl_int_set(dst[0], src[0]);
+       isl_seq_cpy(dst+1, src+1, isl_min(dst_map->nparam, src_map->nparam));
+       isl_seq_cpy(dst+1+dst_map->nparam+in_off,
+                   src+1+src_map->nparam,
+                   isl_min(dst_map->n_in-in_off, src_map->n_in));
+       isl_seq_cpy(dst+1+dst_map->nparam+dst_map->n_in+out_off,
+                   src+1+src_map->nparam+src_map->n_in,
+                   isl_min(dst_map->n_out-out_off, src_map->n_out));
+       isl_seq_cpy(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out+div_off,
+                   src+1+src_map->nparam+src_map->n_in+src_map->n_out,
+                   isl_min(dst_map->extra-div_off, src_map->extra));
+}
+
+static void copy_div(struct isl_basic_map *dst_map, isl_int *dst,
+                    struct isl_basic_map *src_map, isl_int *src,
+                    unsigned in_off, unsigned out_off, unsigned div_off)
+{
+       isl_int_set(dst[0], src[0]);
+       copy_constraint(dst_map, dst+1, src_map, src+1, in_off, out_off, div_off);
+}
+
+static struct isl_basic_map *add_constraints(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2, unsigned i_pos, unsigned o_pos)
+{
+       int i;
+       unsigned div_off = bmap1->n_div;
+
+       for (i = 0; i < bmap2->n_eq; ++i) {
+               int i1 = isl_basic_map_alloc_equality(ctx, bmap1);
+               if (i1 < 0)
+                       goto error;
+               copy_constraint(bmap1, bmap1->eq[i1], bmap2, bmap2->eq[i],
+                               i_pos, o_pos, div_off);
+       }
+
+       for (i = 0; i < bmap2->n_ineq; ++i) {
+               int i1 = isl_basic_map_alloc_inequality(ctx, bmap1);
+               if (i1 < 0)
+                       goto error;
+               copy_constraint(bmap1, bmap1->ineq[i1], bmap2, bmap2->ineq[i],
+                               i_pos, o_pos, div_off);
+       }
+
+       for (i = 0; i < bmap2->n_div; ++i) {
+               int i1 = isl_basic_map_alloc_div(ctx, bmap1);
+               if (i1 < 0)
+                       goto error;
+               copy_div(bmap1, bmap1->div[i1], bmap2, bmap2->div[i],
+                        i_pos, o_pos, div_off);
+       }
+
+       isl_basic_map_free(ctx, bmap2);
+
+       return bmap1;
+
+error:
+       isl_basic_map_free(ctx, bmap1);
+       isl_basic_map_free(ctx, bmap2);
+       return NULL;
+}
+
+struct isl_basic_map *isl_basic_map_extend(struct isl_ctx *ctx,
+               struct isl_basic_map *base,
+               unsigned nparam, unsigned n_in, unsigned n_out, unsigned extra,
+               unsigned n_eq, unsigned n_ineq)
+{
+       struct isl_basic_map *ext;
+       int dims_ok;
+
+       base = isl_basic_map_cow(ctx, base);
+       if (!base)
+               goto error;
+
+       dims_ok = base && base->nparam == nparam &&
+                 base->n_in == n_in && base->n_out == n_out &&
+                 base->extra >= base->n_div + extra;
+
+       if (dims_ok && n_eq == 0 && n_ineq == 0)
+               return base;
+
+       if (base) {
+               isl_assert(ctx, base->nparam <= nparam, goto error);
+               isl_assert(ctx, base->n_in <= n_in, goto error);
+               isl_assert(ctx, base->n_out <= n_out, goto error);
+               extra += base->extra;
+               n_eq += base->n_eq;
+               n_ineq += base->n_ineq;
+       }
+
+       ext = isl_basic_map_alloc(ctx, nparam, n_in, n_out,
+                                       extra, n_eq, n_ineq);
+       if (!base)
+               return ext;
+       if (!ext)
+               goto error;
+
+       ext = add_constraints(ctx, ext, base, 0, 0);
+
+       return ext;
+
+error:
+       isl_basic_map_free(ctx, base);
+       return NULL;
+}
+
+struct isl_basic_set *isl_basic_set_extend(struct isl_ctx *ctx,
+               struct isl_basic_set *base,
+               unsigned nparam, unsigned dim, unsigned extra,
+               unsigned n_eq, unsigned n_ineq)
+{
+       return (struct isl_basic_set *)
+               isl_basic_map_extend(ctx, (struct isl_basic_map *)base,
+                                       nparam, 0, dim, extra, n_eq, n_ineq);
+}
+
+static struct isl_basic_set *isl_basic_set_cow(struct isl_ctx *ctx,
+               struct isl_basic_set *bset)
+{
+       return (struct isl_basic_set *)
+               isl_basic_map_cow(ctx, (struct isl_basic_map *)bset);
+}
+
+struct isl_basic_map *isl_basic_map_cow(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       if (!bmap)
+               return NULL;
+
+       if (bmap->ref > 1) {
+               bmap->ref--;
+               bmap = isl_basic_map_dup(ctx, bmap);
+       }
+       F_CLR(bmap, ISL_PRIMITIVE_SET_FINAL);
+       return bmap;
+}
+
+static struct isl_set *isl_set_cow(struct isl_ctx *ctx, struct isl_set *set)
+{
+       if (set->ref == 1)
+               return set;
+       set->ref--;
+       return isl_set_dup(ctx, set);
+}
+
+struct isl_map *isl_map_cow(struct isl_ctx *ctx, struct isl_map *map)
+{
+       if (map->ref == 1)
+               return map;
+       map->ref--;
+       return isl_map_dup(ctx, map);
+}
+
+static void swap_vars(struct isl_blk blk, isl_int *a,
+                       unsigned a_len, unsigned b_len)
+{
+       isl_seq_cpy(blk.data, a+a_len, b_len);
+       isl_seq_cpy(blk.data+b_len, a, a_len);
+       isl_seq_cpy(a, blk.data, b_len+a_len);
+}
+
+struct isl_basic_set *isl_basic_set_swap_vars(struct isl_ctx *ctx,
+               struct isl_basic_set *bset, unsigned n)
+{
+       int i;
+       struct isl_blk blk;
+
+       if (!bset)
+               goto error;
+
+       isl_assert(ctx, n <= bset->dim, goto error);
+
+       if (n == bset->dim)
+               return bset;
+
+       bset = isl_basic_set_cow(ctx, bset);
+       if (!bset)
+               return NULL;
+
+       blk = isl_blk_alloc(ctx, bset->dim);
+       if (!blk.data)
+               goto error;
+
+       for (i = 0; i < bset->n_eq; ++i)
+               swap_vars(blk,
+                         bset->eq[i]+1+bset->nparam, n, bset->dim - n);
+
+       for (i = 0; i < bset->n_ineq; ++i)
+               swap_vars(blk,
+                         bset->ineq[i]+1+bset->nparam, n, bset->dim - n);
+
+       for (i = 0; i < bset->n_div; ++i)
+               swap_vars(blk,
+                         bset->div[i]+1+1+bset->nparam, n, bset->dim - n);
+
+       isl_blk_free(ctx, blk);
+
+       return bset;
+
+error:
+       isl_basic_set_free(ctx, bset);
+       return NULL;
+}
+
+struct isl_set *isl_set_swap_vars(struct isl_ctx *ctx,
+               struct isl_set *set, unsigned n)
+{
+       int i;
+       set = isl_set_cow(ctx, set);
+       if (!set)
+               return NULL;
+
+       for (i = 0; i < set->n; ++i) {
+               set->p[i] = isl_basic_set_swap_vars(ctx, set->p[i], n);
+               if (!set->p[i]) {
+                       isl_set_free(ctx, set);
+                       return NULL;
+               }
+       }
+       return set;
+}
+
+static void constraint_drop_vars(isl_int *c, unsigned n, unsigned rem)
+{
+       isl_seq_cpy(c, c + n, rem);
+       isl_seq_clr(c + rem, n);
+}
+
+struct isl_basic_set *isl_basic_set_drop_vars(struct isl_ctx *ctx,
+               struct isl_basic_set *bset, unsigned first, unsigned n)
+{
+       int i;
+
+       if (!bset)
+               goto error;
+
+       isl_assert(ctx, first + n <= bset->dim, goto error);
+
+       if (n == 0)
+               return bset;
+
+       bset = isl_basic_set_cow(ctx, bset);
+       if (!bset)
+               return NULL;
+
+       for (i = 0; i < bset->n_eq; ++i)
+               constraint_drop_vars(bset->eq[i]+1+bset->nparam+first, n,
+                                    (bset->dim-first-n)+bset->extra);
+
+       for (i = 0; i < bset->n_ineq; ++i)
+               constraint_drop_vars(bset->ineq[i]+1+bset->nparam+first, n,
+                                    (bset->dim-first-n)+bset->extra);
+
+       for (i = 0; i < bset->n_div; ++i)
+               constraint_drop_vars(bset->div[i]+1+1+bset->nparam+first, n,
+                                    (bset->dim-first-n)+bset->extra);
+
+       bset->dim -= n;
+       bset->extra += n;
+
+       return isl_basic_set_finalize(ctx, bset);
+error:
+       isl_basic_set_free(ctx, bset);
+       return NULL;
+}
+
+/*
+ * We don't cow, as the div is assumed to be redundant.
+ */
+static struct isl_basic_map *isl_basic_map_drop_div(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned div)
+{
+       int i;
+       unsigned pos;
+
+       if (!bmap)
+               goto error;
+
+       pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
+
+       isl_assert(ctx, div < bmap->n_div, goto error);
+
+       for (i = 0; i < bmap->n_eq; ++i)
+               constraint_drop_vars(bmap->eq[i]+pos, 1, bmap->extra-div-1);
+
+       for (i = 0; i < bmap->n_ineq; ++i)
+               constraint_drop_vars(bmap->ineq[i]+pos, 1, bmap->extra-div-1);
+
+       for (i = 0; i < bmap->n_div; ++i)
+               constraint_drop_vars(bmap->div[i]+1+pos, 1, bmap->extra-div-1);
+
+       if (div != bmap->n_div - 1) {
+               int j;
+               isl_int *t = bmap->div[div];
+
+               for (j = div; j < bmap->n_div - 1; ++j)
+                       bmap->div[j] = bmap->div[j+1];
+
+               bmap->div[bmap->n_div - 1] = t;
+       }
+       isl_basic_map_free_div(ctx, bmap, 1);
+
+       return bmap;
+error:
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+static void eliminate_div(struct isl_ctx *ctx, struct isl_basic_map *bmap,
+                               isl_int *eq, unsigned div)
+{
+       int i;
+       unsigned pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
+       unsigned len;
+       len = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+
+       for (i = 0; i < bmap->n_eq; ++i)
+               if (bmap->eq[i] != eq)
+                       isl_seq_elim(bmap->eq[i], eq, pos, len, NULL);
+
+       for (i = 0; i < bmap->n_ineq; ++i)
+               isl_seq_elim(bmap->ineq[i], eq, pos, len, NULL);
+
+       /* We need to be careful about circular definitions,
+        * so for now we just remove the definitions of other divs that
+        * depend on this div and (possibly) recompute them later.
+        */
+       for (i = 0; i < bmap->n_div; ++i)
+               if (!isl_int_is_zero(bmap->div[i][0]) &&
+                   !isl_int_is_zero(bmap->div[i][1 + pos]))
+                       isl_seq_clr(bmap->div[i], 1 + len);
+
+       isl_basic_map_drop_div(ctx, bmap, div);
+}
+
+struct isl_basic_map *isl_basic_map_set_to_empty(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap)
+{
+       int i = 0;
+       unsigned total;
+       if (!bmap)
+               goto error;
+       total = bmap->nparam + bmap->n_in + bmap->n_out;
+       isl_assert(ctx, bmap->n_ineq + bmap->n_eq > 0, goto error);
+       isl_basic_map_free_div(ctx, bmap, bmap->n_div);
+       isl_basic_map_free_inequality(ctx, bmap, bmap->n_ineq);
+       if (bmap->n_eq > 0)
+               isl_basic_map_free_equality(ctx, bmap, bmap->n_eq-1);
+       else {
+               isl_basic_map_alloc_equality(ctx, bmap);
+               if (i < 0)
+                       goto error;
+       }
+       isl_int_set_si(bmap->eq[i][0], 1);
+       isl_seq_clr(bmap->eq[i]+1, total);
+       return isl_basic_map_finalize(ctx, bmap);
+error:
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+static void swap_equality(struct isl_basic_map *bmap, int a, int b)
+{
+       isl_int *t = bmap->eq[a];
+       bmap->eq[a] = bmap->eq[b];
+       bmap->eq[b] = t;
+}
+
+static void swap_inequality(struct isl_basic_map *bmap, int a, int b)
+{
+       isl_int *t = bmap->ineq[a];
+       bmap->ineq[a] = bmap->ineq[b];
+       bmap->ineq[b] = t;
+}
+
+static void swap_div(struct isl_basic_map *bmap, int a, int b)
+{
+       int i;
+       unsigned off = bmap->nparam + bmap->n_in + bmap->n_out;
+       isl_int *t = bmap->div[a];
+       bmap->div[a] = bmap->div[b];
+       bmap->div[b] = t;
+
+       for (i = 0; i < bmap->n_eq; ++i)
+               isl_int_swap(bmap->eq[i][1+off+a], bmap->eq[i][1+off+b]);
+
+       for (i = 0; i < bmap->n_ineq; ++i)
+               isl_int_swap(bmap->ineq[i][1+off+a], bmap->ineq[i][1+off+b]);
+
+       for (i = 0; i < bmap->n_div; ++i)
+               isl_int_swap(bmap->div[i][1+1+off+a], bmap->div[i][1+1+off+b]);
+}
+
+struct isl_basic_map *isl_basic_map_gauss(struct isl_ctx *ctx,
+       struct isl_basic_map *bmap, int *progress)
+{
+       int k;
+       int done;
+       int last_var;
+       unsigned total_var;
+       unsigned total;
+
+       if (!bmap)
+               return NULL;
+
+       total_var = bmap->nparam + bmap->n_in + bmap->n_out;
+       total = total_var + bmap->n_div;
+
+       last_var = total - 1;
+       for (done = 0; done < bmap->n_eq; ++done) {
+               for (; last_var >= 0; --last_var) {
+                       for (k = done; k < bmap->n_eq; ++k)
+                               if (!isl_int_is_zero(bmap->eq[k][1+last_var]))
+                                       break;
+                       if (k < bmap->n_eq)
+                               break;
+               }
+               if (last_var < 0)
+                       break;
+               if (k != done)
+                       swap_equality(bmap, k, done);
+               if (isl_int_is_neg(bmap->eq[done][1+last_var]))
+                       isl_seq_neg(bmap->eq[done], bmap->eq[done], 1+total);
+
+               for (k = 0; k < bmap->n_eq; ++k) {
+                       if (k == done)
+                               continue;
+                       if (isl_int_is_zero(bmap->eq[k][1+last_var]))
+                               continue;
+                       if (progress)
+                               *progress = 1;
+                       isl_seq_elim(bmap->eq[k], bmap->eq[done],
+                                       1+last_var, 1+total, NULL);
+               }
+
+               for (k = 0; k < bmap->n_ineq; ++k) {
+                       if (isl_int_is_zero(bmap->ineq[k][1+last_var]))
+                               continue;
+                       if (progress)
+                               *progress = 1;
+                       isl_seq_elim(bmap->ineq[k], bmap->eq[done],
+                                       1+last_var, 1+total, NULL);
+               }
+
+               for (k = 0; k < bmap->n_div; ++k) {
+                       if (isl_int_is_zero(bmap->div[k][0]))
+                               continue;
+                       if (isl_int_is_zero(bmap->div[k][1+1+last_var]))
+                               continue;
+                       if (progress)
+                               *progress = 1;
+                       isl_seq_elim(bmap->div[k]+1, bmap->eq[done],
+                                       1+last_var, 1+total, &bmap->div[k][0]);
+               }
+
+               if (last_var >= total_var &&
+                   isl_int_is_zero(bmap->div[last_var - total_var][0])) {
+                       unsigned div = last_var - total_var;
+                       isl_seq_neg(bmap->div[div]+1, bmap->eq[done], 1+total);
+                       isl_int_set_si(bmap->div[div][1+1+last_var], 0);
+                       isl_int_set(bmap->div[div][0],
+                                   bmap->eq[done][1+last_var]);
+               }
+       }
+       return bmap;
+}
+
+static unsigned int round_up(unsigned int v)
+{
+       int old_v = v;
+
+       while (v) {
+               old_v = v;
+               v ^= v & -v;
+       }
+       return old_v << 1;
+}
+
+static struct isl_basic_map *remove_duplicate_constraints(
+       struct isl_ctx *ctx,
+       struct isl_basic_map *bmap, int *progress)
+{
+       unsigned int size;
+       int *index;
+       int k, l, h;
+       int bits;
+       unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+
+       if (bmap->n_ineq <= 1)
+               return bmap;
+
+       size = round_up(4 * bmap->n_ineq / 3 - 1);
+       bits = ffs(size) - 1;
+       index = isl_alloc_array(ctx, int, size);
+       memset(index, 0, size * sizeof(int));
+       if (!index)
+               return bmap;
+
+       index[isl_seq_hash(bmap->ineq[0]+1, total, bits)] = 1;
+       for (k = 1; k < bmap->n_ineq; ++k) {
+               u_int32_t hash = isl_seq_hash(bmap->ineq[k]+1, total, bits);
+               for (h = hash; index[h]; h = (h+1) % size)
+                       if (isl_seq_eq(bmap->ineq[k]+1,
+                                      bmap->ineq[index[h]-1]+1, total))
+                               break;
+               if (!index[h]) {
+                       index[h] = k+1;
+                       continue;
+               }
+               *progress = 1;
+               l = index[h] - 1;
+               if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
+                       swap_inequality(bmap, k, l);
+               isl_basic_map_drop_inequality(ctx, bmap, k);
+               --k;
+       }
+
+       free(index);
+       return bmap;
+}
+
+static struct isl_basic_map *remove_duplicate_divs(struct isl_ctx *ctx,
+       struct isl_basic_map *bmap, int *progress)
+{
+       unsigned int size;
+       int *index;
+       int k, l, h;
+       int bits;
+       struct isl_blk eq;
+       unsigned total_var = bmap->nparam + bmap->n_in + bmap->n_out;
+       unsigned total = total_var + bmap->n_div;
+
+       if (bmap->n_div <= 1)
+               return bmap;
+
+       for (k = bmap->n_div - 1; k >= 0; --k)
+               if (!isl_int_is_zero(bmap->div[k][0]))
+                       break;
+       if (k <= 0)
+               return bmap;
+
+       size = round_up(4 * bmap->n_div / 3 - 1);
+       bits = ffs(size) - 1;
+       index = isl_alloc_array(ctx, int, size);
+       memset(index, 0, size * sizeof(int));
+       if (!index)
+               return bmap;
+       eq = isl_blk_alloc(ctx, 1+total);
+       if (!eq.data)
+               goto out;
+
+       isl_seq_clr(eq.data, 1+total);
+       index[isl_seq_hash(bmap->div[k], 2+total, bits)] = k + 1;
+       for (--k; k >= 0; --k) {
+               u_int32_t hash;
+
+               if (isl_int_is_zero(bmap->div[k][0]))
+                       continue;
+
+               hash = isl_seq_hash(bmap->div[k], 2+total, bits);
+               for (h = hash; index[h]; h = (h+1) % size)
+                       if (isl_seq_eq(bmap->div[k],
+                                      bmap->div[index[h]-1], 2+total))
+                               break;
+               if (index[h]) {
+                       *progress = 1;
+                       l = index[h] - 1;
+                       isl_int_set_si(eq.data[1+total_var+k], -1);
+                       isl_int_set_si(eq.data[1+total_var+l], 1);
+                       eliminate_div(ctx, bmap, eq.data, l);
+                       isl_int_set_si(eq.data[1+total_var+k], 0);
+                       isl_int_set_si(eq.data[1+total_var+l], 0);
+               }
+               index[h] = k+1;
+       }
+
+       isl_blk_free(ctx, eq);
+out:
+       free(index);
+       return bmap;
+}
+
+static struct isl_basic_map *eliminate_divs(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap,
+               int *progress)
+{
+       int d;
+       int i;
+       int modified = 0;
+       unsigned off;
+
+       if (!bmap)
+               return NULL;
+
+       off = 1 + bmap->nparam + bmap->n_in + bmap->n_out;
+
+       for (d = bmap->n_div - 1; d >= 0 ; --d) {
+               for (i = 0; i < bmap->n_eq; ++i) {
+                       if (!isl_int_is_one(bmap->eq[i][off + d]) &&
+                           !isl_int_is_negone(bmap->eq[i][off + d]))
+                               continue;
+                       modified = 1;
+                       *progress = 1;
+                       eliminate_div(ctx, bmap, bmap->eq[i], d);
+                       isl_basic_map_drop_equality(ctx, bmap, i);
+               }
+       }
+       if (modified)
+               return eliminate_divs(ctx, bmap, progress);
+       return bmap;
+}
+
+static struct isl_basic_map *normalize_constraints(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap)
+{
+       int i;
+       isl_int gcd;
+       unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+
+       isl_int_init(gcd);
+       for (i = bmap->n_eq - 1; i >= 0; --i) {
+               isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
+               if (isl_int_is_zero(gcd)) {
+                       if (!isl_int_is_zero(bmap->eq[i][0])) {
+                               bmap = isl_basic_map_set_to_empty(ctx, bmap);
+                               break;
+                       }
+                       isl_basic_map_drop_equality(ctx, bmap, i);
+                       continue;
+               }
+               if (isl_int_is_one(gcd))
+                       continue;
+               if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
+                       bmap = isl_basic_map_set_to_empty(ctx, bmap);
+                       break;
+               }
+               isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
+       }
+
+       for (i = bmap->n_ineq - 1; i >= 0; --i) {
+               isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
+               if (isl_int_is_zero(gcd)) {
+                       if (isl_int_is_neg(bmap->ineq[i][0])) {
+                               bmap = isl_basic_map_set_to_empty(ctx, bmap);
+                               break;
+                       }
+                       isl_basic_map_drop_inequality(ctx, bmap, i);
+                       continue;
+               }
+               if (isl_int_is_one(gcd))
+                       continue;
+               isl_int_cdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
+               isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
+       }
+       isl_int_clear(gcd);
+
+       return bmap;
+}
+
+struct isl_basic_map *isl_basic_map_simplify(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap)
+{
+       int progress = 1;
+       if (!bmap)
+               return NULL;
+       while (progress) {
+               progress = 0;
+               bmap = normalize_constraints(ctx, bmap);
+               bmap = eliminate_divs(ctx, bmap, &progress);
+               bmap = isl_basic_map_gauss(ctx, bmap, &progress);
+               bmap = remove_duplicate_divs(ctx, bmap, &progress);
+               bmap = remove_duplicate_constraints(ctx, bmap, &progress);
+       }
+       return bmap;
+}
+
+struct isl_basic_set *isl_basic_set_simplify(
+               struct isl_ctx *ctx, struct isl_basic_set *bset)
+{
+       return (struct isl_basic_set *)
+               isl_basic_map_simplify(ctx,
+                       (struct isl_basic_map *)bset);
+}
+
+static void dump_term(struct isl_basic_map *bmap,
+                       isl_int c, int pos, FILE *out)
+{
+       unsigned in = bmap->n_in;
+       unsigned dim = bmap->n_in + bmap->n_out;
+       if (!pos)
+               isl_int_print(out, c);
+       else {
+               if (!isl_int_is_one(c))
+                       isl_int_print(out, c);
+               if (pos < 1 + bmap->nparam)
+                       fprintf(out, "p%d", pos - 1);
+               else if (pos < 1 + bmap->nparam + in)
+                       fprintf(out, "i%d", pos - 1 - bmap->nparam);
+               else if (pos < 1 + bmap->nparam + dim)
+                       fprintf(out, "o%d", pos - 1 - bmap->nparam - in);
+               else
+                       fprintf(out, "e%d", pos - 1 - bmap->nparam - dim);
+       }
+}
+
+static void dump_constraint_sign(struct isl_basic_map *bmap, isl_int *c,
+                               int sign, FILE *out)
+{
+       int i;
+       int first;
+       unsigned len = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+       isl_int v;
+
+       isl_int_init(v);
+       for (i = 0, first = 1; i < len; ++i) {
+               if (isl_int_sgn(c[i]) * sign <= 0)
+                       continue;
+               if (!first)
+                       fprintf(out, " + ");
+               first = 0;
+               isl_int_abs(v, c[i]);
+               dump_term(bmap, v, i, out);
+       }
+       isl_int_clear(v);
+       if (first)
+               fprintf(out, "0");
+}
+
+static void dump_constraint(struct isl_basic_map *bmap, isl_int *c,
+                               const char *op, FILE *out, int indent)
+{
+       fprintf(out, "%*s", indent, "");
+
+       dump_constraint_sign(bmap, c, 1, out);
+       fprintf(out, " %s ", op);
+       dump_constraint_sign(bmap, c, -1, out);
+
+       fprintf(out, "\n");
+}
+
+static void dump_constraints(struct isl_basic_map *bmap,
+                               isl_int **c, unsigned n,
+                               const char *op, FILE *out, int indent)
+{
+       int i;
+
+       for (i = 0; i < n; ++i)
+               dump_constraint(bmap, c[i], op, out, indent);
+}
+
+static void dump_affine(struct isl_basic_map *bmap, isl_int *exp, FILE *out)
+{
+       int j;
+       int first = 1;
+       unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+
+       for (j = 0; j < 1 + total; ++j) {
+               if (isl_int_is_zero(exp[j]))
+                       continue;
+               if (!first && isl_int_is_pos(exp[j]))
+                       fprintf(out, "+");
+               dump_term(bmap, exp[j], j, out);
+               first = 0;
+       }
+}
+
+static void dump(struct isl_basic_map *bmap, FILE *out, int indent)
+{
+       int i;
+
+       dump_constraints(bmap, bmap->eq, bmap->n_eq, "=", out, indent);
+       dump_constraints(bmap, bmap->ineq, bmap->n_ineq, ">=", out, indent);
+
+       for (i = 0; i < bmap->n_div; ++i) {
+               fprintf(out, "%*s", indent, "");
+               fprintf(out, "e%d = [(", i);
+               dump_affine(bmap, bmap->div[i]+1, out);
+               fprintf(out, ")/");
+               isl_int_print(out, bmap->div[i][0]);
+               fprintf(out, "]\n");
+       }
+}
+
+void isl_basic_set_dump(struct isl_ctx *ctx, struct isl_basic_set *bset,
+                               FILE *out, int indent)
+{
+       fprintf(out, "%*s", indent, "");
+       fprintf(out, "nparam: %d, dim: %d, extra: %d\n",
+                       bset->nparam, bset->dim, bset->extra);
+       dump((struct isl_basic_map *)bset, out, indent);
+}
+
+void isl_basic_map_dump(struct isl_ctx *ctx, struct isl_basic_map *bmap,
+                               FILE *out, int indent)
+{
+       fprintf(out, "%*s", indent, "");
+       fprintf(out, "ref: %d, nparam: %d, in: %d, out: %d, extra: %d\n",
+               bmap->ref,
+               bmap->nparam, bmap->n_in, bmap->n_out, bmap->extra);
+       dump(bmap, out, indent);
+}
+
+int isl_inequality_negate(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned pos)
+{
+       unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+       isl_assert(ctx, pos < bmap->n_ineq, return -1);
+       isl_seq_neg(bmap->ineq[pos], bmap->ineq[pos], 1 + total);
+       isl_int_sub_ui(bmap->ineq[pos][0], bmap->ineq[pos][0], 1);
+       return 0;
+}
+
+struct isl_set *isl_set_alloc(struct isl_ctx *ctx,
+               unsigned nparam, unsigned dim, int n, unsigned flags)
+{
+       struct isl_set *set;
+
+       set = isl_alloc(ctx, struct isl_set,
+                       sizeof(struct isl_set) +
+                       n * sizeof(struct isl_basic_set *));
+       if (!set)
+               return NULL;
+
+       set->ref = 1;
+       set->size = n;
+       set->n = 0;
+       set->nparam = nparam;
+       set->zero = 0;
+       set->dim = dim;
+       set->flags = flags;
+       return set;
+}
+
+struct isl_set *isl_set_dup(struct isl_ctx *ctx, struct isl_set *set)
+{
+       int i;
+       struct isl_set *dup;
+
+       dup = isl_set_alloc(ctx, set->nparam, set->dim, set->n, set->flags);
+       if (!dup)
+               return NULL;
+       for (i = 0; i < set->n; ++i)
+               dup = isl_set_add(ctx, dup,
+                                 isl_basic_set_dup(ctx, set->p[i]));
+       return dup;
+}
+
+struct isl_set *isl_set_from_basic_set(struct isl_ctx *ctx,
+                               struct isl_basic_set *bset)
+{
+       struct isl_set *set;
+
+       set = isl_set_alloc(ctx, bset->nparam, bset->dim, 1, ISL_MAP_DISJOINT);
+       if (!set) {
+               isl_basic_set_free(ctx, bset);
+               return NULL;
+       }
+       return isl_set_add(ctx, set, bset);
+}
+
+struct isl_map *isl_map_from_basic_map(struct isl_ctx *ctx,
+                               struct isl_basic_map *bmap)
+{
+       struct isl_map *map;
+
+       map = isl_map_alloc(ctx, bmap->nparam, bmap->n_in, bmap->n_out, 1,
+                               ISL_MAP_DISJOINT);
+       if (!map) {
+               isl_basic_map_free(ctx, bmap);
+               return NULL;
+       }
+       return isl_map_add(ctx, map, bmap);
+}
+
+struct isl_set *isl_set_add(struct isl_ctx *ctx, struct isl_set *set,
+                               struct isl_basic_set *bset)
+{
+       if (!bset || !set)
+               goto error;
+       isl_assert(ctx, set->nparam == bset->nparam, goto error);
+       isl_assert(ctx, set->dim == bset->dim, goto error);
+       isl_assert(ctx, set->n < set->size, goto error);
+       set->p[set->n] = bset;
+       set->n++;
+       return set;
+error:
+       if (set)
+               isl_set_free(ctx, set);
+       if (bset)
+               isl_basic_set_free(ctx, bset);
+       return NULL;
+}
+
+void isl_set_free(struct isl_ctx *ctx, struct isl_set *set)
+{
+       int i;
+
+       if (!set)
+               return;
+
+       if (--set->ref > 0)
+               return;
+
+       for (i = 0; i < set->n; ++i)
+               isl_basic_set_free(ctx, set->p[i]);
+       free(set);
+}
+
+void isl_set_dump(struct isl_ctx *ctx, struct isl_set *set, FILE *out,
+                 int indent)
+{
+       int i;
+
+       if (!set) {
+               fprintf(out, "null set\n");
+               return;
+       }
+
+       fprintf(out, "%*s", indent, "");
+       fprintf(out, "ref: %d, n: %d\n", set->ref, set->n);
+       for (i = 0; i < set->n; ++i) {
+               fprintf(out, "%*s", indent, "");
+               fprintf(out, "basic set %d:\n", i);
+               isl_basic_set_dump(ctx, set->p[i], out, indent+4);
+       }
+}
+
+void isl_map_dump(struct isl_ctx *ctx, struct isl_map *map, FILE *out,
+                 int indent)
+{
+       int i;
+
+       fprintf(out, "%*s", indent, "");
+       fprintf(out, "ref: %d, n: %d, nparam: %d, in: %d, out: %d\n",
+                       map->ref, map->n, map->nparam, map->n_in, map->n_out);
+       for (i = 0; i < map->n; ++i) {
+               fprintf(out, "%*s", indent, "");
+               fprintf(out, "basic map %d:\n", i);
+               isl_basic_map_dump(ctx, map->p[i], out, indent+4);
+       }
+}
+
+struct isl_basic_map *isl_basic_map_intersect_domain(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap,
+               struct isl_basic_set *bset)
+{
+       struct isl_basic_map *bmap_domain;
+
+       if (!bmap || !bset)
+               goto error;
+
+       isl_assert(ctx, bset->dim == bmap->n_in, goto error);
+       isl_assert(ctx, bset->nparam == bmap->nparam, goto error);
+
+       bmap = isl_basic_map_extend(ctx, bmap,
+                       bset->nparam, bmap->n_in, bmap->n_out,
+                       bset->n_div, bset->n_eq, bset->n_ineq);
+       if (!bmap)
+               goto error;
+       bmap_domain = isl_basic_map_from_basic_set(ctx, bset,
+                                               bset->dim, 0);
+       bmap = add_constraints(ctx, bmap, bmap_domain, 0, 0);
+
+       return isl_basic_map_finalize(ctx, bmap);
+error:
+       isl_basic_map_free(ctx, bmap);
+       isl_basic_set_free(ctx, bset);
+       return NULL;
+}
+
+struct isl_basic_map *isl_basic_map_intersect_range(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap,
+               struct isl_basic_set *bset)
+{
+       struct isl_basic_map *bmap_range;
+
+       if (!bmap || !bset)
+               goto error;
+
+       isl_assert(ctx, bset->dim == bmap->n_out, goto error);
+       isl_assert(ctx, bset->nparam == bmap->nparam, goto error);
+
+       bmap = isl_basic_map_extend(ctx, bmap,
+                       bset->nparam, bmap->n_in, bmap->n_out,
+                       bset->n_div, bset->n_eq, bset->n_ineq);
+       if (!bmap)
+               goto error;
+       bmap_range = isl_basic_map_from_basic_set(ctx, bset,
+                                               0, bset->dim);
+       bmap = add_constraints(ctx, bmap, bmap_range, 0, 0);
+
+       return isl_basic_map_finalize(ctx, bmap);
+error:
+       isl_basic_map_free(ctx, bmap);
+       isl_basic_set_free(ctx, bset);
+       return NULL;
+}
+
+struct isl_basic_map *isl_basic_map_intersect(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2)
+{
+       isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
+       isl_assert(ctx, bmap1->n_in == bmap2->n_in, goto error);
+       isl_assert(ctx, bmap1->n_out == bmap2->n_out, goto error);
+
+       bmap1 = isl_basic_map_extend(ctx, bmap1,
+                       bmap1->nparam, bmap1->n_in, bmap1->n_out,
+                       bmap2->n_div, bmap2->n_eq, bmap2->n_ineq);
+       if (!bmap1)
+               goto error;
+       bmap1 = add_constraints(ctx, bmap1, bmap2, 0, 0);
+
+       return isl_basic_map_finalize(ctx, bmap1);
+error:
+       isl_basic_map_free(ctx, bmap1);
+       isl_basic_map_free(ctx, bmap2);
+       return NULL;
+}
+
+struct isl_basic_set *isl_basic_set_intersect(
+               struct isl_ctx *ctx, struct isl_basic_set *bset1,
+               struct isl_basic_set *bset2)
+{
+       return (struct isl_basic_set *)
+               isl_basic_map_intersect(ctx,
+                       (struct isl_basic_map *)bset1,
+                       (struct isl_basic_map *)bset2);
+}
+
+struct isl_map *isl_map_intersect(struct isl_ctx *ctx, struct isl_map *map1,
+               struct isl_map *map2)
+{
+       unsigned flags = 0;
+       struct isl_map *result;
+       int i, j;
+
+       if (!map1 || !map2)
+               goto error;
+
+       if (F_ISSET(map1, ISL_MAP_DISJOINT) &&
+           F_ISSET(map2, ISL_MAP_DISJOINT))
+               FL_SET(flags, ISL_MAP_DISJOINT);
+
+       result = isl_map_alloc(ctx, map1->nparam, map1->n_in, map1->n_out,
+                               map1->n * map2->n, flags);
+       if (!result)
+               goto error;
+       for (i = 0; i < map1->n; ++i)
+               for (j = 0; j < map2->n; ++j) {
+                       struct isl_basic_map *part;
+                       part = isl_basic_map_intersect(ctx,
+                                   isl_basic_map_copy(ctx, map1->p[i]),
+                                   isl_basic_map_copy(ctx, map2->p[j]));
+                       if (isl_basic_map_is_empty(ctx, part))
+                               isl_basic_map_free(ctx, part);
+                       else
+                               result = isl_map_add(ctx, result, part);
+                       if (!result)
+                               goto error;
+               }
+       isl_map_free(ctx, map1);
+       isl_map_free(ctx, map2);
+       return result;
+error:
+       isl_map_free(ctx, map1);
+       isl_map_free(ctx, map2);
+       return NULL;
+}
+
+struct isl_basic_map *isl_basic_map_reverse(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       struct isl_basic_set *bset;
+       unsigned in;
+
+       if (!bmap)
+               return NULL;
+       bmap = isl_basic_map_cow(ctx, bmap);
+       if (!bmap)
+               return NULL;
+       in = bmap->n_in;
+       bset = isl_basic_set_from_basic_map(ctx, bmap);
+       bset = isl_basic_set_swap_vars(ctx, bset, in);
+       return isl_basic_map_from_basic_set(ctx, bset, bset->dim-in, in);
+}
+
+/* Turn final n dimensions into existentially quantified variables.
+ */
+struct isl_basic_set *isl_basic_set_project_out(
+               struct isl_ctx *ctx, struct isl_basic_set *bset,
+               unsigned n, unsigned flags)
+{
+       int i;
+       size_t row_size;
+       isl_int **new_div;
+       isl_int *old;
+
+       if (!bset)
+               return NULL;
+
+       isl_assert(ctx, n <= bset->dim, goto error);
+
+       if (n == 0)
+               return bset;
+
+       bset = isl_basic_set_cow(ctx, bset);
+
+       row_size = 1 + bset->nparam + bset->dim + bset->extra;
+       old = bset->block2.data;
+       bset->block2 = isl_blk_extend(ctx, bset->block2,
+                                       (bset->extra + n) * (1 + row_size));
+       if (!bset->block2.data)
+               goto error;
+       new_div = isl_alloc_array(ctx, isl_int *, bset->extra + n);
+       if (!new_div)
+               goto error;
+       for (i = 0; i < n; ++i) {
+               new_div[i] = bset->block2.data +
+                               (bset->extra + i) * (1 + row_size);
+               isl_seq_clr(new_div[i], 1 + row_size);
+       }
+       for (i = 0; i < bset->extra; ++i)
+               new_div[n + i] = bset->block2.data + (bset->div[i] - old);
+       free(bset->div);
+       bset->div = new_div;
+       bset->n_div += n;
+       bset->extra += n;
+       bset->dim -= n;
+       bset = isl_basic_set_simplify(ctx, bset);
+       return isl_basic_set_finalize(ctx, bset);
+error:
+       isl_basic_set_free(ctx, bset);
+       return NULL;
+}
+
+struct isl_basic_map *isl_basic_map_apply_range(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2)
+{
+       struct isl_basic_set *bset;
+       unsigned n_in, n_out;
+
+       if (!bmap1 || !bmap2)
+               goto error;
+
+       isl_assert(ctx, bmap1->n_out == bmap2->n_in, goto error);
+       isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
+
+       n_in = bmap1->n_in;
+       n_out = bmap2->n_out;
+
+       bmap2 = isl_basic_map_reverse(ctx, bmap2);
+       if (!bmap2)
+               goto error;
+       bmap1 = isl_basic_map_extend(ctx, bmap1, bmap1->nparam,
+                       bmap1->n_in + bmap2->n_in, bmap1->n_out,
+                       bmap2->extra, bmap2->n_eq, bmap2->n_ineq);
+       if (!bmap1)
+               goto error;
+       bmap1 = add_constraints(ctx, bmap1, bmap2, bmap1->n_in - bmap2->n_in, 0);
+       bset = isl_basic_set_from_basic_map(ctx, bmap1);
+       bset = isl_basic_set_project_out(ctx, bset,
+                                               bset->dim - (n_in + n_out), 0);
+       return isl_basic_map_from_basic_set(ctx, bset, n_in, n_out);
+error:
+       isl_basic_map_free(ctx, bmap1);
+       isl_basic_map_free(ctx, bmap2);
+       return NULL;
+}
+
+struct isl_basic_map *isl_basic_map_apply_domain(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2)
+{
+       if (!bmap1 || !bmap2)
+               goto error;
+
+       isl_assert(ctx, bmap1->n_in == bmap2->n_in, goto error);
+       isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
+
+       bmap1 = isl_basic_map_reverse(ctx, bmap1);
+       bmap1 = isl_basic_map_apply_range(ctx, bmap1, bmap2);
+       return isl_basic_map_reverse(ctx, bmap1);
+error:
+       isl_basic_map_free(ctx, bmap1);
+       isl_basic_map_free(ctx, bmap2);
+       return NULL;
+}
+
+static struct isl_basic_map *var_equal(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned pos)
+{
+       int i;
+       i = isl_basic_map_alloc_equality(ctx, bmap);
+       if (i < 0)
+               goto error;
+       isl_seq_clr(bmap->eq[i],
+                   1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
+       isl_int_set_si(bmap->eq[i][1+bmap->nparam+pos], -1);
+       isl_int_set_si(bmap->eq[i][1+bmap->nparam+bmap->n_in+pos], 1);
+       return isl_basic_map_finalize(ctx, bmap);
+error:
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+static struct isl_basic_map *var_more(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned pos)
+{
+       int i;
+       i = isl_basic_map_alloc_inequality(ctx, bmap);
+       if (i < 0)
+               goto error;
+       isl_seq_clr(bmap->ineq[i],
+                   1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
+       isl_int_set_si(bmap->ineq[i][0], -1);
+       isl_int_set_si(bmap->ineq[i][1+bmap->nparam+pos], -1);
+       isl_int_set_si(bmap->ineq[i][1+bmap->nparam+bmap->n_in+pos], 1);
+       return isl_basic_map_finalize(ctx, bmap);
+error:
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+static struct isl_basic_map *var_less(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned pos)
+{
+       int i;
+       i = isl_basic_map_alloc_inequality(ctx, bmap);
+       if (i < 0)
+               goto error;
+       isl_seq_clr(bmap->ineq[i],
+                   1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
+       isl_int_set_si(bmap->ineq[i][0], -1);
+       isl_int_set_si(bmap->ineq[i][1+bmap->nparam+pos], 1);
+       isl_int_set_si(bmap->ineq[i][1+bmap->nparam+bmap->n_in+pos], -1);
+       return isl_basic_map_finalize(ctx, bmap);
+error:
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+struct isl_basic_map *isl_basic_map_equal(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out, unsigned n_equal)
+{
+       int i;
+       struct isl_basic_map *bmap;
+       bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, n_equal, 0);
+       if (!bmap)
+               return NULL;
+       for (i = 0; i < n_equal && bmap; ++i)
+               bmap = var_equal(ctx, bmap, i);
+       return isl_basic_map_finalize(ctx, bmap);
+}
+
+struct isl_basic_map *isl_basic_map_less_at(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out, unsigned pos)
+{
+       int i;
+       struct isl_basic_map *bmap;
+       bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, pos, 1);
+       if (!bmap)
+               return NULL;
+       for (i = 0; i < pos && bmap; ++i)
+               bmap = var_equal(ctx, bmap, i);
+       if (bmap)
+               bmap = var_less(ctx, bmap, pos);
+       return isl_basic_map_finalize(ctx, bmap);
+}
+
+struct isl_basic_map *isl_basic_map_more_at(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out, unsigned pos)
+{
+       int i;
+       struct isl_basic_map *bmap;
+       bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, pos, 1);
+       if (!bmap)
+               return NULL;
+       for (i = 0; i < pos && bmap; ++i)
+               bmap = var_equal(ctx, bmap, i);
+       if (bmap)
+               bmap = var_more(ctx, bmap, pos);
+       return isl_basic_map_finalize(ctx, bmap);
+}
+
+struct isl_basic_map *isl_basic_map_from_basic_set(
+               struct isl_ctx *ctx, struct isl_basic_set *bset,
+               unsigned n_in, unsigned n_out)
+{
+       struct isl_basic_map *bmap;
+
+       bset = isl_basic_set_cow(ctx, bset);
+       if (!bset)
+               return NULL;
+
+       isl_assert(ctx, bset->dim == n_in + n_out, goto error);
+       bmap = (struct isl_basic_map *) bset;
+       bmap->n_in = n_in;
+       bmap->n_out = n_out;
+       return isl_basic_map_finalize(ctx, bmap);
+error:
+       isl_basic_set_free(ctx, bset);
+       return NULL;
+}
+
+struct isl_basic_set *isl_basic_set_from_basic_map(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap)
+{
+       if (!bmap)
+               goto error;
+       if (bmap->n_in == 0)
+               return (struct isl_basic_set *)bmap;
+       bmap = isl_basic_map_cow(ctx, bmap);
+       if (!bmap)
+               goto error;
+       bmap->n_out += bmap->n_in;
+       bmap->n_in = 0;
+       bmap = isl_basic_map_finalize(ctx, bmap);
+       return (struct isl_basic_set *)bmap;
+error:
+       return NULL;
+}
+
+struct isl_basic_set *isl_basic_map_domain(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       struct isl_basic_set *domain;
+       unsigned n_out;
+       if (!bmap)
+               return NULL;
+       n_out = bmap->n_out;
+       domain = isl_basic_set_from_basic_map(ctx, bmap);
+       return isl_basic_set_project_out(ctx, domain, n_out, 0);
+}
+
+struct isl_basic_set *isl_basic_map_range(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       return isl_basic_map_domain(ctx,
+                       isl_basic_map_reverse(ctx, bmap));
+}
+
+struct isl_set *isl_map_range(struct isl_ctx *ctx, struct isl_map *map)
+{
+       int i;
+       struct isl_set *set;
+
+       if (!map)
+               goto error;
+       map = isl_map_cow(ctx, map);
+       if (!map)
+               goto error;
+
+       set = (struct isl_set *) map;
+       set->zero = 0;
+       for (i = 0; i < map->n; ++i) {
+               set->p[i] = isl_basic_map_range(ctx, map->p[i]);
+               if (!set->p[i])
+                       goto error;
+       }
+       F_CLR(set, ISL_MAP_DISJOINT);
+       return set;
+error:
+       isl_map_free(ctx, map);
+       return NULL;
+}
+
+struct isl_map *isl_map_from_set(
+               struct isl_ctx *ctx, struct isl_set *set,
+               unsigned n_in, unsigned n_out)
+{
+       int i;
+       struct isl_map *map = NULL;
+
+       if (!set)
+               return NULL;
+       isl_assert(ctx, set->dim == n_in + n_out, goto error);
+       set = isl_set_cow(ctx, set);
+       if (!set)
+               return NULL;
+       map = (struct isl_map *)set;
+       for (i = 0; i < set->n; ++i) {
+               map->p[i] = isl_basic_map_from_basic_set(ctx,
+                               set->p[i], n_in, n_out);
+               if (!map->p[i])
+                       goto error;
+       }
+       map->n_in = n_in;
+       map->n_out = n_out;
+       return map;
+error:
+       isl_set_free(ctx, set);
+       return NULL;
+}
+
+struct isl_map *isl_map_alloc(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out, int n,
+               unsigned flags)
+{
+       struct isl_map *map;
+
+       map = isl_alloc(ctx, struct isl_map,
+                       sizeof(struct isl_map) +
+                       n * sizeof(struct isl_basic_map *));
+       if (!map)
+               return NULL;
+
+       map->ref = 1;
+       map->size = n;
+       map->n = 0;
+       map->nparam = nparam;
+       map->n_in = in;
+       map->n_out = out;
+       map->flags = flags;
+       return map;
+}
+
+struct isl_basic_map *isl_basic_map_empty(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out)
+{
+       struct isl_basic_map *bmap;
+       bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, 1, 0);
+       bmap = isl_basic_map_set_to_empty(ctx, bmap);
+       return bmap;
+}
+
+struct isl_map *isl_map_empty(struct isl_ctx *ctx,
+               unsigned nparam, unsigned in, unsigned out)
+{
+       return isl_map_alloc(ctx, nparam, in, out, 0, ISL_MAP_DISJOINT);
+}
+
+struct isl_set *isl_set_empty(struct isl_ctx *ctx,
+               unsigned nparam, unsigned dim)
+{
+       return isl_set_alloc(ctx, nparam, dim, 0, ISL_MAP_DISJOINT);
+}
+
+struct isl_map *isl_map_dup(struct isl_ctx *ctx, struct isl_map *map)
+{
+       int i;
+       struct isl_map *dup;
+
+       dup = isl_map_alloc(ctx, map->nparam, map->n_in, map->n_out, map->n,
+                               map->flags);
+       for (i = 0; i < map->n; ++i)
+               dup = isl_map_add(ctx, dup,
+                                 isl_basic_map_dup(ctx, map->p[i]));
+       return dup;
+}
+
+struct isl_map *isl_map_add(struct isl_ctx *ctx, struct isl_map *map,
+                               struct isl_basic_map *bmap)
+{
+       if (!bmap || !map)
+               goto error;
+       isl_assert(ctx, map->nparam == bmap->nparam, goto error);
+       isl_assert(ctx, map->n_in == bmap->n_in, goto error);
+       isl_assert(ctx, map->n_out == bmap->n_out, goto error);
+       isl_assert(ctx, map->n < map->size, goto error);
+       map->p[map->n] = bmap;
+       map->n++;
+       return map;
+error:
+       if (map)
+               isl_map_free(ctx, map);
+       if (bmap)
+               isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+void isl_map_free(struct isl_ctx *ctx, struct isl_map *map)
+{
+       int i;
+
+       if (!map)
+               return;
+
+       if (--map->ref > 0)
+               return;
+
+       for (i = 0; i < map->n; ++i)
+               isl_basic_map_free(ctx, map->p[i]);
+       free(map);
+}
+
+struct isl_map *isl_map_extend(struct isl_ctx *ctx, struct isl_map *base,
+               unsigned nparam, unsigned n_in, unsigned n_out)
+{
+       int i;
+
+       base = isl_map_cow(ctx, base);
+       if (!base)
+               return NULL;
+
+       isl_assert(ctx, base->nparam <= nparam, goto error);
+       isl_assert(ctx, base->n_in <= n_in, goto error);
+       isl_assert(ctx, base->n_out <= n_out, goto error);
+       base->nparam = nparam;
+       base->n_in = n_in;
+       base->n_out = n_out;
+       for (i = 0; i < base->n; ++i) {
+               base->p[i] = isl_basic_map_extend(ctx, base->p[i],
+                               nparam, n_in, n_out, 0, 0, 0);
+               if (!base->p[i])
+                       goto error;
+       }
+       return base;
+error:
+       isl_map_free(ctx, base);
+       return NULL;
+}
+
+struct isl_basic_map *isl_basic_map_fix_input_si(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap,
+               unsigned input, int value)
+{
+       int j;
+
+       bmap = isl_basic_map_cow(ctx, bmap);
+       if (!bmap)
+               return NULL;
+       isl_assert(ctx, input < bmap->n_in, goto error);
+
+       bmap = isl_basic_map_extend(ctx, bmap,
+                       bmap->nparam, bmap->n_in, bmap->n_out, 0, 1, 0);
+       j = isl_basic_map_alloc_equality(ctx, bmap);
+       if (j < 0)
+               goto error;
+       isl_int_set_si(bmap->eq[j][1+bmap->nparam+input], -1);
+       isl_int_set_si(bmap->eq[j][0], value);
+       return isl_basic_map_finalize(ctx, bmap);
+error:
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+struct isl_map *isl_map_fix_input_si(struct isl_ctx *ctx, struct isl_map *map,
+               unsigned input, int value)
+{
+       int i;
+
+       map = isl_map_cow(ctx, map);
+       if (!map)
+               return NULL;
+
+       isl_assert(ctx, input < map->n_in, goto error);
+       for (i = 0; i < map->n; ++i) {
+               map->p[i] = isl_basic_map_fix_input_si(ctx, map->p[i],
+                                                               input, value);
+               if (!map->p[i])
+                       goto error;
+       }
+       return map;
+error:
+       isl_map_free(ctx, map);
+       return NULL;
+}
+
+struct isl_map *isl_map_reverse(struct isl_ctx *ctx, struct isl_map *map)
+{
+       int i;
+       unsigned t;
+
+       map = isl_map_cow(ctx, map);
+       if (!map)
+               return NULL;
+
+       t = map->n_in;
+       map->n_in = map->n_out;
+       map->n_out = t;
+       for (i = 0; i < map->n; ++i) {
+               map->p[i] = isl_basic_map_reverse(ctx, map->p[i]);
+               if (!map->p[i])
+                       goto error;
+       }
+       return map;
+error:
+       isl_map_free(ctx, map);
+       return NULL;
+}
+
+struct isl_map *isl_basic_map_lexmax(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, struct isl_basic_set *dom,
+               struct isl_set **empty)
+{
+#ifdef ISL_PIPLIB
+       return pip_isl_basic_map_lexmax(ctx, bmap, dom, empty);
+#else
+       isl_basic_map_free(ctx, bmap);
+       isl_basic_set_free(ctx, dom);
+       return NULL;
+#endif
+}
+
+struct isl_map *isl_basic_map_lexmin(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, struct isl_basic_set *dom,
+               struct isl_set **empty)
+{
+#ifdef ISL_PIPLIB
+       return pip_isl_basic_map_lexmin(ctx, bmap, dom, empty);
+#else
+       isl_basic_map_free(ctx, bmap);
+       isl_basic_set_free(ctx, dom);
+       return NULL;
+#endif
+}
+
+struct isl_set *isl_basic_set_lexmin(struct isl_ctx *ctx,
+               struct isl_basic_set *bset)
+{
+       struct isl_basic_map *bmap = NULL;
+       struct isl_basic_set *dom = NULL;
+       struct isl_map *min;
+
+       if (!bset)
+               goto error;
+       bmap = isl_basic_map_from_basic_set(ctx, bset, 0, bset->dim);
+       if (!bmap)
+               goto error;
+       dom = isl_basic_set_alloc(ctx, bmap->nparam, 0, 0, 0, 0);
+       if (!dom)
+               goto error;
+       min = isl_basic_map_lexmin(ctx, bmap, dom, NULL);
+       return isl_map_range(ctx, min);
+error:
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+struct isl_map *isl_basic_map_compute_divs(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       if (bmap->n_div == 0)
+               return isl_map_from_basic_map(ctx, bmap);
+#ifdef ISL_PIPLIB
+       return pip_isl_basic_map_compute_divs(ctx, bmap);
+#else
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+#endif
+}
+
+struct isl_map *isl_map_compute_divs(struct isl_ctx *ctx, struct isl_map *map)
+{
+       int i;
+       struct isl_map *res;
+
+       if (!map)
+               return NULL;
+       if (map->n == 0)
+               return map;
+       res = isl_basic_map_compute_divs(ctx,
+               isl_basic_map_copy(ctx, map->p[0]));
+       for (i = 1 ; i < map->n; ++i) {
+               struct isl_map *r2;
+               r2 = isl_basic_map_compute_divs(ctx,
+                       isl_basic_map_copy(ctx, map->p[i]));
+               if (F_ISSET(map, ISL_MAP_DISJOINT))
+                       res = isl_map_union_disjoint(ctx, res, r2);
+               else
+                       res = isl_map_union(ctx, res, r2);
+       }
+       isl_map_free(ctx, map);
+
+       return res;
+}
+
+struct isl_set *isl_basic_set_compute_divs(struct isl_ctx *ctx,
+               struct isl_basic_set *bset)
+{
+       return (struct isl_set *)
+               isl_basic_map_compute_divs(ctx,
+                       (struct isl_basic_map *)bset);
+}
+
+struct isl_set *isl_map_domain(struct isl_ctx *ctx, struct isl_map *map)
+{
+       int i;
+       struct isl_set *set;
+
+       if (!map)
+               goto error;
+
+       map = isl_map_cow(ctx, map);
+       if (!map)
+               return NULL;
+
+       set = (struct isl_set *)map;
+       set->dim = map->n_in;
+       set->zero = 0;
+       for (i = 0; i < map->n; ++i) {
+               set->p[i] = isl_basic_map_domain(ctx, map->p[i]);
+               if (!set->p[i])
+                       goto error;
+       }
+       F_CLR(set, ISL_MAP_DISJOINT);
+       return set;
+error:
+       isl_map_free(ctx, map);
+       return NULL;
+}
+
+struct isl_map *isl_map_union_disjoint(struct isl_ctx *ctx,
+                       struct isl_map *map1, struct isl_map *map2)
+{
+       int i;
+       unsigned flags = 0;
+       struct isl_map *map = NULL;
+
+       if (!map1 || !map2)
+               goto error;
+
+       if (map1->n == 0) {
+               isl_map_free(ctx, map1);
+               return map2;
+       }
+       if (map2->n == 0) {
+               isl_map_free(ctx, map2);
+               return map1;
+       }
+
+       isl_assert(ctx, map1->nparam == map2->nparam, goto error);
+       isl_assert(ctx, map1->n_in == map2->n_in, goto error);
+       isl_assert(ctx, map1->n_out == map2->n_out, goto error);
+
+       if (F_ISSET(map1, ISL_MAP_DISJOINT) &&
+           F_ISSET(map2, ISL_MAP_DISJOINT))
+               FL_SET(flags, ISL_MAP_DISJOINT);
+
+       map = isl_map_alloc(ctx, map1->nparam, map1->n_in, map1->n_out,
+                               map1->n + map2->n, flags);
+       if (!map)
+               goto error;
+       for (i = 0; i < map1->n; ++i) {
+               map = isl_map_add(ctx, map,
+                                 isl_basic_map_copy(ctx, map1->p[i]));
+               if (!map)
+                       goto error;
+       }
+       for (i = 0; i < map2->n; ++i) {
+               map = isl_map_add(ctx, map,
+                                 isl_basic_map_copy(ctx, map2->p[i]));
+               if (!map)
+                       goto error;
+       }
+       isl_map_free(ctx, map1);
+       isl_map_free(ctx, map2);
+       return map;
+error:
+       isl_map_free(ctx, map);
+       isl_map_free(ctx, map1);
+       isl_map_free(ctx, map2);
+       return NULL;
+}
+
+struct isl_map *isl_map_union(struct isl_ctx *ctx,
+                       struct isl_map *map1, struct isl_map *map2)
+{
+       map1 = isl_map_union_disjoint(ctx, map1, map2);
+       if (!map1)
+               return NULL;
+       if (map1->n > 1)
+               F_CLR(map1, ISL_MAP_DISJOINT);
+       return map1;
+}
+
+struct isl_set *isl_set_union_disjoint(struct isl_ctx *ctx,
+                       struct isl_set *set1, struct isl_set *set2)
+{
+       return (struct isl_set *)
+               isl_map_union_disjoint(ctx,
+                       (struct isl_map *)set1, (struct isl_map *)set2);
+}
+
+struct isl_set *isl_set_union(struct isl_ctx *ctx,
+                       struct isl_set *set1, struct isl_set *set2)
+{
+       return (struct isl_set *)
+               isl_map_union(ctx,
+                       (struct isl_map *)set1, (struct isl_map *)set2);
+}
+
+struct isl_map *isl_map_intersect_range(
+               struct isl_ctx *ctx, struct isl_map *map,
+               struct isl_set *set)
+{
+       unsigned flags = 0;
+       struct isl_map *result;
+       int i, j;
+
+       if (!map || !set)
+               goto error;
+
+       if (F_ISSET(map, ISL_MAP_DISJOINT) &&
+           F_ISSET(set, ISL_MAP_DISJOINT))
+               FL_SET(flags, ISL_MAP_DISJOINT);
+
+       result = isl_map_alloc(ctx, map->nparam, map->n_in, map->n_out,
+                               map->n * set->n, flags);
+       if (!result)
+               goto error;
+       for (i = 0; i < map->n; ++i)
+               for (j = 0; j < set->n; ++j) {
+                       result = isl_map_add(ctx, result,
+                           isl_basic_map_intersect_range(ctx,
+                               isl_basic_map_copy(ctx, map->p[i]),
+                               isl_basic_set_copy(ctx, set->p[j])));
+                       if (!result)
+                               goto error;
+               }
+       isl_map_free(ctx, map);
+       isl_set_free(ctx, set);
+       return result;
+error:
+       isl_map_free(ctx, map);
+       isl_set_free(ctx, set);
+       return NULL;
+}
+
+struct isl_map *isl_map_intersect_domain(
+               struct isl_ctx *ctx, struct isl_map *map,
+               struct isl_set *set)
+{
+       return isl_map_reverse(ctx,
+               isl_map_intersect_range(ctx, isl_map_reverse(ctx, map), set));
+}
+
+struct isl_map *isl_map_apply_domain(
+               struct isl_ctx *ctx, struct isl_map *map1,
+               struct isl_map *map2)
+{
+       if (!map1 || !map2)
+               goto error;
+       map1 = isl_map_reverse(ctx, map1);
+       map1 = isl_map_apply_range(ctx, map1, map2);
+       return isl_map_reverse(ctx, map1);
+error:
+       isl_map_free(ctx, map1);
+       isl_map_free(ctx, map2);
+       return NULL;
+}
+
+struct isl_map *isl_map_apply_range(
+               struct isl_ctx *ctx, struct isl_map *map1,
+               struct isl_map *map2)
+{
+       struct isl_map *result;
+       int i, j;
+
+       if (!map1 || !map2)
+               goto error;
+
+       isl_assert(ctx, map1->nparam == map2->nparam, goto error);
+       isl_assert(ctx, map1->n_out == map2->n_in, goto error);
+
+       result = isl_map_alloc(ctx, map1->nparam, map1->n_in, map2->n_out,
+                               map1->n * map2->n, 0);
+       if (!result)
+               goto error;
+       for (i = 0; i < map1->n; ++i)
+               for (j = 0; j < map2->n; ++j) {
+                       result = isl_map_add(ctx, result,
+                           isl_basic_map_apply_range(ctx,
+                               isl_basic_map_copy(ctx, map1->p[i]),
+                               isl_basic_map_copy(ctx, map2->p[j])));
+                       if (!result)
+                               goto error;
+               }
+       isl_map_free(ctx, map1);
+       isl_map_free(ctx, map2);
+       return result;
+error:
+       isl_map_free(ctx, map1);
+       isl_map_free(ctx, map2);
+       return NULL;
+}
+
+/*
+ * returns range - domain
+ */
+struct isl_basic_set *isl_basic_map_deltas(struct isl_ctx *ctx,
+                               struct isl_basic_map *bmap)
+{
+       struct isl_basic_set *bset;
+       unsigned dim;
+       int i;
+
+       if (!bmap)
+               goto error;
+       isl_assert(ctx, bmap->n_in == bmap->n_out, goto error);
+       dim = bmap->n_in;
+       bset = isl_basic_set_from_basic_map(ctx, bmap);
+       bset = isl_basic_set_extend(ctx, bset, bset->nparam, 3*dim, 0,
+                                       dim, 0);
+       bset = isl_basic_set_swap_vars(ctx, bset, 2*dim);
+       for (i = 0; i < dim; ++i) {
+               int j = isl_basic_map_alloc_equality(ctx,
+                                           (struct isl_basic_map *)bset);
+               if (j < 0)
+                       goto error;
+               isl_int_set_si(bset->eq[j][1+bset->nparam+i], 1);
+               isl_int_set_si(bset->eq[j][1+bset->nparam+dim+i], 1);
+               isl_int_set_si(bset->eq[j][1+bset->nparam+2*dim+i], -1);
+       }
+       return isl_basic_set_project_out(ctx, bset, 2*dim, 0);
+error:
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+static int div_is_redundant(struct isl_ctx *ctx, struct isl_basic_map *bmap,
+                               int div)
+{
+       int i;
+       unsigned pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
+
+       for (i = 0; i < bmap->n_eq; ++i)
+               if (!isl_int_is_zero(bmap->eq[i][pos]))
+                       return 0;
+
+       for (i = 0; i < bmap->n_ineq; ++i)
+               if (!isl_int_is_zero(bmap->ineq[i][pos]))
+                       return 0;
+
+       for (i = 0; i < bmap->n_div; ++i)
+               if (!isl_int_is_zero(bmap->div[i][1+pos]))
+                       return 0;
+
+       return 1;
+}
+
+/*
+ * Remove divs that don't occur in any of the constraints or other divs.
+ * These can arise when dropping some of the variables in a quast
+ * returned by piplib.
+ */
+static struct isl_basic_map *remove_redundant_divs(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       int i;
+
+       for (i = bmap->n_div-1; i >= 0; --i) {
+               if (!div_is_redundant(ctx, bmap, i))
+                       continue;
+               bmap = isl_basic_map_drop_div(ctx, bmap, i);
+       }
+       return bmap;
+}
+
+struct isl_basic_map *isl_basic_map_finalize(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       bmap = remove_redundant_divs(ctx, bmap);
+       if (!bmap)
+               return NULL;
+       F_SET(bmap, ISL_PRIMITIVE_SET_FINAL);
+       return bmap;
+}
+
+struct isl_basic_set *isl_basic_set_finalize(struct isl_ctx *ctx,
+               struct isl_basic_set *bset)
+{
+       return (struct isl_basic_set *)
+               isl_basic_map_finalize(ctx,
+                       (struct isl_basic_map *)bset);
+}
+
+struct isl_set *isl_set_finalize(struct isl_ctx *ctx, struct isl_set *set)
+{
+       int i;
+
+       if (!set)
+               return NULL;
+       for (i = 0; i < set->n; ++i) {
+               set->p[i] = isl_basic_set_finalize(ctx, set->p[i]);
+               if (!set->p[i])
+                       goto error;
+       }
+       return set;
+error:
+       isl_set_free(ctx, set);
+       return NULL;
+}
+
+struct isl_map *isl_map_finalize(struct isl_ctx *ctx, struct isl_map *map)
+{
+       int i;
+
+       if (!map)
+               return NULL;
+       for (i = 0; i < map->n; ++i) {
+               map->p[i] = isl_basic_map_finalize(ctx, map->p[i]);
+               if (!map->p[i])
+                       goto error;
+       }
+       return map;
+error:
+       isl_map_free(ctx, map);
+       return NULL;
+}
+
+struct isl_basic_map *isl_basic_map_identity(struct isl_ctx *ctx,
+               unsigned nparam, unsigned dim)
+{
+       struct isl_basic_map *bmap;
+       int i;
+
+       bmap = isl_basic_map_alloc(ctx, nparam, dim, dim, 0, dim, 0);
+       if (!bmap)
+               goto error;
+
+       for (i = 0; i < dim; ++i) {
+               int j = isl_basic_map_alloc_equality(ctx, bmap);
+               if (j < 0)
+                       goto error;
+               isl_int_set_si(bmap->eq[j][1+nparam+i], 1);
+               isl_int_set_si(bmap->eq[j][1+nparam+dim+i], -1);
+       }
+       return isl_basic_map_finalize(ctx, bmap);
+error:
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+struct isl_map *isl_map_identity(struct isl_ctx *ctx,
+               unsigned nparam, unsigned dim)
+{
+       struct isl_map *map = isl_map_alloc(ctx, nparam, dim, dim, 1,
+                                               ISL_MAP_DISJOINT);
+       if (!map)
+               goto error;
+       map = isl_map_add(ctx, map,
+               isl_basic_map_identity(ctx, nparam, dim));
+       return map;
+error:
+       isl_map_free(ctx, map);
+       return NULL;
+}
+
+int isl_set_is_equal(struct isl_ctx *ctx,
+               struct isl_set *set1, struct isl_set *set2)
+{
+       return isl_map_is_equal(ctx,
+                       (struct isl_map *)set1, (struct isl_map *)set2);
+}
+
+int isl_set_is_subset(struct isl_ctx *ctx,
+               struct isl_set *set1, struct isl_set *set2)
+{
+       return isl_map_is_subset(ctx,
+                       (struct isl_map *)set1, (struct isl_map *)set2);
+}
+
+int isl_basic_map_is_subset(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2)
+{
+       int is_subset;
+       struct isl_map *map1;
+       struct isl_map *map2;
+
+       if (!bmap1 || !bmap2)
+               return -1;
+
+       map1 = isl_map_from_basic_map(ctx,
+                       isl_basic_map_copy(ctx, bmap1));
+       map2 = isl_map_from_basic_map(ctx,
+                       isl_basic_map_copy(ctx, bmap2));
+
+       is_subset = isl_map_is_subset(ctx, map1, map2);
+
+       isl_map_free(ctx, map1);
+       isl_map_free(ctx, map2);
+
+       return is_subset;
+}
+
+int isl_map_is_empty(struct isl_ctx *ctx, struct isl_map *map)
+{
+       int i;
+       int is_empty;
+
+       if (!map)
+               return -1;
+       for (i = 0; i < map->n; ++i) {
+               is_empty = isl_basic_map_is_empty(ctx, map->p[i]);
+               if (is_empty < 0)
+                       return -1;
+               if (!is_empty)
+                       return 0;
+       }
+       return 1;
+}
+
+int isl_set_is_empty(struct isl_ctx *ctx, struct isl_set *set)
+{
+       return isl_map_is_empty(ctx, (struct isl_map *)set);
+}
+
+int isl_map_is_subset(struct isl_ctx *ctx, struct isl_map *map1,
+               struct isl_map *map2)
+{
+       int i;
+       int is_subset = 0;
+       struct isl_map *diff;
+
+       if (!map1 || !map2)
+               return -1;
+
+       if (isl_map_is_empty(ctx, map1))
+               return 1;
+
+       if (isl_map_is_empty(ctx, map2))
+               return 0;
+
+       diff = isl_map_subtract(ctx, isl_map_copy(ctx, map1),
+                                    isl_map_copy(ctx, map2));
+       if (!diff)
+               return -1;
+
+       is_subset = isl_map_is_empty(ctx, diff);
+       isl_map_free(ctx, diff);
+
+       return is_subset;
+}
+
+int isl_map_is_equal(struct isl_ctx *ctx,
+               struct isl_map *map1, struct isl_map *map2)
+{
+       int is_subset;
+
+       if (!map1 || !map2)
+               return -1;
+       is_subset = isl_map_is_subset(ctx, map1, map2);
+       if (is_subset != 1)
+               return is_subset;
+       is_subset = isl_map_is_subset(ctx, map2, map1);
+       return is_subset;
+}
+
+int isl_basic_map_is_strict_subset(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2)
+{
+       int is_subset;
+
+       if (!bmap1 || !bmap2)
+               return -1;
+       is_subset = isl_basic_map_is_subset(ctx, bmap1, bmap2);
+       if (is_subset != 1)
+               return is_subset;
+       is_subset = isl_basic_map_is_subset(ctx, bmap2, bmap1);
+       if (is_subset == -1)
+               return is_subset;
+       return !is_subset;
+}
+
+int isl_basic_map_is_empty(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       struct isl_basic_set *bset = NULL;
+       struct isl_set *min = NULL;
+       int empty;
+
+       if (!bmap)
+               return -1;
+
+       bset = isl_basic_set_from_basic_map(ctx,
+                       isl_basic_map_copy(ctx, bmap));
+       if (!bset)
+               return -1;
+       min = isl_basic_set_lexmin(ctx, bset);
+       if (!min)
+               return -1;
+       empty = min->n == 0;
+       isl_set_free(ctx, min);
+       return empty;
+}
+
+struct isl_map *isl_basic_map_union(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap1,
+               struct isl_basic_map *bmap2)
+{
+       struct isl_map *map;
+       if (!bmap1 || !bmap2)
+               return NULL;
+
+       isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
+       isl_assert(ctx, bmap1->n_in == bmap2->n_in, goto error);
+       isl_assert(ctx, bmap1->n_out == bmap2->n_out, goto error);
+
+       map = isl_map_alloc(ctx, bmap1->nparam,
+                               bmap1->n_in, bmap1->n_out, 2, 0);
+       if (!map)
+               goto error;
+       map = isl_map_add(ctx, map, bmap1);
+       map = isl_map_add(ctx, map, bmap2);
+       return map;
+error:
+       isl_basic_map_free(ctx, bmap1);
+       isl_basic_map_free(ctx, bmap2);
+       return NULL;
+}
+
+struct isl_set *isl_basic_set_union(
+               struct isl_ctx *ctx, struct isl_basic_set *bset1,
+               struct isl_basic_set *bset2)
+{
+       return (struct isl_set *)isl_basic_map_union(ctx,
+                                           (struct isl_basic_map *)bset1,
+                                           (struct isl_basic_map *)bset2);
+}
+
+/* Order divs such that any div only depends on previous divs */
+static struct isl_basic_map *order_divs(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       int i;
+       unsigned off = bmap->nparam + bmap->n_in + bmap->n_out;
+
+       for (i = 0; i < bmap->n_div; ++i) {
+               int pos;
+               pos = isl_seq_first_non_zero(bmap->div[i]+1+1+off+i,
+                                                           bmap->n_div-i);
+               if (pos == -1)
+                       continue;
+               swap_div(bmap, i, pos);
+               --i;
+       }
+       return bmap;
+}
+
+static int find_div(struct isl_basic_map *dst,
+                       struct isl_basic_map *src, unsigned div)
+{
+       int i;
+
+       unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
+
+       for (i = 0; i < dst->n_div; ++i)
+               if (isl_seq_eq(dst->div[i], src->div[div], 1+1+total) &&
+                   isl_seq_first_non_zero(dst->div[i]+1+1+total,
+                                               dst->n_div - src->n_div) == -1)
+                       return i;
+       return -1;
+}
+
+struct isl_basic_map *isl_basic_map_align_divs(struct isl_ctx *ctx,
+               struct isl_basic_map *dst, struct isl_basic_map *src)
+{
+       int i;
+       unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
+
+       if (src->n_div == 0)
+               return dst;
+
+       src = order_divs(ctx, src);
+       dst = isl_basic_map_extend(ctx, dst, dst->nparam, dst->n_in,
+                       dst->n_out, src->n_div, 0, 0);
+       if (!dst)
+               return NULL;
+       for (i = 0; i < src->n_div; ++i) {
+               int j = find_div(dst, src, i);
+               if (j < 0) {
+                       j = isl_basic_map_alloc_div(ctx, dst);
+                       if (j < 0)
+                               goto error;
+                       isl_seq_cpy(dst->div[j], src->div[i], 1+1+total);
+               }
+               if (j != i)
+                       swap_div(dst, i, j);
+       }
+       return dst;
+error:
+       isl_basic_map_free(ctx, dst);
+       return NULL;
+}
+
+static struct isl_map *add_cut_constraint(struct isl_ctx *ctx,
+               struct isl_map *dst,
+               struct isl_basic_map *src, isl_int *c,
+               unsigned len, unsigned extra, int oppose)
+{
+       struct isl_basic_map *copy = NULL;
+       int is_empty;
+       int k;
+
+       copy = isl_basic_map_copy(ctx, src);
+       copy = isl_basic_map_cow(ctx, copy);
+       if (!copy)
+               goto error;
+       copy = isl_basic_map_extend(ctx, copy,
+               copy->nparam, copy->n_in, copy->n_out, 0, 0, 1);
+       k = isl_basic_map_alloc_inequality(ctx, copy);
+       if (k < 0)
+               goto error;
+       if (oppose)
+               isl_seq_neg(copy->ineq[k], c, len);
+       else
+               isl_seq_cpy(copy->ineq[k], c, len);
+       isl_seq_clr(copy->ineq[k]+len, extra);
+       isl_inequality_negate(ctx, copy, k);
+       is_empty = isl_basic_map_is_empty(ctx, copy);
+       if (is_empty < 0)
+               goto error;
+       if (!is_empty)
+               dst = isl_map_add(ctx, dst, copy);
+       else
+               isl_basic_map_free(ctx, copy);
+       return dst;
+error:
+       isl_basic_map_free(ctx, copy);
+       isl_map_free(ctx, dst);
+       return NULL;
+}
+
+static struct isl_map *subtract(struct isl_ctx *ctx, struct isl_map *map,
+               struct isl_basic_map *bmap)
+{
+       int i, j, k;
+       unsigned flags = 0;
+       struct isl_map *rest = NULL;
+       unsigned max;
+       unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+
+       assert(bmap);
+
+       if (!map)
+               goto error;
+
+       if (F_ISSET(map, ISL_MAP_DISJOINT))
+               FL_SET(flags, ISL_MAP_DISJOINT);
+
+       max = map->n * (2 * bmap->n_eq + bmap->n_ineq);
+       rest = isl_map_alloc(ctx, map->nparam, map->n_in, map->n_out,
+                               max, flags);
+       if (!rest)
+               goto error;
+
+       for (i = 0; i < map->n; ++i) {
+               map->p[i] = isl_basic_map_align_divs(ctx, map->p[i], bmap);
+               if (!map->p[i])
+                       goto error;
+       }
+
+       for (j = 0; j < map->n; ++j)
+               map->p[j] = isl_basic_map_cow(ctx, map->p[j]);
+
+       for (i = 0; i < bmap->n_eq; ++i) {
+               for (j = 0; j < map->n; ++j) {
+                       rest = add_cut_constraint(ctx, rest,
+                               map->p[j], bmap->eq[i],
+                               1+total, map->p[j]->n_div - bmap->n_div, 0);
+                       if (!rest)
+                               goto error;
+
+                       rest = add_cut_constraint(ctx, rest,
+                               map->p[j], bmap->eq[i],
+                               1+total, map->p[j]->n_div - bmap->n_div, 1);
+                       if (!rest)
+                               goto error;
+
+                       map->p[j] = isl_basic_map_extend(ctx, map->p[j],
+                               map->p[j]->nparam, map->p[j]->n_in,
+                               map->p[j]->n_out, 0, 1, 0);
+                       if (!map->p[j])
+                               goto error;
+                       k = isl_basic_map_alloc_equality(ctx, map->p[j]);
+                       if (k < 0)
+                               goto error;
+                       isl_seq_cpy(map->p[j]->eq[k], bmap->eq[i], 1+total);
+                       isl_seq_clr(map->p[j]->eq[k]+1+total,
+                                       map->p[j]->n_div - bmap->n_div);
+               }
+       }
+
+       for (i = 0; i < bmap->n_ineq; ++i) {
+               for (j = 0; j < map->n; ++j) {
+                       rest = add_cut_constraint(ctx, rest,
+                               map->p[j], bmap->ineq[i],
+                               1+total, map->p[j]->n_div - bmap->n_div, 0);
+                       if (!rest)
+                               goto error;
+
+                       map->p[j] = isl_basic_map_extend(ctx, map->p[j],
+                               map->p[j]->nparam, map->p[j]->n_in,
+                               map->p[j]->n_out, 0, 0, 1);
+                       if (!map->p[j])
+                               goto error;
+                       k = isl_basic_map_alloc_inequality(ctx, map->p[j]);
+                       if (k < 0)
+                               goto error;
+                       isl_seq_cpy(map->p[j]->ineq[k], bmap->ineq[i], 1+total);
+                       isl_seq_clr(map->p[j]->ineq[k]+1+total,
+                                       map->p[j]->n_div - bmap->n_div);
+               }
+       }
+
+       isl_map_free(ctx, map);
+       return rest;
+error:
+       isl_map_free(ctx, map);
+       isl_map_free(ctx, rest);
+       return NULL;
+}
+
+struct isl_map *isl_map_subtract(struct isl_ctx *ctx, struct isl_map *map1,
+               struct isl_map *map2)
+{
+       int i;
+       if (!map1 || !map2)
+               goto error;
+
+       isl_assert(ctx, map1->nparam == map2->nparam, goto error);
+       isl_assert(ctx, map1->n_in == map2->n_in, goto error);
+       isl_assert(ctx, map1->n_out == map2->n_out, goto error);
+
+       if (isl_map_is_empty(ctx, map2)) {
+               isl_map_free(ctx, map2);
+               return map1;
+       }
+
+       map1 = isl_map_compute_divs(ctx, map1);
+       map2 = isl_map_compute_divs(ctx, map2);
+       if (!map1 || !map2)
+               goto error;
+
+       for (i = 0; map1 && i < map2->n; ++i)
+               map1 = subtract(ctx, map1, map2->p[i]);
+
+       isl_map_free(ctx, map2);
+       return map1;
+error:
+       isl_map_free(ctx, map1);
+       isl_map_free(ctx, map2);
+       return NULL;
+}
+
+struct isl_set *isl_set_subtract(struct isl_ctx *ctx, struct isl_set *set1,
+               struct isl_set *set2)
+{
+       return (struct isl_set *)
+               isl_map_subtract(ctx,
+                       (struct isl_map *)set1, (struct isl_map *)set2);
+}
+
+struct isl_set *isl_set_apply(struct isl_ctx *ctx, struct isl_set *set,
+               struct isl_map *map)
+{
+       isl_assert(ctx, set->dim == map->n_in, goto error);
+       map = isl_map_intersect_domain(ctx, map, set);
+       set = isl_map_range(ctx, map);
+       return set;
+error:
+       isl_set_free(ctx, set);
+       isl_map_free(ctx, map);
+       return NULL;
+}
diff --git a/isl_map_affine_hull.c b/isl_map_affine_hull.c
new file mode 100644 (file)
index 0000000..aec15f6
--- /dev/null
@@ -0,0 +1,271 @@
+#include "isl_ctx.h"
+#include "isl_seq.h"
+#include "isl_set.h"
+#include "isl_lp.h"
+#include "isl_map.h"
+#include "isl_map_private.h"
+
+static void inequality_to_equality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned pos)
+{
+       isl_int *t;
+
+       t = bmap->ineq[pos];
+       bmap->ineq[pos] = bmap->ineq[0];
+       bmap->ineq[0] = bmap->eq[bmap->n_eq];
+       bmap->eq[bmap->n_eq] = t;
+       bmap->n_eq++;
+       bmap->n_ineq--;
+       bmap->ineq++;
+}
+
+struct isl_basic_map *isl_basic_map_affine_hull(struct isl_ctx *ctx,
+                                               struct isl_basic_map *bmap)
+{
+       int i;
+       isl_int opt;
+
+       bmap = isl_basic_map_cow(ctx, bmap);
+       if (!bmap)
+               return NULL;
+
+       isl_int_init(opt);
+       for (i = 0; i < bmap->n_ineq; ++i) {
+               enum isl_lp_result res;
+               res = isl_solve_lp(bmap, 1, bmap->ineq[i]+1, ctx->one, &opt);
+               if (res == isl_lp_unbounded)
+                       continue;
+               if (res == isl_lp_error)
+                       goto error;
+               if (res == isl_lp_empty) {
+                       bmap = isl_basic_map_set_to_empty(ctx, bmap);
+                       break;
+               }
+               isl_int_add(opt, opt, bmap->ineq[i][0]);
+               if (isl_int_is_zero(opt)) {
+                       inequality_to_equality(ctx, bmap, i);
+                       --i;
+               }
+       }
+       isl_basic_map_free_inequality(ctx, bmap, bmap->n_ineq);
+       isl_int_clear(opt);
+       return isl_basic_map_finalize(ctx, bmap);
+error:
+       isl_int_clear(opt);
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+struct isl_basic_set *isl_basic_set_affine_hull(struct isl_ctx *ctx,
+                                               struct isl_basic_set *bset)
+{
+       return (struct isl_basic_set *)
+               isl_basic_map_affine_hull(ctx, (struct isl_basic_map *)bset);
+}
+
+/* Make eq[row][col] of both bmaps equal so we can add the row
+ * add the column to the common matrix.
+ * Note that because of the echelon form, the columns of row row
+ * after column col are zero.
+ */
+static void set_common_multiple(
+       struct isl_basic_map *bmap1, struct isl_basic_map *bmap2,
+       unsigned row, unsigned col)
+{
+       isl_int m, c;
+
+       if (isl_int_eq(bmap1->eq[row][col], bmap2->eq[row][col]))
+               return;
+
+       isl_int_init(c);
+       isl_int_init(m);
+       isl_int_lcm(m, bmap1->eq[row][col], bmap2->eq[row][col]);
+       isl_int_divexact(c, m, bmap1->eq[row][col]);
+       isl_seq_scale(bmap1->eq[row], bmap1->eq[row], c, col+1);
+       isl_int_divexact(c, m, bmap2->eq[row][col]);
+       isl_seq_scale(bmap2->eq[row], bmap2->eq[row], c, col+1);
+       isl_int_clear(c);
+       isl_int_clear(m);
+}
+
+/* Delete a given equality, moving all the following equalities one up.
+ */
+static void delete_row(struct isl_basic_map *bmap, unsigned row)
+{
+       isl_int *t;
+       int r;
+
+       t = bmap->eq[row];
+       bmap->n_eq--;
+       for (r = row; r < bmap->n_eq; ++r)
+               bmap->eq[r] = bmap->eq[r+1];
+       bmap->eq[bmap->n_eq] = t;
+}
+
+/* Make first row entries in column col of bmap1 identical to
+ * those of bmap2, using the fact that entry bmap1->eq[row][col]=a
+ * is non-zero.  Initially, these elements of bmap1 are all zero.
+ * For each row i < row, we set
+ *             A[i] = a * A[i] + B[i][col] * A[row]
+ *             B[i] = a * B[i]
+ * so that
+ *             A[i][col] = B[i][col] = a * old(B[i][col])
+ */
+static void construct_column(
+       struct isl_basic_map *bmap1, struct isl_basic_map *bmap2,
+       unsigned row, unsigned col)
+{
+       int r;
+       isl_int a;
+       isl_int b;
+       unsigned total;
+
+       isl_int_init(a);
+       isl_int_init(b);
+       total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
+       for (r = 0; r < row; ++r) {
+               if (isl_int_is_zero(bmap2->eq[r][col]))
+                       continue;
+               isl_int_gcd(b, bmap2->eq[r][col], bmap1->eq[row][col]);
+               isl_int_divexact(a, bmap1->eq[row][col], b);
+               isl_int_divexact(b, bmap2->eq[r][col], b);
+               isl_seq_combine(bmap1->eq[r], a, bmap1->eq[r],
+                                             b, bmap1->eq[row], total);
+               isl_seq_scale(bmap2->eq[r], bmap2->eq[r], a, total);
+       }
+       isl_int_clear(a);
+       isl_int_clear(b);
+       delete_row(bmap1, row);
+}
+
+/* Make first row entries in column col of bmap1 identical to
+ * those of bmap2, using only these entries of the two matrices.
+ * Let t be the last row with different entries.
+ * For each row i < t, we set
+ *     A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
+ *     B[i] = (A[t][col]-B[t][col]) * B[i] + (B[i][col]-A[i][col) * B[t]
+ * so that
+ *     A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
+ */
+static int transform_column(
+       struct isl_basic_map *bmap1, struct isl_basic_map *bmap2,
+       unsigned row, unsigned col)
+{
+       int i, t;
+       isl_int a, b, g;
+       unsigned total;
+
+       for (t = row-1; t >= 0; --t)
+               if (isl_int_ne(bmap1->eq[t][col], bmap2->eq[t][col]))
+                       break;
+       if (t < 0)
+               return 0;
+
+       total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
+       isl_int_init(a);
+       isl_int_init(b);
+       isl_int_init(g);
+       isl_int_sub(b, bmap1->eq[t][col], bmap2->eq[t][col]);
+       for (i = 0; i < t; ++i) {
+               isl_int_sub(a, bmap2->eq[t][col], bmap1->eq[t][col]);
+               isl_int_gcd(g, a, b);
+               isl_int_divexact(a, a, g);
+               isl_int_divexact(g, b, g);
+               isl_seq_combine(bmap1->eq[i], g, bmap1->eq[i], a, bmap1->eq[t],
+                               total);
+               isl_seq_combine(bmap2->eq[i], g, bmap2->eq[i], a, bmap2->eq[t],
+                               total);
+       }
+       isl_int_clear(a);
+       isl_int_clear(b);
+       isl_int_clear(g);
+       delete_row(bmap1, t);
+       delete_row(bmap2, t);
+       return 1;
+}
+
+/* The implementation is based on Section 5.2 of Michael Karr,
+ * "Affine Relationships Among Variables of a Program",
+ * except that the echelon form we use starts from the last column
+ * and that we are dealing with integer coefficients.
+ */
+static struct isl_basic_map *affine_hull(struct isl_ctx *ctx,
+       struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
+{
+       unsigned total;
+       int col;
+       int row;
+
+       total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
+
+       row = 0;
+       for (col = total-1; col >= 0; --col) {
+               int is_zero1 = row >= bmap1->n_eq ||
+                       isl_int_is_zero(bmap1->eq[row][col]);
+               int is_zero2 = row >= bmap2->n_eq ||
+                       isl_int_is_zero(bmap2->eq[row][col]);
+               if (!is_zero1 && !is_zero2) {
+                       set_common_multiple(bmap1, bmap2, row, col);
+                       ++row;
+               } else if (!is_zero1 && is_zero2) {
+                       construct_column(bmap1, bmap2, row, col);
+               } else if (is_zero1 && !is_zero2) {
+                       construct_column(bmap2, bmap1, row, col);
+               } else {
+                       if (transform_column(bmap1, bmap2, row, col))
+                               --row;
+               }
+       }
+       isl_basic_map_free(ctx, bmap2);
+       return bmap1;
+}
+
+struct isl_basic_map *isl_map_affine_hull(struct isl_ctx *ctx,
+                                               struct isl_map *map)
+{
+       int i;
+       struct isl_basic_map *bmap;
+
+       map = isl_map_compute_divs(ctx, map);
+       map = isl_map_cow(ctx, map);
+       if (!map)
+               return NULL;
+
+       if (map->n == 0) {
+               bmap = isl_basic_map_empty(ctx,
+                                           map->nparam, map->n_in, map->n_out);
+               isl_map_free(ctx, map);
+               return bmap;
+       }
+
+       for (i = 1; i < map->n; ++i)
+               map->p[0] = isl_basic_map_align_divs(ctx, map->p[0], map->p[i]);
+       for (i = 1; i < map->n; ++i)
+               map->p[i] = isl_basic_map_align_divs(ctx, map->p[i], map->p[0]);
+
+       for (i = 0; i < map->n; ++i) {
+               map->p[i] = isl_basic_map_affine_hull(ctx, map->p[i]);
+               map->p[i] = isl_basic_map_gauss(ctx, map->p[i], NULL);
+               if (!map->p[i])
+                       goto error;
+       }
+       while (map->n > 1) {
+               map->p[0] = affine_hull(ctx, map->p[0], map->p[--map->n]);
+               if (!map->p[0])
+                       goto error;
+       }
+       bmap = isl_basic_map_copy(ctx, map->p[0]);
+       isl_map_free(ctx, map);
+       bmap = isl_basic_map_finalize(ctx, bmap);
+       return isl_basic_map_simplify(ctx, bmap);
+error:
+       isl_map_free(ctx, map);
+       return NULL;
+}
+
+struct isl_basic_set *isl_set_affine_hull(struct isl_ctx *ctx,
+                                               struct isl_set *set)
+{
+       return (struct isl_basic_set *)
+               isl_map_affine_hull(ctx, (struct isl_map *)set);
+}
diff --git a/isl_map_piplib.c b/isl_map_piplib.c
new file mode 100644 (file)
index 0000000..8dad265
--- /dev/null
@@ -0,0 +1,490 @@
+#include <piplib/piplibMP.h>
+#include "isl_set.h"
+#include "isl_map.h"
+#include "isl_map_private.h"
+
+static void copy_values_from(isl_int *dst, Entier *src, unsigned n)
+{
+       int i;
+
+       for (i = 0; i < n; ++i)
+               entier_assign(dst[i], src[i]);
+}
+
+static void add_value(isl_int *dst, Entier *src)
+{
+       mpz_add(*dst, *dst, *src);
+}
+
+static void copy_constraint_from(isl_int *dst, PipVector *src,
+               unsigned nparam, unsigned n_in, unsigned n_out,
+               unsigned extra, int *pos)
+{
+       int i;
+
+       copy_values_from(dst, src->the_vector+src->nb_elements-1, 1);
+       copy_values_from(dst+1, src->the_vector, nparam+n_in);
+       isl_seq_clr(dst+1+nparam+n_in, n_out);
+       isl_seq_clr(dst+1+nparam+n_in+n_out, extra);
+       for (i = 0; i + n_in + nparam < src->nb_elements-1; ++i) {
+               int p = pos[i];
+               add_value(&dst[1+nparam+n_in+n_out+p],
+                         &src->the_vector[n_in+nparam+i]);
+       }
+}
+
+static int add_inequality(struct isl_ctx *ctx,
+                  struct isl_basic_map *bmap, int *pos, PipVector *vec)
+{
+       int i = isl_basic_map_alloc_inequality(ctx, bmap);
+       if (i < 0)
+               return -1;
+       copy_constraint_from(bmap->ineq[i], vec,
+           bmap->nparam, bmap->n_in, bmap->n_out, bmap->extra, pos);
+
+       return i;
+}
+
+static int add_equality(struct isl_ctx *ctx,
+                  struct isl_basic_map *bmap, int *pos,
+                  unsigned var, PipVector *vec)
+{
+       int i;
+
+       isl_assert(ctx, var < bmap->n_out, return -1);
+
+       i = isl_basic_map_alloc_equality(ctx, bmap);
+       if (i < 0)
+               return -1;
+       copy_constraint_from(bmap->eq[i], vec,
+           bmap->nparam, bmap->n_in, bmap->n_out, bmap->extra, pos);
+       isl_int_set_si(bmap->eq[i][1+bmap->nparam+bmap->n_in+var], -1);
+
+       return i;
+}
+
+static int find_div(struct isl_ctx *ctx,
+                  struct isl_basic_map *bmap, int *pos, PipNewparm *p)
+{
+       int i, j;
+
+       i = isl_basic_map_alloc_div(ctx, bmap);
+       assert(i != -1);
+
+       copy_constraint_from(bmap->div[i]+1, p->vector,
+           bmap->nparam, bmap->n_in, bmap->n_out, bmap->extra, pos);
+
+       copy_values_from(bmap->div[i], &p->deno, 1);
+       for (j = 0; j < i; ++j)
+               if (isl_seq_eq(bmap->div[i], bmap->div[j],
+                               1+1+bmap->nparam+bmap->n_in+bmap->n_out+j)) {
+                       isl_basic_map_free_div(ctx, bmap, 1);
+                       return j;
+               }
+
+       return i;
+}
+
+/* Count some properties of a quast
+ * - maximal number of new parameters
+ * - maximal depth
+ * - total number of solutions
+ * - total number of empty branches
+ */
+static void quast_count(PipQuast *q, int *maxnew, int depth, int *maxdepth,
+                       int *sol, int *nosol)
+{
+       PipNewparm *p;
+
+       for (p = q->newparm; p; p = p->next)
+               if (p->rank > *maxnew)
+                       *maxnew = p->rank;
+       if (q->condition) {
+               if (++depth > *maxdepth)
+                       *maxdepth = depth;
+               quast_count(q->next_else, maxnew, depth, maxdepth, sol, nosol);
+               quast_count(q->next_then, maxnew, depth, maxdepth, sol, nosol);
+       } else {
+               if (q->list)
+                       ++(*sol);
+               else
+                       ++(*nosol);
+       }
+}
+
+/*
+ * pos: array of length bmap->set.extra, mapping each of the existential
+ *             variables PIP proposes to an existential variable in bmap
+ * bmap: collects the currently active constraints
+ * rest: collects the empty leaves of the quast (if not NULL)
+ */
+struct scan_data {
+       struct isl_ctx                  *ctx;
+       struct isl_basic_map            *bmap;
+       struct isl_set                  **rest;
+       int        *pos;
+};
+
+/*
+ * New existentially quantified variables are places after the existing ones.
+ */
+static struct isl_map *scan_quast_r(struct scan_data *data, PipQuast *q,
+                                   struct isl_map *map)
+{
+       PipNewparm *p;
+       int error;
+       struct isl_basic_map *bmap = data->bmap;
+       unsigned old_n_div = bmap->n_div;
+
+       if (!map)
+               goto error;
+
+       for (p = q->newparm; p; p = p->next) {
+               int pos;
+               PipVector *vec = p->vector;
+               unsigned pip_param = bmap->nparam + bmap->n_in;
+
+               pos = find_div(data->ctx, bmap, data->pos, p);
+               if (pos < 0)
+                       goto error;
+               data->pos[p->rank - pip_param] = pos;
+       }
+
+       if (q->condition) {
+               int j;
+               int pos = add_inequality(data->ctx, bmap, data->pos,
+                                        q->condition);
+               if (pos < 0)
+                       goto error;
+               map = scan_quast_r(data, q->next_then, map);
+
+               if (isl_inequality_negate(data->ctx, bmap, pos))
+                       goto error;
+               map = scan_quast_r(data, q->next_else, map);
+
+               if (isl_basic_map_free_inequality(data->ctx, bmap, 1))
+                       goto error;
+       } else if (q->list) {
+               PipList *l;
+               int j;
+               /* if bmap->n_out is zero, we are only interested in the domains
+                * where a solution exists and not in the actual solution
+                */
+               for (j = 0, l = q->list; j < bmap->n_out && l; ++j, l = l->next)
+                       if (add_equality(data->ctx, bmap, data->pos, j,
+                                               l->vector) < 0)
+                               goto error;
+               map = isl_map_add(data->ctx, map,
+                               isl_basic_map_copy(data->ctx, bmap));
+               if (isl_basic_map_free_equality(data->ctx, bmap,
+                                                   bmap->n_out))
+                       goto error;
+       } else if (map->n && data->rest) {
+               /* not interested in rest if no sol */
+               struct isl_basic_set *bset;
+               bset = isl_basic_set_from_basic_map(data->ctx,
+                               isl_basic_map_copy(data->ctx, bmap));
+               bset = isl_basic_set_drop_vars(data->ctx, bset,
+                               bmap->n_in, bmap->n_out);
+               if (!bset)
+                       goto error;
+               *data->rest = isl_set_add(data->ctx, *data->rest, bset);
+       }
+
+       if (isl_basic_map_free_div(data->ctx, bmap,
+                                          bmap->n_div - old_n_div))
+               goto error;
+       return map;
+error:
+       isl_map_free(data->ctx, map);
+       return NULL;
+}
+
+/*
+ * Returns a map with "context" as domain and as range the first
+ * "keep" variables in the quast lists.
+ */
+struct isl_map *isl_map_from_quast(struct isl_ctx *ctx, PipQuast *q,
+               unsigned keep,
+               struct isl_basic_set *context,
+               struct isl_set **rest)
+{
+       int             pip_param;
+       int             nexist;
+       int             max_depth;
+       int             n_sol, n_nosol;
+       struct scan_data        data;
+       struct isl_map          *map;
+       unsigned                nparam;
+
+       data.ctx = ctx;
+       data.rest = rest;
+       data.bmap = NULL;
+       data.pos = NULL;
+
+       if (!context)
+               goto error;
+
+       nparam = context->nparam;
+       pip_param = nparam + context->dim;
+
+       max_depth = 0;
+       n_sol = 0;
+       n_nosol = 0;
+       nexist = pip_param-1;
+       quast_count(q, &nexist, 0, &max_depth, &n_sol, &n_nosol);
+       nexist -= pip_param-1;
+
+       if (rest) {
+               *rest = isl_set_alloc(data.ctx, nparam, context->dim, n_nosol,
+                                       ISL_MAP_DISJOINT);
+               if (!*rest)
+                       goto error;
+       }
+       map = isl_map_alloc(data.ctx, nparam, context->dim, keep, n_sol,
+                               ISL_MAP_DISJOINT);
+       if (!map)
+               goto error;
+
+       data.bmap = isl_basic_map_from_basic_set(ctx, context,
+                               context->dim, 0);
+       data.bmap = isl_basic_map_extend(data.ctx, data.bmap,
+                       nparam, data.bmap->n_in, keep, nexist, keep, max_depth);
+       if (!data.bmap)
+               goto error;
+
+       if (data.bmap->extra) {
+               int i;
+               data.pos = isl_alloc_array(ctx, int, data.bmap->extra);
+               if (!data.pos)
+                       goto error;
+               for (i = 0; i < data.bmap->n_div; ++i)
+                       data.pos[i] = i;
+       }
+
+       map = scan_quast_r(&data, q, map);
+       map = isl_map_finalize(ctx, map);
+       if (!map)
+               goto error;
+       if (rest) {
+               *rest = isl_set_finalize(ctx, *rest);
+               if (!*rest)
+                       goto error;
+       }
+       isl_basic_map_free(data.ctx, data.bmap);
+       if (data.pos)
+               free(data.pos);
+       return map;
+error:
+       if (data.pos)
+               free(data.pos);
+       isl_basic_map_free(data.ctx, data.bmap);
+       isl_map_free(ctx, map);
+       if (rest) {
+               isl_set_free(ctx, *rest);
+               *rest = NULL;
+       }
+       return NULL;
+}
+
+static void copy_values_to(Entier *dst, isl_int *src, unsigned n)
+{
+       int i;
+
+       for (i = 0; i < n; ++i)
+               entier_assign(dst[i], src[i]);
+}
+
+static void copy_constraint_to(Entier *dst, isl_int *src,
+               unsigned pip_param, unsigned pip_var,
+               unsigned extra_front, unsigned extra_back)
+{
+       copy_values_to(dst+1+extra_front+pip_var+pip_param+extra_back, src, 1);
+       copy_values_to(dst+1+extra_front+pip_var, src+1, pip_param);
+       copy_values_to(dst+1+extra_front, src+1+pip_param, pip_var);
+}
+
+PipMatrix *isl_basic_map_to_pip(struct isl_basic_map *bmap, unsigned pip_param,
+                        unsigned extra_front, unsigned extra_back)
+{
+       int i;
+       unsigned nrow;
+       unsigned ncol;
+       PipMatrix *M;
+       unsigned off;
+       unsigned pip_var = bmap->nparam + bmap->n_in + bmap->n_out
+                               + bmap->n_div - pip_param;
+       unsigned pip_dim = pip_var - bmap->n_div;
+
+       nrow = extra_front + bmap->n_eq + bmap->n_ineq + 2*bmap->n_div;
+       ncol = 1 + extra_front + pip_var + pip_param + extra_back + 1;
+       M = pip_matrix_alloc(nrow, ncol);
+       if (!M)
+               return NULL;
+
+       off = extra_front;
+       for (i = 0; i < bmap->n_eq; ++i) {
+               entier_set_si(M->p[off+i][0], 0);
+               copy_constraint_to(M->p[off+i], bmap->eq[i],
+                                  pip_param, pip_var, extra_front, extra_back);
+       }
+       off += bmap->n_eq;
+       for (i = 0; i < bmap->n_ineq; ++i) {
+               entier_set_si(M->p[off+i][0], 1);
+               copy_constraint_to(M->p[off+i], bmap->ineq[i],
+                                  pip_param, pip_var, extra_front, extra_back);
+       }
+       off += bmap->n_ineq;
+       for (i = 0; i < bmap->n_div; ++i) {
+               unsigned total = bmap->n_in + bmap->n_out
+                                   + bmap->n_div + bmap->nparam + extra_back;
+               if (isl_int_is_zero(bmap->div[i][0]))
+                       continue;
+               entier_set_si(M->p[off+2*i][0], 1);
+               copy_constraint_to(M->p[off+2*i], bmap->div[i]+1,
+                                  pip_param, pip_var, extra_front, extra_back);
+               copy_values_to(M->p[off+2*i]+1+extra_front+pip_dim+i,
+                               bmap->div[i], 1);
+               entier_oppose(M->p[off+2*i][1+extra_front+pip_dim+i],
+                             M->p[off+2*i][1+extra_front+pip_dim+i]);
+
+               entier_set_si(M->p[off+2*i+1][0], 1);
+               Vector_Oppose(M->p[off+2*i]+1+extra_front,
+                             M->p[off+2*i+1]+1+extra_front, total+1);
+               entier_addto(M->p[off+2*i+1][1+extra_front+total],
+                            M->p[off+2*i+1][1+extra_front+total],
+                            M->p[off+2*i+1][1+extra_front+pip_dim+i]);
+               entier_decrement(M->p[off+2*i+1][1+extra_front+total],
+                                M->p[off+2*i+1][1+extra_front+total]);
+       }
+       return M;
+}
+
+static struct isl_map *extremum_on(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, struct isl_basic_set *dom,
+               struct isl_set **empty, int max)
+{
+       PipOptions      *options;
+       PipQuast        *sol;
+       struct isl_map  *map;
+       PipMatrix *domain = NULL, *context = NULL;
+
+       isl_assert(ctx, bmap->nparam == dom->nparam, goto error);
+       isl_assert(ctx, bmap->n_in == dom->dim, goto error);
+
+       domain = isl_basic_map_to_pip(bmap, bmap->nparam + bmap->n_in,
+                                       0, dom->n_div);
+       if (!domain)
+               goto error;
+       context = isl_basic_map_to_pip((struct isl_basic_map *)dom, 0, 0, 0);
+       if (!context)
+               goto error;
+
+       options = pip_options_init();
+       options->Simplify = 1;
+       options->Maximize = max;
+       options->Urs_unknowns = -1;
+       options->Urs_parms = -1;
+       sol = pip_solve(domain, context, -1, options);
+
+       if (sol) {
+               struct isl_basic_set *copy;
+               copy = isl_basic_set_copy(ctx, dom);
+               map = isl_map_from_quast(ctx, sol, bmap->n_out, copy, empty);
+       } else {
+               map = isl_map_empty(ctx, bmap->nparam, bmap->n_in, bmap->n_out);
+               if (empty)
+                       *empty = NULL;
+       }
+       if (!map)
+               goto error;
+       if (map->n == 0 && empty) {
+               isl_set_free(ctx, *empty);
+               *empty = isl_set_from_basic_set(ctx, dom);
+       } else
+               isl_basic_set_free(ctx, dom);
+       isl_basic_map_free(ctx, bmap);
+
+       pip_quast_free(sol);
+       pip_options_free(options);
+       pip_matrix_free(domain);
+       pip_matrix_free(context);
+
+       return map;
+error:
+       if (domain)
+               pip_matrix_free(domain);
+       if (context)
+               pip_matrix_free(context);
+       isl_basic_map_free(ctx, bmap);
+       isl_basic_set_free(ctx, dom);
+       return NULL;
+}
+
+struct isl_map *pip_isl_basic_map_lexmax(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, struct isl_basic_set *dom,
+               struct isl_set **empty)
+{
+       return extremum_on(ctx, bmap, dom, empty, 1);
+}
+
+struct isl_map *pip_isl_basic_map_lexmin(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, struct isl_basic_set *dom,
+               struct isl_set **empty)
+{
+       return extremum_on(ctx, bmap, dom, empty, 0);
+}
+
+struct isl_map *pip_isl_basic_map_compute_divs(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap)
+{
+       PipMatrix *domain = NULL, *context = NULL;
+       PipOptions      *options;
+       PipQuast        *sol;
+       struct isl_map  *map;
+       struct isl_set  *set;
+       struct isl_basic_set    *dom;
+       unsigned         n_in;
+       unsigned         n_out;
+
+       if (!bmap)
+               goto error;
+
+       n_in = bmap->n_in;
+       n_out = bmap->n_out;
+
+       domain = isl_basic_map_to_pip(bmap, bmap->nparam + n_in + n_out, 0, 0);
+       if (!domain)
+               goto error;
+       context = pip_matrix_alloc(0, bmap->nparam + n_in + n_out + 2);
+       if (!context)
+               goto error;
+
+       options = pip_options_init();
+       options->Simplify = 1;
+       options->Urs_unknowns = -1;
+       options->Urs_parms = -1;
+       sol = pip_solve(domain, context, -1, options);
+
+       dom = isl_basic_set_alloc(ctx, bmap->nparam, n_in + n_out, 0, 0, 0);
+       map = isl_map_from_quast(ctx, sol, 0, dom, NULL);
+
+       pip_quast_free(sol);
+       pip_options_free(options);
+       pip_matrix_free(domain);
+       pip_matrix_free(context);
+
+       isl_basic_map_free(ctx, bmap);
+
+       set = isl_map_domain(ctx, map);
+
+       return isl_map_from_set(ctx, set, n_in, n_out);
+error:
+       if (domain)
+               pip_matrix_free(domain);
+       if (context)
+               pip_matrix_free(context);
+       isl_basic_set_free(ctx, dom);
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
diff --git a/isl_map_polylib.c b/isl_map_polylib.c
new file mode 100644 (file)
index 0000000..58848ac
--- /dev/null
@@ -0,0 +1,229 @@
+#include "isl_set.h"
+#include "isl_map.h"
+#include "isl_map_polylib.h"
+
+static void copy_values_from(isl_int *dst, Value *src, unsigned n)
+{
+       int i;
+
+       for (i = 0; i < n; ++i)
+               value_assign(dst[i], src[i]);
+}
+
+static void copy_values_to(Value *dst, isl_int *src, unsigned n)
+{
+       int i;
+
+       for (i = 0; i < n; ++i)
+               value_assign(dst[i], src[i]);
+}
+
+static void copy_constraint_from(isl_int *dst, Value *src,
+               unsigned nparam, unsigned dim, unsigned extra)
+{
+       copy_values_from(dst, src+1+dim+extra+nparam, 1);
+       copy_values_from(dst+1, src+1+dim+extra, nparam);
+       copy_values_from(dst+1+nparam, src+1, dim);
+       copy_values_from(dst+1+nparam+dim, src+1+dim, extra);
+}
+
+static void copy_constraint_to(Value *dst, isl_int *src,
+               unsigned nparam, unsigned dim, unsigned extra)
+{
+       copy_values_to(dst+1+dim+extra+nparam, src, 1);
+       copy_values_to(dst+1+dim+extra, src+1, nparam);
+       copy_values_to(dst+1, src+1+nparam, dim);
+       copy_values_to(dst+1+dim, src+1+nparam+dim, extra);
+}
+
+static int add_equality(struct isl_ctx *ctx, struct isl_basic_map *bmap,
+                        Value *constraint)
+{
+       int i = isl_basic_map_alloc_equality(ctx, bmap);
+       if (i < 0)
+               return -1;
+       copy_constraint_from(bmap->eq[i], constraint, bmap->nparam,
+                            bmap->n_in + bmap->n_out, bmap->extra);
+       return 0;
+}
+
+static int add_inequality(struct isl_ctx *ctx, struct isl_basic_map *bmap,
+                        Value *constraint)
+{
+       int i = isl_basic_map_alloc_inequality(ctx, bmap);
+       if (i < 0)
+               return -1;
+       copy_constraint_from(bmap->ineq[i], constraint, bmap->nparam,
+                            bmap->n_in + bmap->n_out, bmap->extra);
+       return 0;
+}
+
+static struct isl_basic_map *copy_constraints(
+                       struct isl_ctx *ctx, struct isl_basic_map *bmap,
+                       Polyhedron *P)
+{
+       int i;
+       unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra;
+
+       for (i = 0; i < P->NbConstraints; ++i) {
+               if (value_zero_p(P->Constraint[i][0])) {
+                       if (add_equality(ctx, bmap, P->Constraint[i]))
+                               goto error;
+               } else {
+                       if (add_inequality(ctx, bmap, P->Constraint[i]))
+                               goto error;
+               }
+       }
+       for (i = 0; i < bmap->extra; ++i) {
+               int j = isl_basic_map_alloc_div(ctx, bmap);
+               if (j == -1)
+                       goto error;
+               isl_seq_clr(bmap->div[j], 1+1+total);
+       }
+       return bmap;
+error:
+       isl_basic_map_free(ctx, bmap);
+       return NULL;
+}
+
+struct isl_basic_set *isl_basic_set_new_from_polylib(
+                       struct isl_ctx *ctx,
+                       Polyhedron *P, unsigned nparam, unsigned dim)
+{
+       return (struct isl_basic_set *)
+               isl_basic_map_new_from_polylib(ctx, P, nparam, 0, dim);
+}
+
+struct isl_basic_map *isl_basic_map_new_from_polylib(
+                       struct isl_ctx *ctx, Polyhedron *P,
+                       unsigned nparam, unsigned in, unsigned out)
+{
+       struct isl_basic_map *bmap;
+       unsigned extra;
+
+       isl_assert(ctx, P, return NULL);
+       isl_assert(ctx, P->Dimension >= nparam + in + out, return NULL);
+
+       extra = P->Dimension - nparam - in - out;
+       bmap = isl_basic_map_alloc(ctx, nparam, in, out, extra,
+                                       P->NbEq, P->NbConstraints - P->NbEq);
+       if (!bmap)
+               return NULL;
+
+       return copy_constraints(ctx, bmap, P);
+}
+
+struct isl_set *isl_set_new_from_polylib(struct isl_ctx *ctx,
+                       Polyhedron *D, unsigned nparam, unsigned dim)
+{
+       struct isl_set *set = NULL;
+       Polyhedron *P;
+       int n = 0;
+
+       for (P = D; P; P = P->next)
+               ++n;
+
+       set = isl_set_alloc(ctx, nparam, dim, n, ISL_MAP_DISJOINT);
+       if (!set)
+               return NULL;
+
+       for (P = D; P; P = P->next)
+               isl_set_add(ctx, set,
+                   isl_basic_set_new_from_polylib(ctx, P, nparam, dim));
+       return set;
+}
+
+struct isl_map *isl_map_new_from_polylib(struct isl_ctx *ctx,
+                       Polyhedron *D,
+                       unsigned nparam, unsigned in, unsigned out)
+{
+       struct isl_map *map = NULL;
+       Polyhedron *P;
+       int n = 0;
+
+       for (P = D; P; P = P->next)
+               ++n;
+
+       map = isl_map_alloc(ctx, nparam, in, out, n, ISL_MAP_DISJOINT);
+       if (!map)
+               return NULL;
+
+       for (P = D; P; P = P->next)
+               isl_map_add(ctx, map,
+                       isl_basic_map_new_from_polylib(ctx, P,
+                                                           nparam, in, out));
+       return map;
+}
+
+Polyhedron *isl_basic_map_to_polylib(struct isl_ctx *ctx,
+                                               struct isl_basic_map *bmap)
+{
+       int i;
+       Matrix *M;
+       Polyhedron *P;
+       unsigned off;
+
+       M = Matrix_Alloc(bmap->n_eq + bmap->n_ineq + 2*bmap->n_div,
+                1 + bmap->n_in + bmap->n_out + bmap->n_div + bmap->nparam + 1);
+       for (i = 0; i < bmap->n_eq; ++i) {
+               value_set_si(M->p[i][0], 0);
+               copy_constraint_to(M->p[i], bmap->eq[i],
+                          bmap->nparam, bmap->n_in + bmap->n_out, bmap->n_div);
+       }
+       off = bmap->n_eq;
+       for (i = 0; i < bmap->n_ineq; ++i) {
+               value_set_si(M->p[off+i][0], 1);
+               copy_constraint_to(M->p[off+i], bmap->ineq[i],
+                          bmap->nparam, bmap->n_in + bmap->n_out, bmap->n_div);
+       }
+       off += bmap->n_ineq;
+       for (i = 0; i < bmap->n_div; ++i) {
+               unsigned total = bmap->n_in+bmap->n_out+bmap->n_div+bmap->nparam;
+               if (isl_int_is_zero(bmap->div[i][0]))
+                       continue;
+               value_set_si(M->p[off+2*i][0], 1);
+               copy_constraint_to(M->p[off+2*i], bmap->div[i]+1,
+                          bmap->nparam, bmap->n_in + bmap->n_out, bmap->n_div);
+               copy_values_to(M->p[off+2*i]+1+bmap->n_in+bmap->n_out+i,
+                               bmap->div[i], 1);
+               value_oppose(M->p[off+2*i][1+bmap->n_in+bmap->n_out+i],
+                            M->p[off+2*i][1+bmap->n_in+bmap->n_out+i]);
+
+               value_set_si(M->p[off+2*i+1][0], 1);
+               Vector_Oppose(M->p[off+2*i]+1, M->p[off+2*i+1]+1, total+1);
+               value_addto(M->p[off+2*i+1][1+total], M->p[off+2*i+1][1+total],
+                           M->p[off+2*i+1][1+bmap->n_in+bmap->n_out+i]);
+               value_decrement(M->p[off+2*i+1][1+total],
+                               M->p[off+2*i+1][1+total]);
+       }
+       P = Constraints2Polyhedron(M, ctx->MaxRays);
+       Matrix_Free(M);
+
+       return P;
+}
+
+Polyhedron *isl_map_to_polylib(struct isl_ctx *ctx, struct isl_map *map)
+{
+       int i;
+       Polyhedron *R = NULL;
+       Polyhedron **next = &R;
+
+       for (i = 0; i < map->n; ++i) {
+               *next = isl_basic_map_to_polylib(ctx, map->p[i]);
+               next = &(*next)->next;
+       }
+
+       return R;
+}
+
+Polyhedron *isl_basic_set_to_polylib(struct isl_ctx *ctx,
+                       struct isl_basic_set *bset)
+{
+       return isl_basic_map_to_polylib(ctx,
+                       (struct isl_basic_map *)bset);
+}
+
+Polyhedron *isl_set_to_polylib(struct isl_ctx *ctx, struct isl_set *set)
+{
+       return isl_map_to_polylib(ctx, (struct isl_map *)set);
+}
diff --git a/isl_map_private.h b/isl_map_private.h
new file mode 100644 (file)
index 0000000..34c4471
--- /dev/null
@@ -0,0 +1,32 @@
+#include "isl_set.h"
+#include "isl_map.h"
+
+int isl_basic_map_alloc_equality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+int isl_basic_map_free_equality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned n);
+int isl_basic_map_alloc_inequality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+int isl_basic_map_free_inequality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned n);
+int isl_basic_map_alloc_div(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+int isl_basic_map_free_div(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned n);
+
+int isl_inequality_negate(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned pos);
+
+struct isl_basic_map *isl_basic_map_cow(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+struct isl_map *isl_map_cow(struct isl_ctx *ctx, struct isl_map *map);
+
+struct isl_basic_map *isl_basic_map_set_to_empty(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap);
+struct isl_map *isl_basic_map_compute_divs(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap);
+struct isl_map *isl_map_compute_divs(struct isl_ctx *ctx, struct isl_map *map);
+struct isl_basic_map *isl_basic_map_align_divs(struct isl_ctx *ctx,
+               struct isl_basic_map *dst, struct isl_basic_map *src);
+struct isl_basic_map *isl_basic_map_gauss(struct isl_ctx *ctx,
+       struct isl_basic_map *bmap, int *progress);
diff --git a/isl_piplib.c b/isl_piplib.c
new file mode 100644 (file)
index 0000000..d389d73
--- /dev/null
@@ -0,0 +1,9 @@
+#include "isl_piplib.h"
+
+void isl_seq_cpy_to_pip(Entier *dst, isl_int *src, unsigned len)
+{
+       int i;
+
+       for (i = 0; i < len; ++i)
+               entier_assign(dst[i], src[i]);
+}
diff --git a/isl_seq.c b/isl_seq.c
new file mode 100644 (file)
index 0000000..621fe85
--- /dev/null
+++ b/isl_seq.c
@@ -0,0 +1,152 @@
+#include "isl_seq.h"
+
+void isl_seq_clr(isl_int *p, unsigned len)
+{
+       int i;
+       for (i = 0; i < len; ++i)
+               isl_int_set_si(p[i], 0);
+}
+
+void isl_seq_neg(isl_int *dst, isl_int *src, unsigned len)
+{
+       int i;
+       for (i = 0; i < len; ++i)
+               isl_int_neg(dst[i], src[i]);
+}
+
+void isl_seq_cpy(isl_int *dst, isl_int *src, unsigned len)
+{
+       int i;
+       for (i = 0; i < len; ++i)
+               isl_int_set(dst[i], src[i]);
+}
+
+void isl_seq_scale(isl_int *dst, isl_int *src, isl_int m, unsigned len)
+{
+       int i;
+       for (i = 0; i < len; ++i)
+               isl_int_mul(dst[i], src[i], m);
+}
+
+void isl_seq_scale_down(isl_int *dst, isl_int *src, isl_int m, unsigned len)
+{
+       int i;
+       for (i = 0; i < len; ++i)
+               isl_int_divexact(dst[i], src[i], m);
+}
+
+void isl_seq_combine(isl_int *dst, isl_int m1, isl_int *src1,
+                       isl_int m2, isl_int *src2, unsigned len)
+{
+       int i;
+       isl_int tmp;
+
+       isl_int_init(tmp);
+       for (i = 0; i < len; ++i) {
+               isl_int_mul(tmp, m1, src1[i]);
+               isl_int_addmul(tmp, m2, src2[i]);
+               isl_int_set(dst[i], tmp);
+       }
+       isl_int_clear(tmp);
+}
+
+/*
+ * Let d = dst[pos] and s = src[pos]
+ * dst is replaced by |s| dst - sgn(s)d src
+ */
+void isl_seq_elim(isl_int *dst, isl_int *src, unsigned pos, unsigned len,
+                 isl_int *m)
+{
+       isl_int a;
+       isl_int b;
+
+       if (isl_int_is_zero(dst[pos]))
+               return;
+
+       isl_int_init(a);
+       isl_int_init(b);
+
+       isl_int_gcd(a, src[pos], dst[pos]);
+       isl_int_divexact(b, dst[pos], a);
+       if (isl_int_is_pos(src[pos]))
+               isl_int_neg(b, b);
+       isl_int_divexact(a, src[pos], a);
+       isl_int_abs(a, a);
+       isl_seq_combine(dst, a, dst, b, src, len);
+
+       if (m)
+               isl_int_mul(*m, *m, a);
+
+       isl_int_clear(a);
+       isl_int_clear(b);
+}
+
+int isl_seq_eq(isl_int *p1, isl_int *p2, unsigned len)
+{
+       int i;
+       for (i = 0; i < len; ++i)
+               if (isl_int_ne(p1[i], p2[i]))
+                       return 0;
+       return 1;
+}
+
+int isl_seq_first_non_zero(isl_int *p, unsigned len)
+{
+       int i;
+
+       for (i = 0; i < len; ++i)
+               if (!isl_int_is_zero(p[i]))
+                       return i;
+       return -1;
+}
+
+int isl_seq_abs_min_non_zero(isl_int *p, unsigned len)
+{
+       int i, min = isl_seq_first_non_zero(p, len);
+       if (min < 0)
+               return -1;
+       for (i = min + 1; i < len; ++i) {
+               if (isl_int_is_zero(p[i]))
+                       continue;
+               if (isl_int_abs_lt(p[i], p[min]))
+                       min = i;
+       }
+       return min;
+}
+
+void isl_seq_gcd(isl_int *p, unsigned len, isl_int *gcd)
+{
+       int i, min = isl_seq_abs_min_non_zero(p, len);
+
+       if (min < 0) {
+               isl_int_set_si(*gcd, 0);
+               return;
+       }
+       isl_int_abs(*gcd, p[min]);
+       for (i = 0; isl_int_cmp_si(*gcd, 1) > 0 && i < len; ++i) {
+               if (i == min)
+                       continue;
+               if (isl_int_is_zero(p[i]))
+                       continue;
+               isl_int_gcd(*gcd, *gcd, p[i]);
+       }
+}
+
+u_int32_t isl_seq_hash(isl_int *p, unsigned len, unsigned bits)
+{
+       int i;
+       u_int32_t hash = 2166136261u;
+
+       for (i = 0; i < len; ++i) {
+               if (isl_int_is_zero(p[i]))
+                       continue;
+               hash *= 16777619;
+               hash ^= (i & 0xFF);
+               hash = isl_int_hash(p[i], hash);
+       }
+       if (bits == 32)
+               return hash;
+       if (bits >= 16)
+               return (hash >> bits) ^ (hash & (((u_int32_t)1 << bits) - 1));
+       return ((hash >> bits) ^ hash) & (((u_int32_t)1 << bits) - 1);
+}
diff --git a/piplib b/piplib
new file mode 160000 (submodule)
index 0000000..7485397
--- /dev/null
+++ b/piplib
@@ -0,0 +1 @@
+Subproject commit 74853979d06817842d46f3512ecff2d631fa0337