/* Most of our interaction with the Omega library is handled below. */

#include "AC.h"
#include "BoundingBox.h"
#include <map>
#include "basic/util.h"

void printrel(Relation R, const string &s, bool simp)
{
  cout << s << "\n\n";
  Relation Q = copy(R);
  if (simp)
    Q.simplify();
  Q.print();
  cout << "\n\n";
}

/* Output text to cout illustrating why Must_Be_Subset(a, b) is false.
   Does not modify a or b. */
void show_subset_failure(Relation a, Relation b)
{
  printrel(b, "This set:");
  printrel(a, "was supposed to be a superset of this set:");
  Relation diff = Difference(copy(a), copy(b));
  assert(diff.is_satisfiable());
  diff = Sample_Solution(diff);
  cout << "Sample element not in the supposed superset:" << endl;
  diff.print();
}

static Relation *new_generic_relation(int l, Variable_ID *&x, F_And *&root)
{
  Relation *r = new Relation(NonCoercible<int> (l));
  x = new Variable_ID [l + 1];
  {
    char v[3] = {'x', ' ', '\0'};
    for (int i = 1; i <= l; i++) {
      v[1] = i + '0';
      r->name_set_var(i, v);
      x[i] = r->set_var(i);
    }
  }
  root = r->add_and();
  return r;
}

static Relation *new_generic_exists_relation(int l, Variable_ID *&x,
					     F_Exists *&e)
{
  Relation *r = new Relation(NonCoercible<int> (l));
  x = new Variable_ID [l + 1];
  {
    char v[3] = {'x', ' ', '\0'};
    for (int i = 1; i <= l; i++) {
      v[1] = i + '0';
      r->name_set_var(i, v);
      x[i] = r->set_var(i);
    }
  }
  e = r->add_exists();
  return r;
}

/* If part of a space has been done, and we have a dep within that
   space, is the dep violated?  That is, is something we should have
   done according to r not in been_done? */
bool dep_violated(Relation const *been_done, Relation const *r)
{
  DBG(printrel(*r, "Checking dep:"));

  // What is needed by what has been done.  Should be a subset of been_done.
  Relation needed = needed_by(been_done, r);
  Relation diff = Difference(copy(needed), copy(*been_done));
  if (diff.is_satisfiable()) {
    DBG({
      Relation v = Sample_Solution(diff);
      printrel(v, "Problem! Sample element of needed - been_done:");
      Relation x = Composition(copy(*r), v);
      printrel(x, "D(sample element):");
      printrel(needed, "needed:");
    });
    return true;
  }
  return false;
}

/* Non-destructive. */
Relation needed_by(Relation const *been_done, Relation const *r)
{
  return Composition(Inverse(copy(*r)),
		     copy(*been_done));  
}

/* Create a set containing the point specified by g.  Compose with r.
   Non-destructive. */
Relation apply_to_generic(Relation r, llist<Free_Var_Decl *> *g)
{
  const int n = r.n_inp();
  Relation s(n);
  F_And *root = s.add_and();
  for (int i = 0; i < n; i++, g = g->tail()) {
    EQ_Handle e = root->add_EQ();
    e.update_coef(s.set_var(i + 1), -1);
    e.update_coef(s.get_local(g->front()), 1);
    e.finalize();
  }
  return Join(s, copy(r));
}

/* If r = { x -> y | f(x, y) } then return { y | f(v, y) }. */
/* Non-destructive. */
Relation eval_relation_at_v(Relation r, intvector const *v)
{
  Relation *q = intvector_to_singleton_set(v);
  return Join(*q, copy(r));
}

/* If r = { x -> y | f(x, y) } then return { y | f(0, y) }. */
/* Non-destructive. */
Relation eval_relation_at_zero(Relation r)
{
  intvector *z = zero_vector(r.n_inp());
  r = eval_relation_at_v(r, z);
  delete z;
  return r;
}

/* Non-destructive. */
intvector *singleton_set_to_intvector(Relation r)
{
  intvector *v;
  bool b = is_singleton_set(r, &v);
  if (!b) {
    printrel(r, "Fatal error in singleton_set_to_intvector(). r:");
    fatal_error("");
  }
  return v;
}

