// Emacs: -*- C++ -*-

//
//	Copyright 1993, Center for Computer Vision and Visualization,
//	University of Florida.  All rights reserved.
//

static char copyright[] = "Copyright 1993, Center for Computer Vision and Visualization,\nUniversity of Florida.  All rights reserved.\n";

static char rcsid[] = "$Id: YoderIPS.c,v 1.28 1994/01/31 16:39:12 thoth Exp $";

// 
// $Log: YoderIPS.c,v $
// Revision 1.28  1994/01/31  16:39:12  thoth
// We now document internal errors when they occur.
//
// Revision 1.27  1993/09/15  12:34:20  thoth
// Pointset behavior controls have moved to IA::
//
// Revision 1.26  93/07/15  11:35:25  thoth
// Add support for iterator virtual method.
// 
// Revision 1.25  93/05/28  15:05:31  thoth
// removed dependency on with-ia-pointset*
// 
// Revision 1.24  93/05/26  16:57:51  thoth
// Copyright Notices
// 
// Revision 1.23  93/05/19  14:45:45  thoth
// LatticeRelation is now in IA::
// 
// Revision 1.22  93/04/29  12:52:59  thoth
// Mondo PointSets are now templatized.
// 
// Revision 1.21  93/04/19  23:27:48  thoth
// Convert BaseIPS to IA_BasePS<IA_IntPoint>
// 
// Revision 1.20  93/04/17  15:19:55  thoth
// IA_IntPoint upgrades (dim, inf, sup)
// 
// Revision 1.19  93/04/16  15:50:25  thoth
// use IA_IntPoint instead of IntPoint.
// 
// Revision 1.18  93/02/12  16:28:02  thoth
// more IA_ prefixes.
// 
// Revision 1.17  93/01/20  11:38:30  thoth
// Yoder hash code is now cached.
// 
// Revision 1.16  92/12/16  14:48:24  thoth
// conversion to IA_ mostly complete
// 
// Revision 1.15  92/12/07  12:05:13  thoth
// new IntPoint index(unsigned) operation for use with "other" iterators.
// 
// Revision 1.14  92/11/15  15:58:29  thoth
// expect a slight performance increase with preallocation of SetStructure.
// 
// Revision 1.13  92/11/11  12:12:00  thoth
// added another restriction utility.
// 
// Revision 1.12  92/11/04  12:02:14  thoth
// new image extension support routines.
// 
// Revision 1.11  92/10/25  20:54:08  thoth
// changed argument passing on set restriction routines.
// 
// Revision 1.10  92/10/22  15:13:05  thoth
// pointset restriction utilities.
// 
// Revision 1.9  92/09/30  10:39:31  thoth
// White and Black holes are now semi-public classes, unfortunately.
// type() system now based on address of static data member.
// WhiteHole is now a family of pointsets.
// contains now accepts const Point<>&.
// infimum & supremum virtuals now inf & sup.
// 
// Revision 1.8  92/09/20  17:37:46  thoth
// fixed spelling error
// 
// Revision 1.7  92/09/20  14:07:37  thoth
// new infimum and supremum methods for point sets.
// new instance of the "comp undefined" bug.
// 
// Revision 1.6  92/09/14  21:53:35  thoth
// Fixed silly bug in comparison routines.
// 
// Revision 1.5  92/09/03  13:36:21  thoth
// hash_self now deals with unsigneds
// 
// Revision 1.4  92/08/26  13:52:56  thoth
// added some comments, hopefully things are a little clearer now.
// meaningful exception types
// sped up index, contains, empty.
// removed old min/max cruft.
// improved the compare operation.
// 
// Revision 1.3  92/08/25  15:32:44  thoth
// new offsetted virtual functions for set translation
// set difference is operator/, not -
// lattice relationship is determinable on Yoders using compare() friend
// 
// Revision 1.2  92/08/24  10:37:47  thoth
// *** empty log message ***
// 
// Revision 1.1  92/08/23  13:30:17  thoth
// Initial revision
// 
// 

#include <stdlib.h>

#include "ia.h"
#include "YoderIPS.h"
#include "MondoPS.h" // for the declaration of class _WhiteHole
#include "SetStructure.h"
#include "YoderIPSIter.h"

#ifdef IA_PRII
#pragma implementation
#endif

char IA_YoderIPS::type_;

// construct a Yoder from a Boxy point set
IA_YoderIPS::IA_YoderIPS(const IA_BoxyIPS& bps)
: IA_BasePS<IA_IntPoint>(bps.dim())
{
    this->slices = new interval[this->nslices=1];
    this->lastbefore = bps.min()[0]-1;

    this->slices[0].last = bps.max()[0];
    this->slices[0].count_last = bps.card();
    if (bps.dim()>1) {
	this->slices[0].shadow = IA_IntPointSet(bps.shadow());
    } else {
	this->slices[0].shadow = IA_WhiteHole(0);
    }
    compute_hash_cache();
}



