Merge pull request #2 from dhuth/master

Moved omega into chill.

1 files changed, 1489 insertions, 0 deletions

diff --git a/omega/omega_lib/src/hull.cc b/omega/omega_lib/src/hull.cc
new file mode 100644
index 0000000..f1b0601
--- /dev/null
+++ b/omega/omega_lib/src/hull.cc
@@ -0,0 +1,1489 @@
+/*****************************************************************************
+ Copyright (C) 1994-2000 University of Maryland
+ Copyright (C) 2008 University of Southern California
+ Copyright (C) 2009-2010 University of Utah
+ All Rights Reserved.
+
+ Purpose:
+   Various hull calculations.
+
+ Notes:
+
+ History:
+  06/15/09  ConvexRepresentation, Chun Chen
+  11/25/09  RectHull, Chun Chen
+*****************************************************************************/
+
+#include <omega.h>
+#include <omega/farkas.h>
+#include <omega/hull.h>
+#include <basic/Bag.h>
+#include <basic/Map.h>
+#include <basic/omega_error.h>
+#include <list>
+#include <vector>
+#include <set>
+
+namespace omega {
+
+int hull_debug = 0; 
+
+Relation ConvexHull(NOT_CONST Relation &R) {
+  Relation S = Approximate(consume_and_regurgitate(R));
+  if (!S.is_upper_bound_satisfiable())
+    return S;
+  if (S.has_single_conjunct())
+    return S;
+  return Farkas(Farkas(S,Basic_Farkas), Convex_Combination_Farkas);
+}
+
+Relation DecoupledConvexHull(NOT_CONST Relation &R) {
+  Relation S = Approximate(consume_and_regurgitate(R));
+  if (!S.is_upper_bound_satisfiable())
+    return S;
+  if (S.has_single_conjunct())
+    return S;
+  return Farkas(Farkas(S,Decoupled_Farkas), Convex_Combination_Farkas);
+}
+
+Relation AffineHull(NOT_CONST Relation &R) {
+  Relation S = Approximate(consume_and_regurgitate(R));
+  if (!S.is_upper_bound_satisfiable())
+    return S;
+  return Farkas(Farkas(S,Basic_Farkas), Affine_Combination_Farkas);
+}
+
+Relation LinearHull(NOT_CONST Relation &R) {
+  Relation S = Approximate(consume_and_regurgitate(R));
+  if (!S.is_upper_bound_satisfiable())
+    return S;
+  return Farkas(Farkas(S,Basic_Farkas), Linear_Combination_Farkas);
+}
+
+Relation ConicHull(NOT_CONST Relation &R) {
+  Relation S = Approximate(consume_and_regurgitate(R));
+  if (!S.is_upper_bound_satisfiable())
+    return S;  
+  return Farkas(Farkas(S,Basic_Farkas), Positive_Combination_Farkas);
+}
+
+
+Relation FastTightHull(NOT_CONST Relation &input_R, NOT_CONST Relation &input_H) {
+  Relation R = Approximate(consume_and_regurgitate(input_R));
+  Relation H = Approximate(consume_and_regurgitate(input_H));
+
+  if (hull_debug) {
+    fprintf(DebugFile,"[ Computing FastTightHull of:\n");
+    R.prefix_print(DebugFile);
+    fprintf(DebugFile,"given known hull of:\n");
+    H.prefix_print(DebugFile);
+  }
+
+  if (!H.has_single_conjunct()) {
+    if (hull_debug) 
+      fprintf(DebugFile, "] bailing out of FastTightHull, known hull not convex\n");
+    return H;
+  }
+
+  if (!H.is_obvious_tautology()) {
+    R = Gist(R,copy(H));
+    R.simplify(1,0);
+  }
+
+  if (R.has_single_conjunct()) {
+    R = Intersection(R,H);
+    if (hull_debug)  {
+      fprintf(DebugFile, "] quick easy answer to FastTightHull\n");
+      R.prefix_print(DebugFile);
+    }
+    return R;
+  }
+  if (R.has_local(coefficient_of_constant_term)) {
+    if (hull_debug) {
+      fprintf(DebugFile, "] Can't handle recursive application of Farkas lemma\n");
+    }
+    return H;
+  }
+  
+  if (hull_debug) {
+    fprintf(DebugFile,"Gist of R given H is:\n");
+    R.prefix_print(DebugFile);
+  }
+
+  if (1) {
+    Set<Variable_ID> vars;
+    int conjuncts = 0;
+    for (DNF_Iterator s(R.query_DNF()); s.live(); s.next()) {
+      conjuncts++;
+      for (Variable_ID_Iterator v(*((*s)->variables())); v.live(); v++) {
+        bool found = false;
+        for (EQ_Iterator eq = (*s)->EQs(); eq.live(); eq.next())
+          if ((*eq).get_coef(*v) != 0) {
+            if (!found) vars.insert(*v);
+            found = true;
+            break;
+          }
+        if (!found)
+          for (GEQ_Iterator geq = (*s)->GEQs(); geq.live(); geq.next())
+            if ((*geq).get_coef(*v) != 0) {
+              if (!found) vars.insert(*v);
+              found = true;
+              break;
+            }
+      }
+    }
+
+     
+    // We now know which variables appear in R
+    if (hull_debug) {
+      fprintf(DebugFile,"Variables we need a better hull on are: ");
+      foreach(v,Variable_ID,vars,
+              fprintf(DebugFile," %s",v->char_name()));
+      fprintf(DebugFile,"\n");
+    }
+    Conjunct *c = H.single_conjunct();
+    int total=0;
+    int copied = 0;
+    for (EQ_Iterator eq = c->EQs(); eq.live(); eq.next()) {
+      total++;
+      foreach(v,Variable_ID,vars,
+              if ((*eq).get_coef(v) != 0) {
+                R.and_with_EQ(*eq);
+                copied++;
+                break; // out of variable loop
+              }
+        );
+    }
+    for (GEQ_Iterator geq = c->GEQs(); geq.live(); geq.next()) {
+      total++;
+      foreach(v,Variable_ID,vars,
+              if ((*geq).get_coef(v) != 0) {
+                R.and_with_GEQ(*geq);
+                copied++;
+                break; // out of variable loop
+              }
+        );
+    }
+    if (copied < total) {
+      R = Approximate(R);
+
+      if (hull_debug) { 
+        fprintf(DebugFile,"Decomposed relation, copied only %d of %d constraints\n",copied,total);
+        fprintf(DebugFile,"Original R:\n");
+        R.prefix_print(DebugFile);
+        fprintf(DebugFile,"Known hull:\n");
+        H.prefix_print(DebugFile);
+        fprintf(DebugFile,"New R:\n");
+        R.prefix_print(DebugFile);
+      }
+    }
+  }
+
+  Relation F = Farkas(copy(R), Basic_Farkas, true);
+  if (hull_debug)  
+    fprintf(DebugFile,"Farkas Difficulty = " coef_fmt "\n", farkasDifficulty);
+  if (farkasDifficulty > 260) {
+    if (hull_debug)  {
+      fprintf(DebugFile, "] bailing out, farkas is way too complex\n");
+      fprintf(DebugFile,"Farkas:\n");
+      F.prefix_print(DebugFile);
+    }
+    return H;
+  }
+  else if (farkasDifficulty > 130) {
+    // Bail out
+    if (hull_debug)  {
+      fprintf(DebugFile, coef_fmt " non-zeros in original farkas\n", farkasDifficulty);
+    }
+    Relation tmp = Farkas(R, Decoupled_Farkas, true);
+    
+    if (hull_debug)  {
+      fprintf(DebugFile, coef_fmt " non-zeros in decoupled farkas\n", farkasDifficulty);
+    }
+    if (farkasDifficulty > 260)  {
+      if (hull_debug)  {
+        fprintf(DebugFile, "] bailing out, farkas is way too complex\n");
+        fprintf(DebugFile,"Farkas:\n");
+        F.prefix_print(DebugFile);
+      }
+      return H;
+    }
+    else {
+      if (farkasDifficulty > 130) 
+        R = Intersection(H, Farkas(tmp, Affine_Combination_Farkas, true));
+      else R = Intersection(H,
+                            Intersection(Farkas(tmp, Convex_Combination_Farkas, true),
+                                         Farkas(F, Affine_Combination_Farkas, true)));
+      if (hull_debug)  {
+        fprintf(DebugFile, "] bailing out, farkas is too complex, using affine hull\n");
+        fprintf(DebugFile,"Farkas:\n");
+        F.prefix_print(DebugFile);
+        fprintf(DebugFile,"Affine hull:\n");
+        R.prefix_print(DebugFile);
+      }
+      return R;
+    }
+  }
+  
+  R = Intersection(H, Farkas(F, Convex_Combination_Farkas, true));
+  if (hull_debug)  {
+    fprintf(DebugFile, "] Result of FastTightHull:\n");
+    R.prefix_print(DebugFile);
+  }
+  return R;
+}
+
+
+
+namespace {
+  bool parallel(const GEQ_Handle &g1, const GEQ_Handle &g2) {
+    for(Constr_Vars_Iter cvi(g1, false); cvi; cvi++) {
+      coef_t c1 = (*cvi).coef;
+      coef_t c2 = g2.get_coef((*cvi).var);
+      if (c1 != c2) return false;
+    }
+    {
+      for(Constr_Vars_Iter cvi(g2, false); cvi; cvi++) {
+        coef_t c1 = g1.get_coef((*cvi).var);
+        coef_t c2 = (*cvi).coef;
+        if (c1 != c2) return false;
+      }
+    }
+    return true;
+  }
+
+
+  bool hull(const EQ_Handle &e, const GEQ_Handle &g, coef_t &hull) {
+    int sign = 0;
+    for(Constr_Vars_Iter cvi(e, false); cvi; cvi++) {
+      coef_t c1 = (*cvi).coef;
+      coef_t c2 = g.get_coef((*cvi).var);
+      if (sign == 0) sign = (c1*c2>=0?1:-1);
+      if (sign*c1 != c2) return false;
+    }
+    assert(sign != 0);
+    {
+      for(Constr_Vars_Iter cvi(g, false); cvi; cvi++) {
+        coef_t c1 = e.get_coef((*cvi).var);
+        coef_t c2 = (*cvi).coef;
+        if (sign*c1 != c2) return false;
+      }
+    }
+    hull = max(sign * e.get_const(), g.get_const());
+    if (hull_debug) {
+      fprintf(DebugFile,"Hull of:\n %s\n", e.print_to_string().c_str());
+      fprintf(DebugFile," %s\n", g.print_to_string().c_str());
+      fprintf(DebugFile,"is " coef_fmt "\n\n",hull);
+    }
+    return true;
+  }
+
+  bool eq(const EQ_Handle &e1, const EQ_Handle &e2) {
+    int sign = 0;
+    for(Constr_Vars_Iter cvi(e1, false); cvi; cvi++) {
+      coef_t c1 = (*cvi).coef;
+      coef_t c2 = e2.get_coef((*cvi).var);
+      if (sign == 0) sign = (c1*c2>=0?1:-1);
+      if (sign*c1 != c2) return false;
+    }
+    assert(sign != 0);
+    {
+      for(Constr_Vars_Iter cvi(e2, false); cvi; cvi++) {
+        coef_t c1 = e1.get_coef((*cvi).var);
+        coef_t c2 = (*cvi).coef;
+        if (sign*c1 != c2) return false;
+      }
+    }
+    return sign * e1.get_const() == e2.get_const();
+  }
+}
+
+
+// This function is deprecated!!!
+Relation QuickHull(Relation &R) {
+  Tuple<Relation> Rs(1);
+  Rs[1] = R;
+  return QuickHull(Rs);
+}
+
+
+// This function is deprecated!!!
+Relation QuickHull(Tuple<Relation> &Rs) {
+  assert(!Rs.empty());
+
+  // if (Rs.size() == 1) return Rs[1];
+
+  Tuple<Relation> l_Rs;
+  for (int i = 1; i <= Rs.size(); i++)
+    for (DNF_Iterator c(Rs[i].query_DNF()); c; c++) {
+      Relation r = Relation(Rs[i], c.curr());
+      l_Rs.append(Approximate(r));
+    }
+    
+  if (l_Rs.size() == 1)
+    return l_Rs[1];
+  
+  Relation result = Relation::True(Rs[1]);
+  result.copy_names(Rs[1]);
+
+  use_ugly_names++; 
+
+  if (hull_debug > 1) 
+    for (int i = 1; i <= l_Rs.size(); i++) {
+      fprintf(DebugFile,"#%d \n",i);
+      l_Rs[i].prefix_print(DebugFile);
+    }
+
+  
+//   Relation R = copy(Rs[1]);
+//   for (int i = 2; i <= Rs.size(); i++) 
+//     R = Union(R,copy(Rs[i]));
+
+// #if 0
+//   if (!R.is_set()) {
+//     if (R.n_inp() == R.n_out()) {
+//       Relation AC = DeltasToRelation(Hull(Deltas(copy(R),
+//                                                  min(R.n_inp(),R.n_out()))),
+//                                      R.n_inp(),R.n_out());
+//       Relation dH = Hull(Domain(copy(R)),false);
+//       Relation rH = Hull(Range(copy(R)),false);
+//       result = Intersection(AC,Cross_Product(dH,rH)); 
+//     }
+//     else {
+//       Relation dH = Hull(Domain(copy(R)),false);
+//       Relation rH = Hull(Range(copy(R)),false);
+//       result = Cross_Product(dH,rH); 
+//       assert(Must_Be_Subset(copy(R),copy(result)));
+//     }
+//   }
+
+// #endif
+ 
+  Conjunct *first;
+  l_Rs[1] = EQs_to_GEQs(l_Rs[1]);
+  first = l_Rs[1].single_conjunct();
+  for (GEQ_Iterator candidate(first->GEQs()); candidate.live(); candidate.next()) {
+    coef_t maxConstantTerm = (*candidate).get_const();
+    bool found = true; 
+    if (hull_debug > 1) {
+      fprintf(DebugFile,"searching for bound on:\n %s\n", (*candidate).print_to_string().c_str());
+    }
+    for (int i = 2; i <= l_Rs.size(); i++) {
+      Conjunct *other = l_Rs[i].single_conjunct();
+      bool found_for_i = false;
+      for (GEQ_Iterator target(other->GEQs()); target.live(); target.next()) {
+        if (hull_debug > 2) {
+          fprintf(DebugFile,"candidate:\n %s\n", (*candidate).print_to_string().c_str());
+          fprintf(DebugFile,"target:\n %s\n", (*target).print_to_string().c_str());
+        }
+        if (parallel(*candidate,*target)) {
+          if (hull_debug > 1)
+            fprintf(DebugFile,"Found bound:\n %s\n", (*target).print_to_string().c_str());
+          maxConstantTerm = max(maxConstantTerm,(*target).get_const());
+          found_for_i = true;
+          break;
+        }
+      }
+      if (!found_for_i) {
+        for (EQ_Iterator target_e(other->EQs()); target_e.live(); target_e.next()) {
+          coef_t h;
+          if (hull(*target_e,*candidate,h)) {
+            if (hull_debug > 1)
+              fprintf(DebugFile,"Found bound of " coef_fmt ":\n %s\n", h, (*target_e).print_to_string().c_str());
+            maxConstantTerm = max(maxConstantTerm,h);
+            found_for_i = true;
+            break;
+          }
+        };
+        if (!found_for_i) {
+          if (hull_debug > 1) {
+            fprintf(DebugFile,"No bound found in:\n");
+            fprintf(DebugFile, "%s", l_Rs[i].print_with_subs_to_string().c_str());
+          }
+          //if nothing found 
+          found = false;
+          break;
+        }
+      }
+    }
+   
+    if (found) {
+      GEQ_Handle  h = result.and_with_GEQ();
+      copy_constraint(h,*candidate);
+      if (hull_debug > 1)
+        fprintf(DebugFile,"Setting constant term to " coef_fmt " in\n %s\n", maxConstantTerm, h.print_to_string().c_str());
+      h.update_const(maxConstantTerm - (*candidate).get_const());
+      if (hull_debug > 1)
+        fprintf(DebugFile,"Updated constraint is\n %s\n", h.print_to_string().c_str());
+    }
+  }
+
+
+  for (EQ_Iterator candidate_eq(first->EQs()); candidate_eq.live(); candidate_eq.next()) {
+    bool found = true;
+    for (int i = 2; i <= l_Rs.size(); i++) {
+      Conjunct *C = l_Rs[i].single_conjunct();
+      bool found_for_i = false;
+
+      for (EQ_Iterator target(C->EQs()); target.live(); target.next()) {
+        if (eq(*candidate_eq,*target)) {
+          found_for_i = true;
+          break;
+        }
+      }
+      if (!found_for_i) {
+        //if nothing found 
+        found = false;
+        break;
+      }
+    }
+   
+    if (found) {
+      EQ_Handle  h = result.and_with_EQ();
+      copy_constraint(h,*candidate_eq);
+      if (hull_debug > 1)
+        fprintf(DebugFile,"Adding eq constraint: %s\n", h.print_to_string().c_str());
+    }
+  }
+
+  use_ugly_names--;
+  if (hull_debug > 1) {
+    fprintf(DebugFile,"quick hull is of:");
+    result.print_with_subs(DebugFile);
+  }
+  return result;
+}
+
+
+// Relation Hull2(Tuple<Relation> &Rs, Tuple<int> &active) {
+//   assert(Rs.size() == active.size() && Rs.size() > 0);
+
+//   Tuple<Relation> l_Rs;
+//   for (int i = 1; i <= Rs.size(); i++)
+//     if (active[i])
+//       l_Rs.append(copy(Rs[i]));
+
+//   if (l_Rs.size() == 0)
+//     return Relation::False(Rs[1]);
+    
+//   try {
+//     Relation r = l_Rs[1];
+//     for (int i = 2; i <= l_Rs.size(); i++) {
+//       r = Union(r, copy(l_Rs[i]));
+//       r.simplify();
+//     }
+
+//     // Relation F = Farkas(r, Basic_Farkas, true);
+//     // if (farkasDifficulty >= 500)
+//     //   throw std::overflow_error("loop convex hull too complicated.");
+//     // F = Farkas(F, Convex_Combination_Farkas, true);
+//     return Farkas(Farkas(r, Basic_Farkas, true), Convex_Combination_Farkas, true);
+//   }
+//   catch (std::overflow_error) {
+//     return QuickHull(l_Rs);
+//   }
+// }
+  
+
+namespace {
+  void printRs(Tuple<Relation> &Rs) {
+    fprintf(DebugFile,"Rs:\n");
+    for (int i = 1; i <= Rs.size(); i++)
+      fprintf(DebugFile,"#%d : %s\n",i,
+              Rs[i].print_with_subs_to_string().c_str());
+  }
+}
+
+Relation BetterHull(Tuple<Relation> &Rs, bool stridesAllowed, bool checkSubsets,
+                    NOT_CONST Relation &input_knownHull = Relation::Null()) {
+  Relation knownHull = consume_and_regurgitate(input_knownHull);
+  static int OMEGA_WHINGE = -1;
+  if (OMEGA_WHINGE < 0) {
+    OMEGA_WHINGE = getenv("OMEGA_WHINGE") ? atoi(getenv("OMEGA_WHINGE")) : 0;
+  }
+  assert(!Rs.empty());
+  if (Rs.size() == 1) {
+    if (stridesAllowed) return Rs[1];
+    else return Approximate(Rs[1]);
+  }
+
+  if (checkSubsets) {
+    Tuple<bool> live(Rs.size());
+    if (hull_debug) {
+      fprintf(DebugFile,"Checking subsets in hull computation:\n");
+      printRs(Rs);
+    }
+    int i;
+    for(i=1;i <=Rs.size(); i++) live[i] = true;
+    for(i=1;i <=Rs.size(); i++) 
+      for(int j=1;j <=Rs.size(); j++) if (i != j && live[j]) {
+          if (hull_debug) fprintf(DebugFile,"checking %d Is_Obvious_Subset %d\n",i,j);
+          if (Is_Obvious_Subset(copy(Rs[i]),copy(Rs[j]))) {
+            if (hull_debug) fprintf(DebugFile,"yes...\n");
+            live[i] = false;
+            break;
+          }
+        }
+    for(i=1;i <=Rs.size(); i++) if (!live[i]) {
+        if (i < Rs.size()) {
+          Rs[i] = Rs[Rs.size()];
+          live[i] = live[Rs.size()];
+        }
+        Rs[Rs.size()] = Relation();
+        Rs.delete_last();
+        i--;
+      }
+  }
+  Relation hull;
+  if (hull_debug) {
+    fprintf(DebugFile,"Better Hull:\n");
+    printRs(Rs);
+    fprintf(DebugFile,"known hull: %s\n", knownHull.print_with_subs_to_string().c_str());
+  }
+  if (knownHull.is_null()) hull = QuickHull(Rs);
+  else hull = Intersection(QuickHull(Rs),knownHull);
+  // for (int i = 1; i <= Rs.size(); i++)
+  //   hull = RectHull(Union(hull, copy(Rs[i])));
+  // hull = Intersection(hull, knownHull);
+  hull.simplify();
+  if (hull_debug) {
+    fprintf(DebugFile,"quick hull: %s\n", hull.print_with_subs_to_string().c_str());
+  }
+
+  Relation orig = Relation::False(Rs[1]);
+  int i;
+  for (i = 1; i <= Rs.size(); i++) 
+    orig = Union(orig,copy(Rs[i]));
+
+  orig.simplify();
+
+  for (i = 1; i <= Rs.size(); i++) {
+    if (!hull.is_obvious_tautology()) Rs[i] = Gist(Rs[i],copy(hull));
+    Rs[i].simplify();
+    if (Rs[i].is_obvious_tautology()) return hull;
+    if (Rs[i].has_single_conjunct()) {
+      Rs[i] = EQs_to_GEQs(Rs[i]);
+      if (hull_debug) {
+        fprintf(DebugFile,"Checking for hull constraints in:\n  %s\n", Rs[i].print_with_subs_to_string().c_str());
+      }
+      Conjunct *c = Rs[i].single_conjunct();
+      for (GEQ_Iterator g(c->GEQs()); g.live(); g.next()) {
+        Relation tmp = Relation::True(Rs[i]);
+        tmp.and_with_GEQ(*g);
+        if (!Difference(copy(orig),tmp).is_upper_bound_satisfiable()) 
+          hull.and_with_GEQ(*g);
+      }
+      for (EQ_Iterator e(c->EQs()); e.live(); e.next()) {
+        Relation tmp = Relation::True(Rs[i]);
+        tmp.and_with_EQ(*e);
+        if (!Difference(copy(orig),tmp).is_upper_bound_satisfiable()) 
+          hull.and_with_EQ(*e);
+      }
+    }
+  }
+
+  hull = FastTightHull(orig,hull);
+  assert(hull.has_single_conjunct());
+
+  if (stridesAllowed) return hull;
+  else return Approximate(hull);
+
+}
+
+
+
+Relation  Hull(NOT_CONST Relation &S, 
+               bool stridesAllowed,
+               int effort,
+               NOT_CONST Relation &knownHull) {
+  Relation R = consume_and_regurgitate(S);
+  R.simplify(1,0);
+  if (!R.is_upper_bound_satisfiable()) return R;
+  Tuple<Relation> Rs;
+  for (DNF_Iterator c(R.query_DNF()); c.live(); ) {
+    Rs.append(Relation(R,c.curr()));
+    c.next();
+  }
+  if (effort == 1)
+    return BetterHull(Rs,stridesAllowed,false,knownHull);
+  else
+    return QuickHull(Rs);
+}
+
+
+
+Relation Hull(Tuple<Relation> &Rs, 
+              Tuple<int> &validMask, 
+              int effort, 
+              bool stridesAllowed,
+              NOT_CONST Relation &knownHull) {
+  // Use relation of index i only when validMask[i] != 0
+  Tuple<Relation> Rs2;
+  for(int i = 1; i <= Rs.size(); i++) {
+    if (validMask[i]) {
+      Rs[i].simplify();
+      for (DNF_Iterator c(Rs[i].query_DNF()); c.live(); ) {
+        Rs2.append(Relation(Rs[i],c.curr()));
+        c.next();
+      }
+    }
+  }
+  assert(effort == 0 || effort == 1);
+  if (effort == 1)
+    return BetterHull(Rs2,stridesAllowed,true,knownHull);
+  else
+    return QuickHull(Rs2);
+}
+
+
+// This function is deprecated!!!
+Relation CheckForConvexRepresentation(NOT_CONST Relation &R_In) {
+  Relation R = consume_and_regurgitate(R_In);
+  Relation h = Hull(copy(R));
+  if (!Difference(copy(h),copy(R)).is_upper_bound_satisfiable())
+    return h;
+  else
+    return R;
+}
+
+// This function is deprecated!!!
+Relation CheckForConvexPairs(NOT_CONST Relation &S) {
+  Relation R = consume_and_regurgitate(S);
+  Relation hull = FastTightHull(copy(R),Relation::True(R));
+  R.simplify(1,0);
+  if (!R.is_upper_bound_satisfiable() || R.number_of_conjuncts() < 2) return R;
+  Tuple<Relation> Rs;
+  for (DNF_Iterator c(R.query_DNF()); c.live(); ) {
+    Rs.append(Relation(R,c.curr()));
+    c.next();
+  }
+
+  bool *dead = new bool[Rs.size()+1];
+  int i;
+  for(i = 1; i<=Rs.size();i++) dead[i] = false;
+
+  for(i = 1; i<=Rs.size();i++)
+    if (!dead[i]) 
+      for(int j = i+1; j<=Rs.size();j++) if (!dead[j]) {
+          if (hull_debug) {
+            fprintf(DebugFile,"Comparing #%d and %d\n",i,j);
+          }
+          Relation U = Union(copy(Rs[i]),copy(Rs[j]));
+          Relation H_ij = FastTightHull(copy(U),copy(hull));
+          if (!Difference(copy(H_ij),U).is_upper_bound_satisfiable()) {
+            Rs[i] = H_ij;
+            dead[j] = true;
+            if (hull_debug)  {
+              fprintf(DebugFile,"Combined them\n");
+            }
+          }
+        }
+  i = 1;
+  while(i<=Rs.size() && dead[i]) i++;
+  assert(i<=Rs.size());
+  R = Rs[i];
+  i++;
+  for(; i<=Rs.size();i++)
+    if (!dead[i]) 
+      R = Union(R,Rs[i]);
+  delete []dead;
+  return R;
+}
+
+//
+// Supporting functions for ConvexRepresentation
+//
+namespace {
+struct Interval {
+  std::list<std::pair<Relation, Relation> >::iterator pos;
+  coef_t lb;
+  coef_t ub;
+  bool change;
+  coef_t modulo;    
+  Interval(std::list<std::pair<Relation, Relation> >::iterator pos_, coef_t lb_, coef_t ub_):
+    pos(pos_), lb(lb_), ub(ub_) {}
+  friend bool operator<(const Interval &a, const Interval &b);
+};
+
+bool operator<(const Interval &a, const Interval &b) {
+  return a.lb < b.lb;
+}
+
+struct Modulo_Interval {
+  coef_t modulo;
+  coef_t start;
+  coef_t size;
+  Modulo_Interval(coef_t modulo_, coef_t start_, coef_t size_):
+    modulo(modulo_), start(start_), size(size_) {}
+  friend bool operator<(const Interval &a, const Interval &b);
+};
+  
+bool operator<(const Modulo_Interval &a, const Modulo_Interval &b) {
+  if (a.modulo == b.modulo) {
+    if (a.start == b.start)
+      return a.size < b.size;
+    else
+      return a.start < b.start;
+  }
+  else
+    return a.modulo < b.modulo;
+}
+
+void merge_intervals(std::list<Interval> &intervals, coef_t modulo, std::list<std::pair<Relation, Relation> > &Rs, std::list<std::pair<Relation, Relation> >::iterator orig) {
+  // normalize  intervals
+  for (std::list<Interval>::iterator i = intervals.begin(); i != intervals.end(); i++) {
+    (*i).modulo = modulo;
+    (*i).change = false;
+    if ((*i).ub - (*i).lb + 1>= modulo) {
+      (*i).lb = 0;
+      (*i).ub = modulo - 1;
+    } 
+    else if ((*i).ub < 0 || (*i).lb >= modulo) {
+      coef_t range = (*i).ub - (*i).lb;
+      (*i).lb = int_mod((*i).lb, modulo);
+      (*i).ub = (*i).lb + range;
+    }
+  }
+
+  intervals.sort();
+
+  // merge neighboring intervals
+  std::list<Interval>::iterator p = intervals.begin();
+  while (p != intervals.end()) {
+    std::list<Interval>::iterator q = p;
+    q++;
+    while (q != intervals.end()) {
+      if ((*p).ub + 1 >= (*q).lb) {
+        Relation hull = ConvexHull(Union(copy((*(*p).pos).first), copy((*(*q).pos).first)));
+        Relation remainder = Difference(Difference(copy(hull), copy((*(*p).pos).first)), copy((*(*q).pos).first));
+        if (!remainder.is_upper_bound_satisfiable()) {
+          if ((*q).pos == orig)
+            std::swap((*p).pos, (*q).pos);
+          (*(*p).pos).first = hull;
+          (*p).ub = max((*p).ub, (*q).ub);
+          (*p).change = true;
+          Rs.erase((*q).pos);
+          q = intervals.erase(q);
+        }
+        else
+          break;
+      }
+      else
+        break;
+    }
+
+    bool p_moved = false;
+    q = p;
+    q++;
+    while (q != intervals.end()) {
+      if ((*q).ub >= modulo && int_mod((*q).ub, modulo) + 1 >= (*p).lb) {
+        Relation hull = ConvexHull(Union(copy((*(*p).pos).first), copy((*(*q).pos).first)));
+        Relation remainder = Difference(Difference(copy(hull), copy((*(*p).pos).first)), copy((*(*q).pos).first));
+        if (!remainder.is_upper_bound_satisfiable()) {
+          if ((*p).pos == orig)
+            std::swap((*p).pos, (*q).pos);
+          (*(*q).pos).first = hull;
+          coef_t t = (*p).ub - int_mod((*q).ub, modulo);
+          if (t > 0)
+            (*q).ub = (*q).ub + t;              
+          (*q).change = true;
+          Rs.erase((*p).pos);
+          p = intervals.erase(p);
+          p_moved = true;
+          break;
+        }
+        else
+          q++;
+      }
+      else
+        q++;
+    }
+
+    if (!p_moved)
+      p++;
+  }
+
+  // merge by reducing the strengh of modulo
+  std::list<Modulo_Interval> modulo_intervals;
+  coef_t max_distance = modulo/2;
+  for (std::list<Interval>::iterator p = intervals.begin(); p != intervals.end(); p++) {
+    if ((*p).lb >= max_distance)
+      break;
+
+    coef_t size = (*p).ub - (*p).lb;
+      
+    std::list<Interval>::iterator q = p;
+    q++;
+    while (q != intervals.end()) {
+      coef_t distance = (*q).lb - (*p).lb;
+      if (distance > max_distance)
+        break;
+
+      if ((*q).ub - (*q).lb != size || int_mod(modulo, distance) != 0) {
+        q++;
+        continue;
+      }
+
+      int num_reduced = 0;
+      coef_t looking_for = int_mod((*p).lb, distance);
+      for (std::list<Interval>::iterator k = intervals.begin(); k != intervals.end(); k++) {
+        if ((*k).lb == looking_for && (*k).ub - (*k).lb == size) {            
+          num_reduced++;
+          looking_for += distance;
+          if (looking_for >= modulo)
+            break;
+        }            
+        else if ((*k).lb <= looking_for && (*k).ub >= looking_for + size) {
+          looking_for += distance;
+          if (looking_for >= modulo)
+            break;
+        }
+        else if ((*k).lb > looking_for)
+          break;
+      }
+
+      if (looking_for >= modulo && num_reduced > 1)
+        modulo_intervals.push_back(Modulo_Interval(distance, int_mod((*p).lb, distance), size));
+
+      q++;
+    }
+  }
+          
+  modulo_intervals.sort();
+
+  // remove redundant reduced-strength intervals
+  std::list<Modulo_Interval>::iterator p2 = modulo_intervals.begin();
+  while (p2 != modulo_intervals.end()) {
+    std::list<Modulo_Interval>::iterator q2 = p2;
+    q2++;
+    while (q2 != modulo_intervals.end()) {
+      if ((*p2).modulo == (*q2).modulo && (*p2).start == (*q2).start)
+        q2 = modulo_intervals.erase(q2);
+      else if (int_mod((*q2).modulo, (*p2).modulo) == 0 &&
+               (*p2).start == int_mod((*q2).start, (*p2).modulo) &&
+               (*p2).size >= (*q2).size)
+        q2 = modulo_intervals.erase(q2);
+      else
+        q2++;
+    }
+    p2++;
+  }
+                
+  // replace original intervals with new reduced-strength ones
+  for (std::list<Modulo_Interval>::iterator i = modulo_intervals.begin(); i != modulo_intervals.end(); i++) {
+    std::vector<Relation *> candidates;
+    int num_replaced = 0;
+    for (std::list<Interval>::iterator j = intervals.begin(); j != intervals.end(); j++)
+      if (int_mod((*j).modulo, (*i).modulo) == 0 &&
+          (*j).ub - (*j).lb >= (*i).size &&
+          (int_mod((*j).lb, (*i).modulo) == (*i).start ||
+           int_mod((*j).ub, (*i).modulo) == (*i).start + (*i).size)) {
+        candidates.push_back(&((*(*j).pos).first));
+        if (int_mod((*j).lb, (*i).modulo) == (*i).start &&
+            (*j).ub - (*j).lb == (*i).size)
+          num_replaced++;
+      }
+    if (num_replaced <= 1)
+      continue;
+      
+    Relation R = copy(*candidates[0]);
+    for (size_t k = 1; k < candidates.size(); k++)
+      R = Union(R, copy(*candidates[k]));
+    Relation hull = ConvexHull(copy(R));
+    Relation remainder = Difference(copy(hull), copy(R));
+    if (!remainder.is_upper_bound_satisfiable()) {
+      std::list<Interval>::iterator replaced_one = intervals.end();
+      for (std::list<Interval>::iterator j = intervals.begin(); j != intervals.end();)
+        if (int_mod((*j).modulo, (*i).modulo) == 0 &&
+            (*j).ub - (*j).lb >= (*i).size &&
+            (int_mod((*j).lb, (*i).modulo) == (*i).start ||
+             int_mod((*j).ub, (*i).modulo) == (*i).start + (*i).size)) {
+          if (int_mod((*j).lb, (*i).modulo) == (*i).start &&
+              (*j).ub - (*j).lb == (*i).size) {
+            if (replaced_one == intervals.end()) {
+              (*(*j).pos).first = hull;
+              (*j).lb = int_mod((*j).lb, (*i).modulo);
+              (*j).ub = int_mod((*j).ub, (*i).modulo);
+              (*j).modulo = (*i).modulo;
+              (*j).change = true;
+              replaced_one = j;
+              j++;
+            }
+            else {
+              if ((*j).pos == orig) {
+                std::swap((*replaced_one).pos, (*j).pos);
+                (*(*replaced_one).pos).first = (*(*j).pos).first;
+              }
+              Rs.erase((*j).pos);
+              j = intervals.erase(j);
+            }
+          }
+          else {
+            if (int_mod((*j).lb, (*i).modulo) == (*i).start)
+              (*j).lb = (*j).lb + (*i).size + 1;
+            else
+              (*j).ub = (*j).ub - (*i).size - 1;
+            (*j).change = true;
+            j++;
+          }
+        }
+        else
+          j++;
+    }
+  }                
+}    
+} // namespace
+
+
+//
+// Simplify a union of sets/relations to a minimal (may not be
+// optimal) number of convex regions.  It intends to replace
+// CheckForConvexRepresentation and CheckForConvexPairs functions.
+//
+Relation ConvexRepresentation(NOT_CONST Relation &R) {
+  Relation l_R = copy(R);
+  if (!l_R.is_upper_bound_satisfiable() || l_R.number_of_conjuncts() < 2)
+    return R;
+
+  // separate each conjunct into smooth convex region and holes
+  std::list<std::pair<Relation, Relation> > Rs; // pair(smooth region, hole condition)
+  for (DNF_Iterator c(l_R.query_DNF()); c.live(); c++) {
+    Relation r1 = Relation(l_R, c.curr());
+    Relation r2 = Approximate(copy(r1));
+    r1 = Gist(r1, copy(r2));
+    Rs.push_back(std::make_pair(r2, r1));
+  }
+
+  try {
+    bool change = true;
+    while (change) {
+      change = false;
+
+      std::list<std::pair<Relation, Relation> >::iterator i = Rs.begin();
+      while (i != Rs.end()) {
+        // find regions with identical hole conditions to merge
+        {
+          std::list<std::pair<Relation, Relation> >::iterator j = i;
+          j++;
+          while (j != Rs.end()) {
+            if (!Difference(copy((*i).second), copy((*j).second)).is_upper_bound_satisfiable() &&
+                !Difference(copy((*j).second), copy((*i).second)).is_upper_bound_satisfiable()) {
+              if (Must_Be_Subset(copy((*j).first), copy((*i).first))) {
+                j = Rs.erase(j);
+              }
+              else if (Must_Be_Subset(copy((*i).first), copy((*j).first))) {
+                (*i).first = (*j).first;
+                j = Rs.erase(j);
+                change = true;
+              }
+              else {
+                Relation r;
+                bool already_use_recthull = false;
+                try {
+                  // chun's debug
+                  // throw std::runtime_error("dfdf");
+                
+                  r = ConvexHull(Union(copy((*i).first), copy((*j).first)));
+                }
+                catch (const std::overflow_error &e) {
+                  r = RectHull(Union(copy((*i).first), copy((*j).first)));
+                  already_use_recthull = true;
+                }
+                retry_recthull:
+                Relation r2 = Difference(Difference(copy(r), copy((*i).first)), copy((*j).first));
+                if (!r2.is_upper_bound_satisfiable()) { // convex hull is tight
+                  (*i).first = r;
+                  j = Rs.erase(j);
+                  change = true;
+                }
+                else {
+                  if (!already_use_recthull) {
+                    r = RectHull(Union(copy((*i).first), copy((*j).first)));
+                    already_use_recthull = true;
+                    goto retry_recthull;
+                  }
+                  else
+                    j++;
+                }
+              }
+            }
+            else
+              j++;
+          }
+        }
+
+        // find identical smooth regions as candidates for hole merge
+        std::list<std::list<std::pair<Relation, Relation> >::iterator> s;
+        for (std::list<std::pair<Relation, Relation> >::iterator j = Rs.begin(); j != Rs.end(); j++) 
+          if (j != i) {
+            if (!Intersection(Difference(copy((*i).first), copy((*j).first)), copy((*j).second)).is_upper_bound_satisfiable() &&
+                !Intersection(Difference(copy((*j).first), copy((*i).first)), copy((*i).second)).is_upper_bound_satisfiable())
+              s.push_back(j);
+          }
+      
+        if (s.size() != 0) {
+          // convert hole condition c1*x1+c2*x2+... = c*alpha+d to a pair of inequalities
+          (*i).second = EQs_to_GEQs((*i).second, false);
+
+          // find potential wildcards that can be used for hole conditions
+          std::set<Variable_ID> nonsingle_wild;
+          for (EQ_Iterator ei((*i).second.single_conjunct()); ei; ei++)
+            if ((*ei).has_wildcards())
+              for (Constr_Vars_Iter cvi(*ei, true); cvi; cvi++)
+                nonsingle_wild.insert(cvi.curr_var());
+          for (GEQ_Iterator gei((*i).second.single_conjunct()); gei; gei++)
+            if ((*gei).has_wildcards()) {
+              Constr_Vars_Iter cvi(*gei, true);
+              Constr_Vars_Iter cvi2 = cvi;
+              cvi2++;
+              if (cvi2) {
+                nonsingle_wild.insert(cvi.curr_var());
+                for (; cvi2; cvi2++)
+                  nonsingle_wild.insert(cvi2.curr_var());
+              }
+            }
+
+          // find hole condition in c*alpha+d1<=c1*x1+c2*x2+...<=c*alpha+d2 format
+          for (GEQ_Iterator gei((*i).second.single_conjunct()); gei; gei++)
+            if ((*gei).has_wildcards()) {
+              coef_t c;
+              Variable_ID v;
+              {
+                Constr_Vars_Iter cvi(*gei, true);
+                v = cvi.curr_var();
+                c = cvi.curr_coef();
+                if (c < 0 || nonsingle_wild.find(v) != nonsingle_wild.end())
+                  continue;
+              }
+          
+              coef_t lb = posInfinity;
+              for (GEQ_Iterator gei2((*i).second.single_conjunct()); gei2; gei2++) {
+                if (!(*gei2 == *gei) && (*gei2).get_coef(v) != 0) {
+                  if (lb != posInfinity) {
+                    nonsingle_wild.insert(v);
+                    break;
+                  }
+                
+                  bool match = true;
+                  for (Constr_Vars_Iter cvi2(*gei); cvi2; cvi2++)
+                    if (cvi2.curr_coef() != -((*gei2).get_coef(cvi2.curr_var()))) {
+                      match = false;
+                      break;
+                    }
+                  if (match)
+                    for (Constr_Vars_Iter cvi2(*gei2); cvi2; cvi2++)
+                      if (cvi2.curr_coef() != -((*gei).get_coef(cvi2.curr_var()))) {
+                        match = false;
+                        break;
+                      }
+                  if (!match) {
+                    nonsingle_wild.insert(v);
+                    break;
+                  }
+
+                  lb = -(*gei2).get_const();
+                }
+              }
+
+              if (nonsingle_wild.find(v) != nonsingle_wild.end())
+                continue;
+          
+              Relation stride_cond = Relation::True((*i).second);
+              F_Exists *f_exists = stride_cond.and_with_and()->add_exists();
+              Variable_ID e = f_exists->declare();
+              F_And *f_root = f_exists->add_and();
+              GEQ_Handle h1 = f_root->add_GEQ();
+              GEQ_Handle h2 = f_root->add_GEQ();
+              for (Constr_Vars_Iter cvi2(*gei); cvi2; cvi2++) {
+                Variable_ID v = cvi2.curr_var();
+                switch (v->kind()) {
+                case Wildcard_Var: 
+                  h1.update_coef(e, cvi2.curr_coef());
+                  h2.update_coef(e, -cvi2.curr_coef());
+                  break;
+                case Global_Var: {
+                  Global_Var_ID g = v->get_global_var();
+                  Variable_ID v2;
+                  if (g->arity() == 0)
+                    v2 = stride_cond.get_local(g);
+                  else
+                    v2 = stride_cond.get_local(g, v->function_of());
+                  h1.update_coef(v2, cvi2.curr_coef());
+                  h2.update_coef(v2, -cvi2.curr_coef());
+                  break;
+                }
+                default:
+                  h1.update_coef(v, cvi2.curr_coef());
+                  h2.update_coef(v, -cvi2.curr_coef());
+                }
+              }
+              h1.update_const((*gei).get_const());
+              h2.update_const(-lb);
+
+              stride_cond.simplify();
+              Relation other_cond = Gist(copy((*i).second), copy(stride_cond));
+
+              // find regions with potential mergeable stride condition with this one
+              std::list<Interval> intervals;
+              intervals.push_back(Interval(i, lb, (*gei).get_const()));
+            
+              for (std::list<std::list<std::pair<Relation, Relation> >::iterator>::iterator j = s.begin(); j != s.end(); j++)
+                if (Must_Be_Subset(copy((**j).second), copy(other_cond))) {
+                  Relation stride_cond2 = Gist(copy((**j).second), copy(other_cond));
+
+                  // interval can be removed
+                  if (stride_cond2.is_obvious_tautology()) {
+                    intervals.push_back(Interval(*j, 0, c-1));
+                    continue;
+                  }
+
+                  stride_cond2 = EQs_to_GEQs(stride_cond2, false);
+                  coef_t lb, ub;
+                  GEQ_Iterator gei2(stride_cond2.single_conjunct());
+                  coef_t sign = 0;
+                  for (Constr_Vars_Iter cvi(*gei2, true); cvi; cvi++)
+                    if (sign != 0) {
+                      sign = 0;
+                      break;
+                    }
+                    else if (cvi.curr_coef() == c)
+                      sign = 1;
+                    else if (cvi.curr_coef() == -c)
+                      sign = -1;
+                    else {
+                      sign = 0;
+                      break;
+                    }                
+                  if (sign == 0)
+                    continue;
+
+                  bool match = true;
+                  for (Constr_Vars_Iter cvi(*gei2); cvi; cvi++) {
+                    Variable_ID v = cvi.curr_var();
+                    if (v->kind() == Wildcard_Var)
+                      continue;
+                    else if (v->kind() == Global_Var) {
+                      Global_Var_ID g = v->get_global_var();
+                      if (g->arity() == 0)
+                        v = (*i).second.get_local(g);
+                      else
+                        v = (*i).second.get_local(g, v->function_of());
+                    }
+                    
+                    if (cvi.curr_coef() != sign * (*gei).get_coef(v)) {
+                      match = false;
+                      break;
+                    }
+                  }
+                  if (!match)
+                    continue;
+                
+                  for (Constr_Vars_Iter cvi(*gei); cvi; cvi++) {
+                    Variable_ID v = cvi.curr_var();
+                    if (v->kind() == Wildcard_Var)
+                      continue;
+                    else if (v->kind() == Global_Var) {
+                      Global_Var_ID g = v->get_global_var();
+                      if (g->arity() == 0)
+                        v = stride_cond2.get_local(g);
+                      else
+                        v = stride_cond2.get_local(g, v->function_of());
+                    }
+                    
+                    if (cvi.curr_coef() != sign * (*gei2).get_coef(v)) {
+                      match = false;
+                      break;
+                    }
+                  }
+                  if (!match)
+                    continue;
+                  if (sign > 0)
+                    ub = (*gei2).get_const();
+                  else
+                    lb = -(*gei2).get_const();                
+                
+                  gei2++;
+                  if (!gei2)
+                    continue;
+
+                  coef_t sign2 = 0;
+                  for (Constr_Vars_Iter cvi(*gei2, true); cvi; cvi++)
+                    if (sign2 != 0) {
+                      sign2 = 0;
+                      break;
+                    }
+                    else if (cvi.curr_coef() == c)
+                      sign2 = 1;
+                    else if (cvi.curr_coef() == -c)
+                      sign2 = -1;
+                    else {
+                      sign2 = 0;
+                      break;
+                    }
+                  if (sign2 != -sign)
+                    continue;
+
+                  for (Constr_Vars_Iter cvi(*gei2); cvi; cvi++) {
+                    Variable_ID v = cvi.curr_var();
+                    if (v->kind() == Wildcard_Var)
+                      continue;
+                    else if (v->kind() == Global_Var) {
+                      Global_Var_ID g = v->get_global_var();
+                      if (g->arity() == 0)
+                        v = (*i).second.get_local(g);
+                      else
+                        v = (*i).second.get_local(g, v->function_of());
+                    }
+                    
+                    if (cvi.curr_coef() != sign2 * (*gei).get_coef(v)) {
+                      match = false;
+                      break;
+                    }
+                  }
+                  if (!match)
+                    continue;
+                
+                  for (Constr_Vars_Iter cvi(*gei); cvi; cvi++) {
+                    Variable_ID v = cvi.curr_var();
+                    if (v->kind() == Wildcard_Var)
+                      continue;
+                    else if (v->kind() == Global_Var) {
+                      Global_Var_ID g = v->get_global_var();
+                      if (g->arity() == 0)
+                        v = stride_cond2.get_local(g);
+                      else
+                        v = stride_cond2.get_local(g, v->function_of());
+                    }
+                    
+                    if (cvi.curr_coef() != sign2 * (*gei2).get_coef(v)) {
+                      match = false;
+                      break;
+                    }
+                  }
+                  if (!match)
+                    continue;
+                  if (sign2 > 0)
+                    ub = (*gei2).get_const();
+                  else
+                    lb = -(*gei2).get_const();
+                
+                  gei2++;
+                  if (gei2)
+                    continue;
+                    
+                  intervals.push_back(Interval(*j, lb, ub));                
+                }
+
+              merge_intervals(intervals, c, Rs, i);
+
+              // make current region the last one being updated
+              bool invalid = false;
+              for (std::list<Interval>::iterator ii = intervals.begin(); ii != intervals.end(); ii++)
+                if ((*ii).change && (*ii).pos == i) {
+                  invalid = true;
+                  intervals.push_back(*ii);
+                  intervals.erase(ii);
+                  break;
+                }
+
+              // update hole condition for each region
+              for (std::list<Interval>::iterator ii = intervals.begin(); ii != intervals.end(); ii++)
+                if ((*ii).change) {
+                  change = true;
+
+                  if ((*ii).ub - (*ii).lb + 1 >= (*ii).modulo)
+                    (*(*ii).pos).second = copy(other_cond);
+                  else {
+                    Relation stride_cond = Relation::True((*i).second);
+                    F_Exists *f_exists = stride_cond.and_with_and()->add_exists();
+                    Variable_ID e = f_exists->declare();
+                    F_And *f_root = f_exists->add_and();
+                    GEQ_Handle h1 = f_root->add_GEQ();
+                    GEQ_Handle h2 = f_root->add_GEQ();
+                    for (Constr_Vars_Iter cvi2(*gei); cvi2; cvi2++) {
+                      Variable_ID v = cvi2.curr_var();
+                      switch (v->kind()) {
+                      case Wildcard_Var: 
+                        h1.update_coef(e, (*ii).modulo);
+                        h2.update_coef(e, -(*ii).modulo);
+                        break;
+                      case Global_Var: {
+                        Global_Var_ID g = v->get_global_var();
+                        Variable_ID v2;
+                        if (g->arity() == 0)
+                          v2 = stride_cond.get_local(g);
+                        else
+                          v2 = stride_cond.get_local(g, v->function_of());
+                        h1.update_coef(v2, cvi2.curr_coef());
+                        h2.update_coef(v2, -cvi2.curr_coef());
+                        break;
+                      }
+                      default:
+                        h1.update_coef(v, cvi2.curr_coef());
+                        h2.update_coef(v, -cvi2.curr_coef());
+                      }
+                    }
+                    h1.update_const((*ii).ub);
+                    h2.update_const(-(*ii).lb);
+
+                    (*(*ii).pos).second = Intersection(copy(other_cond), stride_cond);
+                    (*(*ii).pos).second.simplify();
+                  }
+                }
+            
+              if (invalid)
+                break;
+            }
+        }      
+        i++;
+      }
+    }
+  }
+  catch (const presburger_error &e) {
+    throw e;
+  }
+
+  Relation R2 = Relation::False(l_R);
+  for (std::list<std::pair<Relation, Relation> >::iterator i = Rs.begin(); i != Rs.end(); i++)
+    R2 = Union(R2, Intersection((*i).first, (*i).second));
+  R2.simplify(0, 1);
+  
+  return R2;
+}
+
+
+//
+// Use gist and value range to calculate a quick rectangular hull. It
+// intends to replace all hull calculations (QuickHull, BetterHull,
+// FastTightHull) beyond the method of ConvexHull (dual
+// representations). In the future, it will support max(...)-like
+// upper bound. So RectHull complements ConvexHull in two ways: first
+// for relations that ConvexHull gets too complicated, second for
+// relations where different conjuncts have different symbolic upper
+// bounds.
+// 
+Relation RectHull(NOT_CONST Relation &Rel) {
+  Relation R = Approximate(consume_and_regurgitate(Rel));
+  if (!R.is_upper_bound_satisfiable())
+    return R;
+  if (R.has_single_conjunct())
+    return R;
+
+  std::vector<std::string> input_names(R.n_inp());
+  for (int i = 1; i <= R.n_inp(); i++)
+    input_names[i-1] = R.input_var(i)->name();
+  std::vector<std::string> output_names(R.n_out());
+  for (int i = 1; i <= R.n_out(); i++)
+    output_names[i-1] = R.output_var(i)->name();
+  
+  DNF_Iterator c(R.query_DNF());
+  Relation r = Relation(R, c.curr());
+  c++;
+  std::vector<std::pair<coef_t, coef_t> > bounds1(R.n_inp());
+  std::vector<std::pair<coef_t, coef_t> > bounds2(R.n_out());
+  {
+    Relation t = Project_Sym(copy(r));
+    t.simplify();
+    for (int i = 1; i <= R.n_inp(); i++) {
+      Tuple<Variable_ID> v;
+      for (int j = 1; j <= R.n_inp(); j++)
+        if (j != i)
+          v.append(r.input_var(j));
+      for (int j = 1; j <= R.n_out(); j++)
+        v.append(r.output_var(j));
+      Relation t2 = Project(copy(t), v);
+      t2.query_variable_bounds(t2.input_var(i), bounds1[i-1].first, bounds1[i-1].second);
+    }
+    for (int i = 1; i <= R.n_out(); i++) {
+      Tuple<Variable_ID> v;
+      for (int j = 1; j <= R.n_out(); j++)
+        if (j != i)
+          v.append(r.output_var(j));
+      for (int j = 1; j <= R.n_inp(); j++)
+        v.append(r.input_var(j));
+      Relation t2 = Project(copy(t), v);
+      t2.query_variable_bounds(t2.output_var(i), bounds2[i-1].first, bounds2[i-1].second);
+    }
+  }
+    
+  while (c.live()) {
+    Relation r2 = Relation(R, c.curr());
+    c++;
+    Relation x = Gist(copy(r), Gist(copy(r), copy(r2), 1), 1);
+    if (Difference(copy(r2), copy(x)).is_upper_bound_satisfiable())
+      x = Relation::True(R);
+    Relation y = Gist(copy(r2), Gist(copy(r2), copy(r), 1), 1);
+    if (Difference(copy(r), copy(y)).is_upper_bound_satisfiable())
+      y = Relation::True(R);
+    r = Intersection(x, y);
+    
+    {
+      Relation t = Project_Sym(copy(r2));
+      t.simplify();
+      for (int i = 1; i <= R.n_inp(); i++) {
+        Tuple<Variable_ID> v;
+        for (int j = 1; j <= R.n_inp(); j++)
+          if (j != i)
+            v.append(r2.input_var(j));
+        for (int j = 1; j <= R.n_out(); j++)
+          v.append(r2.output_var(j));
+        Relation t2 = Project(copy(t), v);
+        coef_t lbound, ubound;
+        t2.query_variable_bounds(t2.input_var(i), lbound, ubound);
+        bounds1[i-1].first = min(bounds1[i-1].first, lbound);
+        bounds1[i-1].second = max(bounds1[i-1].second, ubound);
+      }
+      for (int i = 1; i <= R.n_out(); i++) {
+        Tuple<Variable_ID> v;
+        for (int j = 1; j <= R.n_out(); j++)
+          if (j != i)
+            v.append(r2.output_var(j));
+        for (int j = 1; j <= R.n_inp(); j++)
+          v.append(r2.input_var(j));
+        Relation t2 = Project(copy(t), v);
+        coef_t lbound, ubound;
+        t2.query_variable_bounds(t2.output_var(i), lbound, ubound);
+        bounds2[i-1].first = min(bounds2[i-1].first, lbound);
+        bounds2[i-1].second = max(bounds2[i-1].second, ubound);
+      }
+    }
+
+    Relation r3(R.n_inp(), R.n_out());
+    F_And *f_root = r3.add_and();
+    for (int i = 1; i <= R.n_inp(); i++) {
+      if (bounds1[i-1].first != -posInfinity) {
+        GEQ_Handle h = f_root->add_GEQ();
+        h.update_coef(r3.input_var(i), 1);
+        h.update_const(-bounds1[i-1].first);
+      }
+      if (bounds1[i-1].second != posInfinity) {
+        GEQ_Handle h = f_root->add_GEQ();
+        h.update_coef(r3.input_var(i), -1);
+        h.update_const(bounds1[i-1].second);
+      }
+    }
+    for (int i = 1; i <= R.n_out(); i++) {
+      if (bounds2[i-1].first != -posInfinity) {
+        GEQ_Handle h = f_root->add_GEQ();
+        h.update_coef(r3.output_var(i), 1);
+        h.update_const(-bounds2[i-1].first);
+      }
+      if (bounds2[i-1].second != posInfinity) {
+        GEQ_Handle h = f_root->add_GEQ();
+        h.update_coef(r3.output_var(i), -1);
+        h.update_const(bounds2[i-1].second);
+      }
+    }
+    r = Intersection(r, r3);
+    r.simplify();
+  }
+
+  for (int i = 1; i <= r.n_inp(); i++)
+    r.name_input_var(i, input_names[i-1]);
+  for (int i = 1; i <= r.n_out(); i++)
+    r.name_output_var(i, output_names[i-1]);
+  r.setup_names();
+  return r;
+}
+
+} // namespace