/* Non-destructive. */
string singleton_set_to_string(Relation r, bool force)
{
  intvector *v;
  bool b = is_singleton_set(r, &v);
  if (!b) {
    if (force) {
      printrel(r, "Fatal error in singleton_set_to_intvector(). r:");
      fatal_error("");
    }
    return "UNKNOWN";
  }
  string s = v->to_string();
  delete v;
  return s;
}
Relation *intvector_to_singleton_set(intvector const *n)
{
  int l = n->size();
  Variable_ID *x;
  F_And *root;
  Relation *r = new_generic_relation(l, x, root);
  for (int i = 1; i <= l; i++) {
    EQ_Handle h = root->add_EQ();
    h.update_coef(x[i], 1);
    h.update_const(-(*n)[i - 1]);
  }
  return r;
}

Relation box_to_set(intvector const *lo, intvector const *hi)
{
  int l = lo->size();
  Variable_ID *x;
  F_And *root;
  Relation *r = new_generic_relation(l, x, root);
  for (int i = 1; i <= l; i++) {
    GEQ_Handle g = root->add_GEQ();
    g.update_coef(x[i], 1);
    g.update_const(-(*lo)[i - 1]);
    GEQ_Handle h = root->add_GEQ();
    h.update_coef(x[i], -1);
    h.update_const((*hi)[i - 1]);
  }
  return *r;
}

/* Does not modify its args. */
bool set_contains_intvector(Relation s, intvector const *v)
{
  return Must_Be_Subset(*intvector_to_singleton_set(v), copy(s));
}

/* Lines and halflines */

/* {x | exists(integer t) x = z + tn} */
Relation *wholeline(intvector const *n, intvector const *z)
{
  map<string, Relation *> memo;
  string s = n->to_string() + z->to_string();
  if (memo[s])
    return memo[s];

  int l = n->size();
  Variable_ID *x;
  F_Exists *e;
  Relation *r = memo[s] = new_generic_exists_relation(l, x, e);
  Variable_ID t = e->declare("t");
  F_And *root = e->add_and();
  
  for (int i = 1; i <= l; i++) {
    /* xi - t * ni - zi = 0 */ 
    EQ_Handle c = root->add_EQ();
    c.update_coef(x[i], 1);
    c.update_coef(t, -(*n)[i - 1]);
    c.update_const(-(*z)[i - 1]); 
  }

  r->finalize();
  // printrel(*r, "wholeline is:");
  return r;
}

/* Try to find the smallest q, lo <= q <= hi such that
   r is a subset of halfline(n, z, q).  Return false if none exists.
   Set q and return true if successful. */
bool find_q_in_range_to_cover_halfline(Relation const *r, int lo, int hi,
				       intvector const *n,
				       intvector const *z, int &q)
{
  // printrel(*r, "find_q_in_range_to_cover: ");
  if (Must_Be_Subset(copy(*r), copy(*halfline(n, z, lo))))
    q = lo;
  else if (hi <= lo || !Must_Be_Subset(copy(*r), copy(*halfline(n, z, hi))))
    return false;
  else {
    while (lo + 1 < hi) {
      int mid = (lo + hi) / 2; 
      if (Must_Be_Subset(copy(*r), copy(*halfline(n, z, mid))))
	hi = mid;
      else
	lo = mid;
    }
    q = hi;
  }
  return true;
}

/* Try to find the smallest q, 0 <= q <= MAXINT/2 such that
   r is a subset of halfline(n, z, q).  Return false if none exists.
   Set q and return true if successful. */
bool find_nonneg_q_to_cover_halfline(Relation const *r, intvector const *n,
				     intvector const *z, int &q)
{
  return find_q_in_range_to_cover_halfline(r, 0, MAXINT / 2,
					   n, z, q);
}

