#include <algorithm>
#include <iostream>
#include <iterator>
#include <list>
#include <boost/numeric/interval.hpp>
#include <boost/numeric/interval/io.hpp>

template< class InputIterator1, class InputIterator2, class OutputIterator >
OutputIterator interval_union
   (InputIterator1 first1, InputIterator1 last1,
    InputIterator2 first2, InputIterator2 last2,
    OutputIterator result)
{
  // If a container is empty, return the other one.
  if (first1 == last1) return std::copy(first2, last2, result);
  if (first2 == last2) return std::copy(first1, last1, result);

  // Keep a copy of the current elements, possibly augmented.
  typedef typename std::iterator_traits<InputIterator1>::value_type I;
  I i1 = *first1, i2 = *first2;

  // Iterate as long as both containers have elements.
  for (;;)
  {
    // If intervals do not overlap, output the leftmost one.
    // Otherwise, augment the rightmost one with the leftmost one.
    // In both cases, discard the leftmost one.
    if (upper(i1) < upper(i2))
    {
      if (overlap(i1, i2)) i2 = hull(i1, i2);
      else *(result++) = i1;
      if (++first1 == last1) break;
      i1 = *first1;
    }
    else
    {
      if (overlap(i1, i2)) i1 = hull(i1, i2);
      else *(result++) = i2;
      if (++first2 == last2) break;
      i2 = *first2;
    }
  }

  // Copy the remaining elements.
  if (first1 == last1)
  {
    *(result++) = i2;
    return std::copy(++first2, last2, result);
  }
  else
  {
    *(result++) = i1;
    return std::copy(++first1, last1, result);
  }
}

template< class InputIterator1, class InputIterator2, class OutputIterator >
OutputIterator interval_intersection
   (InputIterator1 first1, InputIterator1 last1,
    InputIterator2 first2, InputIterator2 last2,
    OutputIterator result)
{
  typedef typename std::iterator_traits<InputIterator1>::value_type I;

  // Iterate as long as both containers have elements.
  while (first1 != last1 && first2 != last2)
  {
    I const &i1 = *first1, &i2 = *first2;
    // Output the intersection of the two leftmost elements.
    if (overlap(i1, i2)) *(result++) = intersect(i1, i2);
    // Discard the leftmost one.
    if (upper(i1) < upper(i2)) ++first1;
    else ++first2;
  }

  return result;
}

int main()
{
  typedef boost::numeric::interval< double > I;
  std::list< I > l1, l2, l3, l4;
  l1.push_back(I(0,2)); l1.push_back(I(3,5));
  l2.push_back(I(1,4)); l2.push_back(I(6,7));

  interval_union(l1.begin(), l1.end(), l2.begin(), l2.end(),
                 std::back_inserter(l3));

  interval_intersection(l1.begin(), l1.end(), l2.begin(), l2.end(),
                        std::back_inserter(l4));

  std::cout << "Union: ";
  std::copy(l3.begin(), l3.end(),
            std::ostream_iterator<I>(std::cout, " "));
  std::cout << std::endl;

  std::cout << "Intersection: ";
  std::copy(l4.begin(), l4.end(),
            std::ostream_iterator<I>(std::cout, " "));
  std::cout << std::endl;

  return 0;
}