// aip must be a pointer to the beginning of an array of length len of
// IntPoints.  All IntPoints must have the same dimension.
// aip must be sorted into lexicographical order.
// aip may contain duplicates.
// len must be greater than 0
IA_YoderIPS::IA_YoderIPS(const IA_IntPoint *aip, unsigned len, unsigned offset)
: IA_BasePS<IA_IntPoint>( (len>0) ? aip[0].dim()-offset : 0 )
{
    // cout << "constructing Yoder from " << len << " points "
    // << "of dimen "<<dimen<<'\n' << flush;
    nslices=0;
    if (len==0) {
	slices = 0;
	lastbefore = 0;
	// cout << "maybe not\n" << flush;
	compute_hash_cache();
	return;
    }

    unsigned estimate = 1 + aip[len-1][offset] - aip[0][offset];
    if (estimate> len*2)
	estimate = len*2+1;
    // cout << "allocating "<< estimate <<" intervals\n";
    slices = new interval[estimate];
    lastbefore = aip[0][offset] - 1;
    
    unsigned i=0, j;
    
    if (dimen>1) {
	while (i<len) {
	    j=i;
	    register int	coord=aip[j][offset];
	    for (; i<len && aip[i][offset]==coord; i++)
		;
	    slices[nslices].last = coord;
	    slices[nslices].shadow = IA_YoderIPS(aip+j, i-j, offset+1);
	    nslices++;
	    if (i<len && aip[i][offset] > coord+1) {
		slices[nslices].last = aip[i][offset]-1;
		slices[nslices].shadow = IA_BlackHole(dimen-1);
		// cout << "Black hole dimen " << dimen-1;
		// slices[nslices].shadow.output(cout);
		nslices++;
	    }
	}
	this->canonicalize();
    } else {
	while (i<len) {
	    j=i;
	    for (; i+1<len && aip[i+1][offset] <= aip[i][offset]+1; i++)
		;
	    // cout << '['<<j<<','<<i<<']' << flush;
	    
	    slices[nslices].last = aip[i][offset];
	    slices[nslices].shadow = IA_WhiteHole(0);
	    nslices++;
	    if (i+1<len) {
		slices[nslices].last = aip[i+1][offset] - 1;
		slices[nslices].shadow = IA_BlackHole(0);
		nslices++;
	    }
	    i++;
	}
	this->compute_count_last();
	this->compute_hash_cache();
    }
    // cout << "done\n" << flush;
}



IA_YoderIPS::~IA_YoderIPS()
{
    if (this->slices)
	delete[] this->slices;
}



IA_YoderIPS::IA_YoderIPS(const IA_YoderIPS& ips)
:IA_BasePS<IA_IntPoint>(ips.dimen)
{
    this->lastbefore = ips.lastbefore;
    this->nslices = ips.nslices;

    if (this->nslices)
	this->slices = new interval[this->nslices];
    else
	this->slices = 0;

    for (unsigned i=0; i<this->nslices; i++) {
	this->slices[i].last       = ips.slices[i].last;
	this->slices[i].count_last = ips.slices[i].count_last;
	this->slices[i].shadow     = ips.slices[i].shadow;
    }
    this->hash_cache = ips.hash_cache;
}



IA_YoderIPS& IA_YoderIPS::operator=(const IA_YoderIPS &ips)
{
    if (this->nslices != ips.nslices) {
	delete[] this->slices;
	this->nslices = ips.nslices;
	if (this->nslices)
	    this->slices = new interval[this->nslices];
	else
	    this->slices = 0;
    }
    for (unsigned i=0; i<this->nslices; i++) {
	this->slices[i] = ips.slices[i];
    }
    this->dimen = ips.dimen;
    this->lastbefore = ips.lastbefore;

    compute_hash_cache();
    return *this;
}


//
//
//


int IA_YoderIPS:: equal(IA_BasePS<IA_IntPoint>* ips) const
{
    return ips->type() == type() &&
	*this == *(IA_YoderIPS*)ips;
}

unsigned IA_YoderIPS:: hash_self(unsigned modulus) const
{
    return hash_cache%modulus;
}

int IA_YoderIPS::operator==(const IA_YoderIPS& ips) const
{
    if (this->dimen != ips.dimen ||
	this->lastbefore != ips.lastbefore ||
	this->nslices != ips.nslices)
	return 0;

    for (unsigned i=0; i < this->nslices; i++) {
	if (this->slices[i].last != ips.slices[i].last ||
	    this->slices[i].count_last != ips.slices[i].count_last ||
	    this->slices[i].shadow != ips.slices[i].shadow)
	    return 0;
    }
    return 1;
}

//
//


// compute the index of the cylinder that contains the coord
unsigned IA_YoderIPS::cylinder_index(int coord) const
{
    if (this->nslices<1)
	return ~0;
    if (coord<=this->lastbefore ||
	coord>this->slices[this->nslices-1].last)
	return ~0;
 
    unsigned bottom, middle, top;
    bottom = 0;
    top = this->nslices -1;

    while (bottom<top) {
	middle = (bottom+top)/2;
	if (coord <= this->slices[middle].last) {
	    top = middle;
	} else {
	    bottom = middle+1;
	}
    }

    return bottom;
}

int IA_YoderIPS::contains(const IA_IntPoint &ip) const
{
    unsigned	temp;
    temp = cylinder_index(ip[0]);

    return (temp != ~0) &&
	this->slices[temp].shadow.contains(cdr(ip));
}

unsigned IA_YoderIPS::index(const IA_IntPoint &ip) const
{
    unsigned i = cylinder_index(ip[0]);
    if (i==~0) {
	if (IA::throw_on_index_uncontained)
	    ia_throw(IA::PSET_INDEX_UNCONTAINED, __FILE__,__LINE__);
	return IA::index_of_uncontained;
    }

    unsigned before = i ? slices[i-1].count_last : 0;
    unsigned within = slices[i].shadow.card()
	* (ip[0] - (1+(i ? slices[i-1].last : lastbefore)));
    // The next statement will throw if the point is not contained.
    // This is exactly as it should be.
    unsigned recurse = (dimen>1) ? slices[i].shadow.index(cdr(ip)) : 0;

    return before + within + recurse;
}