/* {x | exists(integer t) t < q && x = z + tn} */
Relation *halfline(intvector const *n, intvector const *z, int q)
{
  map<string, Relation *> memo;
  string s = n->to_string() + z->to_string() + i2s(q);
  if (memo[s])
    return memo[s];

  int l = n->size();
  Variable_ID *x;
  F_Exists *e;
  Relation *r = memo[s] = new_generic_exists_relation(l, x, e);
  Variable_ID t = e->declare("t");
  F_And *root = e->add_and();
  
  /* t < q (i.e., t <= q - 1, i.e., q - t - 1 >= 0) */
  GEQ_Handle tmax = root->add_GEQ();
  tmax.update_coef(t, -1);
  tmax.update_const(q - 1);

  for (int i = 1; i <= l; i++) {
    /* xi - t * ni - zi = 0 */ 
    EQ_Handle c = root->add_EQ();
    c.update_coef(x[i], 1);
    c.update_coef(t, -(*n)[i - 1]);
    c.update_const(-(*z)[i - 1]); 
  }

  r->finalize();
  // printrel(*r, "halfline is:");
  return r;
}

/* {x | exists(integer t) t >= q && x = z + tn} */
Relation *antihalfline(intvector const *n, intvector const *z, int q)
{
  map<string, Relation *> memo;
  string s = n->to_string() + z->to_string() + i2s(q);
  if (memo[s])
    return memo[s];

  int l = n->size();
  Variable_ID *x;
  F_Exists *e;
  Relation *r = memo[s] = new_generic_exists_relation(l, x, e);
  Variable_ID t = e->declare("t");
  F_And *root = e->add_and();
  
  /* t >= q (i.e., t - q >= 0) */
  GEQ_Handle tmin = root->add_GEQ();
  tmin.update_coef(t, 1);
  tmin.update_const(-q);

  for (int i = 1; i <= l; i++) {
    /* xi - t * ni - zi = 0 */ 
    EQ_Handle c = root->add_EQ();
    c.update_coef(x[i], 1);
    c.update_coef(t, -(*n)[i - 1]);
    c.update_const(-(*z)[i - 1]); 
  }

  r->finalize();
  // printrel(*r, "antihalfline is:");
  return r;
}

/* Spaces and halfspaces */

/* Equivalent to preceding([0, ..., 0]). */
Relation lexicographically_negative(int n)
{
  Relation r = Relation(NonCoercible<int> (n));
  F_Or *root = r.add_or();
  for (int i = 1; i <= n; i++) {
    F_And *c = root->add_and();
    for (int j = 1; j <= i; j++)
      if (i == j) {
	GEQ_Handle h = c->add_GEQ();
	h.update_coef(r.set_var(i), -1);
	h.update_const(-1);
      } else {
	EQ_Handle h = c->add_EQ();
	h.update_coef(r.set_var(j), 1);
      }
  }
  return r;
} 

/* {v | v dot n < q } */
Relation *halfspace(intvector const *n, int q)
{
  map<string, Relation *> memo;
  string s = n->to_string() + i2s(q);
  if (memo[s])
    return memo[s];

  int l = n->size();
  Variable_ID *x;
  F_And *root;
  Relation *r = memo[s] = new_generic_relation(l, x, root);
  GEQ_Handle h = root->add_GEQ();
  for (int i = 1; i <= l; i++)
    h.update_coef(x[i], -(*n)[i - 1]);
  h.update_const(-1 + q);

  r->finalize();
  // printrel(*r, "halfspace is:");
  return r;
}

/* {v | v dot n >= q } */
Relation *antihalfspace(intvector const *n, int q)
{
  map<string, Relation *> memo;
  string s = n->to_string() + i2s(q);
  if (memo[s])
    return memo[s];

  int l = n->size();
  Variable_ID *x;
  F_And *root;
  Relation *r = memo[s] = new_generic_relation(l, x, root);
  GEQ_Handle h = root->add_GEQ();
  for (int i = 1; i <= l; i++)
    h.update_coef(x[i], (*n)[i - 1]);
  h.update_const(-q);

  r->finalize();
  // printrel(*r, "halfspace is:");
  return r;
}

/* {v | v dot n < k }, where k is a Free_Var_Decl. */
Relation *generic_halfspace(intvector const *n)
{
  static Free_Var_Decl *k = NULL;
  map<string, Relation *> memo;
  string s = n->to_string();
  if (memo[s])
    return memo[s];

  if (k == NULL)
    k = new Free_Var_Decl("k");
  int l = n->size();
  Variable_ID *x;
  F_And *root;
  Relation *r = memo[s] = new_generic_relation(l, x, root);
  GEQ_Handle h = root->add_GEQ();
  for (int i = 1; i <= l; i++)
    h.update_coef(x[i], -(*n)[i - 1]);
  h.update_coef(r->get_local(k), 1);
  h.update_const(-1);

  r->finalize();
  // printrel(*r, "Generic halfspace is:");
  return r;
}

