//
// ===========
// Rect Domain
// ===========
//
// This represents a domain that looks like an N-dim rectangle,
// with sides parallel to the dimensional axes.  This is optimized
// for a fast "foreach" at the expense of set operations.
//

#include "domains.h"

#define EMPTYREP(rd) do {        \
   rd.p0 = ZERO_POINT;           \
   rd.p1 = ZERO_POINT;           \
   rd.loopStride = NEGONE_POINT; \
  } while (false)

//Is rectdomain obj singular or empty in dimension i?
#define stride_i_zero(obj,i)  (obj.p0.get(i)>=(obj.p1.get(i)-1))

public immutable class tiRectDomain {

  // Constructors
  public tiRectDomain() {
    EMPTYREP(this);
  }

  public tiRectDomain(Point<ARITY> set_p0,Point<ARITY> set_p1) {
    if (!(set_p0 < set_p1)) {
        EMPTYREP(this);
        return;
      }

      p0=set_p0;
      p1=set_p1;
      loopStride=ONES_POINT;
  }

  public tiRectDomain(Point<ARITY> set_p0,
			     Point<ARITY> set_p1,
			     Point<ARITY> set_stride) {
    // check stride is positive

      
    // normalize 
    if (set_p0 < set_p1) { 
      // check stride is positive only for nonempty rectdomains
      CHECK(set_stride>=ZERO_POINT,"stride may not contain a negative Point<N> value.")
      // p0     
      // nothing necessary 
      // p1     
      for (int i = 0; i < set_p1.arity(); i++) {  
	if (set_stride.get(i) == 0) { // singular in i 
	  set_p1 = set_p1.set(i, set_p0.get(i) + 1); 
	  set_stride = set_stride.set(i, 1);
	} else { 
	  int tmp =  tiDomainUtils.align_down((set_p1.get(i) -
					       set_p0.get(i)) - 1,
					      set_stride.get(i)) + 1;
	  set_p1 = set_p1.set(i, set_p0.get(i) + tmp); 
	  if (set_stride.get(i) >= tmp) // singular in i 
	    set_stride = set_stride.set(i, 1);
	} 
      }  
    } else {
      // empty set
      EMPTYREP(this);
      return;
    }
    
    p0 = set_p0;
    p1 = set_p1;
    loopStride = set_stride;
  }

  // THE FOLLOWING TWO CONSTRUCTORS ARE DANGEROUS
//   UNCHECKED ASSERTIONS:
//   (set_p0[i] < set_p1[i] for all i) || ( set_loopStride.get(0)<0 )
//   ^---nonempty                      or  ^----proper empty rep (only for ctor 2)
//   THE FIRST CONSTRUCTOR WILL NOT CONSTRUCT AN EMPTY RECTDOMAIN. DON'T MAKE IT TRY.
//   for the second constructor, if the points would make an empty rectdomain, then caller must ensure set_loopStride == (new tiRectDomain()).stride())
//                               if the RD is to be nonempty, caller must ensure that set_loopStride>=Point.all(1)
//   for the second constructor, the stride must be valid given set_p0 and set_p1 (i.e., set_p1 must be aligned WRT set_p0 and stride)
//                               that is, set_p1[i]=set_p0[i]+1+k*set_loopStride[i] for some k>=0
//                               and if k==0, then set_loopStride[i] must be 1

  inline tiRectDomain(Point<ARITY> set_p0, Point<ARITY> set_p1, boolean bUnused) {
    ASSERTMSG((set_p0<set_p1),"!p0<p1 in tiRD(p,p,b)");
    p0=set_p0;
    p1=set_p1;
    loopStride=ONES_POINT;
  }
  //SEE ABOVE FOR UNCHECKED ASSERTIONS
  inline tiRectDomain(Point<ARITY> set_p0, Point<ARITY> set_p1, Point<ARITY> set_loopStride, boolean bUnused) {
    ASSERT(((set_p0<set_p1)||(set_loopStride.get(0)<0)));
    p0=set_p0;
    p1=set_p1;
    loopStride=set_loopStride;
  }

  public inline tiRectDomain(tiRectDomain set_rd) {
    p0 = set_rd.p0;
    p1 = set_rd.p1;
    loopStride = set_rd.loopStride;
  }
  public tiRectDomain(tiDomain d) {
    if (d.isEmpty()) {
      EMPTYREP(this);
    } else {
      tiRectDomain[] arr=d.isRectangularEx();
      CHECKF((arr!=null),"not a RectDomain in cast.",{EMPTYREP(this); return;})
      tiRectDomain rd=arr[0]; 
      p0 = rd.p0;
      p1 = rd.p1;
      loopStride = rd.loopStride;
    }
  }

  // Set Relations
  public tiDomain op+(tiRectDomain rd) {
    // this + d

    // do union
    if (this <= rd) {
      return rd.promote();
    } else if (rd <= this) {
      return this.promote();
    } else {
      // Can assume this and d are not subsets of each other.  Note that this
      // assumption implies that neither are null sets.

      // information to create result domain
      int expand_left = 0;
      int expand_right = 0;
      int set_stride = -1;
      boolean has_change_dim = false;
      int change_dim;

      // "non encompassing" is where d does not encompass this
      // a "hand" of 0 is left, 1 is middle, 2 is right
      int num_nonencomp_dim = 0;
      int last_nonencomp_dim;
      int last_nonencomp_dim_hand;
      // non overlap is where d does not overlap this
      // a "hand" of 0 is left, 2 is right
      int num_nonoverlap_dim = 0;
      int last_nonoverlap_dim;
      int last_nonoverlap_dim_hand;
      // a dimension where d does not equal this
      int num_unequal_dim = 0;
      int last_unequal_dim;

      //Replacement for the removed stride parameters
      Point<ARITY> stride=getStride();
      Point<ARITY> rdstride=rd.getStride();

      // collect data
      for (int i = 0; i < ARITY; i++) {
	if (p1.get(i) <= rd.p0.get(i)) {
	  // cases:
	  // +- this -+
	  //            +- d -+
	  num_nonoverlap_dim++;
	  last_nonoverlap_dim = i;
	  last_nonoverlap_dim_hand = 2;
	  num_unequal_dim++;
	  last_unequal_dim = i;
	} else if (p0.get(i) >= rd.p1.get(i)) {
	  // cases:
	  //         +- this -+
	  // +- d -+
	  num_nonoverlap_dim++;
	  last_nonoverlap_dim = i;
	  last_nonoverlap_dim_hand = 0;
	  num_unequal_dim++;
	  last_unequal_dim = i;
	} else if ((p0.get(i) >= rd.p0.get(i)) &&
		   (p1.get(i) <= rd.p1.get(i))) {
	  // cases:
	  //   +- this -+
	  // +---- d -----+
	  // - or -
	  // +- this -+
	  // +---- d ----+
	  // - or -
	  //    +- this -+
	  // +---- d ----+
	  // - or -
	  // +- this -+
	  // +-- d ---+
	  if (!((p0.get(i) == rd.p0.get(i)) &&
		(p1.get(i) == rd.p1.get(i)))) {
	    num_unequal_dim++;
	    last_unequal_dim = i;
	  }
	} else {
	  num_unequal_dim++;
	  last_unequal_dim = i;
	  num_nonencomp_dim++;
	  last_nonencomp_dim = i;
	  if ((p0.get(i) < rd.p0.get(i)) && (p1.get(i) > rd.p1.get(i))) {
	    // cases:
	    // +-- this --+
	    //   +- d -+
	    last_nonencomp_dim_hand = 1;
	  } else if ((p0.get(i) < rd.p0.get(i)) &&
		     (p1.get(i) <= rd.p1.get(i))) {
	    // cases:
	    // +- this -+
	    //        +- d -+
	    // - or -
	    // +- this -+
	    //    +- d -+
	    last_nonencomp_dim_hand = 2;
	  } else {
	    // cases:
	    //    +- this -+
	    // +- d -+
	    // - or -
	    // +- this -+
	    // +- d -+
	    last_nonencomp_dim_hand = 0;
	  }
	}
      }
      
      //this block is stride checked
      if ((num_nonoverlap_dim == 1) && (num_unequal_dim == 1)) {
	boolean stridecompatible=true;
	for (int ar=0;ar<ARITY && stridecompatible;ar++) {
	  if (ar!=last_nonoverlap_dim)
	    stridecompatible = stridecompatible && (stride.get(ar)==rdstride.get(ar));
	}
	if (stridecompatible) {
	  if (stride.get(last_nonoverlap_dim) == 0) {
	    if (rdstride.get(last_nonoverlap_dim) == 0) {
	      if (rd.p0.get(last_nonoverlap_dim) < p0.get(last_nonoverlap_dim)) {
		// initiate stride left (stride regeneration)
		expand_left = (p0.get(last_nonoverlap_dim) -
			       rd.p0.get(last_nonoverlap_dim));
		set_stride = expand_left;
	      } else {
		// initiate stride right (stride regeneration)
		expand_right = (rd.p0.get(last_nonoverlap_dim) -
				p0.get(last_nonoverlap_dim));
		set_stride = expand_right;
	      }
	      has_change_dim = true;
	      change_dim = last_nonoverlap_dim;
	    } else {
	    if (p0.get(last_nonoverlap_dim) ==
		(rd.p0.get(last_nonoverlap_dim) -
		 rdstride.get(last_nonoverlap_dim))) {
	      // right append to stride (stride regeneration)
	      expand_right = (rd.p1.get(last_nonoverlap_dim) -
			      1 -
			      rd.p0.get(last_nonoverlap_dim) +
			      rdstride.get(last_nonoverlap_dim));
	      set_stride = rdstride.get(last_nonoverlap_dim);
	      has_change_dim = true;
	      change_dim = last_nonoverlap_dim;
	    } else if (p0.get(last_nonoverlap_dim) ==
		       (rd.p1.get(last_nonoverlap_dim) +
			rdstride.get(last_nonoverlap_dim) - 1)) {
	      // left append to stride (stride regeneration)
	      expand_left = (rd.p1.get(last_nonoverlap_dim) -
			     1 -
			     rd.p0.get(last_nonoverlap_dim) +
			     rdstride.get(last_nonoverlap_dim));
	      set_stride = rdstride.get(last_nonoverlap_dim);
	      has_change_dim = true;
	      change_dim = last_nonoverlap_dim;
	    }
	  }
	} else {
	  if ((rdstride.get(last_nonoverlap_dim) ==
	       stride.get(last_nonoverlap_dim)) ||
	      (rdstride.get(last_nonoverlap_dim) == 0)) {
	    if (last_nonoverlap_dim_hand == 0) {
	      if ((rd.p1.get(last_nonoverlap_dim) - 1) ==
		  (p0.get(last_nonoverlap_dim) -
		   stride.get(last_nonoverlap_dim))) {
		// abut left
		expand_left = (p0.get(last_nonoverlap_dim) -
			       rd.p0.get(last_nonoverlap_dim));
		has_change_dim = true;
		change_dim = last_nonoverlap_dim;
	      }
	    } else if (last_nonoverlap_dim_hand == 2) {
	      if (rd.p0.get(last_nonoverlap_dim) ==
		  (p1.get(last_nonoverlap_dim) -
		   1 +
		   stride.get(last_nonoverlap_dim))) {
		// abut right
		expand_right = (rd.p1.get(last_nonoverlap_dim) -
				p1.get(last_nonoverlap_dim));
		has_change_dim = true;
		change_dim = last_nonoverlap_dim;
	      }
	    }
	  }
	}
	}
      }
      //this block is stride checked
      if ((!has_change_dim) &&
	  (num_nonencomp_dim == 1) &&
	  (num_unequal_dim == 1)) {
	if ((stride.get(last_nonencomp_dim) != 0) &&
	    (rdstride.get(last_nonencomp_dim) ==
	     stride.get(last_nonencomp_dim))) {
	  if (((rd.p0.get(last_nonencomp_dim) -
		p0.get(last_nonencomp_dim)) %
	       stride.get(last_nonencomp_dim)) == 0) {
	     boolean stridecompatible=true;
	    for (int ar=0;ar<ARITY;ar++) {
	      if (ar!=last_nonencomp_dim)
		stridecompatible = stridecompatible && (stride.get(ar)==rdstride.get(ar));
	    }
	    if (stridecompatible) {
	      if (last_nonencomp_dim_hand == 0) {
		// expand left equal stride
		expand_left = (p0.get(last_nonencomp_dim) -
			       rd.p0.get(last_nonencomp_dim));
		has_change_dim = true;
		change_dim = last_nonencomp_dim;
	      } else if (last_nonencomp_dim_hand == 2) {
		// expand right equal stride
		expand_right = (rd.p1.get(last_nonencomp_dim) -
				p1.get(last_nonencomp_dim));
		has_change_dim = true;
		change_dim = last_nonencomp_dim;
	      }
	    }
	  }
	}
      }
      //this block is stride checked
      if ((!has_change_dim) &&
	  (num_nonencomp_dim == 1) &&
	  (num_unequal_dim == 1)) {
	if (stride.get(last_nonencomp_dim) > 0) {
	  if ((rdstride.get(last_nonencomp_dim) == 0)
	      ? false
	      : (((stride.get(last_nonencomp_dim) %
		   rdstride.get(last_nonencomp_dim)) == 0) &&
		 (rdstride.get(last_nonencomp_dim) <
		  stride.get(last_nonencomp_dim)))) {
	    boolean stridecompatible=true;
	    for (int ar=0;ar<ARITY;ar++) {
	      if (ar!=last_nonencomp_dim)
		stridecompatible = stridecompatible && (stride.get(ar)==rdstride.get(ar));
	    }
	    if (stridecompatible) {
	      if (last_nonencomp_dim_hand == 0) {
		if ((rd.p0.get(last_nonencomp_dim) <=
		     p0.get(last_nonencomp_dim)) &&
		    (rd.p1.get(last_nonencomp_dim) ==
		   (p1.get(last_nonencomp_dim) -
		    rdstride.get(last_nonencomp_dim)))) {
		  // overlapped left interpolated points (stride reduction)
		  expand_left = (p0.get(last_nonencomp_dim) -
				 rd.p0.get(last_nonencomp_dim));
		  set_stride = rdstride.get(last_nonencomp_dim);
		  has_change_dim = true;
		  change_dim = last_nonencomp_dim;
		}
	      } else if (last_nonencomp_dim_hand == 2) {
		if ((rd.p1.get(last_nonencomp_dim) >=
		     p1.get(last_nonencomp_dim)) &&
		    (rd.p0.get(last_nonencomp_dim) ==
		     (p0.get(last_nonencomp_dim) +
		      rdstride.get(last_nonencomp_dim)))) {
		  // overlapped right interpolated points (stride reduction)
		  expand_right = (rd.p1.get(last_nonencomp_dim) -
				  p1.get(last_nonencomp_dim));
		  set_stride = rdstride.get(last_nonencomp_dim);
		  has_change_dim = true;
		  change_dim = last_nonencomp_dim;
		}
	      } else {
		if ((rd.p0.get(last_nonencomp_dim) ==
		     (p0.get(last_nonencomp_dim) +
		      rdstride.get(last_nonencomp_dim))) &&
		    (rd.p1.get(last_nonencomp_dim) ==
		     (p1.get(last_nonencomp_dim) -
		      rdstride.get(last_nonencomp_dim)))) {
		  // overlapped middle interpolated points (stride reduction)
		  set_stride = rdstride.get(last_nonencomp_dim);
		  has_change_dim = true;
		  change_dim = last_nonencomp_dim;
		}
	      }
	    }
	  }
	}
      }
      //this block is stride checked
      if ((!has_change_dim) &&
	  (num_unequal_dim == 1)) {
	if (stride.get(last_unequal_dim) > 0) {
	  if ((stride.get(last_unequal_dim) % 2) == 0) {
	    int check_stride = stride.get(last_unequal_dim) / 2;
	    if ((rdstride.get(last_unequal_dim) ==
		stride.get(last_unequal_dim)) ||
		(rdstride.get(last_unequal_dim) == 0)) {
	      boolean stridecompatible=true;
	      for (int ar=0;ar<ARITY;ar++) {
		if (ar!=last_unequal_dim)
		  stridecompatible = stridecompatible && (stride.get(ar)==rdstride.get(ar));
	      }
	      if (stridecompatible) {
		if ((rd.p0.get(last_unequal_dim) ==
		     (p0.get(last_unequal_dim) -
		      check_stride)) &&
		    (rd.p1.get(last_unequal_dim) ==
		     (p1.get(last_unequal_dim) -
		      check_stride))) {
		  // non-overlapped left interpolated points (stride reduction)
		  expand_left = check_stride;
		  set_stride = check_stride;
		  has_change_dim = true;
		  change_dim = last_unequal_dim;
		} else if ((rd.p0.get(last_unequal_dim) ==
			    (p0.get(last_unequal_dim) +
			     check_stride)) &&
			   (rd.p1.get(last_unequal_dim) ==
			    (p1.get(last_unequal_dim) +
			     check_stride))) {
		  // non-overlapped right interpolated points (stride reduction)
		  expand_right = check_stride;
		  set_stride = check_stride;
		  has_change_dim = true;
		  change_dim = last_unequal_dim;
		} else if ((rd.p0.get(last_unequal_dim) ==
			    (p0.get(last_unequal_dim) -
			     check_stride)) &&
			   (rd.p1.get(last_unequal_dim) ==
			    (p1.get(last_unequal_dim) +
			     check_stride))) {
		  // non-overlapped left/right interpolated points (stride reduction)
		  expand_left = check_stride;
		  expand_right = check_stride;
		  set_stride = check_stride;
		  has_change_dim = true;
		  change_dim = last_unequal_dim;
		} else if ((rd.p0.get(last_unequal_dim) ==
			    (p0.get(last_unequal_dim) +
			     check_stride)) &&
			   (rd.p1.get(last_unequal_dim) ==
			    (p1.get(last_unequal_dim) -
			     check_stride))) {
		  // non-overlapped middle interpolated points (stride reduction)
		  set_stride = check_stride;
		  has_change_dim = true;
		  change_dim = last_unequal_dim;
		}
	      }
	    }
	  }
	}
      }
      //this block is stride checked
      if ((!has_change_dim) &&
	  (num_unequal_dim == 1)) {
	if ((stride.get(last_unequal_dim) != 0) &&
	    (rdstride.get(last_unequal_dim) != 0)) {
	  if (((rdstride.get(last_unequal_dim) %
		stride.get(last_unequal_dim)) == 0) &&
	      (rdstride.get(last_unequal_dim) >
		stride.get(last_unequal_dim))) {
	    boolean stridecompatible=true;
	    for (int ar=0;ar<ARITY;ar++) {
	      if (ar!=last_unequal_dim)
		stridecompatible = stridecompatible && (stride.get(ar)==rdstride.get(ar));
	    }
	    if (stridecompatible) {
	      //System.out.println("compatible");
	    if ((rd.p0.get(last_unequal_dim) ==
		 (p0.get(last_unequal_dim) - stride.get(last_unequal_dim))) &&
		(rd.p1.get(last_unequal_dim) <= p1.get(last_unequal_dim))) {
	      // expand left greater stride
	      expand_left = stride.get(last_unequal_dim);
	      has_change_dim = true;
	      change_dim = last_unequal_dim;
	    } else if ((rd.p0.get(last_unequal_dim) >=
			p0.get(last_unequal_dim)) &&
		       (rd.p1.get(last_unequal_dim) ==
			(p1.get(last_unequal_dim) +
			 stride.get(last_unequal_dim)))) {
	      // expand right greater stride
	      expand_right = stride.get(last_unequal_dim);
	      has_change_dim = true;
	      change_dim = last_unequal_dim;
	    } else if ((rd.p0.get(last_unequal_dim) ==
			(p0.get(last_unequal_dim) -
			 stride.get(last_unequal_dim))) &&
		       (rd.p1.get(last_unequal_dim) ==
			(p1.get(last_unequal_dim) +
			 stride.get(last_unequal_dim)))) {
	      // expand left/right greater stride
	      expand_left = stride.get(last_unequal_dim);
	      expand_right = stride.get(last_unequal_dim);
	      has_change_dim = true;
	      change_dim = last_unequal_dim;
	    }
	  }
	  }
	}
      }
      //this block is stride checked
      if ((!has_change_dim) &&
	  (num_unequal_dim == 1)) {
	boolean stridecompatible=true;
	for (int ar=0;ar<ARITY;ar++) {
	  if (ar!=last_unequal_dim)
	    stridecompatible = stridecompatible && (stride.get(ar)==rdstride.get(ar));
	}
	if (stridecompatible) {
	  if ((stride.get(last_unequal_dim) == 0) &&
	      ((rdstride.get(last_unequal_dim) % 2) == 0)) {
	    int check_stride = rdstride.get(last_unequal_dim) / 2;
	    if ((rd.p0.get(last_unequal_dim) ==
		 (p0.get(last_unequal_dim) - check_stride)) &&
		(rd.p1.get(last_unequal_dim) ==
		 (p1.get(last_unequal_dim) + check_stride))) {
	      // initiate stride left/right (stride regeneration)
	      expand_left = check_stride;
	      expand_right = check_stride;
	      set_stride = check_stride;
	      has_change_dim = true;
	      change_dim = last_unequal_dim;
	    }
	    
	  }
	}
      }

      // execute changes
      if (has_change_dim) {
	Point<ARITY> new_p0 = p0;
	Point<ARITY> new_p1 = p1;
	Point<ARITY> new_stride = stride;

	if (set_stride >= 0) {
	  new_stride = new_stride.set(change_dim, set_stride);
	}
	if (expand_left > 0) {
	  new_p0 =
	    new_p0.set(change_dim,
		       (new_p0.get(change_dim) -
			tiDomainUtils.
		          align_up(expand_left,
				   (new_stride.get(change_dim) != 0)
				     ? new_stride.get(change_dim)
				     : 1)));
	}
	if (expand_right > 0) {
	  new_p1 = new_p1.set(change_dim,
			      new_p1.get(change_dim) + expand_right);
	}
	return (new tiRectDomain(new_p0, new_p1, new_stride)).promote();
      } else {
	return promote() + rd;
      }
    }
  }
  public tiDomain op-(tiRectDomain rd) {
    // this - d

    // do difference
    // information to create result domain
    int shrink_left = 0;
    int shrink_right = 0;
    int mult_stride = 1;
    boolean has_change_dim = false;
    int change_dim;

    tiRectDomain di = this * rd;

    if (di.isEmpty()) {
      // sets disjoint
      return this.promote();
    } else if (di == this) {
      // null set
      return (new tiRectDomain()).promote();
    }

    // Replace the removed stride variable
    Point<ARITY> stride=getStride();
    Point<ARITY> distride=di.getStride();

    // "non encompassing" is where di does not encompass this
    // a "hand" of 0 is left, 1 is middle, 2 is right
    int num_nonencomp_dim = 0;
    int last_nonencomp_dim;
    int last_nonencomp_dim_hand;
    // when a dimension has the same limits, set if stride of di is 2x of this
    int num_equal_dim_stride_di_2x_this = 0;
    int last_equal_dim_stride_di_2x_this = 0;

    // collect data
    for (int i = 0; i < ARITY; i++) {
      if ((p1.get(i) <= di.p0.get(i)) || (p0.get(i) >= di.p1.get(i))) {
	// cases:
	// +- this -+
	//            +- di -+
	// - or -
	//          +- this -+
	// +- di -+
      } else if ((p0.get(i) >= di.p0.get(i)) && (p1.get(i) <= di.p1.get(i))) {
	// cases:
	//   +- this -+
	// +---- di ----+
	// - or -
	// +- this -+
	// +---- di ----+
	// - or -
	//     +- this -+
	// +---- di ----+
	// - or -
	// +- this -+
	// +-- di --+
	if ((p0.get(i) == di.p0.get(i)) && (p1.get(i) == di.p1.get(i))) {
	  if (stride.get(i) != 0) {
	    if (distride.get(i) == (2*(stride.get(i)))) {
	      num_equal_dim_stride_di_2x_this++;
	      last_equal_dim_stride_di_2x_this = i;
	    }
	  }
	}
      } else {
	num_nonencomp_dim++;
	last_nonencomp_dim = i;
	if ((p0.get(i) < di.p0.get(i)) && (p1.get(i) > di.p1.get(i))) {
	  // cases:
	  // +-- this --+
	  //   +- di -+
	  last_nonencomp_dim_hand = 1;
	} else if ((p0.get(i) < di.p0.get(i)) && (p1.get(i) <= di.p1.get(i))) {
	  // cases:
	  // +- this -+
	  //        +- di -+
	  // - or -
	  // +- this -+
	  //   +- di -+
	  last_nonencomp_dim_hand = 2;
	} else {
	  // cases:
	  //    +- this -+
	  // +- di -+
	  // - or -
	  // +- this -+
	  // +- di -+
	  last_nonencomp_dim_hand = 0;
	}
      }
    }

    // interpret collected data
    // this block is stride checked
    if (num_nonencomp_dim == 1) {
      if ((stride.set(last_nonencomp_dim,distride.get(last_nonencomp_dim)))==distride) { //stride compatibility check
	if ((last_nonencomp_dim_hand == 0) &&
	    ((stride.get(last_nonencomp_dim) ==
	      distride.get(last_nonencomp_dim)) ||
	     distride.get(last_nonencomp_dim) == 0)) {
	  // slice off left side
	  shrink_left = (di.p1.get(last_nonencomp_dim) -
			 p0.get(last_nonencomp_dim));
	  has_change_dim = true;
	  change_dim = last_nonencomp_dim;
	} else if ((last_nonencomp_dim_hand == 2) &&
		   ((stride.get(last_nonencomp_dim) ==
		     distride.get(last_nonencomp_dim)) ||
		    distride.get(last_nonencomp_dim) == 0)) {
	  // slice off right side
	  shrink_right = (p1.get(last_nonencomp_dim) -
			  di.p0.get(last_nonencomp_dim));
	  has_change_dim = true;
	  change_dim = last_nonencomp_dim;
	} else if ((last_nonencomp_dim_hand == 1) &&
		   ((stride.get(last_nonencomp_dim) == 0)
		    ? false
		    : (new tiRectDomain
		       (p0.set(last_nonencomp_dim,
			       (p0.get(last_nonencomp_dim) +
				stride.get(last_nonencomp_dim))),
			p1.set(last_nonencomp_dim,
			       (p1.get(last_nonencomp_dim) -
				stride.get(last_nonencomp_dim))),
			stride)) == (di))) {
	  // middle elimination (stride magnification)
	  mult_stride = ((p1.get(last_nonencomp_dim) -
			  p0.get(last_nonencomp_dim) - 1) /
			 stride.get(last_nonencomp_dim));
	  has_change_dim = true;
	  change_dim = last_nonencomp_dim;
	}
      }
    }
    // the relevant subblock here is stride-checked
    if (!has_change_dim) {
      if (num_nonencomp_dim <= 1) {
	boolean has_check_dim = false;
	if (num_nonencomp_dim == 1) {
	  change_dim = last_nonencomp_dim;
	  has_check_dim = true;
	} else {
	  if (num_equal_dim_stride_di_2x_this == 1) {
	    change_dim = last_equal_dim_stride_di_2x_this;
	    has_check_dim = true;
	  }
	}
	if (has_check_dim) {
	  if ((stride.get(change_dim) != 0) &&
	      (distride.get(change_dim) != 0)) {
	    if ((stride.get(change_dim)*2) == distride.get(change_dim)) {
	      if ((di.p1.get(change_dim) == p1.get(change_dim)) ||
		  (di.p1.get(change_dim) ==
		   (p1.get(change_dim) - stride.get(change_dim)))) {
		if ((stride.set(change_dim,distride.get(change_dim)))==distride) { //stride compatibility check
		  if (p0.get(change_dim) == di.p0.get(change_dim)) {
		    // shifted interpolated points (stride magnification)
		    shrink_left = stride.get(change_dim);
		    mult_stride = 2;
		    has_change_dim = true;
		  } else if ((p0.get(change_dim) + stride.get(change_dim)) ==
			     di.p0.get(change_dim)) {
		    // interpolated points (stride magnification)
		    mult_stride = 2;
		    has_change_dim = true;
		  }
		}
	      }
	    }
	  }
	}
      }
    }

    // execute changes
    if (has_change_dim) {
      Point<ARITY> new_p0 = p0;
      Point<ARITY> new_p1 = p1;
      Point<ARITY> new_stride = stride;

      if (shrink_left > 0) {
	new_p0 =
	  new_p0.set(change_dim,
		     (new_p0.get(change_dim) +
		      tiDomainUtils.
		        align_up(shrink_left,
				 (stride.get(change_dim) != 0)
				   ? stride.get(change_dim)
				   : 1)));
      }
      if (shrink_right > 0) {
	new_p1 = new_p1.set(change_dim,
			    new_p1.get(change_dim) - shrink_right);
      }
      if (mult_stride > 1) {
	new_stride = new_stride.set(change_dim,
				    (new_stride.get(change_dim) *
				     mult_stride));
      }
      return (new tiRectDomain(new_p0, new_p1, new_stride)).promote();
    } else {
      return promote() - rd;
    }
  }
  public inline tiRectDomain op*=(tiRectDomain d) { return this * d; }
  public tiRectDomain op*(tiRectDomain rd) {
    // this * d
    // Intersection of a single rectangle with a single rectangle
    
    if ((loopStride == ONES_POINT) && (rd.loopStride == ONES_POINT))
      return unitStrideIntersect(rd);
    Point<ARITY> new_p0;
    Point<ARITY> new_p1;
    Point<ARITY> new_stride;
    Point<ARITY> stride=getStride();
    Point<ARITY> rdstride=rd.getStride();
    if (p0 == rd.p0) {
      // p0 equal
      new_p0 = p0;
      new_p1 = p1.lowerBound(rd.p1);
      new_stride = stride.getLcm(rdstride);
    } else {
      // p0 not equal
      boolean detect_null_set = false;

      for (int i = 0; (i < ARITY) && (!detect_null_set); i++) {
	int start_p0 = p0.get(i);
	int start_d_p0 = rd.p0.get(i);
	int max_start = tiDomainUtils.max(start_p0, start_d_p0);

	if (rdstride.get(i) != 0) {
	  start_d_p0 = (start_d_p0 +
			tiDomainUtils.align_up((max_start - start_d_p0),
					       rdstride.get(i)));
	}
	if (stride.get(i) != 0) {
	  start_p0 = (start_p0 +
		      tiDomainUtils.
		        align_up((tiDomainUtils.max(max_start, start_d_p0) -
				  start_p0),
				 stride.get(i)));
	}

	if ((stride.get(i) == 0) || (rdstride.get(i) == 0)) {
	  if (start_p0 != start_d_p0) {
	    detect_null_set = true;
	  } else {
	    new_p0 = new_p0.set(i, start_p0);
	  }
	} else {
	  boolean found = false;
	  for (int k = 0; (k < stride.get(i)) && (!found); k++) {
	    if (((rdstride.get(i) * k) % stride.get(i)) ==
		((start_p0 - start_d_p0) % stride.get(i))) {
	      found = true;
	      new_p0 = new_p0.set(i, start_d_p0 + k * rdstride.get(i));
	    }
	  }
	  if (!found) {
	    detect_null_set = true;
	  }
	}
      }

      if (detect_null_set) {
	// null set
	return new tiRectDomain();
      } else {
	new_p1 = p1.lowerBound(rd.p1);
	new_stride = stride.getLcm(rdstride);
      }
    }

    return new tiRectDomain(new_p0, new_p1, new_stride);
  }

  private inline tiRectDomain unitStrideIntersect(tiRectDomain other)
    {
      Point<ARITY> new_p0=p0,new_p1=p1;
      int lbCor,ubCor; //lower and upper bounds
      for (int i=0;i<ARITY;i++)
	{
	  lbCor=tiDomainUtils.max(p0.get(i),other.p0.get(i));
	  ubCor=tiDomainUtils.min(p1.get(i),other.p1.get(i));
	  if (ubCor <= lbCor)
	      //EMPTY SET
	      return new tiRectDomain();
	  new_p0=new_p0.set(i,lbCor);
	  new_p1=new_p1.set(i,ubCor);
	}

      return new tiRectDomain(new_p0,new_p1,true); //the new ctor
    }

  //DEPRECATED - use isEmpty
  public inline boolean isNull() {
    // this is the null set
    return (loopStride.get(0) < 0);
  }
  //DEPRECATED - use isNotEmpty
  public inline boolean isNotNull() {
    return !isNull();
  }
  public inline boolean isEmpty() {
    // this is empty
    // WARNING: the code-generator explicitly assumes this representation for empty RD
    // see code-call.cc:421
    return (loopStride.get(0) < 0);
  }
  public inline boolean isNotEmpty() {
    return !isEmpty();
  }

  //isADR() returns true iff the domain is a unit-stride rectangle or union of same
  //      (ADR stands for "all dense rects")
  // Since singular rectangles are still unit stride, we can use loopStride
  public boolean isADR() {
    return isEmpty() || loopStride == ONES_POINT;
  }
  public boolean single op==(tiRectDomain single rd) {
    return this.equals(rd);
  }
  public inline boolean single equals(tiRectDomain single rd) {
    // this == d
    return (boolean single)
           ((isEmpty() && rd.isEmpty()) || (p0 == rd.p0 &&
	    p1 == rd.p1 &&
	    loopStride == rd.loopStride));
  }
  public inline boolean single op!=(tiRectDomain single rd) {
    return !(this == rd);
  }
  public boolean op<=(tiRectDomain rd) { // isSubset
    // this <= d
    
    //The empty set is a subset of any set
    if (isEmpty()) return true;
    
    //The subset must be contained within the bounding box of rd
    if (!(p0 >= rd.p0) || !(p1 <= rd.p1)) return false;
    
    //If we are within the bounding box, and are unit stride, then we are a subset
    if ((loopStride == ONES_POINT) && (rd.loopStride == ONES_POINT)) return true;
    
    //If our strides are offset such that the bottom corners don't line up, we are not a subset
    if (!rd.contains(p0)) return false;
    
    for (int i=0; i < ARITY ; i++) {
      int stride_i= (stride_i_zero(this,i)) ? 0 : loopStride.get(i);
      if ( stride_i_zero(rd,i)) {
	if ( stride_i != 0 )
	  //If the other set is singular in dimension i and we are not, we must not be a subset
	  return false;
      } else if ( (stride_i % rd.loopStride.get(i)) != 0 )
	//If our stride is not divisible by the other's stride, we are not a subset
	return false;
    }
    return true;
  }
  public boolean op<(tiRectDomain rd) { // isStrictSubset
    // this < d
    return (this <= rd) && !(this == rd);
  }
  public boolean op>=(tiRectDomain rd) { // isSuperset
    return rd <= this;
  }
  public boolean op>(tiRectDomain rd) { // isStrictSuperset
    return rd < this;
  }

  // Object relations
  public boolean equals(Object d) {
    return ((d instanceof tiDomain) && (this.equals((tiDomain)d)));
  }

  // Domain relations
  public tiDomain op+(tiDomain d) {
    return promote() + d;
  }
  public tiDomain op*(tiDomain d) {
    return promote() * d;
  }
  public tiDomain op-(tiDomain d) {
    return promote() - d;
  }
  public inline boolean equals(tiDomain d) {
    return (promote()).equals(d);
  }
  public boolean op<=(tiDomain d) { // isSubset
    return promote() <= d;
  }
  public boolean op<(tiDomain d) { // isStrictSubset
    return promote() < d;
  }
  public boolean op>=(tiDomain d) { // isSuperset
    return promote() >= d;
  }
  public boolean op>(tiDomain d) { // isStrictSuperset
    return promote() > d;
  }

  // RectDomain Resizing
  private tiRectDomain resize(int k, int dir, int s) {
    if (isEmpty()) return this;
    int direction;
    Point<ARITY> new_p0, new_p1, new_stride;

    CHECK(s > 0,"s <= 0 in accrete/shrink.")
    CHECK(Math.abs(dir) > 0 && Math.abs(dir) <= ARITY,"bad dir in accrete/shrink.")

    new_p0 = p0;
    new_p1 = p1;
    new_stride = loopStride;
    direction = tiDomainUtils.direction_value(dir);

    if ((s == loopStride.get(direction)) || (stride_i_zero(this,direction)) ||
        ((k+1)*s >= loopStride.get(direction) && loopStride.get(direction) % s == 0)) { // PR668
      if (tiDomainUtils.direction_backward(dir)) {
	new_p0 = p0.set(direction, p0.get(direction) - k*s);
      } else {
	new_p1 = p1.set(direction, p1.get(direction) + k*s);
      }
      new_stride = loopStride.set(direction, s);

      return new tiRectDomain(new_p0, new_p1, new_stride);
    } else {
      // !!! error - result not a Rect domain
      System.out.println("Runtime: (Error) not a RectDomain in accrete/shrink.");
      return new tiRectDomain();
    }
  }

  public tiRectDomain accrete(int k, int dir, int s) {
    CHECK(k >= 0,"k < 0 in accrete.")
    return resize(k, dir, s);
  }
  public tiRectDomain accrete(int k, int dir) {
    CHECK(k >= 0,"k < 0 in accrete.")
    return resize(k, dir, 1);
  }
  public tiRectDomain accrete(int k, Point<ARITY> p) {
    CHECK(k >= 0,"k < 0 in accrete.")
    tiRectDomain result = this;
    for (int i = 1; i <= ARITY; i++) {
      result = result.resize(k, i, p.get(tiDomainUtils.direction_value(i)));
      result = result.resize(k, -i, p.get(tiDomainUtils.direction_value(-i)));
    }
    return result;
  }
  public tiRectDomain accrete(int k) {
    CHECK(k >= 0,"k < 0 in accrete.")
    return accrete(k, ONES_POINT);
  }

  public tiRectDomain shrink(int k, int dir) {
    int dir_v=tiDomainUtils.direction_value(dir);
    return resize(-k,
		  dir,
		  loopStride.get(dir_v));
  }
  public tiRectDomain shrink(int k) {
    tiRectDomain result = this;
    for (int i = 1; i <= ARITY; i++) {
      result = result.resize(-k,
			     i,
			     loopStride.get(tiDomainUtils.direction_value(i)));
      result = result.resize(-k,
			     -i,
			     loopStride.get(tiDomainUtils.direction_value(-i)));
    }
    return result;
  }
  public tiRectDomain border(int k, int dir, int s) {
    if (isEmpty()) return this;
    int direction;
    Point<ARITY> new_p0, new_p1;

    new_p0 = p0;
    new_p1 = p1;
    direction = tiDomainUtils.direction_value(dir);
    CHECK((k >= loopStride.get(direction)),"border requires k >= this.stride()[|dir|]")
    if (tiDomainUtils.direction_backward(dir)) {
      new_p1 = new_p1.set(direction, new_p0.get(direction) - s + k);
      new_p0 = new_p0.set(direction, new_p0.get(direction) - s);
    } else {
      new_p0 = new_p0.set(direction, new_p1.get(direction) + s - k);
      new_p1 = new_p1.set(direction, new_p1.get(direction) + s);
    }

    return new tiRectDomain(new_p0, new_p1, loopStride.set(direction, 1), false);
  }
  public tiRectDomain border(int k, int dir) {
    return border(k, dir, 1);
  }
  public tiRectDomain border(int dir) {
    return border(1, dir, 1);
  }

  // Point Relations
  public tiRectDomain op+(Point<ARITY> p) {
    // this + p
    //translation will not affect validity
    ASSERTMSG((new tiRectDomain(p0+p,p1+p,loopStride))==(new tiRectDomain(p0+p,p1+p,loopStride,true)),"op+(p)");
    return new tiRectDomain(p0 + p, p1 + p, loopStride,true);
  }
  public tiRectDomain op-(Point<ARITY> p) {
    // this - p
    //translation will not affect validity
    ASSERT((new tiRectDomain(p0-p,p1-p,loopStride))==(new tiRectDomain(p0-p,p1-p,loopStride,true)));
    return new tiRectDomain(p0 - p, p1 - p, loopStride,true);
  }
  public tiRectDomain op*(Point<ARITY> p) {
    // this * p
    if (isEmpty()) return this;

    Point<ARITY> new_p0 = p0 * p;
    Point<ARITY> new_p1 = (p1 - ONES_POINT) * p;
    Point<ARITY> new_stride = getStride() * p;

    for (int i = 0; i < ARITY; i++) {
      if (p.get(i) < 0) {
	int saved_p0 = new_p0.get(i);
	new_p0 = new_p0.set(i, new_p1.get(i));
	new_p1 = new_p1.set(i, saved_p0);
	new_stride = new_stride.set(i, -new_stride.get(i));
      }
      //Bench to see if this is faster than just calling the normalizing ctor
      if (new_stride.get(i) == 0)
	new_stride = new_stride.set(i,1);
    }
    new_p1 = new_p1 + ONES_POINT;

    return new tiRectDomain(new_p0, new_p1, new_stride,true);
  }
  public tiRectDomain op/(Point<ARITY> p) {
    // this / p
    if (isEmpty()) return this;

    // TEST THIS
    Point<ARITY> new_p0 = p0 / p;
    Point<ARITY> new_p1 = (p1 - ONES_POINT)/p;
    Point<ARITY> new_stride = loopStride / p;

    for (int i = 0; i < ARITY; i++) {
      if (loopStride.get(i) < p.get(i)) {
	new_stride = new_stride.set(i, 1);
      }
      if (p.get(i) < 0) {
	int saved_p0 = new_p0.get(i);
	new_p0 = new_p0.set(i, new_p1.get(i));
	new_p1 = new_p1.set(i, saved_p0);
	new_stride = new_stride.set(i, -new_stride.get(i));
      }
    }
    new_p1 = new_p1 + ONES_POINT;
    //Not confident in proper behavior here to use the new ctor
    return new tiRectDomain(new_p0, new_p1, new_stride);
  }
  /* op-assign operators */
  public inline tiRectDomain op+=(Point<ARITY> p) { return this + p; }
  public inline tiRectDomain op-=(Point<ARITY> p) { return this - p; }
  public inline tiRectDomain op*=(Point<ARITY> p) { return this * p; }
  public inline tiRectDomain op/=(Point<ARITY> p) { return this / p; }

  public tiRectDomain permute(Point<ARITY> p) {
    //TEST THIS - DO WE NEED GETSTRIDE FOR THE CTOR?
    //I don't think permutation affects validity
    ASSERT((new tiRectDomain(p0.permute(p),p1.permute(p),loopStride.permute(p)))==(new tiRectDomain(p0.permute(p),p1.permute(p),loopStride.permute(p),true)));
    return new tiRectDomain(p0.permute(p),
				   p1.permute(p),
				   loopStride.permute(p),true);
  }
  public tiRectDomain translate(Point<ARITY> p) {
    // this + p
    ASSERT((new tiRectDomain(p0+p,p1+p,loopStride))==(new tiRectDomain(p0+p,p1+p,loopStride,true)));
    return new tiRectDomain(p0 + p, p1 + p, loopStride,true);
  }

  // Shape information
  public boolean contains(Point<ARITY> p) {
    // this contains p
    if ((p >= p0) && (p < p1)) {
      for (int i = 0; i < ARITY; i++) {
	int st = loopStride.get(i);
	if (stride_i_zero(this,i)) {
	  if (p.get(i) != p0.get(i))
	    return false;
	} else if (st != 1) { //skip st==1
	   if (((p.get(i) - p0.get(i)) % st) != 0) {
	    return false;
	  }
	}
      }
      return true;
    }
    return false;
  }
  public inline boolean isRectangular() {
    return true;
  }
  public int size() {
    if (isEmpty()) return 0;
    int side_mult = 1;
    for (int i = 0; i < ARITY; i++) {
       side_mult *= ((p1.get(i) - p0.get(i) - 1)/loopStride.get(i))+1;
    }
    return side_mult;
  }
  public inline Point<ARITY> min() {
    return p0;
  }
  public Point<ARITY> lwb() {
    return p0;
  }
  public inline Point<ARITY> upb() {
    return p1;
  }
  public Point<ARITY> max() {
    return p1 - (ONES_POINT);
  }
  public inline Point<ARITY> getStride() {
      if (isEmpty()) return ZERO_POINT;
      Point<ARITY> retval=loopStride;
      for (int i=0;i<ARITY;i++)
	if (stride_i_zero(this,i))
	  retval=retval.set(i,0);

      return retval;

  }
  public inline Point<ARITY> stride() {
    CHECK(isNotEmpty(),"called stride() on an empty RectDomain")
     return loopStride;
  }
  
  public inline tiRectDomain boundingBox() {
    if (isEmpty()) return this; // The special ctor doesn't catch empty and set stride correctly
    return new tiRectDomain(min(),upb(),true);
  }
  public int getNumSidePoints(int dir) {
    if (stride_i_zero(this,dir)) {
      if (p1.get(dir) <= p0.get(dir)) {
        return 0;
      } else {
        return 1;
      }
    } else {
      return ((p1.get(dir) - 1 - p0.get(dir)) / loopStride.get(dir)) + 1;
    }
  }

  // Dynamic optimizations
  inline private tiDomain promote() {
    // Promote the object to be represented in a more general type.
    // For now, always promote to a MultiRectDomainA.
    return tiMultiRectADomain.toDomain(this);
  }

  // Convert to a string
  public String local toString() {
    if (isEmpty()) {
      return "[" + ZERO_POINT.toString() + ":" + NEGONE_POINT.toString() + ":" + ONES_POINT.toString() + "]";
    } else {
      return "[" + min().toString() + ":" + max().toString() + ":" + stride().toString() + "]";
    }
  }

#if ARITY > 1
  public tiRectDomainM1 slice(int k) {
    k--; // make it zero-based
    if (loopStride == ONES_POINT) {
      Point<ARITYM1> new_p0, new_p1;
      for (int i=0, j=0; i < ARITY; i++) {
        if (i != k) {
          new_p0 = new_p0.set(j, p0.get(i));
          new_p1 = new_p1.set(j, p1.get(i));
          j++;
        }
      }
      //If the dimensions were fine before, they'll be fine now too
      ASSERT(((new tiRectDomainM1(new_p0,new_p1))==(new tiRectDomainM1(new_p0,new_p1,true))));
      return new tiRectDomainM1(new_p0, new_p1,true);
    } else {
      Point<ARITYM1> new_p0, new_p1, new_stride;
      //Point<ARITY> stride=getStride();
      for (int i=0, j=0; i < ARITY; i++) {
        if (i != k) {
          new_p0 = new_p0.set(j, p0.get(i));
          new_p1 = new_p1.set(j, p1.get(i));
          new_stride = new_stride.set(j, loopStride.get(i));
          j++;
        }
      }
      //If the dimensions were valid before, they'll be fine now too
      ASSERT(((new tiRectDomainM1(new_p0,new_p1,new_stride))==(new tiRectDomainM1(new_p0,new_p1,new_stride,true))));
      return new tiRectDomainM1(new_p0, new_p1, new_stride,true);
    }
  }
 #endif
  // Debugging
  /*
  public void print(String desc) {
    Point<ARITY> p = min();
    System.out.print(desc);
    System.out.print("[[");
    for (int i = 0; i < p.arity(); i++) {
      if (i != 0) System.out.print(",");
      System.out.print(p.get(i));
    }
    System.out.print("] : ");
    Point<ARITY> p2 = upb();
    System.out.print("[");
    for (int i = 0; i < p2.arity(); i++) {
      if (i != 0) System.out.print(",");
      System.out.print(p2.get(i));
    }
    System.out.print("] : ");
    Point<ARITY> p3 = getStride();
    System.out.print("[");
    for (int i = 0; i < p3.arity(); i++) {
      if (i != 0) System.out.print(",");
      System.out.print(p3.get(i));
    }
    System.out.println("]]");
  }
  */

  // State
  private Point<ARITY> p0;
  private Point<ARITY> p1;
  private Point<ARITY> loopStride;

  // debugging helpers
  
  private boolean trivialStride() {
    for (int dim = 1; dim <= ARITY; ++dim)
      if (loopStride[dim] != 1)
	return false;
    return true;
  }

  private void println() {
    System.out.println(toString());
  }
}