IA_IntPoint IA_YoderIPS::index(unsigned idx) const
{
    unsigned bottom = 0;
    unsigned top = nslices-1;

    if (idx>=slices[top].count_last) {
	if (IA::throw_on_index_uncontained)
	    ia_throw(IA::PSET_INDEX_UNCONTAINED, __FILE__,__LINE__);
	// the index was outside the card
	return IA_IntPoint();
    }

    unsigned middle;

    while (bottom<top) {
	middle = (bottom+top)/2;
	// cout << bottom << ".." << top << "["<<slices[middle].count_last<<"\n";
	if (idx<slices[middle].count_last) {
	    top = middle;
	    while (top>bottom && slices[top].shadow.empty())
		top--;
	} else {
	    bottom = middle+1;
	    while (bottom < top && slices[bottom].shadow.empty())
		bottom++;
	}
    }

    idx -= (bottom>0) ? slices[bottom-1].count_last : 0;

    unsigned scard = slices[bottom].shadow.card();

    int	car = idx/scard + last_before(bottom)+1;
    unsigned frag = idx%scard;
    // cout << idx<<","<<bottom<<","<<scard <<","<<car<<","<<frag<<endl;
    return concat(car, slices[bottom].shadow.index(frag));
}

//
//

// lexically minimum point in the pointset
IA_IntPoint IA_YoderIPS::min() const
{
    if (this->nslices<1) {
	ia_throw(IA::PSET_INDEX_UNCONTAINED, __FILE__, __LINE__);
	return IA_IntPoint();
    }
    if (this->dimen>1)
	return concat(this->lastbefore+1 , this->slices[0].shadow.min());
    else
	return IA_IntPoint(this->lastbefore+1);
}

// the lexically maxmimum point in the pointset
IA_IntPoint IA_YoderIPS::max() const
{
    if (this->nslices<1) {
	ia_throw(IA::PSET_INDEX_UNCONTAINED, __FILE__, __LINE__);
	return IA_IntPoint();
    }
    if (this->dimen>1)
	return concat(this->slices[this->nslices -1].last,
		      this->slices[this->nslices -1].shadow.max());
    else
	return IA_IntPoint(this->slices[this->nslices -1].last);
}

// the infimum of a pointset, the coordinate-by-coordinate minimum
IA_IntPoint IA_YoderIPS::inf() const
{
    if (this->nslices<1) {
	ia_throw(IA::PSET_INDEX_UNCONTAINED, __FILE__, __LINE__);
	return IA_IntPoint();
    }

    if (dimen <= 1)
	return lastbefore+1;
    else {
	IA_IntPoint	temp(slices[0].shadow.inf());
	for (unsigned i=1; i<nslices; i++)
	    if (!slices[i].shadow.empty())
		temp = ::inf(temp, slices[i].shadow.inf());
	return concat(lastbefore+1, temp);
    }
}

// the supremum of a pointset, the coordinate-by-coordinate maximum
IA_IntPoint IA_YoderIPS::sup() const
{
    if (this->nslices<1) {
	ia_throw(IA::PSET_INDEX_UNCONTAINED, __FILE__, __LINE__);
	return IA_IntPoint();
    }

    if (dimen <= 1)
	return slices[nslices-1].last;
    else {
	IA_IntPoint	temp(slices[0].shadow.sup());
	for (unsigned i=1; i<nslices; i++) {
	    if (!slices[i].shadow.empty())
		temp = ::sup(temp, slices[i].shadow.sup());
	}
	return concat(slices[nslices-1].last, temp);
    }
}

// is the IA_YoderIPS empty?
int IA_YoderIPS::empty() const
{
    // consider the canonical Yoder which can't have empty pointsets
    // at the beginning or end of the cylinder array
    return this->nslices < 1;
}

int IA_YoderIPS::isboxy() const
{
    // the canonical Yoder with more than one cylinder would not be boxy
    if (this->nslices>1)
	return 0;

    // an empty Yoder is boxy
    if (this->nslices<1)
	return 1;

    // we only have one cylinder.  It must be boxy
    if  (this->slices[0].shadow.type() == IA_BoxyIPS::s_type() ||
	 this->slices[0].shadow.type() == IA_WhiteHole_<IA_IntPoint>::s_type()) {
	// white holes are boxy
	return 1;
    } else {
	// A black hole can't occur (canonical form) and
	// if the shadow is a Yoder and boxy, then IntPointSet is
	// not doing its job.
	return 0;
    }
}

IA_BasePSIter<IA_IntPoint> * IA_YoderIPS::iterator() const
{
    // cast away const
    return new IA_YoderIPSIter((IA_YoderIPS*)this);
}

//
//

// Part of canonicalization that compresses adjacent shadow cylinders
// that have the same shadow into one larger shadow cylinder.
// It also eliminates empty pointsets at the beginning and end.
void IA_YoderIPS::compress_duplicates()
{
    unsigned low=0, high=0, scan;
    
    ////////////////
    // compress duplicates
    while (high<this->nslices && this->slices[high].shadow.empty() )
	high++;
    if (high>0)
	this->lastbefore = this->slices[high-1].last;
    // cout << "squashed " << high << "empty cylinders\n";

    while (high<this->nslices) {
	for (scan=high+1;
	     scan<this->nslices &&
	     this->slices[high].shadow == this->slices[scan].shadow;
	     scan++)
	    ;

	if (low<high) {
	    this->slices[low].shadow = this->slices[high].shadow;
	}
	this->slices[low].last = this->slices[scan-1].last;
	//
	low++;
	high = scan;
    }

    while (low>0 && this->slices[low-1].shadow.empty())
	low--;

    this->nslices = low;
}

// Part of canonicalization that computes the count_last field of each
// shadow cylinder
void IA_YoderIPS::compute_count_last()
{
  unsigned	total=0;
  ////////////////
  // scan through and set the count_last fields
      
  total = slices[0].count_last =
    (slices[0].last - lastbefore) *
      slices[0].shadow.card();
  
  for (unsigned i=1; i<nslices; i++) {
    total += (slices[i].last - slices[i-1].last) *
      slices[i].shadow.card();
    slices[i].count_last = total;
  }
}