/* {v | v dot n < k }, where k is a new Free_Var_Decl. unused. */
Relation *generic_halfspace_(intvector const *n)
{
  static int kcount = 0;
  map<string, Relation *> memo;
  string s = n->to_string();
  if (memo[s])
    return memo[s];

  int l = n->size();
  char kname[20];
  sprintf(kname, "k%d", kcount++);
  Free_Var_Decl *k = new Free_Var_Decl(kname);
  Variable_ID *x;
  F_And *root;
  Relation *r = memo[s] = new_generic_relation(l, x, root);
  GEQ_Handle h = root->add_GEQ();
  for (int i = 1; i <= l; i++)
    h.update_coef(x[i], -(*n)[i - 1]);
  h.update_coef(r->get_local(k), 1);
  h.update_const(-1);

  r->finalize();
  // printrel(*r, "Generic halfspace is:");
  return r;
}

/* Return an n-dimensional empty set. */
Relation *emptyset(int n)
{
  return new Relation(Relation::False(n));
}

Relation *wholespace(int n)
{
  return new Relation(Relation::True(n));
}

/* Given S is done, and D is a dependence, what appears ready to execute? */
Relation ready_wrt_dep(Relation S, Relation D)
{
  Relation r = Complement(Composition(copy(D), Complement(copy(S))));
  r.simplify();
  return r;
}

/* z is a set.
   Let s = {x | exists(p lexicographically positive) s.t. x + p is in z}.
   Return s (or s union z if include_z is true). */
/* non-destructive */
Relation preceding(Relation z, bool include_z)
{
  int n = z.n_set();
  /* Construct relation {x -> y | x precedes y lexicographically}. */
  Relation r = Relation(n, n);
  F_Or *root = r.add_or();
  for (int i = 1; i <= n; i++) {
    F_And *c = root->add_and();
    for (int j = 1; j <= i; j++)
      if (i == j) {
	GEQ_Handle h = c->add_GEQ();
	h.update_coef(r.input_var(j), -1);
	h.update_coef(r.output_var(j), 1);
	h.update_const(-1);
      } else {
	EQ_Handle h = c->add_EQ();
	h.update_coef(r.input_var(j), -1);
	h.update_coef(r.output_var(j), 1);
      }
  }
  r = Join(copy(z), Inverse(r));
  return (include_z ? Union(r, copy(z)) : r);
}

/* z is a set.
   Let s = {(x1, ..., xn) | exists(p > 0) s.t.
			   (x1, ..., xn) + (0, ..., 0, p) is in z}.
   Return s (or s union z if include_z is true). */
/* non-destructive */
Relation halflines(Relation z, bool include_z)
{
  int n = z.n_set();
  Relation r = Relation(n, n);
  F_And *c = r.add_and();
  for (int i = 1; i <= n; i++)
    if (i == n) {
      GEQ_Handle h = c->add_GEQ();
      h.update_coef(r.input_var(i), -1);
      h.update_coef(r.output_var(i), 1);
      h.update_const(-1);
    } else {
      EQ_Handle h = c->add_EQ();
      h.update_coef(r.input_var(i), -1);
      h.update_coef(r.output_var(i), 1);
    }
  r = Join(copy(z), Inverse(r));
  return (include_z ? Union(r, copy(z)) : r);
}

/* Return the lexicographically greatest element of the set z. */
/* destructive */
Relation last(Relation z)
{
  Relation prec = preceding(z, false);
  return Difference(z, prec);
}

/* destructive */
Relation keep_last_per_bundle(Relation z)
{
  Relation prec = halflines(z, false);
  return Difference(z, prec);
}

