IntervalListAlgo.h

// IntervalListAlgo.h
//
// Author David Barrett-Lennard
// (C)opyright Cedanet Pty Ltd 2007

#pragma once
#ifndef Ceda_cxUtils_IntervalListAlgo_H
#define Ceda_cxUtils_IntervalListAlgo_H

#include "IntervalList.h"
#include "Range.h"
#include <cstddef>
#include <algorithm>

namespace ceda
{
/*
Use to collect the output of algorithms such as IntersectIntervalLists, UnionIntervalLists,
SubtractIntervalLists and EraseIntervalLists that write the intervals in increasing order and 
satisfying the requirement that there are nonzero gaps between intervals.
*/

template <typename T>
class IntervalListWriter
{
    cxNotCloneable(IntervalListWriter)
public:
    typedef typename IntervalList<T>::Interval Interval;

    IntervalListWriter(IntervalList<T>& w) : m_w(w) {}
    
    IntervalListWriter& operator<<(Interval i)
    {
        m_w.safe_push_back(i);
        return *this;
    }

private:
    IntervalList<T>& m_w;
};


/*
r1, r2 represent interval sets, expressed as iterator-ranges for ordered lists of half
open intervals.  

Writes half open intervals [t1,t2) from (r1 union r2) to 'out' in left to right order, 
either calling Add1, Add2 or Add3 depending on whether [t1,t2) is part of r1\r2, 
r2\r1 or the intersection of r1,r2.

r1              [---------------)            [---------------)       [---)
r2                      [--------------)           [---)   [----)
                      
r1 union r2     [------)[-------)[-----)     [----)[---)[-)[-)[-)    [---)
                   |        |       |
                   |        |       |
                 r1\r2    r1*r2   r2\r1
*/

template <typename Range1, typename Range2, typename Out>
void ScanTwoIntervalLists(Range1 r1, Range2 r2, Out& out)
{
    // Set to false for slightly higher performance if it can be assumed that consecutive 
    // intervals in r1 are separated by nonzero gaps, and similarly for r2.
    const bool allowGaps = true;

    typedef typename Out::Interval Interval;
    typedef typename Interval::value_type T;

    if (!r1) goto end1;
    if (!r2) goto end2;
    for(;;)
    {
        if (r1->end() <= r2->begin())
        {
            /*
                r1 comes before r2 (possibly adjacent)

                r1:      [-----)
                r2:                 [-----)
            */
            out.Add1(r1->begin(), r1->end());
            ++r1;
            if (!r1) goto end1;
        }
        else if (r2->end() <= r1->begin())
        {
            /*
                r2 comes before r1 (possibly adjacent)

                r1:                 [-----)
                r2:      [-----)
            */
            out.Add2(r2->begin(), r2->end());
            ++r2;
            if (!r2) goto end2;
        }
        else
        {
            T t1;
            if (r1->begin() < r2->begin())
            {
                out.Add1(r1->begin(),r2->begin());
                t1 = r2->begin();
            }
            else 
            {
                if (r2->begin() < r1->begin())
                {
                    out.Add2(r2->begin(),r1->begin());
                }
                t1 = r1->begin();
            }

            if (r2->end() < r1->end()) goto stepw2;
            
            for(;;)
            {
                /*
                Step through r1 coming before r2->end()
                
                    W1:   [----)  [---)  [---)              [-------------)      
                    W2:     [-------------------------------)
                                                          r2->end()
                */
                for(;;)
                {
                    cxAssert(r1->end() <= r2->end());
                    
                    if (t1 < r1->end()) out.Add3(t1,r1->end());
                    t1 = r1->end();
                    
                    if (!++r1)
                    {
                        if (t1 < r2->end()) out.Add2(t1,r2->end());
                        ++r2;
                        goto end1;   // Will add remainder of W2 if any
                    }
                    
                    if (r1->end() > r2->end()) break;
                    
                    if (!allowGaps || t1 < r1->begin())
                    {
                        cxAssert(t1 < r1->begin());
                        out.Add2(t1,r1->begin());
                    }
                    t1 = r1->begin();
                } 
                
                if (r1->begin() > r2->end())
                {
                    /* r1 comes after r2 with a nonzero gap.  
                       Therefore r2->end() marks the end of h

                        W1:                                        [------) 
                        W2:     [-------------------------------)
                                                              r2->end()
                    */
                    if (t1 < r2->end()) out.Add2(t1,r2->end());
                    
                    if (!++r2) goto end2;
                    break;  // Break out of inner loop, continue with outer loop
                }
                else
                {
                    if (!allowGaps || t1 < r1->begin())
                    {
                        cxAssert(t1 < r1->begin());
                        out.Add2(t1,r1->begin());
                    }
                    t1 = r1->begin();
                }
                
            stepw2:
                /*
                Step through r2 coming before r1->end()
                
                    W2:   [----)  [---)  [---)   [------) 
                    W1:     [-------------------------------)
                                                          r1->end()
                */
                for(;;)
                {
                    cxAssert(r2->end() <= r1->end());
                    
                    if (t1 < r2->end()) out.Add3(t1,r2->end());
                    t1 = r2->end();
                    
                    if (!++r2)
                    {
                        if (t1 < r1->end()) out.Add1(t1,r1->end());
                        ++r1;
                        goto end2;   // Will add remainder of W1 if any
                    }
                    
                    if (r2->end() > r1->end()) break;
                    
                    if (!allowGaps || t1 < r2->begin())
                    {
                        cxAssert(t1 < r2->begin());
                        out.Add1(t1,r2->begin());
                    }
                    t1 = r2->begin();
                } 
                
                if (r2->begin() > r1->end())
                {
                    /* r2 comes after r1 with a nonzero gap.  

                        W2:                                        [------) 
                        W1:     [-------------------------------)
                                                              r1->end()
                    */
                    if (t1 < r1->end()) out.Add1(t1,r1->end());
                    
                    if (!++r1) goto end1;
                    break;  // Break out of inner loop, continue with outer loop
                }
                else
                {
                    if (!allowGaps || t1 < r2->begin())
                    {
                        cxAssert(t1 < r2->begin());
                        out.Add1(t1,r2->begin());
                    }
                    t1 = r2->begin();
                }
            }
        }
    }

end1:
    cxAssert(!r1);
    while(r2) 
    {
        out.Add2(r2->begin(), r2->end());
        ++r2;
    }
    return;
    
end2:
    cxAssert(!r2);
    while(r1) 
    {
        out.Add1(r1->begin(), r1->end());
        ++r1;
    }
}

/*
Like ScanTwoIntervalLists except 

    -   for each call to Add1(e2,t1,t2) for an interval [t1,t2) in r1\r2, the summed extent e2 of 
        the intervals in r2 to the left of position t1 is provided.
    
    -   for each call to Add2(e1,t1,t2) for an interval [t1,t2) in r2\r1, the summed extent e1 of 
        the intervals in r1 to the left of position t1 is provided.
        
    -   for each call to Add3(e1,e2,t1,t2) for an interval [t1,t2) in the intersection of r1 and r2,
        the summed extent e1 of the intervals in r1 and the summed extent e2 of the intervals in r2 
        to the left of position t1 is provided.
*/

template <typename Range1, typename Range2, typename Out>
void ScanTwoIntervalListsWithAccumulatedExtents(Range1 r1, Range2 r2, Out& out)
{
    // Set to false for slightly higher performance if it can be assumed that consecutive 
    // intervals in r1 are separated by nonzero gaps, and similarly for r2.
    const bool allowGaps = true;

    typedef typename Out::Interval Interval;
    typedef typename Interval::value_type T;

    T extent1 = 0;
    T extent2 = 0;
    
    if (!r1) goto end1;
    if (!r2) goto end2;
    for(;;)
    {
        if (r1->end() <= r2->begin())
        {
            /*
                r1 comes before r2 (possibly adjacent)

                r1:      [-----)
                r2:                 [-----)
            */
            out.Add1(extent2, r1->begin(), r1->end());
            extent1 += r1->size();
            ++r1;
            if (!r1) goto end1;
        }
        else if (r2->end() <= r1->begin())
        {
            /*
                r2 comes before r1 (possibly adjacent)

                r1:                 [-----)
                r2:      [-----)
            */
            out.Add2(extent1, r2->begin(), r2->end());
            extent2 += r2->size();
            ++r2;
            if (!r2) goto end2;
        }
        else
        {
            T t1;
            if (r1->begin() < r2->begin())
            {
                out.Add1(extent2, r1->begin(), r2->begin());
                t1 = r2->begin();
            }
            else 
            {
                if (r2->begin() < r1->begin())
                {
                    out.Add2(extent1, r2->begin(), r1->begin());
                }
                t1 = r1->begin();
            }

            if (r2->end() < r1->end()) goto stepw2;
            
            for(;;)
            {
                /*
                Step through r1 coming before r2->end()
                
                    W1:   [----)  [---)  [---)              [-------------)      
                    W2:     [-------------------------------)
                                                          r2->end()
                */
                for(;;)
                {
                    cxAssert(r1->end() <= r2->end());
                    
                    if (t1 < r1->end()) out.Add3(extent1, extent2, t1, r1->end());
                    t1 = r1->end();
                    
                    extent1 += r1->size();
                    if (!++r1)
                    {
                        if (t1 < r2->end()) out.Add2(extent1, t1,r2->end());
                        extent2 += r2->size();
                        ++r2;
                        goto end1;   // Will add remainder of W2 if any
                    }
                    
                    if (r1->end() > r2->end()) break;
                    
                    if (!allowGaps || t1 < r1->begin())
                    {
                        cxAssert(t1 < r1->begin());
                        out.Add2(extent1, t1,r1->begin());
                    }
                    t1 = r1->begin();
                } 
                
                if (r1->begin() > r2->end())
                {
                    /* r1 comes after r2 with a nonzero gap.  
                       Therefore r2->end() marks the end of h

                        W1:                                        [------) 
                        W2:     [-------------------------------)
                                                              r2->end()
                    */
                    if (t1 < r2->end()) out.Add2(extent1, t1,r2->end());
                    
                    extent2 += r2->size();
                    if (!++r2) goto end2;
                    break;  // Break out of inner loop, continue with outer loop
                }
                else
                {
                    if (!allowGaps || t1 < r1->begin())
                    {
                        cxAssert(t1 < r1->begin());
                        out.Add2(extent1, t1, r1->begin());
                    }
                    t1 = r1->begin();
                }
                
            stepw2:
                /*
                Step through r2 coming before r1->end()
                
                    W2:   [----)  [---)  [---)   [------) 
                    W1:     [-------------------------------)
                                                          r1->end()
                */
                for(;;)
                {
                    cxAssert(r2->end() <= r1->end());
                    
                    if (t1 < r2->end()) out.Add3(extent1, extent2, t1, r2->end());
                    t1 = r2->end();
                    
                    extent2 += r2->size();
                    if (!++r2)
                    {
                        if (t1 < r1->end()) out.Add1(extent2, t1, r1->end());
                        extent1 += r1->size();
                        ++r1;
                        goto end2;   // Will add remainder of W1 if any
                    }
                    
                    if (r2->end() > r1->end()) break;
                    
                    if (!allowGaps || t1 < r2->begin())
                    {
                        cxAssert(t1 < r2->begin());
                        out.Add1(extent2, t1, r2->begin());
                    }
                    t1 = r2->begin();
                } 
                
                if (r2->begin() > r1->end())
                {
                    /* r2 comes after r1 with a nonzero gap.  

                        W2:                                        [------) 
                        W1:     [-------------------------------)
                                                              r1->end()
                    */
                    if (t1 < r1->end()) out.Add1(extent2, t1, r1->end());
                    
                    extent1 += r1->size();
                    if (!++r1) goto end1;
                    break;  // Break out of inner loop, continue with outer loop
                }
                else
                {
                    if (!allowGaps || t1 < r2->begin())
                    {
                        cxAssert(t1 < r2->begin());
                        out.Add1(extent2, t1, r2->begin());
                    }
                    t1 = r2->begin();
                }
            }
        }
    }

end1:
    cxAssert(!r1);
    while(r2) 
    {
        out.Add2(extent1, r2->begin(), r2->end());
        extent2 += r2->size();
        ++r2;
    }
    return;
    
end2:
    cxAssert(!r2);
    while(r1) 
    {
        out.Add1(extent2, r1->begin(), r1->end());
        extent1 += r1->size();
        ++r1;
    }
}


// r1, r2 represent interval sets, expressed as iterator-ranges for ordered lists of half
// open intervals.  Writes the set theoretic intersection of the two interval sets to 'out' 
// in the manner consistent with an IntervalListWriter<T>.
template <typename Range1, typename Range2, typename Out>
void IntersectIntervalLists(Range1 r1, Range2 r2, Out& out)
{
    /*
    The intersection of two interval sets W1,W2 can be achieved by linear scans through both 
    interval sets, giving an O(n1+n2) algorithm where n1,n2 are the sizes of the W1,W2.
    
    The algorithm checks whether the next interval from one list comes before the next interval
    from the other list.  In that case we increment the iterators in the manner of merge sort.
    This continues until we find a pair of intervals that overlap.  The intersection of two 
    intervals is always a single contiguous interval, and this cannot be adjacent to other 
    intervals in W so these interval intersections can be calculated and pushed onto the result 
    W independently.
   
                       next interval 
                          from W1
                             |
                             |
    W1:                    [----)         [---) [----)    [---------) [--) [-----)
                    
    W2:       [---) [----)          [---)                      [---------------)
                |
                |
          next interval
            from W2
            
     W:                                                        [----) [--) [---)         
    */

    if (!r1 || !r2) return;  // Intersection is empty
    for(;;)
    {
        if (r1->end() <= r2->begin())
        {
            // r1 before r2
            ++r1;
            if (!r1) return;
        }
        else if (r2->end() <= r1->begin())
        {
            // r2 before r1
            ++r2;
            if (!r2) return;
        }
        else
        {
            // Intervals *r1 and *r2 intersect.  Write their nonempty intersection
            out << Intersect(*r1,*r2);

            // Step r1 or r2 depending on which interval is lagging behind
            if (r1->end() < r2->end()) { ++r1; if (!r1) return; }
            else                       { ++r2; if (!r2) return; }
        }
    }
}

// r1, r2 represent interval sets, expressed as iterator-ranges for ordered lists of half
// open intervals.  Writes the set theoretic union of the two interval sets to 'out' 
// in the manner consistent with an IntervalListWriter<T>.
template <typename Range1, typename Range2, typename Out>
void UnionIntervalLists(Range1 r1, Range2 r2, Out& out)
{
    if (!r1 || !r2) goto end;
    for(;;)
    {
        if (r1->end() < r2->begin())
        {
            /*
                r1 comes before r2 with a nonzero gap

                W1:      [-----)
                W2:                 [-----)
                               |    |
                            nonzero gap
            */
            out << *r1;
            if (!++r1) break;
        }
        else if (r2->end() < r1->begin())
        {
            /*
                r2 comes before r1 with a nonzero gap

                W1:                 [-----)
                W2:      [-----)
                               |    |
                            nonzero gap
            */
            out << *r2;
            if (!++r2) break;
        }
        else
        {
            /*
            There is no nonzero gap between r1 and r2.  Therefore they are either adjacent or else 
            they overlap.  Either way they will coalesce in the union.
            Scan forwards through intervals finding a single contiguous interval h to be added to 
            the union W.  We stop when we reach the end of W1 or W2 or else we have a gap.
            E.g
            
            W1:    [-----r1-----)  [--) [--)  [---------------------)           [---)
            
            W2:      [-r2-)  [------------------) [--) [--)   [---)      [---)
            
            W:     [------------------------------------------------)
                                         h                          |    |
                                                              gap causes h to finish
                                                              
            Most generally there is a need to alternate between two loops in a leap frog effect.  
            In one loop we step through r2's until they go past r1->end(), and in the other loop we 
            step through r1's until they go past r2->end().
            */
            
            typename Out::Interval h;
            h.setbegin(std::min(r1->begin(), r2->begin()));
            
            if (r2->end() < r1->end()) goto stepw2;
            
            for(;;)
            {
                /*
                Step through r1 coming before r2->end()
                
                    W1:   [----)  [---)  [---)       [------) 
                    W2:     [-------------------------------)
                                                          r2->end()
                */
                do
                {
                    cxAssert(r1->end() <= r2->end());
                    if (!++r1)
                    {
                        // Reached end of W1, so r2->end() marks the end of h
                        h.setend(r2->end());
                        out << h;
                        ++r2;
                        goto end;   // Will add remainder of W2 if any
                    }
                } while(r1->end() <= r2->end());
                if (r1->begin() > r2->end())
                {
                    /* r1 comes after r2 with a nonzero gap.  Therefore r2->end() marks the end of h

                        W1:                                        [------) 
                        W2:     [-------------------------------)
                                                              r2->end()
                    */
                    h.setend(r2->end());
                    out << h;
                    if (!++r2) goto end;
                    break;  // Break out of inner loop, continue with outer loop
                }
            stepw2:
                /*
                Step through r2 coming before r1->end()
                
                    W2:   [----)  [---)  [---)   [------) 
                    W1:     [-------------------------------)
                                                          r1->end()
                */
                do 
                {
                    cxAssert(r2->end() <= r1->end());
                    if (!++r2)
                    {
                        h.setend(r1->end());
                        out << h;
                        ++r1;
                        goto end;   // Will add remainder of W1 if any
                    }
                } while(r2->end() <= r1->end());
                if (r2->begin() > r1->end())
                {
                    /* r2 comes after r1 with a nonzero gap.  Therefore r1->end() marks the end of h

                        W2:                                        [------) 
                        W1:     [-------------------------------)
                                                              r1->end()
                    */
                    h.setend(r1->end());
                    out << h;
                    if (!++r1) goto end;
                    break;  // Break out of inner loop, continue with outer loop
                }
            }
        }
    }
end:    
    cxAssert(!r1 || !r2);
    while(r1) out << *r1++;
    while(r2) out << *r2++;
}

// r1, r2 represent interval sets, expressed as iterator-ranges for ordered lists of half
// open intervals.  Writes the set theoretic difference r1\r2 of the two interval sets to 'out' 
// in the manner consistent with an IntervalListWriter<T>.
template <typename Range1, typename Range2, typename Out>
void SubtractIntervalLists(Range1 r1, Range2 r2, Out& out)
{
    using Interval = typename Out::Interval;
    using T = typename Interval::value_type;
    
    /*
    W = W1 \ W2
    
    If intervals don't overlap then we can scan through the intervals of W1, W2, stepping
    iterator w1 or w2 depending on which interval comes next - in the manner of merge sort.  
    Intervals from W1 are added to W whereas intervals from W2 are otherwise ignored.
    
        W1:     [------) [---)    [---)        [---)
        W2:                              [---)
        W:      [------) [---)    [---)        [---)
    
    If an interval of W1 intersects one or more intervals from W2 then we tend to keep 
    the gaps between the intervals of W2.
    
        W1:     [-----------)
        W2:   [---) [--)  [----)
        W:        [-)  [--)
    */
    
    if (!r1 || !r2) goto end;
    for(;;)
    {
        if (r1->end() <= r2->begin())
        {
            /*
                r1 comes before r2 (possibly adjacent)

                W1:      [-----)
                W2:                 [-----)
            */
            out << *r1;
            if (!++r1) return;
        }
        else if (r2->end() <= r1->begin())
        {
            /*
                r2 comes before r1 (possibly adjacent)

                W1:                 [-----)
                W2:      [-----)
            */
            if (!++r2) break;
        }
        else
        {
            /*
            r1 and r2 overlap - i.e. they have a non-empty intersection.
            We need to work out what parts of r1 to keep (if any).

            W1:     [-----------)
            W2:   [---) [--)  [----)
            W:        [-)  [--)
            */
            
            if (r2->begin() > r1->begin())
            {
                /*
                    There is a prefix of W1 to be added to W
                    
                    W1:     [-----------------------)
                    W2:         [---)...
                    W:      [---)    ...
                */
                out << Interval(r1->begin(),r2->begin());
            }
            
            while(r2->end() < r1->end())
            {
                // There are necessarily one or more gaps to be added to W
                T gapStart = r2->end();
                if (!++r2)
                {
                    /*
                        W1:     [-----------r1-----------)
                        W2:   [------r2-------)     
                                              |          |
                                                  gap
                    */
                    out << Interval(gapStart,r1->end());
                    ++r1;
                    goto end;
                }
                
                if (r2->begin() < r1->end())
                {
                    /*
                        W1:     [-----------r1---------------)
                        W2:   [---------------)     [--r2
                                              |     |
                                                gap
                    */
                    out << Interval(gapStart,r2->begin());
                }
                else
                {
                    /*
                        W1:     [-----------r1---------------)
                        W2:   [---------------)                 [--r2
                                              |              |
                                                     gap
                    */
                    out << Interval(gapStart,r1->end());
                    break;
                }
            }

            // No more of r1 to add to W
            if (!++r1) return;
        }
    }
    
end:
    cxAssert(!r1 || !r2);
    while(r1) out << *r1++;
}

// Calculates r1\r2.  Similar to SubtractIntervalLists() above, except the way the intervals are 
// output to 'out' is different.
// Rather than provide them using out << i where i is an interval,
// for each interval [t1,t2) in the output, this function calls out.Add(e2, t1, t2) where e2 is the summed 
// extent of all intervals from r2 to the left of position t1.
template <typename Range1, typename Range2, typename Out>
void SubtractIntervalLists_2(Range1 r1, Range2 r2, Out& out)
{
    using T = decltype(r1->begin());
    T extent2 = 0;
    
    if (!r1 || !r2) goto end;
    for(;;)
    {
        if (r1->end() <= r2->begin())
        {
            /*
                r1 comes before r2 (possibly adjacent)

                W1:      [-----)
                W2:                 [-----)
            */
            out.Add(extent2, r1->begin(), r1->end());
            if (!++r1) return;
        }
        else if (r2->end() <= r1->begin())
        {
            /*
                r2 comes before r1 (possibly adjacent)

                W1:                 [-----)
                W2:      [-----)
            */
            extent2 += r2->size();
            if (!++r2) break;
        }
        else
        {
            /*
            r1 and r2 overlap - i.e. they have a non-empty intersection.
            We need to work out what parts of r1 to keep (if any).

            W1:     [-----------)
            W2:   [---) [--)  [----)
            W:        [-)  [--)
            */
            
            if (r2->begin() > r1->begin())
            {
                /*
                    There is a prefix of W1 to be added to W
                    
                    W1:     [-----------------------)
                    W2:         [---)...
                    W:      [---)    ...
                */
                out.Add(extent2,r1->begin(),r2->begin());
            }
            
            while(r2->end() < r1->end())
            {
                // There are necessarily one or more gaps to be added to W
                T gapStart = r2->end();
                extent2 += r2->size();
                if (!++r2)
                {
                    /*
                        W1:     [-----------r1-----------)
                        W2:   [------r2-------)     
                                              |          |
                                                  gap
                    */
                    out.Add(extent2,gapStart,r1->end());
                    ++r1;
                    goto end;
                }
                
                if (r2->begin() < r1->end())
                {
                    /*
                        W1:     [-----------r1---------------)
                        W2:   [---------------)     [--r2
                                              |     |
                                                gap
                    */
                    out.Add(extent2,gapStart,r2->begin());
                }
                else
                {
                    /*
                        W1:     [-----------r1---------------)
                        W2:   [---------------)                 [--r2
                                              |              |
                                                     gap
                    */
                    out.Add(extent2,gapStart,r1->end());
                    break;
                }
            }

            // No more of r1 to add to W
            if (!++r1) return;
        }
    }
    
end:
    cxAssert(!r1 || !r2);
    while(r1) 
    {
        out.Add(extent2, r1->begin(), r1->end());
        ++r1;
    }
}

// r1, r2 represent interval sets, expressed as iterator-ranges for ordered lists of half
// open intervals.  Writes the result of erasing the intervals of r2 from r1 to 'out' in 
// the manner consistent with an IntervalListWriter<T>.
template <typename Range1, typename Range2, typename Out>
void EraseIntervalLists(Range1 r1, Range2 r2, Out& out)
{
    typedef typename Out::Interval Interval;
    typedef typename Interval::value_type T;

    /*
    W = W1 - W2
    
    Very curiously the looping strategy resembles the code to find U = W1 union W2.
    
                W1:    [-----r1-----)  [--) [--)  [---------------------)           [---)
                
                W2:      [-r2-)  [------------------) [--) [--)   [---)      [---)
                
                U:     [------------------------------------------------)
                                             h                          |    |
                                                                  gap causes h to finish
            
    Our idea is to apply different shifts to h.begin() and h.end(),  which shrinks the size
    of h according to what is erased by W2.
    */

    // Accumulated shift from W2
    T shift = 0;

    if (!r1 || !r2) goto end;
    for(;;)
    {
        if (r1->end() < r2->begin())
        {
            /*
                r1 comes before r2 with a nonzero gap

                W1:      [-----)
                W2:                 [-----)
                               |    |
                            nonzero gap
            */
            out << *r1 + shift;
            if (!++r1) break;
        }
        else if (r2->end() < r1->begin())
        {
            /*
                r2 comes before r1 with a nonzero gap

                W1:                 [-----)
                W2:      [-----)
                               |    |
                            nonzero gap
            */
            shift -= r2->size();
            if (!++r2) break;
        }
        else
        {
            /*
            There is no nonzero gap between r1 and r2.  Therefore they are either adjacent or else 
            they overlap.  Either way they will coalesce in the union.
            Scan forwards through intervals finding a single contiguous interval h to be added to 
            the union.  We stop when we reach the end of W1 or W2 or else we have a gap.
            E.g
            
            W1:    [-----r1-----)  [--) [--)  [---------------------)           [---)
            
            W2:      [-r2-)  [------------------) [--) [--)   [---)      [---)
            
            U:     [------------------------------------------------)
                                         h                          |    |
                                                              gap causes h to finish
                                                              
            Most generally there is a need to alternate between two loops in a leap frog effect.  
            In one loop we step through r2's until they go past r1->end(), and in the other loop we 
            step through r1's until they go past r2->end().
            */
            
            typename Out::Interval h;
            h.setbegin(std::min(r1->begin(), r2->begin()) + shift);
            
            if (r2->end() < r1->end()) goto stepw2;
            
            for(;;)
            {
                /*
                Step through r1 coming before r2->end()
                
                    W1:   [----)  [---)  [---)       [------) 
                    W2:     [-------------------------------)
                                                          r2->end()
                */
                do
                {
                    cxAssert(r1->end() <= r2->end());
                    if (!++r1)
                    {
                        // Reached end of W1, so r2->end() marks the end of h
                        shift -= r2->size();
                        h.setend(r2->end() + shift);
                        if (!h.empty()) out << h;
                        ++r2;
                        goto end;   // Will add remainder of W2 if any
                    }
                } while(r1->end() <= r2->end());
                if (r1->begin() > r2->end())
                {
                    /* r1 comes after r2 with a nonzero gap.  Therefore r2->end() marks the end of h

                        W1:                                        [------) 
                        W2:     [-------------------------------)
                                                              r2->end()
                    */
                    shift -= r2->size();
                    h.setend(r2->end() + shift);
                    if (!h.empty()) out << h;
                    if (!++r2) goto end;
                    break;  // Break out of inner loop, continue with outer loop
                }
            stepw2:
                /*
                Step through r2 coming before r1->end()
                
                    W2:   [----)  [---)  [---)   [------) 
                    W1:     [-------------------------------)
                                                          r1->end()
                */
                do 
                {
                    cxAssert(r2->end() <= r1->end());
                    shift -= r2->size();
                    if (!++r2)
                    {
                        h.setend(r1->end() + shift);
                        out << h;
                        ++r1;
                        goto end;   // Will add remainder of W1 if any
                    }
                } while(r2->end() <= r1->end());
                if (r2->begin() > r1->end())
                {
                    /* r2 comes after r1 with a nonzero gap.  Therefore r1->end() marks the end of h

                        W2:                                        [------) 
                        W1:     [-------------------------------)
                                                              r1->end()
                    */
                    h.setend(r1->end() + shift);
                    out << h;
                    if (!++r1) goto end;
                    break;  // Break out of inner loop, continue with outer loop
                }
            }
        }
    }
end:    
    cxAssert(!r1 || !r2);
    while(r1) out << *r1++ + shift;
}

// out1 = r1-r2 and out2 = r2-r1
template <typename Range1, typename Range2, typename Out1, typename Out2>
void SymmetricEraseIntervalLists(Range1 r1, Range2 r2, Out1& out1, Out2& out2)
{
    typedef typename Out1::Interval Interval1;
    typedef typename Out2::Interval Interval2;
    typedef typename Interval1::value_type T;

    /*
    W = W1 - W2
    
    Very curiously the looping strategy resembles the code to find U = W1 union W2.
    
                W1:    [-----r1-----)  [--) [--)  [---------------------)           [---)
                
                W2:      [-r2-)  [------------------) [--) [--)   [---)      [---)
                
                U:     [------------------------------------------------)
                                             h                          |    |
                                                                  gap causes h to finish
            
    Our idea is to apply different shifts to h.begin() and h.end(),  which shrinks the size
    of h according to what is erased by W2.
    */

    T shift1 = 0;   // Shift to be applied to W1
    T shift2 = 0;   // Shift to be applied to W2

    if (r1 && r2)
    {
        for(;;)
        {
            if (r1->end() < r2->begin())
            {
                /*
                    r1 comes before r2 with a nonzero gap

                    W1:      [-----)
                    W2:                 [-----)
                                   |    |
                                nonzero gap
                */
                out1 << *r1 + shift1;
                shift2 -= r1->size();
                if (!++r1) break;
            }
            else if (r2->end() < r1->begin())
            {
                /*
                    r2 comes before r1 with a nonzero gap

                    W1:                 [-----)
                    W2:      [-----)
                                   |    |
                                nonzero gap
                */
                out2 << *r2 + shift2;
                shift1 -= r2->size();
                if (!++r2) break;
            }
            else
            {
                /*
                There is no nonzero gap between r1 and r2.  Therefore they are either adjacent or else 
                they overlap.  Either way they will coalesce in the union.
                Scan forwards through intervals finding a single contiguous interval h to be added to 
                the union.  We stop when we reach the end of W1 or W2 or else we have a gap.
                E.g
                
                W1:    [-----r1-----)  [--) [--)  [---------------------)           [---)
                
                W2:      [-r2-)  [------------------) [--) [--)   [---)      [---)
                
                U:     [------------------------------------------------)
                                             h                          |    |
                                                                  gap causes h to finish
                                                                  
                Most generally there is a need to alternate between two loops in a leap frog effect.  
                In one loop we step through r2's until they go past r1->end(), and in the other loop we 
                step through r1's until they go past r2->end().
                */
                
                const T hstart = std::min(r1->begin(), r2->begin());
                
                Interval1 h1;
                h1.setbegin(hstart + shift1);

                Interval2 h2;
                h2.setbegin(hstart + shift2);
                
                bool finished = false;
                T hend;

                /*
                Want efficient loop that 
                1.  sets finished = true if need to exit
                2.  sets hend to end position
                3.  Is ready to execute
                
                        h2.setend(hend + shift2);
                        if (!h2.empty()) out2 << h2;
                        h1.setend(hend + shift1);
                        if (!h1.empty()) out1 << h1;
                */
                
                if (r2->end() < r1->end()) goto stepw2;
                
                for(;;)
                {
                    /*
                    Step through r1 coming before r2->end()
                    
                        W1:   [----)  [---)  [---)       [------) 
                        W2:     [-------------------------------)
                                                              r2->end()
                    */
                    do
                    {
                        cxAssert(r1->end() <= r2->end());
                        shift2 -= r1->size();
                        if (!++r1)
                        {
                            finished = true;
                            goto r2_marks_end;
                        }
                    } while(r1->end() <= r2->end());
                    
                    if (r1->begin() > r2->end())
                    {
                        /* r1 comes after r2 with a nonzero gap.  Therefore r2->end() marks the end of h1, h2

                            W1:                                        [------) 
                            W2:     [-------------------------------)
                                                                  r2->end()
                        */
                    r2_marks_end:    
                        shift1 -= r2->size();
                        hend = r2->end();
                        if (!++r2) finished = true;
                        break;
                    }
                stepw2:
                    /*
                    Step through r2 coming before r1->end()
                    
                        W2:   [----)  [---)  [---)   [------) 
                        W1:     [-------------------------------)
                                                              r1->end()
                    */
                    do 
                    {
                        cxAssert(r2->end() <= r1->end());
                        shift1 -= r2->size();
                        if (!++r2)
                        {
                            finished = true;
                            goto r1_marks_end;
                        }
                    } while(r2->end() <= r1->end());
                    
                    if (r2->begin() > r1->end())
                    {
                        /* r2 comes after r1 with a nonzero gap.  Therefore r1->end() marks the end of h1

                            W2:                                        [------) 
                            W1:     [-------------------------------)
                                                                  r1->end()
                        */
                    r1_marks_end:
                        shift2 -= r1->size();
                        hend = r1->end();
                        if (!++r1) finished = true;
                        break;
                    }
                }
                
                h2.setend(hend + shift2);
                if (!h2.empty()) out2 << h2;
                
                h1.setend(hend + shift1);
                if (!h1.empty()) out1 << h1;
                
                if (finished) break;
            }
        }
    }
    cxAssert(!r1 || !r2);
    while(r1) out1 << *r1++ + shift1;
    while(r2) out2 << *r2++ + shift2;
}


// W = W1 intersect W2
template <typename T>
void Intersect(IntervalList<T>& W, const IntervalList<T>& W1, const IntervalList<T>& W2)
{
    // In-place not supported
    cxAssert(&W != &W1);
    cxAssert(&W != &W2);

    W.clear();
    IntersectIntervalLists(
        MakeRangeFromIterators(W1),
        MakeRangeFromIterators(W2),
        CEDA_REFERENCE_TEMPORARY(IntervalListWriter<T>,IntervalListWriter<T>(W)));
}

// W = W1 union W2
template <typename T>
void Union(IntervalList<T>& W, const IntervalList<T>& W1, const IntervalList<T>& W2)
{
    // In-place not supported
    cxAssert(&W != &W1);
    cxAssert(&W != &W2);
    
    W.clear();
    UnionIntervalLists(
        MakeRangeFromIterators(W1),
        MakeRangeFromIterators(W2),
        CEDA_REFERENCE_TEMPORARY(IntervalListWriter<T>, IntervalListWriter<T>(W)));
}

// W = W1 \ W2
template <typename T>
void Subtract(IntervalList<T>& W, const IntervalList<T>& W1, const IntervalList<T>& W2)
{
    // In-place not supported
    cxAssert(&W != &W1);
    cxAssert(&W != &W2);
    
    W.clear();
    SubtractIntervalLists(
        MakeRangeFromIterators(W1),
        MakeRangeFromIterators(W2),
        CEDA_REFERENCE_TEMPORARY(IntervalListWriter<T>,IntervalListWriter<T>(W)));
}

/*
Erase W2 from W1.  An erase W1-W2 is different from a set difference W1\W2 because it 
also causes intervals of W1 to be shifted to the left by the number of elements erased by 
W2.  The shift increases going left to right.
*/
template <typename T>
void Erase(IntervalList<T>& W, const IntervalList<T>& W1, const IntervalList<T>& W2)
{
    // In-place not supported
    cxAssert(&W != &W1);
    cxAssert(&W != &W2);
    
    W.clear();
    EraseIntervalLists(
        MakeRangeFromIterators(W1),
        MakeRangeFromIterators(W2),
        CEDA_REFERENCE_TEMPORARY(IntervalListWriter<T>,IntervalListWriter<T>(W)));
}

template <typename T>
void SymmetricErase(IntervalList<T>& E1, IntervalList<T>& E2, const IntervalList<T>& W1, const IntervalList<T>& W2)
{
    // In-place not supported
    cxAssert(&E1 != &W1);
    cxAssert(&E1 != &W2);
    cxAssert(&E2 != &W1);
    cxAssert(&E2 != &W2);
    
    E1.clear();
    E2.clear();
    SymmetricEraseIntervalLists(
        MakeRangeFromIterators(W1),
        MakeRangeFromIterators(W2),
        CEDA_REFERENCE_TEMPORARY(IntervalListWriter<T>,IntervalListWriter<T>(E1)),
        CEDA_REFERENCE_TEMPORARY(IntervalListWriter<T>,IntervalListWriter<T>(E2)));
}


/*
The functions below can be efficient assuming Named RVO (RVO = Return Value Optimisation)
Google RVO if you've never heard of it before

Note the following:

*   NRVO first appeared in VC2005

*   The functions below comply with the requirements for NVRO.  Therefore the optimsation
    will be made (in release builds)
    
*   The optimisation is not made in /Od (debug builds).

*   No copy or assignment will occur in the following case.  i.e. it is optimal

        IntervalList<T> W = Intersect(W1,W2);

*   There is no risk of aliasing issues.  E.g.

        W1 = Intersect(W1,W2);
    
    will cause creation of a temporary, and an assignment to W1 from the temporary.  This is 
    exactly what we want for correctness.
    
*   Unfortunately the following example is sub-optimal even with NVRO

        IntervalList<T> W3;
        ...
        W3 = Intersect(W1,W2);
    
    i.e. an assignment will be made.  In that case it is better to use
    
        IntervalList<T> W3;
        ...
        Intersect(W3,W1,W2);    // avoids assignment
*/

// Returns W1 intersect W2.  
template <typename T>
IntervalList<T> Intersect(const IntervalList<T>& W1, const IntervalList<T>& W2)
{
    IntervalList<T> W;
    Intersect(W,W1,W2);
    return W;
}

// Returns W1 union W2.  
template <typename T>
IntervalList<T> Union(const IntervalList<T>& W1, const IntervalList<T>& W2)
{
    IntervalList<T> W;
    Union(W,W1,W2);
    return W;
}

// Returns W1 \ W2.  
template <typename T>
IntervalList<T> Subtract(const IntervalList<T>& W1, const IntervalList<T>& W2)
{
    IntervalList<T> W;
    Subtract(W,W1,W2);
    return W;
}

// Returns W1 - W2.  
template <typename T>
IntervalList<T> Erase(const IntervalList<T>& W1, const IntervalList<T>& W2)
{
    IntervalList<T> W;
    Erase(W,W1,W2);
    return W;
}

} // namespace ceda

#endif // include guard