void IA_YoderIPS::compute_hash_cache()
{
    hash_cache = this->dimen;
    hash_cache += this->lastbefore;
    for (unsigned i=0; i<this->nslices; i++) {
	hash_cache += this->slices[i].last;
	hash_cache += this->slices[i].shadow.hash_self(~0);
    }
}

// complete canonicalization of a Yoder IPS.  This must be performed
// before a Yoder is handed to the user.
void IA_YoderIPS::canonicalize()
{
  compress_duplicates();
  compute_count_last();
  compute_hash_cache();
}

// Copy a block of intervals from one location to another
void IA_YoderIPS::interval_BlT(interval *source, interval*dest, unsigned count)
{
    for (unsigned i=0; i<count; i++) {
	dest[i].last = source[i].last;
	dest[i].shadow = source[i].shadow;
    }
}

//
//

static void squirgle(IA_YoderIPS::interval* ivl, IA_IntPoint p)
{
    ivl->last += p[0];
    ivl->shadow += cdr(p);
}

// this constructs and returns a copy of the Yoder translated
// by `ip'.
IA_BasePS<IA_IntPoint> * IA_YoderIPS::offsetted(const IA_IntPoint&ip) const
{
    // this could be more efficient
    IA_YoderIPS	*rval = new IA_YoderIPS(*this);
    rval->lastbefore += ip[0];
    for (unsigned i=0; i<rval->nslices; i++) {
#if 0
	// koddam "comp undefined" bug
	squirgle(&rval->slices[i], ip);
#else
	rval->slices[i].last += ip[0];
	rval->slices[i].shadow += cdr(ip);
#endif
    }
    return rval;
}

//
//
//


// This static function unions two non-overlapping Yoder pointsets into
// dest.  It requires that first.slices[first.nslices-1].last <=
// second.lastbefore.
void IA_YoderIPS::union_nonoverlapping(
    const IA_YoderIPS &first, const IA_YoderIPS &second, IA_YoderIPS &dest)

{
    dest.lastbefore = first.lastbefore;
    
    interval_BlT(first.slices, dest.slices, first.nslices);
    dest.nslices = first.nslices;
    
    if (dest.slices[dest.nslices-1].last < second.lastbefore) {
	// there is a gab between the two
	dest.slices[dest.nslices].last = second.lastbefore;
	dest.slices[dest.nslices].shadow = IA_BlackHole(second.dimen -1);
	dest.nslices++;
    }
    
    interval_BlT(second.slices, dest.slices+dest.nslices, second.nslices);
    
    dest.nslices += second.nslices;
}

// Given two pointsets that don't begin at the same first coord,
// this static function copies the beginning part that only one pointset
// covers into dest.  The pointset with the "leader" portion must be
// `first' and `f_idx' is an index into the cylinder array that will be
// adjusted to point to the first cylinder of `first' that overlaps with
// `second'.
// It requires (first.lastbefore < second.lastbefore)
// dest.lastbefore is set to first.lastbefore
// `dest' must have a large enough slices array allocated before the call.
// normal usage has *f_idx==0 and dest.nslices==0
void IA_YoderIPS:: copy_nonoverlapping_leader(
    const IA_YoderIPS &first, unsigned *f_idx, const IA_YoderIPS &second,
    IA_YoderIPS &dest)
{
    // cout << "copying non-overlapping leader from " << first.lastbefore <<
    // " until " << second.lastbefore << "\n";
    dest.lastbefore = first.lastbefore;
    while (first.slices[*f_idx].last <= second.lastbefore) {
	dest.slices[dest.nslices] = first.slices[*f_idx];
	dest.nslices++;
	(*f_idx)++;
    }
    // cout << "ate " << *f_idx << " slices";

    // the simple *f_idx==0 case caused me the MOST difficulty
    if ( *f_idx==0 || first.slices[*f_idx-1].last != second.lastbefore) {
	dest.slices[dest.nslices].last = second.lastbefore;
	dest.slices[dest.nslices].shadow =
	    first.slices[*f_idx].shadow;
	dest.nslices++;
	// cout << " and filled a hole";
    }
    // cout << '\n';
}

// This function takes the cylinders from source.slices[*s_idx] to the
// end of the source's slices array and copies them onto the end of dest.
// *s_idx is set to source.nslices by this operation.
// This operation is usually performed after a call to
// merge_overlapping_cylinders (below).
void IA_YoderIPS:: copy_nonoverlapping_trailer(
    const IA_YoderIPS &source, unsigned *s_idx, IA_YoderIPS &dest)
{
    while (*s_idx < source.nslices) {
	dest.slices[dest.nslices].last = source.slices[*s_idx].last;
	dest.slices[dest.nslices].shadow =
	    source.slices[*s_idx].shadow;
	(*s_idx)++;
	dest.nslices++;
    }
}

// Given two pointsets that don't begin at the same first coord,
// this static function "burns" the nonoverlapping beginning.
// `first' must contain the nonoverlapping leader portion and f_idx
// is incremented to the first cylinder that overlaps with second.
// It requires (first.lastbefore < second.lastbefore)
// dest.lastbefore is set to second.lastbefore
// normal usage has *f_idx==0 and dest.nslices==0
void IA_YoderIPS:: burn_nonoverlapping_leader(
    const IA_YoderIPS &first, unsigned *f_idx, const IA_YoderIPS &second,
    IA_YoderIPS &dest)
{
    dest.lastbefore = second.lastbefore;
    while (first.slices[*f_idx].last <= second.lastbefore) {
	(*f_idx)++;
    }
}