/* Translate s by the sum over i of a_i * v_i. */
/* s is a set of n dimensional points; there are m v's each of dimension n. */
Relation translate_by_weighted_vectors(Relation s, Free_Var_Decl **a,
				       llist<intvector *> *v)
{
  size_t n = s.n_set(), m = v->size();
  assert(n == v->front()->size());
  Relation r(n, n);
  F_And *root = r.add_and();
  for (int i = 0; i < (int) n; i++) {
    EQ_Handle e = root->add_EQ();
    e.update_coef(r.input_var(i + 1), 1);
    e.update_coef(r.output_var(i + 1), -1);
    for (int j = 0; j < (int) m; j++)
      e.update_coef(r.get_local(a[j]), (*((*v)[j]))[i]);
    e.finalize();
  }
  Relation result = Composition(r, copy(s));
  // printrel(s, "Translated this set generically:");
  // printrel(result, "And got:");
  return result;
}

/* Given a set s and a vector v, return the set { v + x | x is in s }. */
/* Does not modify its args.  Modification of returned result will not
   modify args. */
Relation translate_set(Relation s, intvector const *v)
{
  return translate_relation(s, v);
}

Relation translate_range(Relation r, intvector const *v)
{
  return translate_relation(r, v);
}

Relation translate_domain(Relation r, intvector const *v)
{
  return Inverse(translate_range(Inverse(copy(r)), v));
}

/* Non-destructive. */
Relation translate_relation(Relation s, intvector const *v, int n)
{
  if (n == -1)
    n = (s.is_set() ? s.n_set() : s.n_out());
  Relation r(n, n);
  F_And *root = r.add_and();
  for (int j = 0; j < n; j++) {
    EQ_Handle e = root->add_EQ();
    e.update_coef(r.input_var(j + 1), 1);
    e.update_coef(r.output_var(j + 1), -1);
    e.update_const((*v)[j]);
  }
  Relation result = Composition(r, copy(s));
  // printrel(s, "Translated this set by " + v->to_string() + ":");
  // printrel(result, "And got:");
  return result;
}


/* Return the identity relation: { (x1, ..., xn) -> (x1, ..., xn) }. */
Relation identity_relation(int n)
{
  Relation r(n, n);
  F_And *root = r.add_and();
  for (int j = 0; j < n; j++) {
    EQ_Handle e = root->add_EQ();
    e.update_coef(r.input_var(j + 1), 1);
    e.update_coef(r.output_var(j + 1), -1);
  }
  return r;
}

/* Apply { (x1, ..., xn) -> (xn, ..., x1) } to set s. */
Relation dreverse(Relation s)
{
  int n = s.n_set();
  if (n == 1)
    return s;
  Relation r(n, n);
  F_And *root = r.add_and();
  for (int j = 1; j <= n; j++) {
    EQ_Handle e = root->add_EQ();
    e.update_coef(r.input_var(j), 1);
    e.update_coef(r.output_var(n + 1 - j), -1);
  }
  return Join(s, r);
}

/* Given a set s = { (x1, ...., xn), (y1, ...., yn), ... }, return the set
   { (x1, x2, ..., -xi, ... xn), (y1, y2, ..., -yi, ... yn), ... }. */
/* i counts from 0. */
Relation negate_ith_var(Relation s, int i)
{
  int n = s.n_set();
  Relation r(n, n);
  F_And *root = r.add_and();
  for (int j = 0; j < n; j++) {
    EQ_Handle e = root->add_EQ();
    e.update_coef(r.input_var(j + 1), 1);
    e.update_coef(r.output_var(j + 1), (i == j) ? 1 : -1);
  }
  Relation result = Composition(r, copy(s));
  // printrel(s, "Negated dimension " + int2string(i) + " of:");
  // printrel(result, "And got:");
  return result;
}

inline Relation restrict_ith_dimension_to_be_below_z(Relation s, int i, int z)
{
  return Intersection(*halfspace(unit_vector(i, s.n_set()), z), copy(s));
}

inline Relation restrict_ith_dimension_to_be_at_or_above_z(Relation s, int i, int z)
{
  return Intersection(*antihalfspace(unit_vector(i, s.n_set()), z), copy(s));
}

inline Relation restrict_ith_dimension_to_be_above_z(Relation s, int i, int z)
{
  return restrict_ith_dimension_to_be_at_or_above_z(s, i, z + 1);
}

/* Try to find a lower bound for set s in the i'th dimension.  If one
   cannot be found then return false.  Otherwise set result and return true. */