// Given two pointsets that have overlapping portions, this function
// "merges" them into dest. *f_idx is the index of a cylinder of `first'
// that overlaps with the (*s_idx)th cylinder of `second'.  The new
// slices of dest are formed by taking the overlapping portions of the
// slices of first and second and applying `f' to their shadows to
// form the shadow of the new cylinder.  *f_idx and *s_idx are incremented
// in a sort of leapfrog until one of them reaches the end of their
// respective cylinder arrays.
// dest must have a large enough `slices' array allocated before the call.
// The new cylinders are appended to the end of any cylinders that have
// been placed in dest already (dest.nslices).
void IA_YoderIPS:: merge_overlapping_cylinders
  (const IA_YoderIPS &first, unsigned *f_idx,
   const IA_YoderIPS &second, unsigned *s_idx,
   IA_YoderIPS &dest, IA_IntPointSet(*f)(const IA_IntPointSet&,const IA_IntPointSet&))
{
  while ( *f_idx < first.nslices && *s_idx < second.nslices ) {
      // cout << "merging shadows... lasts:" << first.slices[*f_idx].last <<
      // ' ' << second.slices[*s_idx].last << '\n';
    dest.slices[dest.nslices].shadow =
      f(first.slices[*f_idx].shadow, second.slices[*s_idx].shadow);

    if (first.slices[*f_idx].last > second.slices[*s_idx].last) {
      dest.slices[dest.nslices].last = second.slices[*s_idx].last;
      (*s_idx)++;
    } else if (first.slices[*f_idx].last < second.slices[*s_idx].last) {
      dest.slices[dest.nslices].last = first.slices[*f_idx].last;
      (*f_idx)++;
    } else { /* == */
      dest.slices[dest.nslices].last = second.slices[*s_idx].last;
      (*s_idx)++; (*f_idx)++;
    }

    
    dest.nslices++;
  }
}

//
//
//

static IA_IntPointSet tmp_or(const IA_IntPointSet &lhs, const IA_IntPointSet &rhs)
{
  return lhs | rhs;
}

// set union
IA_YoderIPS IA_YoderIPS::operator|(const IA_YoderIPS &that) const
{
    IA_YoderIPS	rval(this->dimen);

    if (this->dimen != that.dimen) {
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return rval;
    }

    // this size is an upper bound.  It may be too loose by 1.
    rval.slices = new interval[ this->nslices+that.nslices+1 ];
    // rval.nslices is already 0

    if (this->slices[this->nslices-1].last <= that.lastbefore) {
      // *this does not overlap that at all
	union_nonoverlapping(*this, that, rval);
    } else if (that.slices[that.nslices-1].last <= this->lastbefore) {
      // *this does not overlap that at all
	union_nonoverlapping(that, *this, rval);
      
    } else {
      unsigned lidx=0, ridx=0;
      
      ////////////////
      // union the two Yoder PointSets
	  
      if (this->lastbefore < that.lastbefore) {
	copy_nonoverlapping_leader(*this, &lidx, that, rval);
      } else if (this->lastbefore > that.lastbefore) {
	copy_nonoverlapping_leader(that, &ridx, *this, rval);
      } else /* == */
	rval.lastbefore = that.lastbefore;
      
      merge_overlapping_cylinders (*this,&lidx, that,&ridx, rval, tmp_or);
      
      copy_nonoverlapping_trailer(*this, &lidx, rval);
      copy_nonoverlapping_trailer(that, &ridx, rval);
    }

    rval.canonicalize();

    return rval;
}

//

static IA_IntPointSet tmp_and(const IA_IntPointSet &lhs, const IA_IntPointSet &rhs)
{
  return lhs & rhs;
}

// set intersection
IA_YoderIPS IA_YoderIPS::operator&(const IA_YoderIPS &that) const
{
    IA_YoderIPS	rval(this->dimen);

    if (this->dimen != that.dimen) {
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return rval;
    }

    // this size is an upper bound.  It may be too loose by 1.
    rval.slices = new interval[ this->nslices+that.nslices+1 ];
    // rval.nslices is already 0

    if (this->slices[this->nslices-1].last <= that.lastbefore ||
	that.slices[that.nslices-1].last <= this->lastbefore) {
	return rval;
    } else {
      unsigned lidx=0, ridx=0;
      
      ////////////////
      // intersect the two Yoder PointSets
	  
      if (this->lastbefore < that.lastbefore) {
	burn_nonoverlapping_leader(*this, &lidx, that, rval);
      } else if (this->lastbefore > that.lastbefore) {
	burn_nonoverlapping_leader(that, &ridx, *this, rval);
      } else /* == */
	rval.lastbefore = that.lastbefore;
      
      merge_overlapping_cylinders (*this,&lidx, that,&ridx, rval, tmp_and);

      // burn the nonoverlapping trailers by ignoring them
    }

    rval.canonicalize();

    return rval;
}



// set restriction stuff, dangerous

void IA_YoderIPS:: burn_nonoverlapping_leader(
    const IA_YoderIPS &first, unsigned *f_idx, const IA_YoderIPS &second,
    IA_YoderIPS &dest, IA_SetStructure *ss_result, IA_SetStructure sub)
{
    unsigned	count=0;
    dest.lastbefore = second.lastbefore;
    while (first.slices[*f_idx].last <= second.lastbefore) {
	count += (first.slices[*f_idx].last - first.last_before(*f_idx)) *
	    (first.slices[*f_idx].shadow.card());
	(*f_idx)++;
    }
    if (*f_idx==0 || first.slices[*f_idx-1].last < second.lastbefore)
	count += (second.lastbefore - first.last_before(*f_idx)) *
	    (first.slices[*f_idx].shadow.card());
    ss_result->add_interval(count, sub);
}


void IA_YoderIPS:: merge_overlapping_cylinders
  (const IA_YoderIPS &first, unsigned *f_idx,
   const IA_YoderIPS &second, unsigned *s_idx,
   IA_YoderIPS &dest, IA_SetStructure * ss_result,
   IA_IntPointSet(*f)(const IA_IntPointSet&,const IA_IntPointSet&, IA_SetStructure *)
   )
{
    while ( *f_idx < first.nslices && *s_idx < second.nslices ) {
	// cout << "merging shadows... lasts:" << first.slices[*f_idx].last <<
	// ' ' << second.slices[*s_idx].last << '\n';
	struct interval	&ivl=dest.slices[dest.nslices];
	IA_SetStructure	ss;
	ivl.shadow = f(first.slices[*f_idx].shadow,
		       second.slices[*s_idx].shadow, &ss);
	
	if (first.slices[*f_idx].last > second.slices[*s_idx].last) {
	    ivl.last = second.slices[*s_idx].last;
	    (*s_idx)++;
	} else if (first.slices[*f_idx].last < second.slices[*s_idx].last) {
	    ivl.last = first.slices[*f_idx].last;
	    (*f_idx)++;
	} else { /* == */
	    ivl.last = second.slices[*s_idx].last;
	    (*s_idx)++; (*f_idx)++;
	}
	
	ss_result->add_interval(ivl.last - dest.last_before(dest.nslices), ss);
	
	dest.nslices++;
    }
}

void IA_YoderIPS:: burn_nonoverlapping_trailer
(const IA_YoderIPS &first, unsigned *f_idx,
 const IA_YoderIPS &dest, IA_SetStructure *ss_result, IA_SetStructure sub)
{
    if (*f_idx>=first.nslices)
	return;

    unsigned count=0;
    count += (first.slices[*f_idx].last - dest.slices[dest.nslices-1].last) *
	first.slices[*f_idx].shadow.card();
    (*f_idx)++;
    while (*f_idx<first.nslices) {
	count += (first.slices[*f_idx].last - first.last_before(*f_idx)) *
	    first.slices[*f_idx].shadow.card();
	(*f_idx)++;
    }

    ss_result->add_interval(count, sub);
}


static IA_IntPointSet tmp_and1(const IA_IntPointSet &lhs, const IA_IntPointSet &rhs,
			    IA_SetStructure *result)
{
    return intersect_with_structure(lhs, rhs, result);
}


IA_YoderIPS
IA_YoderIPS::intersect_with_structure (const IA_YoderIPS &that, IA_SetStructure *result)
    const
{
    IA_YoderIPS	rval(this->dimen);

    if (this->dimen != that.dimen) {
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return rval;
    }

    // this size is an upper bound.  It may be too loose by 1.
    rval.slices = new interval[ this->nslices+that.nslices+1 ];
    // rval.nslices is already 0

    if (this->slices[this->nslices-1].last <= that.lastbefore ||
	that.slices[that.nslices-1].last <= this->lastbefore) {
	result->add_interval(this->card(), IA_SetStructure::FIRST_ONLY);
	return rval;
    } else {
      unsigned lidx=0, ridx=0;
      
      ////////////////
      // intersect the two Yoder PointSets
	  
      if (this->lastbefore < that.lastbefore) {
	burn_nonoverlapping_leader(*this, &lidx, that, rval,
				   result, IA_SetStructure::FIRST_ONLY);
      } else if (this->lastbefore > that.lastbefore) {
	burn_nonoverlapping_leader(that, &ridx, *this, rval);
      } else /* == */
	rval.lastbefore = that.lastbefore;
      
      merge_overlapping_cylinders (*this,&lidx, that,&ridx, rval,
				   result, tmp_and1);

      burn_nonoverlapping_trailer(*this, &lidx, rval,
				  result, IA_SetStructure::FIRST_ONLY);
    }

    rval.canonicalize();

    return rval;
}

static IA_IntPointSet tmp_and2(const IA_IntPointSet &lhs, const IA_IntPointSet &rhs,
			    IA_SetStructure *result)
{
    return intersect_with_dualstruct(lhs, rhs, result);
}

IA_YoderIPS
IA_YoderIPS::intersect_with_dualstruct (const IA_YoderIPS &that, IA_SetStructure *result)
    const
{
    IA_YoderIPS	rval(this->dimen);

    if (this->dimen != that.dimen) {
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return rval;
    }

    // this size is an upper bound.  It may be too loose by 1.
    rval.slices = new interval[ this->nslices+that.nslices+1 ];
    // rval.nslices is already 0
    result->ensure_space( this->nslices+that.nslices+1 );

    if (this->slices[this->nslices-1].last <= that.lastbefore) {
	result->add_interval(this->card(), IA_SetStructure::FIRST_ONLY);
	result->add_interval(this->card(), IA_SetStructure::SECOND_ONLY);
    } else if ( that.slices[that.nslices-1].last <= this->lastbefore) {
	result->add_interval(this->card(), IA_SetStructure::SECOND_ONLY);
	result->add_interval(this->card(), IA_SetStructure::FIRST_ONLY);
    } else {
      unsigned lidx=0, ridx=0;
      
      ////////////////
      // intersect the two Yoder PointSets
	  
      if (this->lastbefore < that.lastbefore) {
	burn_nonoverlapping_leader(*this, &lidx, that, rval,
				   result, IA_SetStructure::FIRST_ONLY);
      } else if (this->lastbefore > that.lastbefore) {
	burn_nonoverlapping_leader(that, &ridx, *this, rval,
				   result, IA_SetStructure::SECOND_ONLY);
      } else /* == */
	rval.lastbefore = that.lastbefore;
      
      merge_overlapping_cylinders (*this,&lidx, that,&ridx, rval,
				   result, tmp_and2);

      burn_nonoverlapping_trailer(*this, &lidx, rval,
				  result, IA_SetStructure::FIRST_ONLY);
      burn_nonoverlapping_trailer(that, &ridx, rval,
				  result, IA_SetStructure::SECOND_ONLY);
    }

    rval.canonicalize();

    return rval;
}