/* i counts from 0. */
bool measure_lower_bound(Relation s, int i, int &result)
{
  if (!s.is_satisfiable())
    return false;
  int hi = MAXINT / 2, lo = -hi;
  /* The following line is a hack to try to get around the fact that the
     Omega library can yield wrong answers when dealing with large
     integers.  For example, Intersection({[x] Exists ( alpha : x = 4alpha) },
     {x < -1073741823}) returns the empty set. */
  if (!Must_Be_Subset(copy(s), restrict_ith_dimension_to_be_below_z(restrict_ith_dimension_to_be_above_z(Relation::True(s.n_set()), i, lo/10), i, hi/10)))
    return false;
  if (restrict_ith_dimension_to_be_below_z(s, i, lo).is_satisfiable() ||
      !restrict_ith_dimension_to_be_below_z(s, i, hi).is_satisfiable())
    return false;
  do {
    int mid = (lo + hi) / 2;
    if (restrict_ith_dimension_to_be_below_z(s, i, mid).is_satisfiable())
      hi = mid;
    else
      lo = mid;
  } while (hi - lo > 1);
  result = lo;
  return true;
}

/* Analogous to measure_lower_bound().  i counts from 0. */
bool measure_upper_bound(Relation s, int i, int &result)
{
  if (measure_lower_bound(negate_ith_var(s, i), i, result)) {
    result = -result;
    return true;
  } else
    return false;
}

/* Sets b and returns true if successful. */
/* Non-destructive. */
bool find_bounding_box(Relation s, BoundingBox &b)
{
  int n = s.n_set();
  intvector *lo = zero_vector(n), *hi = zero_vector(n);
  for (int i = 0; i < n; i++) 
    if (!measure_lower_bound(s, i, (*lo)[i]) ||
	!measure_upper_bound(s, i, (*hi)[i]))
      return false;
  b = BoundingBox(lo, hi);
  // DBG(printrel(s, "Found bounding box " + b.to_string() + " for set:"));
  return true;
}

llist<intvector *> *set_bounded_by_box_to_llist(Relation s, BoundingBox *b)
{
  // cout << "set_bounded_by_box_to_llist " << b->to_string() << endl;
  if (Intersection(copy(s), box_to_set(b->lwb(), b->upb())).is_satisfiable()) {
    BoundingBox *o = b->calve();
    if (o == NULL) {
      llist<intvector *> *result = NULL;
      for (BoxIter i = b->start(); !i.is_done(); i.next())
	if (set_contains_intvector(s, *i))
	  result = cons((*i)->copy(), result);
      return result;
    } else
      return extend(set_bounded_by_box_to_llist(s, b),
		    set_bounded_by_box_to_llist(s, o));
  } else
    return NULL;
}

/* Return a set that contains all elts of s except those listed in l. 
   Destructive: may modify s (but not l).  s need not contain elements of l. */
Relation remove_from_set(Relation s, llist<intvector *> *l)
{
  return Difference(s, llist_to_set(l, s.n_set()));
}

/* n is the arity of the vectors.  If 0, use n = l->front->size(). */
Relation llist_to_set(llist<intvector *> *l, int n)
{
  assert(l != NULL || n > 0);
  if (n == 0)
    n = l->front()->size();
  Relation r = Relation::False(n);
  while (l != NULL) {
    r = Union(r, *intvector_to_singleton_set(l->front()));
    l = l->tail();
  }
  return r;
}

/* Returns whether s is finite.  If not, result is set to NULL. If so,
   result is set to a list of the points in s. */
bool set_to_llist(Relation s, llist<intvector *> *&result)
{
  result = NULL;
  if (!s.is_lower_bound_satisfiable())
    return true;
  BoundingBox b;
  if (!find_bounding_box(s, b))
    return false;
  result = set_bounded_by_box_to_llist(s, &b);
  return true;
}

/* Inefficient. */
int size_of_finite_set(Relation s)
{
  llist<intvector *> *l;
  assert(set_to_llist(s, l));
  int n = l->size();
  for (ListIterator<intvector *> i(l); !i.isDone(); i.next())
    delete *i;
  free_all(l);
  return n;
}

/* Is the set s finite?  (inefficient implementation) */
bool is_finite(Relation s)
{
  BoundingBox b;
  return (!s.is_satisfiable() || find_bounding_box(s, b));
}

/* If s is a singleton set containing one fixed integer tuple then set
   *v (if v not null) to its one element and return true. */
bool is_singleton_set(Relation s, intvector **v)
{
  BoundingBox b;
  if (find_bounding_box(s, b) &&
      b.lwb()->equals(b.upb())) {
    if (v != NULL)
      *v = b.lwb();
    return true;
  }
  return false;
}

/* True iff r is a subset of s and v.v. */
bool same_relation(Relation r, Relation s)
{
  return Must_Be_Subset(copy(r), copy(s)) && Must_Be_Subset(copy(s), copy(r));
}

/* Does there exist a p such that r1(p) equals r2(p)?  Does not modify args. */
bool exists_range_overlap(Relation r1, Relation r2)
{
  // true iff r(p) = p for some p
  Relation r = Composition(Inverse(copy(r1)), copy(r2));
  Relation result = Intersection(r, identity_relation(r1.n_inp()));
  return result.is_satisfiable();
}

void test_omint()
{
  {
    intvector *l = new intvector(cons(3, cons(7)));
    intvector *h = new intvector(cons(5, cons(9)));
    intvector *o = new intvector(cons(1, cons(2)));
    Relation box = box_to_set(l, h);
    printrel(box, "box = box_to_set(" + l->to_string() + ", " + h->to_string() +
	     ") = ");
    BoundingBox b;
    assert(find_bounding_box(box, b));
    cout << "find_bounding_box yields " << b.to_string() << endl;
    printrel(translate_set(box, o),
	     "Translation of box by " + o->to_string() + ":");
    
    
    intvector *f = new intvector(cons(-10));
    printrel(preceding(*intvector_to_singleton_set(f)),
	     "preceding(" + f->to_string() + "):");
    printrel(preceding(*intvector_to_singleton_set(l)),
	     "preceding(" + l->to_string() + "):");
    printrel(preceding(*intvector_to_singleton_set(h), true),
	     "preceding(" + h->to_string() + ", true):");

    printrel(translate_set(preceding(*intvector_to_singleton_set(l)),
			     new intvector(cons(2, cons(1)))),
	     "translate_set(preceding(" + l->to_string() + "), [2, 1]):");
    printrel(translate_domain(identity_relation(2),
			      new intvector(cons(2, cons(1)))),
	     "translate_domain(I, [2, 1]):");
    printrel(translate_range(identity_relation(2),
			      new intvector(cons(2, cons(1)))),
	     "translate_range(I, [2, 1]):");
    assert(same_relation
	   (lexicographically_negative(3),
	    preceding(*intvector_to_singleton_set
		      (new intvector(cons(0, cons(0, cons(0))))))));
  }
  
  intvector *a = new intvector(cons(1, cons(2, cons(3))));
  intvector *b = new intvector(cons(1, cons(0, cons(0))));
  cout << "trying to create halfspace {v | v dot " << a->to_string() <<
    " < 3}" << endl;
  halfspace(a, 3);
  cout << "trying to create halfspace {v | v dot " << b->to_string() <<
    " < 5}" << endl;
  halfspace(b, 5);
  cout << "trying to create halfspace {v | v dot " << b->to_string() <<
    " >= 5}" << endl;
  antihalfspace(b, 5);
  cout << "trying to create wholeline {x | exists(t): x = " <<
    b->to_string() << " + t " << a->to_string() << "}" << endl;
  wholeline(a, b);
  cout << "trying to create halfline {x | exists(t < 0): x = " <<
    b->to_string() << " + t " << a->to_string() << "}" << endl;
  halfline(a, b);
  cout << "trying to create antihalfline {x | exists(t >= 0): x = " <<
    b->to_string() << " + t " << a->to_string() << "}" << endl;
  antihalfline(a, b);
  /*
    printrel(*emptyset(4), "Empty set in 4D:");
    printrel(*wholespace(4), "Whole space (4D):");
  */

  Relation s = Union(*intvector_to_singleton_set(a),
		     *intvector_to_singleton_set(b));
  
  printrel(s, "s:");
  llist<intvector *> *x;
  assert(set_to_llist(s, x));
  cout << "s contains:";
  for (; x != NULL; x = x->tail())
    cout << " " << x->front()->to_string();
  cout << endl;
}