//
//
//

void IA_YoderIPS:: copy_nonoverlapping_leader(
    const IA_YoderIPS &first, unsigned *f_idx, const IA_YoderIPS &second,
    IA_YoderIPS &dest, IA_SetStructure *result, IA_SetStructure sub)
{
    unsigned	count=0;
    // cout << "copying non-overlapping leader from " << first.lastbefore <<
    // " until " << second.lastbefore << "\n";
    dest.lastbefore = first.lastbefore;
    while (first.slices[*f_idx].last <= second.lastbefore) {
	dest.slices[dest.nslices] = first.slices[*f_idx];

	count += (first.last_before(*f_idx+1) - first.last_before(*f_idx)) *
	    first.slices[*f_idx].shadow.card();

	dest.nslices++;
	(*f_idx)++;
    }
    // cout << "ate " << *f_idx << " slices";

    // the simple *f_idx==0 case caused me the MOST difficulty
    if ( *f_idx==0 || first.slices[*f_idx-1].last != second.lastbefore) {
	dest.slices[dest.nslices].last = second.lastbefore;
	dest.slices[dest.nslices].shadow =
	    first.slices[*f_idx].shadow;

	count += (dest.last_before(dest.nslices+1) - dest.last_before(dest.nslices)) * dest.slices[dest.nslices].shadow.card();

	dest.nslices++;
	// cout << " and filled a hole";
    }
    // cout << '\n';

    result->add_interval(count, sub);
}



void IA_YoderIPS:: copy_nonoverlapping_trailer(
    const IA_YoderIPS &source, unsigned *s_idx, IA_YoderIPS &dest,
    IA_SetStructure *result, IA_SetStructure sub)
{
    unsigned	count=0;

    while (*s_idx < source.nslices) {
	dest.slices[dest.nslices].last = source.slices[*s_idx].last;
	dest.slices[dest.nslices].shadow =
	    source.slices[*s_idx].shadow;

	count += (dest.last_before(dest.nslices+1) - dest.last_before(dest.nslices)) * dest.slices[dest.nslices].shadow.card();

	(*s_idx)++;
	dest.nslices++;
    }

    result->add_interval(count, sub);
}

static IA_IntPointSet tmp_or(const IA_IntPointSet &lhs, const IA_IntPointSet &rhs,
			    IA_SetStructure *result)
{
    return union_with_structure(lhs, rhs, result);
}

IA_YoderIPS
IA_YoderIPS::union_with_structure (const IA_YoderIPS &that, IA_SetStructure *result)
    const
{
    IA_YoderIPS	rval(this->dimen);

    if (this->dimen != that.dimen) {
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return rval;
    }

    // this size is an upper bound.  It may be too loose by 1.
    rval.slices = new interval[ this->nslices+that.nslices+1 ];
    // rval.nslices is already 0

    if (this->slices[this->nslices-1].last <= that.lastbefore) {
	union_nonoverlapping(*this, that, rval);
	result->add_interval(this->card(), IA_SetStructure::FIRST_ONLY);
	result->add_interval(that.card(), IA_SetStructure::SECOND_ONLY);
    } else if ( that.slices[that.nslices-1].last <= this->lastbefore) {
	union_nonoverlapping(that, *this, rval);
	result->add_interval(that.card(), IA_SetStructure::SECOND_ONLY);
	result->add_interval(this->card(), IA_SetStructure::FIRST_ONLY);
    } else {
      unsigned lidx=0, ridx=0;
      
      ////////////////
      // union the two Yoder PointSets
	  
      if (this->lastbefore < that.lastbefore) {
	copy_nonoverlapping_leader(*this, &lidx, that, rval,
				   result, IA_SetStructure::FIRST_ONLY);
      } else if (this->lastbefore > that.lastbefore) {
	copy_nonoverlapping_leader(that, &ridx, *this, rval,
				   result, IA_SetStructure::SECOND_ONLY);
      } else /* == */
	rval.lastbefore = that.lastbefore;
      
      merge_overlapping_cylinders (*this,&lidx, that,&ridx, rval,
				   result, tmp_or);

      copy_nonoverlapping_trailer(*this, &lidx, rval,
				  result, IA_SetStructure::FIRST_ONLY);
      copy_nonoverlapping_trailer(that, &ridx, rval,
				  result, IA_SetStructure::SECOND_ONLY);
    }

    rval.canonicalize();

    return rval;
}

//

static IA_IntPointSet tmp_xor(const IA_IntPointSet &lhs, const IA_IntPointSet &rhs)
{
  return lhs ^ rhs;
}

// set symmetric difference
IA_YoderIPS IA_YoderIPS::operator^(const IA_YoderIPS &that) const
{
    IA_YoderIPS	rval(this->dimen);

    if (this->dimen != that.dimen) {
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return rval;
    }

    // this size is an upper bound.  It may be too loose by 1.
    rval.slices = new interval[ this->nslices+that.nslices+1 ];
    // rval.nslices is already 0

    if (this->slices[this->nslices-1].last <= that.lastbefore) {
      // *this does not overlap that at all
	union_nonoverlapping(*this, that, rval);
    } else if (that.slices[that.nslices-1].last <= this->lastbefore) {
      // *this does not overlap that at all
	union_nonoverlapping(that, *this, rval);
      
    } else {
      unsigned lidx=0, ridx=0;
      
      ////////////////
      // union the two Yoder PointSets
      //cout << "lastbefores : " << this->lastbefore << ' ' <<
      // that.lastbefore << '\n';
      if (this->lastbefore < that.lastbefore) {
	copy_nonoverlapping_leader(*this, &lidx, that, rval);
      } else if (this->lastbefore > that.lastbefore) {
	copy_nonoverlapping_leader(that, &ridx, *this, rval);
      } else /* == */
	rval.lastbefore = that.lastbefore;
      
      merge_overlapping_cylinders (*this,&lidx, that,&ridx, rval, tmp_xor);
      
      copy_nonoverlapping_trailer(*this, &lidx, rval);
      copy_nonoverlapping_trailer(that, &ridx, rval);
    }

    rval.canonicalize();

    return rval;
}

//

static IA_IntPointSet tmp_minus(const IA_IntPointSet &lhs, const IA_IntPointSet &rhs)
{
  return lhs / rhs;
}

// set asymmetric difference
IA_YoderIPS IA_YoderIPS::operator/(const IA_YoderIPS &that) const
{
    IA_YoderIPS	rval(this->dimen);

    if (this->dimen != that.dimen) {
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return rval;
    }

    // this size is an upper bound.  It may be too loose by 1.
    rval.slices = new interval[ this->nslices+that.nslices+1 ];
    // rval.nslices is already 0

    if (this->slices[this->nslices-1].last <= that.lastbefore ||
	that.slices[that.nslices-1].last <= this->lastbefore) {
      return *this;
    } else {
      unsigned lidx=0, ridx=0;
      
      ////////////////
      // union the two Yoder PointSets
	  
      if (this->lastbefore < that.lastbefore) {
	copy_nonoverlapping_leader(*this, &lidx, that, rval);
      } else if (this->lastbefore > that.lastbefore) {
	burn_nonoverlapping_leader(that, &ridx, *this, rval);
      } else /* == */
	rval.lastbefore = that.lastbefore;
      
      merge_overlapping_cylinders (*this,&lidx, that,&ridx, rval, tmp_minus);
      
      copy_nonoverlapping_trailer(*this, &lidx, rval);
    }

    rval.canonicalize();

    return rval;
}

//
//
//

// This function can return
// NO_REL, P_SUBSET, P_SUPERSET, or EQUAL
// It will abort if internal inconsistencies would cause it
// to return otherwise.
IA::LatticeRelation
compare(const IA_YoderIPS &lhs, const IA_YoderIPS &rhs)
{
    // the arguments must be canonical

    IA::LatticeRelation	rval=IA::EQUAL, temp;
    unsigned	lidx=0, ridx=0;

    if (lhs.lastbefore<rhs.lastbefore) {
	rval = IA::P_SUPERSET;
	while (lhs.slices[lidx].last <= rhs.lastbefore)
	    lidx++;
    } else if (lhs.lastbefore>rhs.lastbefore) {
	rval = IA::P_SUBSET;
	while (rhs.slices[ridx].last <= lhs.lastbefore)
	    ridx++;
    }

    while (lidx<lhs.nslices && ridx<rhs.nslices) {
	temp = compare(lhs.slices[lidx].shadow, rhs.slices[ridx].shadow);
	if (rval==IA::P_SUBSET) {
	    if (temp==IA::P_SUPERSET ||
		temp==IA::NO_REL)
		return IA::NO_REL;
	} else if (rval==IA::EQUAL) {
	    if (temp==IA::NO_REL)
		return temp;
	    else if (temp==IA::P_SUBSET ||
		     temp==IA::EQUAL ||
		     temp==IA::P_SUPERSET)
		rval=temp;
	    else {
		cout << "weird comparison result: " << rval;
		lhs.slices[lidx].shadow.output(cout);
		cout << " with\n";
		rhs.slices[ridx].shadow.output(cout);
		cout << ".\n";

		IA::internal_error(__FILE__,__LINE__);
	    }
	} else if (rval==IA::P_SUPERSET) {
	    if (temp==IA::P_SUBSET ||
		temp==IA::NO_REL)
		return IA::NO_REL;
	} else {
	    // rval should have been completely screened!
	    // this is impossible
	    IA::internal_error(__FILE__,__LINE__);
	}

	if (lhs.slices[lidx].last > rhs.slices[ridx].last)
	    ridx++;
	else if (lhs.slices[lidx].last < rhs.slices[ridx].last)
	    lidx++;
	else { /* == */
	    ridx++;
	    lidx++;
	}	    
    }

    switch (rval) {
    case IA::P_SUBSET:
	if (lidx < lhs.nslices)
	    return IA::NO_REL;
	else
	    return rval;
	break;
    case IA::EQUAL:
	if (lidx < lhs.nslices)
	    return IA::P_SUPERSET;
	else if (ridx < rhs.nslices)
	    return IA::P_SUBSET;
	else
	    return rval;
	// NOTREACHED
	break;
    case IA::P_SUPERSET:
	if (ridx < rhs.nslices)
	    return IA::NO_REL;
	else
	    return rval;
	break;
    default:
	IA::internal_error(__FILE__,__LINE__);
	break;
    }
}

//
//

void IA_YoderIPS::output(ostream &o, unsigned indent) const
{
  for (unsigned i=0; i<indent; i++)
    o << ' ';
  o << lastbefore << '\n';
  for (i=0; i<nslices; i++) {
    slices[i].shadow.output(o,indent+2);
    for (unsigned j=0; j<indent; j++)
      o << ' ';
    o << slices[i].last << '\n';
  }
}
