/* This is the header file of the Parma Polyhedra Library.
   Copyright (C) 2001-2004 Roberto Bagnara <bagnara@cs.unipr.it>

This file is part of the Parma Polyhedra Library (PPL).

The PPL is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.

The PPL is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
USA.

For the most up-to-date information see the Parma Polyhedra Library
site: http://www.cs.unipr.it/ppl/ . */

#ifndef PPL_ppl_hh
#define PPL_ppl_hh 1

#ifdef NDEBUG
# define PPL_SAVE_NDEBUG 1
# undef NDEBUG
#endif

/* Declaration of macros and functions providing version  -*- C++ -*-
   and licensing information.
   Copyright (C) 2001-2004 Roberto Bagnara <bagnara@cs.unipr.it>

This file is part of the Parma Polyhedra Library (PPL).

The PPL is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.

The PPL is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
USA.

For the most up-to-date information see the Parma Polyhedra Library
site: http://www.cs.unipr.it/ppl/ . */

#ifndef PPL_version_hh
#define PPL_version_hh 1

/*! \defgroup PPL_defines Library Defines */
//@{

//! The major number of the PPL version.
#define PPL_VERSION_MAJOR 0

//! The minor number of the PPL version.
#define PPL_VERSION_MINOR 6

//! The revision number of the PPL version.
#define PPL_VERSION_REVISION 1

//! \brief
//! The beta number of the PPL version.  This is zero for official
//! releases and nonzero for development snapshots.
#define PPL_VERSION_BETA 0

//! A string containing the PPL version.
/*!
  Let <CODE>M</CODE> and <CODE>m</CODE> denote the numbers associated
  to PPL_VERSION_MAJOR and PPL_VERSION_MINOR, respectively.  The
  format of PPL_VERSION is <CODE>M "." m</CODE> if both
  PPL_VERSION_REVISION (<CODE>r</CODE>) and PPL_VERSION_BETA
  (<CODE>b</CODE>)are zero, <CODE>M "." m "pre" b</CODE> if
  PPL_VERSION_REVISION is zero and PPL_VERSION_BETA is not zero,
  <CODE>M "." m "." r</CODE> if PPL_VERSION_REVISION is not zero and
  PPL_VERSION_BETA is zero, <CODE>M "." m "." r "pre" b</CODE> if
  neither PPL_VERSION_REVISION nor PPL_VERSION_BETA are zero.
*/
#define PPL_VERSION "0.6.1"

//@}

namespace Parma_Polyhedra_Library {

//! Returns the major number of the PPL version.
unsigned
version_major();

//! Returns the minor number of the PPL version.
unsigned
version_minor();

//! Returns the revision number of the PPL version.
unsigned
version_revision();

//! Returns the beta number of the PPL version.
unsigned
version_beta();

//! Returns a character string containing the PPL version.
const char* version();

//! Returns a character string containing the PPL banner.
/*!
  The banner provides information about the PPL version, the licensing,
  the lack of any warranty whatsoever, the C++ compiler used to build
  the library, where to report bugs and where to look for further
  information.
*/
const char* banner();

} // namespace Parma_Polyhedra_Library

#endif // !defined(PPL_version_hh)
/* config.h.  Generated by configure.  */
/* config.h.in.  Generated from configure.ac by autoheader.  */


#ifndef PPL_ppl_config_h
#define PPL_ppl_config_h 1


/* Not zero if the C++ compiler supports __attribute__ ((weak)). */
#define CXX_SUPPORTS_ATTRIBUTE_WEAK 1

/* Not zero if the C++ compiler supports flexible arrays. */
#define CXX_SUPPORTS_FLEXIBLE_ARRAYS 1

/* Enable more assertions when defined. */
/* #undef EXTRA_ROW_DEBUG */

/* Not zero if GMP has been compiled with support for exceptions. */
#define GMP_SUPPORTS_EXCEPTIONS 1

/* Define to 1 if you have the declaration of `RLIMIT_AS', and to 0 if you
   don't. */
#define HAVE_DECL_RLIMIT_AS 1

/* Define to 1 if you have the declaration of `RLIMIT_DATA', and to 0 if you
   don't. */
#define HAVE_DECL_RLIMIT_DATA 1

/* Define to 1 if you have the declaration of `RLIMIT_RSS', and to 0 if you
   don't. */
#define HAVE_DECL_RLIMIT_RSS 1

/* Define to 1 if you have the declaration of `RLIMIT_VMEM', and to 0 if you
   don't. */
#define HAVE_DECL_RLIMIT_VMEM 1

/* Define to 1 if you have the <dlfcn.h> header file. */
#define HAVE_DLFCN_H 1

/* Define to 1 if you have the <getopt.h> header file. */
#define HAVE_GETOPT_H 1

/* Define to 1 if you have the <inttypes.h> header file. */
#define HAVE_INTTYPES_H 1

/* Define to 1 if you have the <memory.h> header file. */
#define HAVE_MEMORY_H 1

/* Define to 1 if you have the <signal.h> header file. */
#define HAVE_SIGNAL_H 1

/* Define to 1 if you have the <stdint.h> header file. */
#define HAVE_STDINT_H 1

/* Define to 1 if you have the <stdlib.h> header file. */
#define HAVE_STDLIB_H 1

/* Define to 1 if you have the <strings.h> header file. */
#define HAVE_STRINGS_H 1

/* Define to 1 if you have the <string.h> header file. */
#define HAVE_STRING_H 1

/* Define to 1 if you have the <sys/resource.h> header file. */
#define HAVE_SYS_RESOURCE_H 1

/* Define to 1 if you have the <sys/stat.h> header file. */
#define HAVE_SYS_STAT_H 1

/* Define to 1 if you have the <sys/time.h> header file. */
#define HAVE_SYS_TIME_H 1

/* Define to 1 if you have the <sys/types.h> header file. */
#define HAVE_SYS_TYPES_H 1

/* Define to 1 if you have the <unistd.h> header file. */
#define HAVE_UNISTD_H 1

/* Assertions are disabled when this is defined. */
#define NDEBUG 1

/* Name of package */
#define PACKAGE "ppl"

/* Define to the address where bug reports for this package should be sent. */
#define PACKAGE_BUGREPORT "ppl-devel@cs.unipr.it"

/* Define to the full name of this package. */
#define PACKAGE_NAME "the Parma Polyhedra Library"

/* Define to the full name and version of this package. */
#define PACKAGE_STRING "the Parma Polyhedra Library 0.6.1"

/* Define to the one symbol short name of this package. */
#define PACKAGE_TARNAME "ppl"

/* Define to the version of this package. */
#define PACKAGE_VERSION "0.6.1"

/* Size of GMP's mp_limb_t. */
#define SIZEOF_MP_LIMB_T 8

/* Define to 1 if you have the ANSI C header files. */
#define STDC_HEADERS 1

/* Version number of package */
#define VERSION "0.6.1"

/* Define to empty if `const' does not conform to ANSI C. */
/* #undef const */

/* Define to `__inline__' or `__inline' if that's what the C compiler
   calls it, or to nothing if 'inline' is not supported under any name.  */
#ifndef __cplusplus
/* #undef inline */
#endif


#endif // !defined(PPL_ppl_config_h)

// This file intentionally contains only this comment line.
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <deque>
#include <exception>
#include <gmp.h>
#include <gmpxx.h>
#include <iosfwd>
#include <iostream>
#include <iterator>
#include <limits>
#include <list>
#include <map>
#include <set>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
namespace Parma_Polyhedra_Library {


typedef mpz_class Integer;

}
namespace Parma_Polyhedra_Library {

class ERational;
class Boundary;
class LBoundary;
class UBoundary;
class Interval;

}
namespace Parma_Polyhedra_Library {

class BoundingBox;

}
namespace Parma_Polyhedra_Library {

class Constraint;

}
namespace Parma_Polyhedra_Library {

class Generator;

}
namespace Parma_Polyhedra_Library {

class Init;

}
namespace Parma_Polyhedra_Library {

class SatRow;

}
namespace Parma_Polyhedra_Library {

class SatMatrix;

}
namespace Parma_Polyhedra_Library {

class Row;

}
namespace Parma_Polyhedra_Library {

class Matrix;

}
namespace Parma_Polyhedra_Library {

class Variable;

}
namespace Parma_Polyhedra_Library {

class LinExpression;

}
namespace Parma_Polyhedra_Library {

class ConSys;

}
namespace Parma_Polyhedra_Library {

class GenSys;

}
namespace Parma_Polyhedra_Library {

class Poly_Con_Relation;

}
namespace Parma_Polyhedra_Library {

class Poly_Gen_Relation;

}
namespace Parma_Polyhedra_Library {

class BHRZ03_Certificate;

}
namespace Parma_Polyhedra_Library {

class H79_Certificate;

}
namespace Parma_Polyhedra_Library {

class Polyhedron;
namespace IO_Operators {
}

}
namespace std {
}
namespace Parma_Polyhedra_Library {

class C_Polyhedron;

}
namespace Parma_Polyhedra_Library {

class NNC_Polyhedron;

}
namespace Parma_Polyhedra_Library {


typedef size_t dimension_type;


inline dimension_type
not_a_dimension() {
  return std::numeric_limits<dimension_type>::max();
}


enum Relation_Symbol {

  LESS_THAN,

  LESS_THAN_OR_EQUAL,

  EQUAL,

  GREATER_THAN_OR_EQUAL,

  GREATER_THAN
};


enum Complexity_Class {

  POLYNOMIAL,

  SIMPLEX,

  ANY
};




extern Integer* tmp_Integer;
inline dimension_type
compute_capacity(dimension_type requested_size) {

  return 2*(requested_size + 1);


}







class Throwable {
public:

  virtual void throw_me() const = 0;
};
extern const Throwable* volatile abandon_expensive_computations;







inline void
maybe_abandon() {
  if (const Throwable* p = abandon_expensive_computations)
    p->throw_me();
}





struct From_Bounding_Box {
};





const Integer& operator*(const Constraint& c, const Generator& g);






const Integer& reduced_scalar_product(const Constraint& c,
          const Generator& g);

}
namespace Parma_Polyhedra_Library {

inline void negate(Integer& x);
inline void gcd_assign(Integer& x, const Integer& y);
inline void gcd_assign(Integer& x, const Integer& y, const Integer& z);
inline void lcm_assign(Integer& x, const Integer& y);
inline void lcm_assign(Integer& x, const Integer& y, const Integer& z);
inline void exact_div_assign(Integer& x, const Integer& y);
inline void exact_div_assign(Integer& x, const Integer& y, const Integer& z);
inline void sqrt_assign(Integer& x);
inline void sqrt_assign(Integer& x, const Integer& y);
inline int cmp(const Integer& x, const Integer& y);
inline const Integer& Integer_zero();
inline const Integer& Integer_one();

}


namespace std {


void swap(Parma_Polyhedra_Library::Integer& x,
   Parma_Polyhedra_Library::Integer& y);

}
namespace Parma_Polyhedra_Library {

inline void
negate(Integer& x) {
  mpz_neg(x.get_mpz_t(), x.get_mpz_t());
}

inline void
gcd_assign(Integer& x, const Integer& y) {
  mpz_gcd(x.get_mpz_t(), x.get_mpz_t(), y.get_mpz_t());
}

inline void
gcd_assign(Integer& x, const Integer& y, const Integer& z) {
  mpz_gcd(x.get_mpz_t(), y.get_mpz_t(), z.get_mpz_t());
}

inline void
lcm_assign(Integer& x, const Integer& y) {
  mpz_lcm(x.get_mpz_t(), x.get_mpz_t(), y.get_mpz_t());
}

inline void
lcm_assign(Integer& x, const Integer& y, const Integer& z) {
  mpz_lcm(x.get_mpz_t(), y.get_mpz_t(), z.get_mpz_t());
}

inline void
exact_div_assign(Integer& x, const Integer& y) {
  mpz_divexact(x.get_mpz_t(), x.get_mpz_t(), y.get_mpz_t());
}

inline void
exact_div_assign(Integer& x, const Integer& y, const Integer& z) {
  mpz_divexact(x.get_mpz_t(), y.get_mpz_t(), z.get_mpz_t());
}

inline void
sqrt_assign(Integer& x) {
  mpz_sqrt(x.get_mpz_t(), x.get_mpz_t());
}

inline void
sqrt_assign(Integer& x, const Integer& y) {
  mpz_sqrt(x.get_mpz_t(), y.get_mpz_t());
}

inline int
cmp(const Integer& x, const Integer& y) {
  return mpz_cmp(x.get_mpz_t(), y.get_mpz_t());
}

inline const Integer&
Integer_zero() {
  static Integer z(0);
  return z;
}

inline const Integer&
Integer_one() {
  static Integer o(1);
  return o;
}

}

inline void
std::swap(Parma_Polyhedra_Library::Integer& x,
   Parma_Polyhedra_Library::Integer& y) {
  mpz_swap(x.get_mpz_t(), y.get_mpz_t());
}
namespace Parma_Polyhedra_Library {





bool operator==(const ERational& x, const ERational& y);





bool operator!=(const ERational& x, const ERational& y);







bool operator<=(const ERational& x, const ERational& y);







bool operator>=(const ERational& x, const ERational& y);





bool operator<(const ERational& x, const ERational& y);





bool operator>(const ERational& x, const ERational& y);

namespace IO_Operators {





std::ostream& operator<<(std::ostream& s, const ERational& x);

}

}
class Parma_Polyhedra_Library::ERational {
public:




  ERational(const Integer& num, const Integer& den);







  explicit ERational(char sign);


  ERational(const ERational& y);


  ERational& operator=(const ERational& y);





  int direction_of_infinity() const;





  const Integer& numerator() const;





  const Integer& denominator() const;

  friend bool
  Parma_Polyhedra_Library::operator==(const ERational& x, const ERational& y);
  friend bool
  Parma_Polyhedra_Library::operator!=(const ERational& x, const ERational& y);
  friend bool
  Parma_Polyhedra_Library::operator<=(const ERational& x, const ERational& y);
  friend bool
  Parma_Polyhedra_Library::operator>=(const ERational& x, const ERational& y);
  friend bool
  Parma_Polyhedra_Library::operator<(const ERational& x, const ERational& y);
  friend bool
  Parma_Polyhedra_Library::operator>(const ERational& x, const ERational& y);
  friend std::ostream&
  Parma_Polyhedra_Library::IO_Operators::operator<<(std::ostream& s,
          const ERational& x);

private:

  int e;


  mpq_class v;
};


namespace Parma_Polyhedra_Library {





bool operator<(const Boundary& x, const Boundary& y);





bool operator>(const Boundary& x, const Boundary& y);

}
class Parma_Polyhedra_Library::Boundary {
protected:

  enum Flag {

    NEG = -1,

    ZERO = 0,

    POS = 1
  };


  ERational value;

  Flag flag;


  Boundary(const ERational& v, Flag f);

  friend bool
  Parma_Polyhedra_Library::operator<(const Boundary& x, const Boundary& y);
  friend bool
  Parma_Polyhedra_Library::operator>(const Boundary& x, const Boundary& y);

public:

  bool is_closed() const;


  const ERational& bound() const;
};





class Parma_Polyhedra_Library::LBoundary : public Boundary {
public:

  enum OpenClosed {

    OPEN = Boundary::POS,

    CLOSED = Boundary::ZERO
  };


  LBoundary(const ERational& v, OpenClosed f);


  bool OK() const;
};




class Parma_Polyhedra_Library::UBoundary : public Boundary {
public:

  enum OpenClosed {

    OPEN = Boundary::NEG,

    CLOSED = Boundary::ZERO
  };


  UBoundary(const ERational& v, OpenClosed f);


  bool OK() const;
};
class Parma_Polyhedra_Library::Interval {
public:

  Interval();


  bool is_empty() const;


  const LBoundary& lower_bound() const;


  const UBoundary& upper_bound() const;




  void raise_lower_bound(LBoundary new_lower);




  void lower_upper_bound(UBoundary new_upper);


  void set_empty();


  bool OK() const;

private:

  LBoundary lower;


  UBoundary upper;
};
namespace Parma_Polyhedra_Library {

inline
ERational::ERational(const Integer& num, const Integer& den)
  : e(0), v(num, den) {
  assert(den != 0);
  v.canonicalize();
}

inline
ERational::ERational(char sign)
  : e(sign == '+' ? 1 : -1) {
  assert(sign == '+' || sign == '-');
}

inline
ERational::ERational(const ERational& y)
  : e(y.e) {
  if (e == 0)
    v = y.v;
}

inline ERational&
ERational::operator=(const ERational& y) {
  e = y.e;
  if (e == 0)
    v = y.v;
  return *this;
}

inline int
ERational::direction_of_infinity() const {
  return e;
}

inline const Integer&
ERational::numerator() const {
  assert(e == 0);
  return v.get_num();
}

inline const Integer&
ERational::denominator() const {
  assert(e == 0);
  return v.get_den();
}


inline bool
operator==(const ERational& x, const ERational& y) {
  return x.e == y.e && (x.e != 0 || x.v == y.v);
}


inline bool
operator!=(const ERational& x, const ERational& y) {
  return !(x == y);
}


inline bool
operator<(const ERational& x, const ERational& y) {
  return x.e < y.e || (x.e == 0 && y.e == 0 && x.v < y.v);
}


inline bool
operator>(const ERational& x, const ERational& y) {
  return y < x;
}


inline bool
operator<=(const ERational& x, const ERational& y) {
  return x < y || x == y;
}


inline bool
operator>=(const ERational& x, const ERational& y) {
  return y <= x;
}

inline
Boundary::Boundary(const ERational& v, Flag f)
  : value(v), flag(f) {
}

inline bool
Boundary::is_closed() const {
  return flag == ZERO;
}

inline const ERational&
Boundary::bound() const {
  return value;
}

inline
LBoundary::LBoundary(const ERational& v, OpenClosed f)
  : Boundary(v, f == CLOSED ? ZERO : POS) {
}

inline
UBoundary::UBoundary(const ERational& v, OpenClosed f)
  : Boundary(v, f == CLOSED ? ZERO : NEG) {
}


inline bool
operator<(const Boundary& x, const Boundary& y) {
  return x.value < y.value ||
    (x.value == y.value && x.flag < y.flag);
}


inline bool
operator>(const Boundary& x, const Boundary& y) {
  return y < x;
}

inline
Interval::Interval()
  : lower(ERational('-'), LBoundary::OPEN),
    upper(ERational('+'), UBoundary::OPEN) {
}

inline bool
Interval::is_empty() const {
  return lower > upper;
}

inline const LBoundary&
Interval::lower_bound() const {
  return lower;
}

inline const UBoundary&
Interval::upper_bound() const {
  return upper;
}

inline void
Interval::raise_lower_bound(LBoundary new_lower) {
  if (new_lower > lower)
    lower = new_lower;
}

inline void
Interval::lower_upper_bound(UBoundary new_upper) {
  if (new_upper < upper)
    upper = new_upper;
}

inline void
Interval::set_empty() {
  lower = LBoundary(ERational('+'), LBoundary::OPEN);
  upper = UBoundary(ERational('-'), UBoundary::OPEN);
  assert(is_empty());
}

}
namespace Parma_Polyhedra_Library {
class BoundingBox {
public:

  BoundingBox(dimension_type num_dimensions);


  dimension_type space_dimension() const;




  const Interval& operator[](dimension_type k) const;


  bool is_empty() const;
  bool get_lower_bound(dimension_type k, bool& closed,
         Integer& n, Integer& d) const;
  bool get_upper_bound(dimension_type k, bool& closed,
         Integer& n, Integer& d) const;


  void set_empty();
  void raise_lower_bound(dimension_type k, bool closed,
    const Integer& n, const Integer& d);
  void lower_upper_bound(dimension_type k, bool closed,
    const Integer& n, const Integer& d);

private:



  std::vector<Interval> vec;



  mutable bool empty;

  mutable bool empty_up_to_date;
};

namespace IO_Operators {





std::ostream& operator<<(std::ostream& s, const BoundingBox& bbox);

}

}
namespace Parma_Polyhedra_Library {

inline
BoundingBox::BoundingBox(dimension_type num_dimensions)
  : vec(num_dimensions), empty(false), empty_up_to_date(true) {
}

inline dimension_type
BoundingBox::space_dimension() const {
  return vec.size();
}

inline const Interval&
BoundingBox::operator[](const dimension_type k) const {
  assert(k < vec.size());
  return vec[k];
}

inline bool
BoundingBox::is_empty() const {
  if (empty_up_to_date)
    return empty;
  else {
    empty_up_to_date = true;
    for (dimension_type k = vec.size(); k-- > 0; )
      if (vec[k].is_empty()) {
 empty = true;
 return true;
      }
    empty = false;
    return false;
  }
}

inline bool
BoundingBox::get_lower_bound(const dimension_type k, bool& closed,
        Integer& n, Integer& d) const {
  assert(k < vec.size());
  const LBoundary& lb = vec[k].lower_bound();
  const ERational& lr = lb.bound();

  if (lr.direction_of_infinity() != 0)
    return false;

  closed = lb.is_closed();
  n = lr.numerator();
  d = lr.denominator();

  return true;
}

inline bool
BoundingBox::get_upper_bound(const dimension_type k, bool& closed,
        Integer& n, Integer& d) const {
  assert(k < vec.size());
  const UBoundary& ub = vec[k].upper_bound();
  const ERational& ur = ub.bound();

  if (ur.direction_of_infinity() != 0)
    return false;

  closed = ub.is_closed();
  n = ur.numerator();
  d = ur.denominator();
  return true;
}

inline void
BoundingBox::set_empty() {
  for (dimension_type k = vec.size(); k-- > 0; )
    vec[k].set_empty();
  empty = empty_up_to_date = true;
}

inline void
BoundingBox::raise_lower_bound(const dimension_type k, const bool closed,
          const Integer& n, const Integer& d) {
  assert(k < vec.size());
  assert(d != 0);
  vec[k].raise_lower_bound(LBoundary(ERational(n, d),
         (closed
          ? LBoundary::CLOSED
          : LBoundary::OPEN)));
  empty_up_to_date = false;
}

inline void
BoundingBox::lower_upper_bound(const dimension_type k, const bool closed,
          const Integer& n, const Integer& d) {
  assert(k < vec.size());
  assert(d != 0);
  vec[k].lower_upper_bound(UBoundary(ERational(n, d),
         (closed
          ? UBoundary::CLOSED
          : UBoundary::OPEN)));
  empty_up_to_date = false;
}

}
namespace Parma_Polyhedra_Library {

namespace IO_Operators {



std::ostream&
operator<<(std::ostream& s, const Variable& v);

}



bool less(Variable v, Variable w);

}
class Parma_Polyhedra_Library::Variable {

public:

  explicit Variable(dimension_type i);


  dimension_type id() const;


  typedef void output_function_type(std::ostream& s, const Variable& v);


  static void set_output_function(output_function_type* p);


  static output_function_type* get_output_function();


  struct Compare {

    bool operator()(Variable x, Variable y) const;
  };

private:

  dimension_type varid;


  friend class Init;

  friend std::ostream&
  Parma_Polyhedra_Library::IO_Operators::operator<<(std::ostream& s,
          const Variable& v);


  static output_function_type* current_output_function;


  static void default_output_function(std::ostream& s, const Variable& v);
};


namespace Parma_Polyhedra_Library {


typedef std::set<Variable, Variable::Compare> Variables_Set;

}
namespace Parma_Polyhedra_Library {

inline
Variable::Variable(dimension_type i)
  : varid(i) {
}

inline dimension_type
Variable::id() const {
  return varid;
}

inline void
Variable::set_output_function(output_function_type* p) {
  current_output_function = p;
}

inline Variable::output_function_type*
Variable::get_output_function() {
  return current_output_function;
}


inline bool
less(const Variable v, const Variable w) {
  return v.id() < w.id();
}

inline bool
Variable::Compare::operator()(const Variable x, const Variable y) const {
  return less(x, y);
}

}
class Parma_Polyhedra_Library::Init {
private:

  static unsigned int count;

public:

  Init();


  ~Init();
};
namespace {

Parma_Polyhedra_Library::Init Parma_Polyhedra_Library_initializer;

}
namespace Parma_Polyhedra_Library {




  enum Topology {
    NECESSARILY_CLOSED = 0,
    NOT_NECESSARILY_CLOSED = 1
  };

}
namespace Parma_Polyhedra_Library {







bool operator==(const SatRow& x, const SatRow& y);





bool operator!=(const SatRow& x, const SatRow& y);
int compare(const SatRow& x, const SatRow& y);





bool subset_or_equal(const SatRow& x, const SatRow& y);





bool strict_subset(const SatRow& x, const SatRow& y);





void set_union(const SatRow& x, const SatRow& y, SatRow& z);

}
class Parma_Polyhedra_Library::SatRow {
private:

  mpz_t vec;


  static unsigned int first_one(mp_limb_t w);


  static unsigned int last_one(mp_limb_t w);

public:

  SatRow();


  SatRow(const SatRow& y);


  ~SatRow();


  SatRow& operator=(const SatRow& y);


  void swap(SatRow& y);


  bool operator[](unsigned int k) const;


  void set(unsigned int k);


  void clear(unsigned int k);


  void clear_from(unsigned int k);


  void clear();

  friend int
  Parma_Polyhedra_Library::compare(const SatRow& x,
       const SatRow& y);

  friend bool
  Parma_Polyhedra_Library::operator==(const SatRow& x, const SatRow& y);

  friend bool
  Parma_Polyhedra_Library::operator!=(const SatRow& x, const SatRow& y);

  friend bool
  Parma_Polyhedra_Library::subset_or_equal(const SatRow& x, const SatRow& y);

  friend bool
  Parma_Polyhedra_Library::strict_subset(const SatRow& x, const SatRow& y);

  friend void
  Parma_Polyhedra_Library::set_union(const SatRow& x,
         const SatRow& y,
         SatRow& z);


  int first() const;




  int next(int position) const;


  int last() const;




  int prev(int position) const;


  unsigned int count_ones() const;


  bool empty() const;


  bool OK() const;
};

namespace std {





void swap(Parma_Polyhedra_Library::SatRow& x,
   Parma_Polyhedra_Library::SatRow& y);





void
iter_swap(std::vector<Parma_Polyhedra_Library::SatRow>::iterator x,
   std::vector<Parma_Polyhedra_Library::SatRow>::iterator y);

}
namespace Parma_Polyhedra_Library {

inline
SatRow::SatRow() {
  mpz_init(vec);
}

inline
SatRow::SatRow(const SatRow& y) {
  mpz_init_set(vec, y.vec);
}

inline
SatRow::~SatRow() {
  mpz_clear(vec);
}

inline SatRow&
SatRow::operator=(const SatRow& y) {
  mpz_set(vec, y.vec);
  return *this;
}

inline bool
SatRow::operator[](const unsigned int k) const {
  return mpz_tstbit(vec, k);
}

inline void
SatRow::set(const unsigned int k) {
  mpz_setbit(vec, k);
}

inline void
SatRow::clear(const unsigned int k) {
  mpz_clrbit(vec, k);
}

inline void
SatRow::clear_from(const unsigned int k) {
  mpz_tdiv_r_2exp(vec, vec, k);
}

inline unsigned int
SatRow::count_ones() const {
  return mpz_popcount(vec);
}

inline bool
SatRow::empty() const {
  return mpz_sgn(vec) == 0;
}

inline void
SatRow::swap(SatRow& y) {
  mpz_swap(vec, y.vec);
}

inline void
SatRow::clear() {
  mpz_set_ui(vec, 0UL);
}


inline bool
operator==(const SatRow& x, const SatRow& y) {
  return mpz_cmp(x.vec, y.vec) == 0;
}


inline bool
operator!=(const SatRow& x, const SatRow& y) {
  return mpz_cmp(x.vec, y.vec) != 0;
}


inline void
set_union(const SatRow& x, const SatRow& y, SatRow& z) {
  mpz_ior(z.vec, x.vec, y.vec);
}

}


namespace std {


inline void
swap(Parma_Polyhedra_Library::SatRow& x,
     Parma_Polyhedra_Library::SatRow& y) {
  x.swap(y);
}


inline void
iter_swap(std::vector<Parma_Polyhedra_Library::SatRow>::iterator x,
   std::vector<Parma_Polyhedra_Library::SatRow>::iterator y) {
  swap(*x, *y);
}

}
class Parma_Polyhedra_Library::SatMatrix {
public:

  SatMatrix();




  SatMatrix(dimension_type n_rows, dimension_type n_columns);


  SatMatrix(const SatMatrix& y);


  ~SatMatrix();


  SatMatrix& operator=(const SatMatrix& y);


  void swap(SatMatrix& y);


  SatRow& operator[](dimension_type k);


  const SatRow& operator[](dimension_type k) const;


  void clear();


  void transpose();


  void transpose_assign(const SatMatrix& y);


  dimension_type num_columns() const;


  dimension_type num_rows() const;


  void sort_rows();
  bool sorted_contains(const SatRow& row) const;


  void add_row(const SatRow& row);


  void rows_erase_to_end(dimension_type first_to_erase);


  void columns_erase_to_end(dimension_type first_to_erase);


  void resize(dimension_type new_n_rows, dimension_type new_n_columns);


  bool OK() const;




  void ascii_dump(std::ostream& s) const;





  bool ascii_load(std::istream& s);






private:

  std::vector<SatRow> rows;


  dimension_type row_size;

  struct RowCompare {
    bool operator()(const SatRow& x, const SatRow& y) const;
  };
};

namespace std {





void swap(Parma_Polyhedra_Library::SatMatrix& x,
   Parma_Polyhedra_Library::SatMatrix& y);

}
namespace Parma_Polyhedra_Library {

inline
SatMatrix::SatMatrix()
  : rows(),
    row_size(0) {
}

inline
SatMatrix::SatMatrix(const dimension_type n_rows,
       const dimension_type n_columns)
  : rows(n_rows),
    row_size(n_columns) {
}

inline
SatMatrix::SatMatrix(const SatMatrix& y)
  : rows(y.rows),
    row_size(y.row_size) {
}

inline
SatMatrix::~SatMatrix() {
}

inline void
SatMatrix::rows_erase_to_end(const dimension_type first_to_erase) {


  assert(first_to_erase <= rows.size());
  if (first_to_erase < rows.size())
    rows.erase(rows.begin() + first_to_erase, rows.end());
  assert(OK());
}

inline void
SatMatrix::columns_erase_to_end(const dimension_type first_to_erase) {


  assert(first_to_erase <= row_size);
  row_size = first_to_erase;
  assert(OK());
}

inline void
SatMatrix::swap(SatMatrix& y) {
  std::swap(row_size, y.row_size);
  std::swap(rows, y.rows);
}

inline SatRow&
SatMatrix::operator[](const dimension_type k) {
  assert(k < rows.size());
  return rows[k];
}

inline const SatRow&
SatMatrix::operator[](const dimension_type k) const {
  assert(k < rows.size());
  return rows[k];
}

inline dimension_type
SatMatrix::num_columns() const {
  return row_size;
}

inline dimension_type
SatMatrix::num_rows() const {
  return rows.size();
}

inline void
SatMatrix::clear() {

  std::vector<SatRow>().swap(rows);
  row_size = 0;
}

inline bool
SatMatrix::RowCompare::
operator()(const SatRow& x, const SatRow& y) const {
  return compare(x, y) < 0;
}

}


namespace std {


inline void
swap(Parma_Polyhedra_Library::SatMatrix& x,
     Parma_Polyhedra_Library::SatMatrix& y) {
  x.swap(y);
}

}
class Parma_Polyhedra_Library::Row {
public:



  class Type;


  enum Kind {
    LINE_OR_EQUALITY = 0,
    RAY_OR_POINT_OR_INEQUALITY = 1
  };


  Row();







  void construct(Type t, dimension_type sz);
  void construct(Type t, dimension_type sz, dimension_type capacity);



  Row(Type t, dimension_type sz);


  Row(Type t, dimension_type sz, dimension_type capacity);


  Row(const Row& y);





  Row(const Row& y, dimension_type capacity);






  Row(const Row& y, dimension_type sz, dimension_type capacity);


  ~Row();


  Row& operator=(const Row& y);


  void swap(Row& y);






  void assign(Row& y);






  void resize_no_copy(dimension_type new_sz);






  void grow_no_copy(dimension_type new_sz);






  void shrink(dimension_type new_sz);




  Type type() const;


  Topology topology() const;



  bool is_necessarily_closed() const;



  bool is_line_or_equality() const;



  bool is_ray_or_point_or_inequality() const;





  void set_necessarily_closed();


  void set_not_necessarily_closed();


  void set_is_line_or_equality();


  void set_is_ray_or_point_or_inequality();



  static dimension_type max_size();


  dimension_type size() const;


  dimension_type space_dimension() const;


  const Integer& inhomogeneous_term() const;


  const Integer& coefficient(dimension_type n) const;




  Integer& operator[](dimension_type k);


  const Integer& operator[](dimension_type k) const;







  void normalize();




  void sign_normalize();







  void strong_normalize();



  bool check_strong_normalized() const;
  void linear_combine(const Row& y, dimension_type k);




  bool all_homogeneous_terms_are_zero() const;


  bool OK(dimension_type row_size, dimension_type row_capacity) const;

private:
  friend class Parma_Polyhedra_Library::LinExpression;
  friend class Parma_Polyhedra_Library::Constraint;
  friend class Parma_Polyhedra_Library::Generator;

  class Impl;


  Impl* impl;
};

namespace Parma_Polyhedra_Library {





const Integer& operator*(const Row& x, const Row& y);







const Integer& reduced_scalar_product(const Row& x, const Row& y);
int compare(const Row& x, const Row& y);






bool operator==(const Row& x, const Row& y);

bool operator!=(const Row& x, const Row& y);

bool operator<=(const Row& x, const Row& y);

bool operator <(const Row& x, const Row& y);

bool operator>=(const Row& x, const Row& y);

bool operator >(const Row& x, const Row& y);




}
class Parma_Polyhedra_Library::Row::Type {
public:

  Type();




  Type(Topology topol, Kind kind);



  Topology topology() const;
  bool is_necessarily_closed() const;
  bool is_line_or_equality() const;
  bool is_ray_or_point_or_inequality() const;

  void set_necessarily_closed();
  void set_not_necessarily_closed();
  void set_is_line_or_equality();
  void set_is_ray_or_point_or_inequality();


private:

  typedef unsigned int flags_t;


  Type(flags_t mask);


  flags_t flags;



  static const flags_t NNC = NOT_NECESSARILY_CLOSED << 0;
  static const flags_t RPI = RAY_OR_POINT_OR_INEQUALITY << 1;



  bool test_all(flags_t mask) const;


  void set(flags_t mask);


  void reset(flags_t mask);
};
class Parma_Polyhedra_Library::Row::Impl {
public:
  void* operator new(size_t fixed_size, dimension_type capacity);


  void operator delete(void* p);



  void operator delete(void* p, dimension_type capacity);



  Impl(Type t, dimension_type sz);


  Impl(const Impl& y);


  Impl(const Impl& y, dimension_type sz);






  ~Impl();







  void resize_no_copy(dimension_type new_sz);


  void grow_no_copy(dimension_type new_sz);


  void shrink(dimension_type new_sz);


  static dimension_type max_size();




  dimension_type size() const;


  void set_size(dimension_type new_sz);


  void bump_size();





  Integer& operator[](dimension_type k);


  const Integer& operator[](dimension_type k) const;


private:
  friend class Parma_Polyhedra_Library::Row;


  dimension_type size_;


  Type type;


  Integer vec_[



  ];


  Impl();


  Impl& operator=(const Impl&);


  void copy_construct(const Impl& y);
};

namespace std {





void swap(Parma_Polyhedra_Library::Row& x,
   Parma_Polyhedra_Library::Row& y);





void iter_swap(std::vector<Parma_Polyhedra_Library::Row>::iterator x,
        std::vector<Parma_Polyhedra_Library::Row>::iterator y);

}
namespace Parma_Polyhedra_Library {

inline void*
Row::Impl::operator new(const size_t fixed_size,
   const dimension_type capacity) {

  return ::operator new(fixed_size + capacity*sizeof(Integer));




}

inline void
Row::Impl::operator delete(void* p) {
  ::operator delete(p);
}

inline void
Row::Impl::operator delete(void* p, dimension_type) {
  ::operator delete(p);
}

inline dimension_type
Row::Impl::max_size() {
  return size_t(-1)/sizeof(Integer);
}

inline dimension_type
Row::Impl::size() const {
  return size_;
}

inline void
Row::Impl::set_size(const dimension_type new_sz) {
  size_ = new_sz;
}

inline void
Row::Impl::bump_size() {
  ++size_;
}

inline void
Row::Impl::resize_no_copy(const dimension_type new_sz) {
  if (new_sz < size())
    shrink(new_sz);
  else
    grow_no_copy(new_sz);
}

inline
Row::Impl::Impl(const Type t, const dimension_type sz)
  : size_(0), type(t) {
  grow_no_copy(sz);
}

inline
Row::Impl::Impl(const Impl& y)
  : size_(0), type(y.type) {
  copy_construct(y);
}

inline
Row::Impl::Impl(const Impl& y, const dimension_type sz)
  : size_(0), type(y.type) {
  copy_construct(y);
  grow_no_copy(sz);
}

inline
Row::Impl::~Impl() {

  shrink(0);



}

inline Integer&
Row::Impl::operator[](const dimension_type k) {
  assert(k < size());
  return vec_[k];
}

inline const Integer&
Row::Impl::operator[](const dimension_type k) const {
  assert(k < size());
  return vec_[k];
}

inline
Row::Type::Type()
  : flags(0) {
}

inline
Row::Type::Type(const Topology topol, const Kind kind)
  : flags(static_cast<flags_t>(topol | (kind << 1))) {
}

inline
Row::Type::Type(const flags_t mask)
  : flags(mask) {
}

inline bool
Row::Type::test_all(const flags_t mask) const {
  return (flags & mask) == mask;
}

inline void
Row::Type::set(const flags_t mask) {
  flags |= mask;
}

inline void
Row::Type::reset(const flags_t mask) {
  flags &= ~mask;
}

inline bool
Row::Type::is_ray_or_point_or_inequality() const {
  return test_all(RPI);
}

inline void
Row::Type::set_is_ray_or_point_or_inequality() {
  set(RPI);
}

inline bool
Row::Type::is_line_or_equality() const {
  return !is_ray_or_point_or_inequality();
}

inline void
Row::Type::set_is_line_or_equality() {
  reset(RPI);
}

inline Topology
Row::Type::topology() const {
  return test_all(NNC) ? NOT_NECESSARILY_CLOSED : NECESSARILY_CLOSED;
}

inline bool
Row::Type::is_necessarily_closed() const {
  return !test_all(NNC);
}

inline void
Row::Type::set_necessarily_closed() {
  reset(NNC);
}

inline void
Row::Type::set_not_necessarily_closed() {
  set(NNC);
}

inline dimension_type
Row::max_size() {
  return Impl::max_size();
}

inline dimension_type
Row::size() const {
  return impl->size();
}

inline Row::Type
Row::type() const {
  return impl->type;
}

inline bool
Row::is_necessarily_closed() const {
  return type().is_necessarily_closed();
}

inline dimension_type
Row::space_dimension() const {
  const dimension_type sz = size();
  return (sz == 0)
    ? 0
    : sz - (is_necessarily_closed() ? 1 : 2);
}
inline
Row::Row()
  : impl(0) {
}

inline void
Row::construct(const Type t,
        const dimension_type sz,

        const

        dimension_type capacity) {
  assert(capacity >= sz);




  impl = new (capacity) Impl(t, sz);



}

inline void
Row::construct(const Type t, const dimension_type sz) {
  construct(t, sz, sz);
}

inline
Row::Row(const Type t,
  const dimension_type sz, const dimension_type capacity) {
  construct(t, sz, capacity);
}

inline
Row::Row(const Type t, const dimension_type sz) {
  construct(t, sz);
}

inline
Row::Row(const Row& y)
  : impl(y.impl
  ? new (compute_capacity(y.size())) Impl(*y.impl)
  : 0) {







}

inline
Row::Row(const Row& y,

  const

  dimension_type capacity) {
  assert(capacity >= y.size());




  impl = y.impl ? new (capacity) Impl(*y.impl) : 0;



}

inline
Row::Row(const Row& y,
  const dimension_type sz,

  const

  dimension_type capacity) {
  assert(capacity >= sz);
  assert(sz >= y.size());




  impl = y.impl ? new (capacity) Impl(*y.impl, sz) : 0;



}

inline
Row::~Row() {
  delete impl;
}

inline void
Row::resize_no_copy(const dimension_type new_sz) {
  assert(impl);



  impl->resize_no_copy(new_sz);
}

inline void
Row::grow_no_copy(const dimension_type new_sz) {
  assert(impl);



  impl->grow_no_copy(new_sz);
}

inline void
Row::shrink(const dimension_type new_sz) {
  assert(impl);
  impl->shrink(new_sz);
}

inline void
Row::swap(Row& y) {
  std::swap(impl, y.impl);



}

inline void
Row::assign(Row& y) {
  impl = y.impl;



}

inline Row&
Row::operator=(const Row& y) {

  Row tmp(y);

  swap(tmp);

  return *this;
}

inline bool
Row::is_line_or_equality() const {
  return type().is_line_or_equality();
}

inline bool
Row::is_ray_or_point_or_inequality() const {
  return type().is_ray_or_point_or_inequality();
}

inline Topology
Row::topology() const {
  return type().topology();
}

inline void
Row::set_is_line_or_equality() {
  impl->type.set_is_line_or_equality();
}

inline void
Row::set_is_ray_or_point_or_inequality() {
  impl->type.set_is_ray_or_point_or_inequality();
}

inline void
Row::set_necessarily_closed() {
  impl->type.set_necessarily_closed();
}

inline void
Row::set_not_necessarily_closed() {
  impl->type.set_not_necessarily_closed();
}

inline Integer&
Row::operator[](const dimension_type k) {
  return (*impl)[k];
}

inline const Integer&
Row::operator[](const dimension_type k) const {
  return (*impl)[k];
}

inline const Integer&
Row::inhomogeneous_term() const {
  return (*this)[0];
}

inline const Integer&
Row::coefficient(const dimension_type k) const {
  return (*this)[k+1];
}

inline void
Row::strong_normalize() {
  normalize();
  sign_normalize();
}


inline bool
operator==(const Row& x, const Row& y) {
  return compare(x, y) == 0;
}


inline bool
operator!=(const Row& x, const Row& y) {
  return compare(x, y) != 0;
}


inline bool
operator<=(const Row& x, const Row& y) {
  return compare(x, y) <= 0;
}


inline bool
operator<(const Row& x, const Row& y) {
  return compare(x, y) < 0;
}


inline bool
operator>=(const Row& x, const Row& y) {
  return compare(x, y) >= 0;
}


inline bool
operator>(const Row& x, const Row& y) {
  return compare(x, y) > 0;
}

}


namespace std {


inline void
swap(Parma_Polyhedra_Library::Row& x, Parma_Polyhedra_Library::Row& y) {
  x.swap(y);
}


inline void
iter_swap(std::vector<Parma_Polyhedra_Library::Row>::iterator x,
   std::vector<Parma_Polyhedra_Library::Row>::iterator y) {
  swap(*x, *y);
}

}
class Parma_Polyhedra_Library::Matrix {
protected:

  static dimension_type max_num_rows();


  static dimension_type max_num_columns();





  Matrix(Topology topol);
  Matrix(Topology topol, dimension_type n_rows, dimension_type n_columns);


  Matrix(const Matrix& y);
  Matrix(Matrix& y, dimension_type first_stolen);


  ~Matrix();


  Matrix& operator=(const Matrix& y);

public:





  class const_iterator {
  private:
    typedef std::vector<Row>::const_iterator Iter;

    Iter i;

  public:
    typedef std::forward_iterator_tag iterator_category;
    typedef std::iterator_traits<Iter>::value_type value_type;
    typedef std::iterator_traits<Iter>::difference_type difference_type;
    typedef std::iterator_traits<Iter>::pointer pointer;
    typedef std::iterator_traits<Iter>::reference reference;


    const_iterator();




    explicit const_iterator(const Iter& b);


    const_iterator(const const_iterator& y);


    const_iterator& operator=(const const_iterator& y);


    reference operator*() const;


    pointer operator->() const;


    const_iterator& operator++();


    const_iterator operator++(int);




    bool operator==(const const_iterator& y) const;




    bool operator!=(const const_iterator& y) const;
  };





  const_iterator begin() const;


  const_iterator end() const;

private:

  std::vector<Row> rows;


  Topology row_topology;


  dimension_type row_size;




  dimension_type row_capacity;


  dimension_type index_first_pending;





  bool sorted;

public:

  void swap(Matrix& y);


  void set_sorted(bool value);


  void set_necessarily_closed();


  void set_not_necessarily_closed();


  void set_rows_topology();


  void unset_pending_rows();


  void set_index_first_pending_row(dimension_type first_pending);
  void grow(dimension_type new_n_rows, dimension_type new_n_columns);




  void remove_columns(dimension_type new_n_columns);
  void resize_no_copy(dimension_type new_n_rows, dimension_type new_n_columns);
  void add_zero_columns(dimension_type n);
  void add_rows_and_columns(dimension_type n);


  void swap_columns(dimension_type i, dimension_type j);
  void permute_columns(const std::vector<dimension_type>& cycles);




  Topology topology() const;




  bool is_necessarily_closed() const;


  bool is_sorted() const;
  dimension_type space_dimension() const;




  dimension_type num_columns() const;


  dimension_type num_rows() const;


  dimension_type first_pending_row() const;


  dimension_type num_pending_rows() const;




  dimension_type num_lines_or_equalities() const;





  bool check_sorted() const;




  Row& operator[](dimension_type k);


  const Row& operator[](dimension_type k) const;



  void normalize();


  void strong_normalize();




  void sort_rows();




  void sort_rows(dimension_type first_row, dimension_type last_row);




  void sort_pending_and_remove_duplicates();


  void add_row(const Row& row);


  void add_pending_row(const Row& row);


  void add_pending_row(Row::Type type);




  void insert(const Row& row);




  void insert_pending(const Row& row);





  void add_rows(const Matrix& y);


  void add_pending_rows(const Matrix& y);
  void merge_rows_assign(const Matrix& y);


  void clear();
  void ascii_dump(std::ostream& s) const;
  bool ascii_load(std::istream& s);




  void erase_to_end(dimension_type first_to_erase);
  void sort_and_remove_with_sat(SatMatrix& sat);
  dimension_type gauss();
  void back_substitute(dimension_type rank);






  void gram_shmidt();
  bool OK(bool check_strong_normalized = true) const;
};

namespace std {





void swap(Parma_Polyhedra_Library::Matrix& x,
   Parma_Polyhedra_Library::Matrix& y);

}


namespace Parma_Polyhedra_Library {





bool operator==(const Matrix& x, const Matrix& y);





bool operator!=(const Matrix& x, const Matrix& y);

}
namespace Parma_Polyhedra_Library {

inline dimension_type
Matrix::max_num_rows() {


  return std::vector<Row>().max_size();
}

inline dimension_type
Matrix::max_num_columns() {
  return Row::max_size();
}

inline
Matrix::const_iterator::const_iterator()
  : i(Iter()) {
}

inline
Matrix::const_iterator::const_iterator(const Iter& b)
  : i(b) {
}

inline
Matrix::const_iterator::const_iterator(const const_iterator& y)
  : i(y.i) {
}

inline Matrix::const_iterator&
Matrix::const_iterator::operator=(const const_iterator& y) {
  i = y.i;
  return *this;
}

inline Matrix::const_iterator::reference
Matrix::const_iterator::operator*() const {
  return *i;
}

inline Matrix::const_iterator::pointer
Matrix::const_iterator::operator->() const {
  return &*i;
}

inline Matrix::const_iterator&
Matrix::const_iterator::operator++() {
  ++i;
  return *this;
}

inline Matrix::const_iterator
Matrix::const_iterator::operator++(int) {
  return const_iterator(i++);
}

inline bool
Matrix::const_iterator::operator==(const const_iterator& y) const {
  return i == y.i;
}

inline bool
Matrix::const_iterator::operator!=(const const_iterator& y) const {
  return !operator==(y);
}

inline Matrix::const_iterator
Matrix::begin() const {
  return const_iterator(rows.begin());
}

inline Matrix::const_iterator
Matrix::end() const {
  return const_iterator(rows.end());
}

inline void
Matrix::swap(Matrix& y) {
  std::swap(rows, y.rows);
  std::swap(row_topology, y.row_topology);
  std::swap(row_size, y.row_size);
  std::swap(row_capacity, y.row_capacity);
  std::swap(index_first_pending, y.index_first_pending);
  std::swap(sorted, y.sorted);
}

inline
Matrix::Matrix(const Topology topol)
  : rows(),
    row_topology(topol),
    row_size(0),
    row_capacity(0),
    index_first_pending(0),
    sorted(true) {
}

inline
Matrix::~Matrix() {
}

inline Row&
Matrix::operator[](const dimension_type k) {
  assert(k < rows.size());
  return rows[k];
}

inline const Row&
Matrix::operator[](const dimension_type k) const {
  assert(k < rows.size());
  return rows[k];
}

inline dimension_type
Matrix::num_rows() const {
  return rows.size();
}

inline dimension_type
Matrix::first_pending_row() const {
  return index_first_pending;
}

inline dimension_type
Matrix::num_pending_rows() const {
  assert(num_rows() >= first_pending_row());
  return num_rows() - first_pending_row();
}

inline void
Matrix::unset_pending_rows() {
  index_first_pending = num_rows();
}

inline void
Matrix::set_index_first_pending_row(const dimension_type first_pending) {
  index_first_pending = first_pending;
}

inline void
Matrix::set_necessarily_closed() {
  row_topology = NECESSARILY_CLOSED;
  if (num_rows() > 0)
    set_rows_topology();
}

inline void
Matrix::set_not_necessarily_closed() {
  row_topology = NOT_NECESSARILY_CLOSED;
  if (num_rows() > 0)
    set_rows_topology();
}

inline bool
Matrix::is_necessarily_closed() const {
  return row_topology == NECESSARILY_CLOSED;
}

inline Topology
Matrix::topology() const {
  return row_topology;
}

inline void
Matrix::set_sorted(const bool value) {
  sorted = value;
}

inline bool
Matrix::is_sorted() const {



  assert(!sorted || check_sorted());
  return sorted;
}


inline dimension_type
Matrix::num_columns() const {
  return row_size;
}

inline dimension_type
Matrix::space_dimension() const {
  const dimension_type n_columns = num_columns();
  return (n_columns == 0)
    ? 0
    : n_columns - (is_necessarily_closed() ? 1 : 2);
}


inline bool
operator!=(const Matrix& x, const Matrix& y) {
  return !(x == y);
}

inline void
Matrix::add_zero_columns(const dimension_type n) {
  assert(n > 0);
  grow(num_rows(), num_columns() + n);
}

inline void
Matrix::erase_to_end(const dimension_type first_to_erase) {
  assert(first_to_erase <= rows.size());
  if (first_to_erase < rows.size())
    rows.erase(rows.begin() + first_to_erase, rows.end());
}

inline void
Matrix::clear() {


  std::vector<Row>().swap(rows);
  row_size = 0;
  row_capacity = 0;
  index_first_pending = 0;
  sorted = true;
}

inline void
Matrix::remove_columns(const dimension_type new_n_columns) {
  assert(new_n_columns < num_columns());



  resize_no_copy(num_rows(), new_n_columns);


  strong_normalize();
}

}

namespace std {


inline void
swap(Parma_Polyhedra_Library::Matrix& x,
     Parma_Polyhedra_Library::Matrix& y) {
  x.swap(y);
}

}
namespace Parma_Polyhedra_Library {

namespace IO_Operators {







std::ostream& operator<<(std::ostream& s, const ConSys& cs);

}


bool operator==(const Polyhedron& x, const Polyhedron& y);

}


namespace std {



void swap(Parma_Polyhedra_Library::ConSys& x,
   Parma_Polyhedra_Library::ConSys& y);

}
class Parma_Polyhedra_Library::ConSys : private Matrix {
public:

  ConSys();


  ConSys(const Constraint& c);


  ConSys(const ConSys& cs);


  ~ConSys();


  ConSys& operator=(const ConSys& y);


  static dimension_type max_space_dimension();


  dimension_type space_dimension() const;




  void clear();




  void insert(const Constraint& c);




  static const ConSys& zero_dim_empty();
  class const_iterator
    : public std::iterator<std::forward_iterator_tag,
      Constraint,
      void,
      const Constraint*,
      const Constraint&> {
  public:

    const_iterator();


    const_iterator(const const_iterator& y);


    ~const_iterator();


    const_iterator& operator=(const const_iterator& y);


    const Constraint& operator*() const;


    const Constraint* operator->() const;


    const_iterator& operator++();


    const_iterator operator++(int);




    bool operator==(const const_iterator& y) const;




    bool operator!=(const const_iterator& y) const;

  private:
    friend class ConSys;


    Parma_Polyhedra_Library::Matrix::const_iterator i;


    const Parma_Polyhedra_Library::Matrix* csp;


    const_iterator(const Parma_Polyhedra_Library::Matrix::const_iterator& iter,
     const ConSys& csys);


    void skip_forward();
  };





  const_iterator begin() const;


  const_iterator end() const;
  bool OK() const;
  void ascii_dump(std::ostream& s) const;
  bool ascii_load(std::istream& s);

private:
  friend class const_iterator;
  friend class Parma_Polyhedra_Library::Polyhedron;

  friend bool
  Parma_Polyhedra_Library::operator==(const Polyhedron& x,
          const Polyhedron& y);

  friend void std::swap(Parma_Polyhedra_Library::ConSys& x,
   Parma_Polyhedra_Library::ConSys& y);


  ConSys(Topology topol);





  ConSys(Topology topol, dimension_type n_rows, dimension_type n_columns);
  ConSys(ConSys& y, dimension_type first_stolen);


  void swap(ConSys& y);
  bool adjust_topology_and_dimension(Topology topol,
         dimension_type num_dimensions);




  bool has_strict_inequalities() const;


  Constraint& operator[](dimension_type k);


  const Constraint& operator[](dimension_type k) const;


  bool satisfies_all_constraints(const Generator& g) const;
  void affine_preimage(dimension_type v,
         const LinExpression& expr,
         const Integer& denominator);


  dimension_type num_equalities() const;


  dimension_type num_inequalities() const;





  void insert_pending(const Constraint& c);
};
namespace Parma_Polyhedra_Library {

namespace IO_Operators {







std::ostream& operator<<(std::ostream& s, const GenSys& gs);

}


bool operator==(const Polyhedron& x, const Polyhedron& y);

}


namespace std {



void swap(Parma_Polyhedra_Library::GenSys& x,
   Parma_Polyhedra_Library::GenSys& y);

}
class Parma_Polyhedra_Library::GenSys : private Matrix {
public:

  GenSys();


  GenSys(const Generator& g);


  GenSys(const GenSys& gs);


  ~GenSys();


  GenSys& operator=(const GenSys& y);


  static dimension_type max_space_dimension();


  dimension_type space_dimension() const;




  void clear();




  void insert(const Generator& g);




  static const GenSys& zero_dim_univ();
  class const_iterator
    : public std::iterator<std::forward_iterator_tag,
    Generator,
    void,
    const Generator*,
                                const Generator&> {
  public:

    const_iterator();


    const_iterator(const const_iterator& y);


    ~const_iterator();


    const_iterator& operator=(const const_iterator& y);


    const Generator& operator*() const;


    const Generator* operator->() const;


    const_iterator& operator++();


    const_iterator operator++(int);




    bool operator==(const const_iterator& y) const;




    bool operator!=(const const_iterator& y) const;

  private:
    friend class GenSys;


    Parma_Polyhedra_Library::Matrix::const_iterator i;


    const Parma_Polyhedra_Library::Matrix* gsp;


    const_iterator(const Parma_Polyhedra_Library::Matrix::const_iterator& iter,
     const GenSys& gsys);




    void skip_forward();
  };





  const_iterator begin() const;


  const_iterator end() const;







  bool OK() const;
  void ascii_dump(std::ostream& s) const;
  bool ascii_load(std::istream& s);

private:
  friend class const_iterator;
  friend class Parma_Polyhedra_Library::Polyhedron;

  friend bool
  Parma_Polyhedra_Library::operator==(const Polyhedron& x,
          const Polyhedron& y);

  friend void std::swap(Parma_Polyhedra_Library::GenSys& x,
   Parma_Polyhedra_Library::GenSys& y);


  GenSys(Topology topol);





  GenSys(Topology topol, dimension_type n_rows, dimension_type n_columns);
  GenSys(GenSys& y, dimension_type first_stolen);


  void swap(GenSys& y);
  bool adjust_topology_and_dimension(Topology topol,
         dimension_type num_dimensions);
  void add_corresponding_points();




  bool has_points() const;
  void add_corresponding_closure_points();




  bool has_closure_points() const;


  Generator& operator[](dimension_type k);


  const Generator& operator[](dimension_type k) const;




  Parma_Polyhedra_Library::Poly_Con_Relation
  relation_with(const Constraint& c) const;


  bool satisfied_by_all_generators(const Constraint& c) const;
  void affine_image(dimension_type v,
      const LinExpression& expr,
      const Integer& denominator);


  dimension_type num_lines() const;


  dimension_type num_rays() const;






  void remove_invalid_lines_and_rays();





  void insert_pending(const Generator& g);
};
namespace Parma_Polyhedra_Library {




LinExpression operator+(const LinExpression& e1, const LinExpression& e2);


LinExpression operator+(const Integer& n, const LinExpression& e);


LinExpression operator+(const LinExpression& e, const Integer& n);



LinExpression operator+(const LinExpression& e);


LinExpression operator-(const LinExpression& e);



LinExpression operator-(const LinExpression& e1, const LinExpression& e2);


LinExpression operator-(const Integer& n, const LinExpression& e);


LinExpression operator-(const LinExpression& e, const Integer& n);



LinExpression operator*(const Integer& n, const LinExpression& e);


LinExpression operator*(const LinExpression& e, const Integer& n);



LinExpression& operator+=(LinExpression& e1, const LinExpression& e2);


LinExpression& operator+=(LinExpression& e, const Variable v);


LinExpression& operator+=(LinExpression& e, const Integer& n);



LinExpression& operator-=(LinExpression& e1, const LinExpression& e2);


LinExpression& operator-=(LinExpression& e, const Variable v);


LinExpression& operator-=(LinExpression& e, const Integer& n);



LinExpression& operator*=(LinExpression& e, const Integer& n);

namespace IO_Operators {



std::ostream& operator<<(std::ostream& s, const LinExpression& e);

}

}


namespace std {



void swap(Parma_Polyhedra_Library::LinExpression& x,
   Parma_Polyhedra_Library::LinExpression& y);

}
class Parma_Polyhedra_Library::LinExpression : private Row {
public:

  LinExpression();


  LinExpression(const LinExpression& e);


  ~LinExpression();




  explicit LinExpression(const Integer& n);




  LinExpression(const Variable v);
  explicit LinExpression(const Constraint& c);
  explicit LinExpression(const Generator& g);


  dimension_type space_dimension() const;


  const Integer& coefficient(Variable v) const;


  const Integer& inhomogeneous_term() const;


  static const LinExpression& zero();

private:
  friend class Parma_Polyhedra_Library::Constraint;
  friend class Parma_Polyhedra_Library::Generator;
  friend class Parma_Polyhedra_Library::Polyhedron;



  friend class Parma_Polyhedra_Library::ConSys;



  friend class Parma_Polyhedra_Library::GenSys;

  friend void std::swap(Parma_Polyhedra_Library::LinExpression& x,
         Parma_Polyhedra_Library::LinExpression& y);


  LinExpression(const LinExpression& e, dimension_type sz);






  LinExpression(dimension_type sz, bool);


  void swap(LinExpression& y);

  friend LinExpression
  Parma_Polyhedra_Library::operator+(const LinExpression& e1,
         const LinExpression& e2);
  friend LinExpression
  Parma_Polyhedra_Library::operator+(const Integer& n,
         const LinExpression& e);
  friend LinExpression
  Parma_Polyhedra_Library::operator+(const LinExpression& e,
         const Integer& n);

  friend LinExpression
  Parma_Polyhedra_Library::operator-(const LinExpression& e);

  friend LinExpression
  Parma_Polyhedra_Library::operator-(const LinExpression& e1,
         const LinExpression& e2);
  friend LinExpression
  Parma_Polyhedra_Library::operator-(const Integer& n,
         const LinExpression& e);
  friend LinExpression
  Parma_Polyhedra_Library::operator-(const LinExpression& e,
         const Integer& n);

  friend LinExpression
  Parma_Polyhedra_Library::operator*(const Integer& n,
         const LinExpression& e);
  friend LinExpression
  Parma_Polyhedra_Library::operator*(const LinExpression& e,
         const Integer& n);

  friend LinExpression&
  Parma_Polyhedra_Library::operator+=(LinExpression& e1,
          const LinExpression& e2);
  friend LinExpression&
  Parma_Polyhedra_Library::operator+=(LinExpression& e,
          const Variable v);
  friend LinExpression&
  Parma_Polyhedra_Library::operator+=(LinExpression& e,
          const Integer& n);

  friend LinExpression&
  Parma_Polyhedra_Library::operator-=(LinExpression& e1,
          const LinExpression& e2);
  friend LinExpression&
  Parma_Polyhedra_Library::operator-=(LinExpression& e,
          const Variable v);
  friend LinExpression&
  Parma_Polyhedra_Library::operator-=(LinExpression& e,
          const Integer& n);

  friend LinExpression&
  Parma_Polyhedra_Library::operator*=(LinExpression& e,
          const Integer& n);

  friend std::ostream&
  Parma_Polyhedra_Library::IO_Operators::operator<<(std::ostream& s,
          const LinExpression& e);
};
namespace Parma_Polyhedra_Library {

inline
LinExpression::LinExpression()
  : Row(Row::Type(), 1) {
}

inline
LinExpression::LinExpression(dimension_type sz, bool)
  : Row(Row::Type(), sz) {
}

inline
LinExpression::LinExpression(const Variable v)
  : Row(Row::Type(), v.id() + 2) {
  (*this)[v.id() + 1] = 1;
}

inline
LinExpression::LinExpression(const LinExpression& e)
  : Row(e) {
}

inline
LinExpression::~LinExpression() {
}

inline
LinExpression::LinExpression(const LinExpression& e, dimension_type sz)
  : Row(e, sz, sz) {
}

inline
LinExpression::LinExpression(const Integer& n)
  : Row(Row::Type(), 1) {
  (*this)[0] = n;
}

inline dimension_type
LinExpression::space_dimension() const {
  return size() - 1;
}

inline const Integer&
LinExpression::coefficient(Variable v) const {
  if (v.id() >= space_dimension())
    return Integer_zero();
  return Row::coefficient(v.id());
}

inline const Integer&
LinExpression::inhomogeneous_term() const {
  return Row::inhomogeneous_term();
}

inline const LinExpression&
LinExpression::zero() {
  static LinExpression z = LinExpression(Integer_zero());
  return z;
}


inline LinExpression
operator+(const LinExpression& e) {
  return e;
}


inline LinExpression
operator+(const LinExpression& e, const Integer& n) {
  return n + e;
}


inline LinExpression
operator-(const LinExpression& e, const Integer& n) {
  return -n + e;
}


inline LinExpression
operator*(const LinExpression& e, const Integer& n) {
  return n * e;
}


inline LinExpression&
operator+=(LinExpression& e, const Integer& n) {
  e[0] += n;
  return e;
}


inline LinExpression&
operator-=(LinExpression& e, const Integer& n) {
  e[0] -= n;
  return e;
}

inline void
LinExpression::swap(LinExpression& y) {
  Row::swap(y);
}

}


namespace std {


inline void
swap(Parma_Polyhedra_Library::LinExpression& x,
     Parma_Polyhedra_Library::LinExpression& y) {
  x.swap(y);
}

}
namespace Parma_Polyhedra_Library {

namespace IO_Operators {



std::ostream& operator<<(std::ostream& s, const Constraint& c);

}





Constraint operator==(const LinExpression& e1, const LinExpression& e2);


Constraint operator==(const LinExpression& e, const Integer& n);


Constraint operator==(const Integer& n, const LinExpression& e);



Constraint operator<=(const LinExpression& e1, const LinExpression& e2);


Constraint operator<=(const LinExpression& e, const Integer& n);


Constraint operator<=(const Integer& n, const LinExpression& e);



Constraint operator>=(const LinExpression& e1, const LinExpression& e2);


Constraint operator>=(const LinExpression& e, const Integer& n);


Constraint operator>=(const Integer& n, const LinExpression& e);



Constraint operator<(const LinExpression& e1, const LinExpression& e2);


Constraint operator<(const LinExpression& e, const Integer& n);


Constraint operator<(const Integer& n, const LinExpression& e);



Constraint operator>(const LinExpression& e1, const LinExpression& e2);


Constraint operator>(const LinExpression& e, const Integer& n);


Constraint operator>(const Integer& n, const LinExpression& e);

}


namespace std {



void swap(Parma_Polyhedra_Library::Constraint& x,
   Parma_Polyhedra_Library::Constraint& y);

}
class Parma_Polyhedra_Library::Constraint : private Row {
public:

  Constraint(const Constraint& c);


  ~Constraint();


  Constraint& operator=(const Constraint& c);


  dimension_type space_dimension() const;


  enum Type {

    EQUALITY,

    NONSTRICT_INEQUALITY,

    STRICT_INEQUALITY
  };


  Type type() const;




  bool is_equality() const;




  bool is_inequality() const;




  bool is_nonstrict_inequality() const;




  bool is_strict_inequality() const;






  const Integer& coefficient(Variable v) const;


  const Integer& inhomogeneous_term() const;


  static const Constraint& zero_dim_false();




  static const Constraint& zero_dim_positivity();


  bool OK() const;

private:
  friend class Parma_Polyhedra_Library::ConSys;
  friend class Parma_Polyhedra_Library::ConSys::const_iterator;
  friend class Parma_Polyhedra_Library::Polyhedron;


  friend class Parma_Polyhedra_Library::GenSys;

  friend const Integer&
  Parma_Polyhedra_Library::operator*(const Constraint& c,
         const Generator& g);
  friend const Integer&
  Parma_Polyhedra_Library::reduced_scalar_product(const Constraint& c,
        const Generator& g);
  friend
  Parma_Polyhedra_Library::LinExpression::LinExpression(const Constraint& c);
  friend void std::swap(Parma_Polyhedra_Library::Constraint& x,
   Parma_Polyhedra_Library::Constraint& y);


  Constraint();




  explicit Constraint(LinExpression& e);




  Constraint(Row::Type t, dimension_type sz);


  void swap(Constraint& y);




  void
  throw_dimension_incompatible(const char* method,
          const char* name_var,
          Variable v) const;

  friend Constraint
  Parma_Polyhedra_Library::operator==(const LinExpression& e1,
          const LinExpression& e2);
  friend Constraint
  Parma_Polyhedra_Library::operator==(const LinExpression& e,
          const Integer& n);
  friend Constraint
  Parma_Polyhedra_Library::operator==(const Integer& n,
          const LinExpression& e);

  friend Constraint
  Parma_Polyhedra_Library::operator>=(const LinExpression& e1,
          const LinExpression& e2);
  friend Constraint
  Parma_Polyhedra_Library::operator>=(const LinExpression& e,
          const Integer& n);
  friend Constraint
  Parma_Polyhedra_Library::operator>=(const Integer& n,
          const LinExpression& e);

  friend Constraint
  Parma_Polyhedra_Library::operator<=(const LinExpression& e1,
          const LinExpression& e2);
  friend Constraint
  Parma_Polyhedra_Library::operator<=(const LinExpression& e,
          const Integer& n);
  friend Constraint
  Parma_Polyhedra_Library::operator<=(const Integer& n,
          const LinExpression& e);

  friend Constraint
  Parma_Polyhedra_Library::operator>(const LinExpression& e1,
         const LinExpression& e2);
  friend Constraint
  Parma_Polyhedra_Library::operator>(const LinExpression& e,
         const Integer& n);
  friend Constraint
  Parma_Polyhedra_Library::operator>(const Integer& n,
         const LinExpression& e);

  friend Constraint
  Parma_Polyhedra_Library::operator<(const LinExpression& e1,
         const LinExpression& e2);
  friend Constraint
  Parma_Polyhedra_Library::operator<(const LinExpression& e,
         const Integer& n);
  friend Constraint
  Parma_Polyhedra_Library::operator<(const Integer& n,
         const LinExpression& e);


  Constraint(const Constraint& c, dimension_type sz);




  static Constraint construct_epsilon_geq_zero();


  static const Constraint& epsilon_geq_zero();




  static const Constraint& epsilon_leq_one();
  bool is_trivial_true() const;
  bool is_trivial_false() const;


  void set_is_equality();


  void set_is_inequality();
};
namespace Parma_Polyhedra_Library {

inline
Constraint::Constraint(LinExpression& e) {
  Row::swap(e);
}

inline
Constraint::Constraint(const Constraint& c)
  : Row(c) {
}

inline
Constraint::Constraint(const Constraint& c, const dimension_type sz)
  : Row(c, sz, sz) {
}

inline
Constraint::Constraint(Row::Type t, const dimension_type sz)
  : Row(t, sz) {
}

inline
Constraint::~Constraint() {
}

inline Constraint&
Constraint::operator=(const Constraint& c) {
  Row::operator=(c);
  return *this;
}

inline dimension_type
Constraint::space_dimension() const {
  return Row::space_dimension();
}

inline bool
Constraint::is_equality() const {
  return is_line_or_equality();
}

inline bool
Constraint::is_inequality() const {
  return is_ray_or_point_or_inequality();
}

inline Constraint::Type
Constraint::type() const {
  if (is_equality())
    return EQUALITY;
  if (is_necessarily_closed())
    return NONSTRICT_INEQUALITY;
  else
    return ((*this)[size() - 1] < 0)
      ? STRICT_INEQUALITY
      : NONSTRICT_INEQUALITY;
}

inline bool
Constraint::is_nonstrict_inequality() const {
  return type() == NONSTRICT_INEQUALITY;
}

inline bool
Constraint::is_strict_inequality() const {
  return type() == STRICT_INEQUALITY;
}

inline void
Constraint::set_is_equality() {
  set_is_line_or_equality();
}

inline void
Constraint::set_is_inequality() {
  set_is_ray_or_point_or_inequality();
}

inline const Integer&
Constraint::coefficient(const Variable v) const {
  const dimension_type v_id = v.id();
  if (v_id >= space_dimension())
    throw_dimension_incompatible("coefficient(v)", "v", v);
  return Row::coefficient(v_id);
}

inline const Integer&
Constraint::inhomogeneous_term() const {
  return Row::inhomogeneous_term();
}


inline Constraint
operator==(const LinExpression& e1, const LinExpression& e2) {
  LinExpression diff = e1 - e2;
  Constraint c(diff);
  c.set_is_equality();

  c.strong_normalize();
  return c;
}


inline Constraint
operator>=(const LinExpression& e1, const LinExpression& e2) {
  LinExpression diff = e1 - e2;
  Constraint c(diff);
  c.set_is_inequality();

  c.normalize();
  return c;
}


inline Constraint
operator>(const LinExpression& e1, const LinExpression& e2) {
  LinExpression diff;


  const dimension_type e1_dim = e1.space_dimension();
  const dimension_type e2_dim = e2.space_dimension();
  if (e1_dim > e2_dim)
    diff -= Variable(e1_dim);
  else
    diff -= Variable(e2_dim);
  diff += e1;
  diff -= e2;

  Constraint c(diff);
  c.set_not_necessarily_closed();
  c.set_is_inequality();
  return c;
}


inline Constraint
operator==(const Integer& n, const LinExpression& e) {
  LinExpression diff = n - e;
  Constraint c(diff);
  c.set_is_equality();

  c.strong_normalize();
  return c;
}


inline Constraint
operator>=(const Integer& n, const LinExpression& e) {
  LinExpression diff = n - e;
  Constraint c(diff);
  c.set_is_inequality();

  c.normalize();
  return c;
}


inline Constraint
operator>(const Integer& n, const LinExpression& e) {
  LinExpression diff;


  diff -= Variable(e.space_dimension());
  diff += n;
  diff -= e;

  Constraint c(diff);
  c.set_not_necessarily_closed();
  c.set_is_inequality();
  return c;
}


inline Constraint
operator==(const LinExpression& e, const Integer& n) {
  LinExpression diff = e - n;
  Constraint c(diff);
  c.set_is_equality();

  c.strong_normalize();
  return c;
}


inline Constraint
operator>=(const LinExpression& e, const Integer& n) {
  LinExpression diff = e - n;
  Constraint c(diff);
  c.set_is_inequality();

  c.normalize();
  return c;
}


inline Constraint
operator>(const LinExpression& e, const Integer& n) {
  LinExpression diff;


  diff -= Variable(e.space_dimension());
  diff += e;
  diff -= n;

  Constraint c(diff);
  c.set_not_necessarily_closed();
  c.set_is_inequality();
  return c;
}


inline Constraint
operator<=(const LinExpression& e1, const LinExpression& e2) {
  return e2 >= e1;
}


inline Constraint
operator<=(const Integer& n, const LinExpression& e) {
  return e >= n;
}


inline Constraint
operator<=(const LinExpression& e, const Integer& n) {
  return n >= e;
}


inline Constraint
operator<(const LinExpression& e1, const LinExpression& e2) {
  return e2 > e1;
}


inline Constraint
operator<(const Integer& n, const LinExpression& e) {
  return e > n;
}


inline Constraint
operator<(const LinExpression& e, const Integer& n) {
  return n > e;
}

inline const Constraint&
Constraint::zero_dim_false() {
  static const Constraint zdf(LinExpression::zero() == Integer_one());
  return zdf;
}

inline const Constraint&
Constraint::zero_dim_positivity() {
  static const Constraint zdp(LinExpression::zero() <= Integer_one());
  return zdp;
}

inline const Constraint&
Constraint::epsilon_geq_zero() {
  static const Constraint eps_geq_zero = construct_epsilon_geq_zero();
  return eps_geq_zero;
}

inline const Constraint&
Constraint::epsilon_leq_one() {
  static const Constraint eps_leq_one(LinExpression::zero() < Integer_one());
  return eps_leq_one;
}

inline void
Constraint::swap(Constraint& y) {
  Row::swap(y);
}

}

namespace std {


inline void
swap(Parma_Polyhedra_Library::Constraint& x,
     Parma_Polyhedra_Library::Constraint& y) {
  x.swap(y);
}

}
namespace Parma_Polyhedra_Library {

inline
ConSys::ConSys()
  : Matrix(NECESSARILY_CLOSED) {
}

inline
ConSys::ConSys(const Constraint& c)
  : Matrix(c.topology()) {
  Matrix::insert(c);
}

inline
ConSys::ConSys(const ConSys& cs)
  : Matrix(cs) {
}

inline
ConSys::ConSys(const Topology topol)
  : Matrix(topol) {
}

inline
ConSys::ConSys(const Topology topol,
        const dimension_type n_rows, const dimension_type n_columns)
  : Matrix(topol, n_rows, n_columns) {
}

inline
ConSys::ConSys(ConSys& y, const dimension_type first_stolen)
  : Matrix(y, first_stolen) {
}

inline
ConSys::~ConSys() {
}

inline ConSys&
ConSys::operator=(const ConSys& y) {
  Matrix::operator=(y);
  return *this;
}

inline Constraint&
ConSys::operator[](const dimension_type k) {
  return static_cast<Constraint&> (Matrix::operator[](k));
}

inline const Constraint&
ConSys::operator[](const dimension_type k) const {
  return static_cast<const Constraint&> (Matrix::operator[](k));
}

inline dimension_type
ConSys::max_space_dimension() {



  return Matrix::max_num_columns() - 2;
}

inline dimension_type
ConSys::space_dimension() const {
  return Matrix::space_dimension();
}

inline void
ConSys::clear() {
  Matrix::clear();
}

inline const ConSys&
ConSys::zero_dim_empty() {
  static const ConSys zdf(Constraint::zero_dim_false());
  return zdf;
}

inline
ConSys::const_iterator::
const_iterator(const Parma_Polyhedra_Library::Matrix::const_iterator& iter,
        const ConSys& csys)
  : i(iter), csp(&csys) {
}

inline
ConSys::const_iterator::const_iterator()
  : i(), csp(0) {
}

inline
ConSys::const_iterator::const_iterator(const const_iterator& y)
  : i(y.i), csp(y.csp) {
}

inline
ConSys::const_iterator::~const_iterator() {
}

inline ConSys::const_iterator&
ConSys::const_iterator::operator=(const const_iterator& y) {
  i = y.i;
  csp = y.csp;
  return *this;
}

inline const Constraint&
ConSys::const_iterator::operator*() const {
  return static_cast<const Constraint&>(*i);
}

inline const Constraint*
ConSys::const_iterator::operator->() const {
  return static_cast<const Constraint*>(i.operator->());
}

inline ConSys::const_iterator&
ConSys::const_iterator::operator++() {
  ++i;
  skip_forward();
  return *this;
}

inline ConSys::const_iterator
ConSys::const_iterator::operator++(int) {
  const const_iterator tmp = *this;
  operator++();
  return tmp;
}

inline bool
ConSys::const_iterator::operator==(const const_iterator& y) const {
  return i == y.i;
}

inline bool
ConSys::const_iterator::operator!=(const const_iterator& y) const {
  return i != y.i;
}

inline ConSys::const_iterator
ConSys::begin() const {
  const_iterator i(Matrix::begin(), *this);
  i.skip_forward();
  return i;
}

inline ConSys::const_iterator
ConSys::end() const {
  const const_iterator i(Matrix::end(), *this);
  return i;
}

inline void
ConSys::swap(ConSys& y) {
  Matrix::swap(y);
}

}


namespace std {


inline void
swap(Parma_Polyhedra_Library::ConSys& x,
     Parma_Polyhedra_Library::ConSys& y) {
  x.swap(y);
}

}
namespace Parma_Polyhedra_Library {

namespace IO_Operators {



std::ostream& operator<<(std::ostream& s, const Generator& g);

}

}


namespace std {



void swap(Parma_Polyhedra_Library::Generator& x,
   Parma_Polyhedra_Library::Generator& y);

}
class Parma_Polyhedra_Library::Generator : private Row {
public:






  static Generator line(const LinExpression& e);







  static Generator ray(const LinExpression& e);
  static Generator point(const LinExpression& e = LinExpression::zero(),
    const Integer& d = Integer_one());
  static Generator
  closure_point(const LinExpression& e = LinExpression::zero(),
  const Integer& d = Integer_one());


  Generator(const Generator& g);


  ~Generator();


  Generator& operator=(const Generator& g);


  dimension_type space_dimension() const;


  enum Type {

    LINE,

    RAY,

    POINT,

    CLOSURE_POINT
  };


  Type type() const;


  bool is_line() const;


  bool is_ray() const;


  bool is_point() const;


  bool is_closure_point() const;







  const Integer& coefficient(Variable v) const;






  const Integer& divisor() const;


  static const Generator& zero_dim_point();




  static const Generator& zero_dim_closure_point();


  bool OK() const;

private:



  Generator(LinExpression& e);


  void swap(Generator& y);




  void
  throw_dimension_incompatible(const char* method,
          const char* name_var,
          Variable v) const;




  void
  throw_invalid_argument(const char* method, const char* reason) const;

  friend class Parma_Polyhedra_Library::GenSys;
  friend class Parma_Polyhedra_Library::GenSys::const_iterator;
  friend class Parma_Polyhedra_Library::Polyhedron;
  friend const Integer&
  Parma_Polyhedra_Library::operator*(const Constraint& c,
         const Generator& g);
  friend const Integer&
  Parma_Polyhedra_Library::reduced_scalar_product(const Constraint& c,
        const Generator& g);
  friend
  Parma_Polyhedra_Library::LinExpression::LinExpression(const Generator& g);



  friend class Parma_Polyhedra_Library::ConSys;

  friend void std::swap(Parma_Polyhedra_Library::Generator& x,
   Parma_Polyhedra_Library::Generator& y);

  friend std::ostream&
  Parma_Polyhedra_Library::IO_Operators::operator<<(std::ostream& s,
          const Generator& g);


  Generator(const Generator& g, dimension_type dimension);


  bool is_ray_or_point() const;


  void set_is_line();


  void set_is_ray_or_point();
  bool is_matching_closure_point(const Generator& p) const;


  Generator();
};


namespace Parma_Polyhedra_Library {



Generator line(const LinExpression& e);



Generator ray(const LinExpression& e);





Generator point(const LinExpression& e = LinExpression::zero(),
  const Integer& d = Integer_one());





Generator closure_point(const LinExpression& e = LinExpression::zero(),
   const Integer& d = Integer_one());

}
namespace Parma_Polyhedra_Library {

inline
Generator::Generator(LinExpression& e) {
  Row::swap(e);
}

inline
Generator::Generator(const Generator& g)
  : Row(g) {
}

inline
Generator::Generator(const Generator& g, dimension_type dimension)
  : Row(g, dimension, dimension) {
}

inline
Generator::~Generator() {
}

inline Generator&
Generator::operator=(const Generator& g) {
  Row::operator=(g);
  return *this;
}

inline dimension_type
Generator::space_dimension() const {
  return Row::space_dimension();
}

inline bool
Generator::is_line() const {
  return is_line_or_equality();
}

inline bool
Generator::is_ray_or_point() const {
  return is_ray_or_point_or_inequality();
}

inline bool
Generator::is_ray() const {
  return is_ray_or_point() && ((*this)[0] == 0);
}

inline Generator::Type
Generator::type() const {
  if (is_line())
    return LINE;
  const Generator& g = *this;
  if (g[0] == 0)
    return RAY;
  if (is_necessarily_closed())
    return POINT;
  else

    return (g[size() - 1] == 0) ? CLOSURE_POINT : POINT;
}

inline bool
Generator::is_point() const {
  return type() == POINT;
}

inline bool
Generator::is_closure_point() const {
  return type() == CLOSURE_POINT;
}

inline void
Generator::set_is_line() {
  set_is_line_or_equality();
}

inline void
Generator::set_is_ray_or_point() {
  set_is_ray_or_point_or_inequality();
}

inline const Integer&
Generator::coefficient(const Variable v) const {
  const dimension_type v_id = v.id();
  if (v_id >= space_dimension())
    throw_dimension_incompatible("coefficient(v)", "v", v);
  return Row::coefficient(v_id);
}

inline const Integer&
Generator::divisor() const {
  const Integer& d = Row::inhomogeneous_term();
  if (!is_ray_or_point() || d == 0)
    throw_invalid_argument("divisor()",
      "*this is neither a point nor a closure point");
  return d;
}


inline const Generator&
Generator::zero_dim_point() {
  static const Generator zdp = point();
  return zdp;
}

inline const Generator&
Generator::zero_dim_closure_point() {
  static const Generator zdcp = closure_point();
  return zdcp;
}


inline Generator
line(const LinExpression& e) {
  return Generator::line(e);
}


inline Generator
ray(const LinExpression& e) {
  return Generator::ray(e);
}


inline Generator
point(const LinExpression& e, const Integer& d) {
  return Generator::point(e, d);
}


inline Generator
closure_point(const LinExpression& e, const Integer& d) {
  return Generator::closure_point(e, d);
}

inline void
Generator::swap(Generator& y) {
  Row::swap(y);
}

}

namespace std {


inline void
swap(Parma_Polyhedra_Library::Generator& x,
     Parma_Polyhedra_Library::Generator& y) {
  x.swap(y);
}

}
namespace Parma_Polyhedra_Library {

inline
GenSys::GenSys()
  : Matrix(NECESSARILY_CLOSED) {
}

inline
GenSys::GenSys(const Generator& g)
  : Matrix(g.topology()) {
  Matrix::insert(g);
}

inline
GenSys::GenSys(const GenSys& gs)
  : Matrix(gs) {
}

inline
GenSys::GenSys(const Topology topol)
  : Matrix(topol) {
}

inline
GenSys::GenSys(const Topology topol,
        const dimension_type n_rows, const dimension_type n_columns)
  : Matrix(topol, n_rows, n_columns) {
}

inline
GenSys::GenSys(GenSys& y, const dimension_type first_stolen)
  : Matrix(y, first_stolen) {
}

inline
GenSys::~GenSys() {
}

inline GenSys&
GenSys::operator=(const GenSys& y) {
  Matrix::operator=(y);
  return *this;
}

inline dimension_type
GenSys::max_space_dimension() {

  return Matrix::max_num_columns() - 1;
}

inline dimension_type
GenSys::space_dimension() const {
  return Matrix::space_dimension();
}

inline void
GenSys::clear() {
  Matrix::clear();
}

inline Generator&
GenSys::operator[](const dimension_type k) {
  return static_cast<Generator&>(Matrix::operator[](k));
}

inline const Generator&
GenSys::operator[](const dimension_type k) const {
  return static_cast<const Generator&>(Matrix::operator[](k));
}

inline
GenSys::const_iterator::const_iterator(const Matrix::const_iterator& iter,
           const GenSys& gsys)
  : i(iter), gsp(&gsys) {
}

inline
GenSys::const_iterator::const_iterator()
  : i(), gsp(0) {
}

inline
GenSys::const_iterator::const_iterator(const const_iterator& y)
  : i(y.i), gsp(y.gsp) {
}

inline
GenSys::const_iterator::~const_iterator() {
}

inline
GenSys::const_iterator&
GenSys::const_iterator::operator=(const const_iterator& y) {
  i = y.i;
  gsp = y.gsp;
  return *this;
}

inline const Generator&
GenSys::const_iterator::operator*() const {
  return static_cast<const Generator&>(*i);
}

inline const Generator*
GenSys::const_iterator::operator->() const {
  return static_cast<const Generator*>(i.operator->());
}

inline GenSys::const_iterator&
GenSys::const_iterator::operator++() {
  ++i;
  if (!gsp->is_necessarily_closed())
    skip_forward();
  return *this;
}

inline GenSys::const_iterator
GenSys::const_iterator::operator++(int) {
  const const_iterator tmp = *this;
  operator++();
  return tmp;
}

inline bool
GenSys::const_iterator::operator==(const const_iterator& y) const {
  return i == y.i;
}

inline bool
GenSys::const_iterator::operator!=(const const_iterator& y) const {
  return i != y.i;
}

inline GenSys::const_iterator
GenSys::begin() const {
  const_iterator i(Matrix::begin(), *this);
  if (!is_necessarily_closed())
    i.skip_forward();
  return i;
}

inline GenSys::const_iterator
GenSys::end() const {
  const const_iterator i(Matrix::end(), *this);
  return i;
}

inline const GenSys&
GenSys::zero_dim_univ() {
  static const GenSys zdu(Generator::zero_dim_point());
  return zdu;
}

inline void
GenSys::swap(GenSys& y) {
  Matrix::swap(y);
}

}


namespace std {


inline void
swap(Parma_Polyhedra_Library::GenSys& x,
     Parma_Polyhedra_Library::GenSys& y) {
  x.swap(y);
}

}
namespace Parma_Polyhedra_Library {





bool operator==(const Poly_Con_Relation& x, const Poly_Con_Relation& y);



bool operator!=(const Poly_Con_Relation& x, const Poly_Con_Relation& y);



Poly_Con_Relation operator&&(const Poly_Con_Relation& x,
        const Poly_Con_Relation& y);





Poly_Con_Relation operator-(const Poly_Con_Relation& x,
       const Poly_Con_Relation& y);

namespace IO_Operators {



std::ostream& operator<<(std::ostream& s, const Poly_Con_Relation& r);

}

}







class Parma_Polyhedra_Library::Poly_Con_Relation {
private:

  typedef unsigned int flags_t;



  static const flags_t NOTHING = 0U;
  static const flags_t IS_DISJOINT = 1U << 0;
  static const flags_t STRICTLY_INTERSECTS = 1U << 1;
  static const flags_t IS_INCLUDED = 1U << 2;
  static const flags_t SATURATES = 1U << 3;



  static const flags_t EVERYTHING
  = IS_DISJOINT
  | STRICTLY_INTERSECTS
  | IS_INCLUDED
  | SATURATES;


  flags_t flags;


  static bool implies(flags_t x, flags_t y);


  Poly_Con_Relation(flags_t mask);


  void ascii_dump(std::ostream& s) const;

  friend bool
  Parma_Polyhedra_Library::operator==(const Poly_Con_Relation& x,
          const Poly_Con_Relation& y);

  friend bool
  Parma_Polyhedra_Library::operator!=(const Poly_Con_Relation& x,
          const Poly_Con_Relation& y);

  friend Poly_Con_Relation
  Parma_Polyhedra_Library::operator&&(const Poly_Con_Relation& x,
          const Poly_Con_Relation& y);

  friend Poly_Con_Relation
  Parma_Polyhedra_Library::operator-(const Poly_Con_Relation& x,
         const Poly_Con_Relation& y);

  friend std::ostream&
  Parma_Polyhedra_Library::
  IO_Operators::operator<<(std::ostream& s, const Poly_Con_Relation& r);

public:





  flags_t get_flags() const;

public:

  static Poly_Con_Relation nothing();




  static Poly_Con_Relation is_disjoint();




  static Poly_Con_Relation strictly_intersects();




  static Poly_Con_Relation is_included();




  static Poly_Con_Relation saturates();


  bool implies(const Poly_Con_Relation& y) const;


  bool OK() const;
};
namespace Parma_Polyhedra_Library {

inline
Poly_Con_Relation::Poly_Con_Relation(flags_t mask)
  : flags(mask) {
}

inline Poly_Con_Relation::flags_t
Poly_Con_Relation::get_flags() const {
  return flags;
}

inline Poly_Con_Relation
Poly_Con_Relation::nothing() {
  return Poly_Con_Relation(NOTHING);
}

inline Poly_Con_Relation
Poly_Con_Relation::is_disjoint() {
  return Poly_Con_Relation(IS_DISJOINT);
}

inline Poly_Con_Relation
Poly_Con_Relation::strictly_intersects() {
  return Poly_Con_Relation(STRICTLY_INTERSECTS);
}

inline Poly_Con_Relation
Poly_Con_Relation::is_included() {
  return Poly_Con_Relation(IS_INCLUDED);
}

inline Poly_Con_Relation
Poly_Con_Relation::saturates() {
  return Poly_Con_Relation(SATURATES);
}

inline bool
Poly_Con_Relation::implies(flags_t x, flags_t y) {
  return (x & y) == y;
}

inline bool
Poly_Con_Relation::implies(const Poly_Con_Relation& y) const {
  return implies(flags, y.flags);
}


inline bool
operator==(const Poly_Con_Relation& x, const Poly_Con_Relation& y) {
  return x.flags == y.flags;
}


inline bool
operator!=(const Poly_Con_Relation& x, const Poly_Con_Relation& y) {
  return x.flags != y.flags;
}


inline Poly_Con_Relation
operator&&(const Poly_Con_Relation& x, const Poly_Con_Relation& y) {
  return Poly_Con_Relation(x.flags | y.flags);
}


inline Poly_Con_Relation
operator-(const Poly_Con_Relation& x, const Poly_Con_Relation& y) {
  return Poly_Con_Relation(x.flags & ~y.flags);
}

}
namespace Parma_Polyhedra_Library {





bool operator==(const Poly_Gen_Relation& x, const Poly_Gen_Relation& y);



bool operator!=(const Poly_Gen_Relation& x, const Poly_Gen_Relation& y);



Poly_Gen_Relation operator&&(const Poly_Gen_Relation& x,
        const Poly_Gen_Relation& y);





Poly_Gen_Relation operator-(const Poly_Gen_Relation& x,
       const Poly_Gen_Relation& y);

namespace IO_Operators {



std::ostream& operator<<(std::ostream& s, const Poly_Gen_Relation& r);

}

}







class Parma_Polyhedra_Library::Poly_Gen_Relation {
private:

  typedef unsigned int flags_t;



  static const flags_t NOTHING = 0U;
  static const flags_t SUBSUMES = 1U << 0;



  static const flags_t EVERYTHING
  = SUBSUMES;


  flags_t flags;


  static bool implies(flags_t x, flags_t y);


  Poly_Gen_Relation(flags_t mask);


  void ascii_dump(std::ostream& s) const;

  friend bool
  Parma_Polyhedra_Library::operator==(const Poly_Gen_Relation& x,
          const Poly_Gen_Relation& y);

  friend bool
  Parma_Polyhedra_Library::operator!=(const Poly_Gen_Relation& x,
          const Poly_Gen_Relation& y);

  friend Poly_Gen_Relation
  Parma_Polyhedra_Library::operator&&(const Poly_Gen_Relation& x,
          const Poly_Gen_Relation& y);

  friend Poly_Gen_Relation
  Parma_Polyhedra_Library::operator-(const Poly_Gen_Relation& x,
         const Poly_Gen_Relation& y);

  friend std::ostream&
  Parma_Polyhedra_Library::
  IO_Operators::operator<<(std::ostream& s, const Poly_Gen_Relation& r);

public:





  flags_t get_flags() const;

public:

  static Poly_Gen_Relation nothing();


  static Poly_Gen_Relation subsumes();


  bool implies(const Poly_Gen_Relation& y) const;


  bool OK() const;
};
namespace Parma_Polyhedra_Library {

inline
Poly_Gen_Relation::Poly_Gen_Relation(flags_t mask)
  : flags(mask) {
}

inline Poly_Gen_Relation::flags_t
Poly_Gen_Relation::get_flags() const {
  return flags;
}

inline Poly_Gen_Relation
Poly_Gen_Relation::nothing() {
  return Poly_Gen_Relation(NOTHING);
}

inline Poly_Gen_Relation
Poly_Gen_Relation::subsumes() {
  return Poly_Gen_Relation(SUBSUMES);
}

inline bool
Poly_Gen_Relation::implies(flags_t x, flags_t y) {
  return (x & y) == y;
}

inline bool
Poly_Gen_Relation::implies(const Poly_Gen_Relation& y) const {
  return implies(flags, y.flags);
}


inline bool
operator==(const Poly_Gen_Relation& x, const Poly_Gen_Relation& y) {
  return x.flags == y.flags;
}


inline bool
operator!=(const Poly_Gen_Relation& x, const Poly_Gen_Relation& y) {
  return x.flags != y.flags;
}


inline Poly_Gen_Relation
operator&&(const Poly_Gen_Relation& x, const Poly_Gen_Relation& y) {
  return Poly_Gen_Relation(x.flags | y.flags);
}


inline Poly_Gen_Relation
operator-(const Poly_Gen_Relation& x, const Poly_Gen_Relation& y) {
  return Poly_Gen_Relation(x.flags & ~y.flags);
}

}
class Parma_Polyhedra_Library::BHRZ03_Certificate {
public:

  BHRZ03_Certificate();


  BHRZ03_Certificate(const Polyhedron& ph);


  BHRZ03_Certificate(const BHRZ03_Certificate& y);


  ~BHRZ03_Certificate();
  int compare(const BHRZ03_Certificate& y) const;


  int compare(const Polyhedron& ph) const;






  bool is_stabilizing(const Polyhedron& ph) const;






  struct Compare {

    bool operator()(const BHRZ03_Certificate& x,
      const BHRZ03_Certificate& y) const;
  };




  bool OK() const;

private:

  dimension_type poly_dim;

  dimension_type lin_space_dim;

  dimension_type num_constraints;



  dimension_type num_points;




  std::vector<dimension_type> num_rays_null_coord;
};
namespace Parma_Polyhedra_Library {

inline
BHRZ03_Certificate::BHRZ03_Certificate()
  : poly_dim(0), lin_space_dim(0), num_constraints(0), num_points(1),
    num_rays_null_coord() {

  assert(OK());
}

inline
BHRZ03_Certificate::BHRZ03_Certificate(const BHRZ03_Certificate& y)
  : poly_dim(y.poly_dim), lin_space_dim(y.lin_space_dim),
    num_constraints(y.num_constraints), num_points(y.num_points),
    num_rays_null_coord(y.num_rays_null_coord) {
}

inline
BHRZ03_Certificate::~BHRZ03_Certificate() {
}

inline bool
BHRZ03_Certificate::is_stabilizing(const Polyhedron& ph) const {
  return compare(ph) == 1;
}

inline bool
BHRZ03_Certificate::Compare::operator()(const BHRZ03_Certificate& x,
     const BHRZ03_Certificate& y) const {


  return (x.compare(y) == 1);
}

}
class Parma_Polyhedra_Library::H79_Certificate {
public:

  H79_Certificate();


  H79_Certificate(const Polyhedron& ph);


  H79_Certificate(const H79_Certificate& y);


  ~H79_Certificate();
  int compare(const H79_Certificate& y) const;


  int compare(const Polyhedron& ph) const;






  struct Compare {

    bool operator()(const H79_Certificate& x,
      const H79_Certificate& y) const;
  };

private:

  dimension_type poly_dim;

  dimension_type num_constraints;
};
namespace Parma_Polyhedra_Library {

inline
H79_Certificate::H79_Certificate()
  : poly_dim(0), num_constraints(0) {

}

inline
H79_Certificate::H79_Certificate(const H79_Certificate& y)
  : poly_dim(y.poly_dim), num_constraints(y.num_constraints) {
}

inline
H79_Certificate::~H79_Certificate() {
}

inline bool
H79_Certificate::Compare::operator()(const H79_Certificate& x,
         const H79_Certificate& y) const {


  return (x.compare(y) == 1);
}

}
namespace Parma_Polyhedra_Library {

namespace IO_Operators {
std::ostream&
operator<<(std::ostream& s, const Polyhedron& ph);

}
bool operator==(const Polyhedron& x, const Polyhedron& y);
bool operator!=(const Polyhedron& x, const Polyhedron& y);

}
class Parma_Polyhedra_Library::Polyhedron {
public:

  static dimension_type max_space_dimension();


  enum Degenerate_Kind {

    UNIVERSE,

    EMPTY
  };

protected:
  Polyhedron(Topology topol,
      dimension_type num_dimensions,
      Degenerate_Kind kind);


  Polyhedron(const Polyhedron& y);
  Polyhedron(Topology topol, const ConSys& cs);
  Polyhedron(Topology topol, ConSys& cs);
  Polyhedron(Topology topol, const GenSys& gs);
  Polyhedron(Topology topol, GenSys& gs);
  template <typename Box>
  Polyhedron(Topology topol, const Box& box);




  Polyhedron& operator=(const Polyhedron& y);

public:




  dimension_type space_dimension() const;


  const ConSys& constraints() const;


  const ConSys& minimized_constraints() const;


  const GenSys& generators() const;


  const GenSys& minimized_generators() const;
  Poly_Con_Relation relation_with(const Constraint& c) const;
  Poly_Gen_Relation relation_with(const Generator& g) const;




  bool is_empty() const;




  bool is_universe() const;




  bool is_topologically_closed() const;







  bool is_disjoint_from(const Polyhedron& y) const;




  bool is_bounded() const;
  bool bounds_from_above(const LinExpression& expr) const;
  bool bounds_from_below(const LinExpression& expr) const;
  bool maximize(const LinExpression& expr,
  Integer& sup_n, Integer& sup_d, bool& maximum) const;
  bool maximize(const LinExpression& expr,
  Integer& sup_n, Integer& sup_d, bool& maximum,
  const Generator** const pppoint) const;
  bool minimize(const LinExpression& expr,
  Integer& inf_n, Integer& inf_d, bool& minimum) const;
  bool minimize(const LinExpression& expr,
  Integer& inf_n, Integer& inf_d, bool& minimum,
  const Generator** const pppoint) const;







  bool contains(const Polyhedron& y) const;







  bool strictly_contains(const Polyhedron& y) const;
  template <typename Box>
  void shrink_bounding_box(Box& box, Complexity_Class complexity = ANY) const;
  bool OK(bool check_not_empty = false) const;
  void add_constraint(const Constraint& c);
  bool add_constraint_and_minimize(const Constraint& c);
  void add_generator(const Generator& g);
  bool add_generator_and_minimize(const Generator& g);
  void add_constraints(const ConSys& cs);
  void add_recycled_constraints(ConSys& cs);
  bool add_constraints_and_minimize(const ConSys& cs);
  bool add_recycled_constraints_and_minimize(ConSys& cs);
  void add_generators(const GenSys& gs);
  void add_recycled_generators(GenSys& gs);
  bool add_generators_and_minimize(const GenSys& gs);
  bool add_recycled_generators_and_minimize(GenSys& gs);
  void intersection_assign(const Polyhedron& y);
  bool intersection_assign_and_minimize(const Polyhedron& y);
  void poly_hull_assign(const Polyhedron& y);
  bool poly_hull_assign_and_minimize(const Polyhedron& y);
  void poly_difference_assign(const Polyhedron& y);
  void affine_image(Variable var,
      const LinExpression& expr,
      const Integer& denominator = Integer_one());
  void affine_preimage(Variable var,
         const LinExpression& expr,
         const Integer& denominator = Integer_one());
  void generalized_affine_image(Variable var,
    const Relation_Symbol relsym,
    const LinExpression& expr,
    const Integer& denominator = Integer_one());
  void generalized_affine_image(const LinExpression& lhs,
    const Relation_Symbol relsym,
    const LinExpression& rhs);
  void time_elapse_assign(const Polyhedron& y);


  void topological_closure_assign();
  void BHRZ03_widening_assign(const Polyhedron& y, unsigned* tp = 0);
  void limited_BHRZ03_extrapolation_assign(const Polyhedron& y,
        const ConSys& cs,
        unsigned* tp = 0);
  void bounded_BHRZ03_extrapolation_assign(const Polyhedron& y,
        const ConSys& cs,
        unsigned* tp = 0);
  void H79_widening_assign(const Polyhedron& y, unsigned* tp = 0);
  void limited_H79_extrapolation_assign(const Polyhedron& y,
     const ConSys& cs,
     unsigned* tp = 0);
  void bounded_H79_extrapolation_assign(const Polyhedron& y,
     const ConSys& cs,
     unsigned* tp = 0);
  void add_dimensions_and_embed(dimension_type m);
  void add_dimensions_and_project(dimension_type m);
  void concatenate_assign(const Polyhedron& y);
  void remove_dimensions(const Variables_Set& to_be_removed);
  void remove_higher_dimensions(dimension_type new_dimension);
  template <typename PartialFunction>
  void map_dimensions(const PartialFunction& pfunc);
  void expand_dimension(Variable var, dimension_type m);
  void fold_dimensions(const Variables_Set& to_be_folded, Variable var);



  friend bool Parma_Polyhedra_Library::operator==(const Polyhedron& x,
        const Polyhedron& y);





  ~Polyhedron();
  void swap(Polyhedron& y);






  void ascii_dump(std::ostream& s) const;







  bool ascii_load(std::istream& s);



private:

  ConSys con_sys;


  GenSys gen_sys;


  SatMatrix sat_c;


  SatMatrix sat_g;
class Status {
public:

  Status();



  bool test_zero_dim_univ() const;
  void reset_zero_dim_univ();
  void set_zero_dim_univ();

  bool test_empty() const;
  void reset_empty();
  void set_empty();

  bool test_c_up_to_date() const;
  void reset_c_up_to_date();
  void set_c_up_to_date();

  bool test_g_up_to_date() const;
  void reset_g_up_to_date();
  void set_g_up_to_date();

  bool test_c_minimized() const;
  void reset_c_minimized();
  void set_c_minimized();

  bool test_g_minimized() const;
  void reset_g_minimized();
  void set_g_minimized();

  bool test_sat_c_up_to_date() const;
  void reset_sat_c_up_to_date();
  void set_sat_c_up_to_date();

  bool test_sat_g_up_to_date() const;
  void reset_sat_g_up_to_date();
  void set_sat_g_up_to_date();

  bool test_c_pending() const;
  void reset_c_pending();
  void set_c_pending();

  bool test_g_pending() const;
  void reset_g_pending();
  void set_g_pending();



  bool OK() const;




  void ascii_dump(std::ostream& s) const;





  bool ascii_load(std::istream& s);

private:

  typedef unsigned int flags_t;



  static const flags_t ZERO_DIM_UNIV = 0U;
  static const flags_t EMPTY = 1U << 0;
  static const flags_t C_UP_TO_DATE = 1U << 1;
  static const flags_t G_UP_TO_DATE = 1U << 2;
  static const flags_t C_MINIMIZED = 1U << 3;
  static const flags_t G_MINIMIZED = 1U << 4;
  static const flags_t SAT_C_UP_TO_DATE = 1U << 5;
  static const flags_t SAT_G_UP_TO_DATE = 1U << 6;
  static const flags_t CS_PENDING = 1U << 7;
  static const flags_t GS_PENDING = 1U << 8;



  flags_t flags;


  Status(flags_t mask);


  bool test_all(flags_t mask) const;


  bool test_any(flags_t mask) const;


  void set(flags_t mask);


  void reset(flags_t mask);
};




  Status status;


  dimension_type space_dim;


  Topology topology() const;




  bool is_necessarily_closed() const;
  bool marked_empty() const;


  bool constraints_are_up_to_date() const;


  bool generators_are_up_to_date() const;






  bool constraints_are_minimized() const;






  bool generators_are_minimized() const;


  bool has_pending_constraints() const;


  bool has_pending_generators() const;




  bool has_something_pending() const;




  bool can_have_something_pending() const;




  bool sat_c_is_up_to_date() const;




  bool sat_g_is_up_to_date() const;
  void set_zero_dim_univ();




  void set_empty();


  void set_constraints_up_to_date();


  void set_generators_up_to_date();


  void set_constraints_minimized();


  void set_generators_minimized();


  void set_constraints_pending();


  void set_generators_pending();


  void set_sat_c_up_to_date();


  void set_sat_g_up_to_date();







  void clear_empty();






  void clear_constraints_up_to_date();






  void clear_generators_up_to_date();


  void clear_constraints_minimized();


  void clear_generators_minimized();


  void clear_pending_constraints();


  void clear_pending_generators();


  void clear_sat_c_up_to_date();


  void clear_sat_g_up_to_date();
  bool process_pending() const;
  bool process_pending_constraints() const;





  void process_pending_generators() const;
  void remove_pending_to_obtain_constraints() const;
  bool remove_pending_to_obtain_generators() const;
  void update_constraints() const;
  bool update_generators() const;
  void update_sat_c() const;
  void update_sat_g() const;
  void obtain_sorted_constraints() const;
  void obtain_sorted_generators() const;
  void obtain_sorted_constraints_with_sat_c() const;
  void obtain_sorted_generators_with_sat_g() const;
  bool minimize() const;







  bool strongly_minimize_constraints() const;







  bool strongly_minimize_generators() const;



  enum Three_Valued_Boolean {
    TVB_TRUE,
    TVB_FALSE,
    TVB_DONT_KNOW
  };


  Three_Valued_Boolean quick_equivalence_test(const Polyhedron& y) const;



  bool is_included_in(const Polyhedron& y) const;
  bool bounds(const LinExpression& expr, bool from_above) const;
  bool max_min(const LinExpression& expr,
        const bool maximize,
        Integer& ext_n, Integer& ext_d, bool& included,
        const Generator** const pppoint = 0) const;







  void select_CH78_constraints(const Polyhedron& y,
          ConSys& cs_selected) const;





  void select_H79_constraints(const Polyhedron& y,
         ConSys& cs_selected,
         ConSys& cs_not_selected) const;

  friend class Parma_Polyhedra_Library::BHRZ03_Certificate;
  friend class Parma_Polyhedra_Library::H79_Certificate;

  bool BHRZ03_combining_constraints(const Polyhedron& y,
        const BHRZ03_Certificate& y_cert,
         const Polyhedron& H79,
        const ConSys& x_minus_H79_con_sys);

  bool BHRZ03_evolving_points(const Polyhedron& y,
         const BHRZ03_Certificate& y_cert,
         const Polyhedron& H79);

  bool BHRZ03_evolving_rays(const Polyhedron& y,
       const BHRZ03_Certificate& y_cert,
       const Polyhedron& H79);




  static void add_low_level_constraints(ConSys& cs);
  static void add_dimensions(Matrix& mat1,
                             Matrix& mat2,
                             SatMatrix& sat1,
                             SatMatrix& sat2,
        dimension_type add_dim);
  static bool is_BHRZ03_stabilizing(const Polyhedron& x, const Polyhedron& y);






  static bool minimize(bool con_to_gen,
         Matrix& source, Matrix& dest, SatMatrix& sat);





  static bool add_and_minimize(bool con_to_gen,
          Matrix& source1, Matrix& dest, SatMatrix& sat,
          const Matrix& source2);





  static bool add_and_minimize(bool con_to_gen,
          Matrix& source, Matrix& dest, SatMatrix& sat);




  static dimension_type conversion(Matrix& source,
       dimension_type start,
       Matrix& dest,
       SatMatrix& sat,
       dimension_type num_lines_or_equalities);





  static int simplify(Matrix& mat, SatMatrix& sat);





  void throw_runtime_error(const char* method) const;
  void throw_invalid_argument(const char* method, const char* reason) const;

  void throw_topology_incompatible(const char* method,
       const char* ph_name,
       const Polyhedron& ph) const;
  void throw_topology_incompatible(const char* method,
       const char* c_name,
       const Constraint& c) const;
  void throw_topology_incompatible(const char* method,
       const char* g_name,
       const Generator& g) const;
  void throw_topology_incompatible(const char* method,
       const char* cs_name,
       const ConSys& cs) const;
  void throw_topology_incompatible(const char* method,
       const char* gs_name,
       const GenSys& gs) const;

  void throw_dimension_incompatible(const char* method,
        const char* ph_name,
        const Polyhedron& ph) const;
  void throw_dimension_incompatible(const char* method,
        const char* var_name,
        const Variable var) const;
  void throw_dimension_incompatible(const char* method,
        const char* row_name,
        const Row& row) const;
  void throw_dimension_incompatible(const char* method,
        const char* sys_name,
        const Matrix& sys) const;
  void throw_dimension_incompatible(const char* method,
        dimension_type required_space_dim) const;

  void throw_invalid_generator(const char* method,
          const char* g_name) const;
  void throw_invalid_generators(const char* method,
    const char* gs_name) const;



};


namespace std {



void swap(Parma_Polyhedra_Library::Polyhedron& x,
   Parma_Polyhedra_Library::Polyhedron& y);

}
namespace Parma_Polyhedra_Library {

inline
Polyhedron::Status::Status(flags_t mask)
  : flags(mask) {
}

inline
Polyhedron::Status::Status()
  : flags(ZERO_DIM_UNIV) {
}

inline bool
Polyhedron::Status::test_all(flags_t mask) const {
  return (flags & mask) == mask;
}

inline bool
Polyhedron::Status::test_any(flags_t mask) const {
  return flags & mask;
}

inline void
Polyhedron::Status::set(flags_t mask) {
  flags |= mask;
}

inline void
Polyhedron::Status::reset(flags_t mask) {
  flags &= ~mask;
}

inline bool
Polyhedron::Status::test_zero_dim_univ() const {
  return flags == ZERO_DIM_UNIV;
}

inline void
Polyhedron::Status::reset_zero_dim_univ() {

  if (flags == ZERO_DIM_UNIV)

    flags = EMPTY;
}

inline void
Polyhedron::Status::set_zero_dim_univ() {

  flags = ZERO_DIM_UNIV;
}

inline bool
Polyhedron::Status::test_empty() const {
  return test_any(EMPTY);
}

inline void
Polyhedron::Status::reset_empty() {
  reset(EMPTY);
}

inline void
Polyhedron::Status::set_empty() {
  flags = EMPTY;
}

inline bool
Polyhedron::Status::test_c_up_to_date() const {
  return test_any(C_UP_TO_DATE);
}

inline void
Polyhedron::Status::reset_c_up_to_date() {
  reset(C_UP_TO_DATE);
}

inline void
Polyhedron::Status::set_c_up_to_date() {
  set(C_UP_TO_DATE);
}

inline bool
Polyhedron::Status::test_g_up_to_date() const {
  return test_any(G_UP_TO_DATE);
}

inline void
Polyhedron::Status::reset_g_up_to_date() {
  reset(G_UP_TO_DATE);
}

inline void
Polyhedron::Status::set_g_up_to_date() {
  set(G_UP_TO_DATE);
}

inline bool
Polyhedron::Status::test_c_minimized() const {
  return test_any(C_MINIMIZED);
}

inline void
Polyhedron::Status::reset_c_minimized() {
  reset(C_MINIMIZED);
}

inline void
Polyhedron::Status::set_c_minimized() {
  set(C_MINIMIZED);
}

inline bool
Polyhedron::Status::test_g_minimized() const {
  return test_any(G_MINIMIZED);
}

inline void
Polyhedron::Status::reset_g_minimized() {
  reset(G_MINIMIZED);
}

inline void
Polyhedron::Status::set_g_minimized() {
  set(G_MINIMIZED);
}


inline bool
Polyhedron::Status::test_c_pending() const {
  return test_any(CS_PENDING);
}

inline void
Polyhedron::Status::reset_c_pending() {
  reset(CS_PENDING);
}

inline void
Polyhedron::Status::set_c_pending() {
  set(CS_PENDING);
}

inline bool
Polyhedron::Status::test_g_pending() const {
  return test_any(GS_PENDING);
}

inline void
Polyhedron::Status::reset_g_pending() {
  reset(GS_PENDING);
}

inline void
Polyhedron::Status::set_g_pending() {
  set(GS_PENDING);
}


inline bool
Polyhedron::Status::test_sat_c_up_to_date() const {
  return test_any(SAT_C_UP_TO_DATE);
}

inline void
Polyhedron::Status::reset_sat_c_up_to_date() {
  reset(SAT_C_UP_TO_DATE);
}

inline void
Polyhedron::Status::set_sat_c_up_to_date() {
  set(SAT_C_UP_TO_DATE);
}

inline bool
Polyhedron::Status::test_sat_g_up_to_date() const {
  return test_any(SAT_G_UP_TO_DATE);
}

inline void
Polyhedron::Status::reset_sat_g_up_to_date() {
  reset(SAT_G_UP_TO_DATE);
}

inline void
Polyhedron::Status::set_sat_g_up_to_date() {
  set(SAT_G_UP_TO_DATE);
}

}
namespace Parma_Polyhedra_Library {

inline Topology
Polyhedron::topology() const {


  return con_sys.topology();
}

inline bool
Polyhedron::is_necessarily_closed() const {


  return con_sys.is_necessarily_closed();
}

inline dimension_type
Polyhedron::space_dimension() const {
  return space_dim;
}

inline
Polyhedron::~Polyhedron() {
}

inline void
Polyhedron::swap(Polyhedron& y) {
  if (topology() != y.topology())
    throw_topology_incompatible("swap(y)", "y", y);
  std::swap(con_sys, y.con_sys);
  std::swap(gen_sys, y.gen_sys);
  std::swap(sat_c, y.sat_c);
  std::swap(sat_g, y.sat_g);
  std::swap(status, y.status);
  std::swap(space_dim, y.space_dim);
}

}


inline void
std::swap(Parma_Polyhedra_Library::Polyhedron& x,
   Parma_Polyhedra_Library::Polyhedron& y) {
  x.swap(y);
}

namespace Parma_Polyhedra_Library {

inline bool
Polyhedron::marked_empty() const {
  return status.test_empty();
}

inline bool
Polyhedron::constraints_are_up_to_date() const {
  return status.test_c_up_to_date();
}

inline bool
Polyhedron::generators_are_up_to_date() const {
  return status.test_g_up_to_date();
}

inline bool
Polyhedron::constraints_are_minimized() const {
  return status.test_c_minimized();
}

inline bool
Polyhedron::generators_are_minimized() const {
  return status.test_g_minimized();
}

inline bool
Polyhedron::sat_c_is_up_to_date() const {
  return status.test_sat_c_up_to_date();
}

inline bool
Polyhedron::sat_g_is_up_to_date() const {
  return status.test_sat_g_up_to_date();
}

inline bool
Polyhedron::has_pending_constraints() const {
  return status.test_c_pending();
}

inline bool
Polyhedron::has_pending_generators() const {
  return status.test_g_pending();
}

inline bool
Polyhedron::has_something_pending() const {
  return status.test_c_pending() || status.test_g_pending();
}

inline bool
Polyhedron::can_have_something_pending() const {
  return constraints_are_minimized()
    && generators_are_minimized()
    && (sat_c_is_up_to_date() || sat_g_is_up_to_date());
}

inline void
Polyhedron::set_constraints_up_to_date() {
  status.set_c_up_to_date();
}

inline void
Polyhedron::set_generators_up_to_date() {
  status.set_g_up_to_date();
}

inline void
Polyhedron::set_constraints_minimized() {
  set_constraints_up_to_date();
  status.set_c_minimized();
}

inline void
Polyhedron::set_generators_minimized() {
  set_generators_up_to_date();
  status.set_g_minimized();
}

inline void
Polyhedron::set_constraints_pending() {
  status.set_c_pending();
}

inline void
Polyhedron::set_generators_pending() {
  status.set_g_pending();
}

inline void
Polyhedron::set_sat_c_up_to_date() {
  status.set_sat_c_up_to_date();
}

inline void
Polyhedron::set_sat_g_up_to_date() {
  status.set_sat_g_up_to_date();
}

inline void
Polyhedron::clear_empty() {
  status.reset_empty();
}

inline void
Polyhedron::clear_constraints_minimized() {
  status.reset_c_minimized();
}

inline void
Polyhedron::clear_generators_minimized() {
  status.reset_g_minimized();
}

inline void
Polyhedron::clear_pending_constraints() {
  status.reset_c_pending();
}

inline void
Polyhedron::clear_pending_generators() {
  status.reset_g_pending();
}

inline void
Polyhedron::clear_sat_c_up_to_date() {
  status.reset_sat_c_up_to_date();

}

inline void
Polyhedron::clear_sat_g_up_to_date() {
  status.reset_sat_g_up_to_date();

}

inline void
Polyhedron::clear_constraints_up_to_date() {
  clear_pending_constraints();
  clear_constraints_minimized();
  clear_sat_c_up_to_date();
  clear_sat_g_up_to_date();
  status.reset_c_up_to_date();

}

inline void
Polyhedron::clear_generators_up_to_date() {
  clear_pending_generators();
  clear_generators_minimized();
  clear_sat_c_up_to_date();
  clear_sat_g_up_to_date();
  status.reset_g_up_to_date();

}

inline bool
Polyhedron::process_pending() const {
  assert(space_dim > 0 && !marked_empty());
  assert(has_something_pending());

  Polyhedron& x = const_cast<Polyhedron&>(*this);

  if (x.has_pending_constraints())
    return x.process_pending_constraints();

  assert(x.has_pending_generators());
  x.process_pending_generators();
  return true;
}

inline void
Polyhedron::add_low_level_constraints(ConSys& cs) {
  if (cs.is_necessarily_closed())

    cs.insert(Constraint::zero_dim_positivity());
  else {

    cs.insert(Constraint::epsilon_leq_one());
    cs.insert(Constraint::epsilon_geq_zero());
  }
}

inline bool
Polyhedron::is_empty() const {
  if (marked_empty())
    return true;



  if (generators_are_up_to_date() && !has_pending_constraints())
    return false;
  return !minimize();
}

inline bool
Polyhedron::bounds_from_above(const LinExpression& expr) const {
  return bounds(expr, true);
}

inline bool
Polyhedron::bounds_from_below(const LinExpression& expr) const {
  return bounds(expr, false);
}

inline bool
Polyhedron::maximize(const LinExpression& expr,
       Integer& sup_n, Integer& sup_d, bool& maximum) const {
  return max_min(expr, true, sup_n, sup_d, maximum);
}

inline bool
Polyhedron::maximize(const LinExpression& expr,
       Integer& sup_n, Integer& sup_d, bool& maximum,
       const Generator** const pppoint) const {
  return max_min(expr, true, sup_n, sup_d, maximum, pppoint);
}

inline bool
Polyhedron::minimize(const LinExpression& expr,
       Integer& inf_n, Integer& inf_d, bool& minimum) const {
  return max_min(expr, false, inf_n, inf_d, minimum);
}

inline bool
Polyhedron::minimize(const LinExpression& expr,
       Integer& inf_n, Integer& inf_d, bool& minimum,
       const Generator** const pppoint) const {
  return max_min(expr, false, inf_n, inf_d, minimum, pppoint);
}


inline bool
operator!=(const Polyhedron& x, const Polyhedron& y) {
  return !(x == y);
}

inline bool
Polyhedron::strictly_contains(const Polyhedron& y) const {
  const Polyhedron& x = *this;
  return x.contains(y) && !y.contains(x);
}

template <typename Box>
Polyhedron::Polyhedron(Topology topol, const Box& box)
  : con_sys(topol),
    gen_sys(topol),
    sat_c(),
    sat_g() {

  space_dim = box.space_dimension();


  if (box.is_empty()) {
    set_empty();
    return;
  }


  if (space_dim == 0) {
    set_zero_dim_univ();
    return;
  }




  con_sys.insert(Variable(space_dim - 1) >= 0);

  for (dimension_type k = space_dim; k-- > 0; ) {

    bool l_closed = false;
    Integer l_n, l_d;
    bool l_bounded = box.get_lower_bound(k, l_closed, l_n, l_d);
    if (l_bounded && topol == NECESSARILY_CLOSED && !l_closed)
      throw_invalid_argument("C_Polyhedron(const Box& box)",
        "box has an open lower bound");

    bool u_closed = false;
    Integer u_n, u_d;
    bool u_bounded = box.get_upper_bound(k, u_closed, u_n, u_d);
    if (u_bounded && topol == NECESSARILY_CLOSED && !u_closed)
      throw_invalid_argument("C_Polyhedron(const Box& box)",
        "box has an open upper bound");


    if (l_bounded && u_bounded
 && l_closed && u_closed
 && l_n == u_n && l_d == u_d) {

      con_sys.insert(l_d * Variable(k) == l_n);
    }
    else {

      if (l_bounded) {
       if (l_closed)

  con_sys.insert(l_d * Variable(k) >= l_n);
       else

  con_sys.insert(l_d * Variable(k) > l_n);
      }

      if (u_bounded) {
       if (u_closed)

  con_sys.insert(u_d * Variable(k) <= u_n);
       else

  con_sys.insert(u_d * Variable(k) < u_n);
      }
    }
  }


  add_low_level_constraints(con_sys);

  dimension_type n_rows = con_sys.num_rows() - 1;
  con_sys[0].swap(con_sys[n_rows]);
  con_sys.set_sorted(false);

  con_sys.set_index_first_pending_row(n_rows);
  con_sys.erase_to_end(n_rows);


  set_constraints_up_to_date();
  assert(OK());
}

template <typename Box>
void
Polyhedron::shrink_bounding_box(Box& box, Complexity_Class complexity) const {
  bool polynomial = (complexity != ANY);
  if ((polynomial && !has_something_pending()
       && constraints_are_minimized()) || !polynomial) {


    if (is_universe())
      return;
  }
  if (polynomial) {
    if (marked_empty() ||
 (generators_are_up_to_date() && gen_sys.num_rows() == 0)) {
      box.set_empty();
      return;
    }
    if (constraints_are_up_to_date()) {
      for (ConSys::const_iterator i = con_sys.begin(); i != con_sys.end(); ++i)
 if ((*i).is_trivial_false()){
   box.set_empty();
   return;
 }
    }
  }
  else


    if (is_empty()) {
      box.set_empty();
      return;
    }

  if (space_dim == 0)
    return;



  std::vector<LBoundary>
    lower_bound(space_dim, LBoundary(ERational('+'), LBoundary::OPEN));

  std::vector<UBoundary>
    upper_bound(space_dim, UBoundary(ERational('-'), UBoundary::OPEN));

  if (!polynomial && has_something_pending())
    process_pending();

  if (polynomial &&
       (!generators_are_up_to_date() || has_pending_constraints())) {

    assert(constraints_are_up_to_date());




    ConSys cs(con_sys);
    if (cs.num_pending_rows() > 0) {
      cs.unset_pending_rows();
      cs.sort_rows();
    }
    else if (!cs.is_sorted())
      cs.sort_rows();

    if (has_pending_constraints() || !constraints_are_minimized())
      cs.back_substitute(cs.gauss());

    const ConSys::const_iterator cs_begin = cs.begin();
    const ConSys::const_iterator cs_end = cs.end();

    for (ConSys::const_iterator i = cs_begin; i != cs_end; ++i) {
      dimension_type varid = space_dim;
      const Constraint& c = *i;


      for (dimension_type j = space_dim; j-- > 0; ) {
 if (c.is_trivial_false()) {
   box.set_empty();
   return;
 }


 if (c.coefficient(Variable(j)) != 0)
   if (varid != space_dim) {
     varid = space_dim;
     break;
   }
   else
     varid = j;
      }
      if (varid != space_dim) {
 const Integer& d = c.coefficient(Variable(varid));
 const Integer& n = c.inhomogeneous_term();





 ERational r(-n, d);
 Constraint::Type c_type = c.type();
 switch (c_type) {
 case Constraint::EQUALITY:
   lower_bound[varid] = LBoundary(r, LBoundary::CLOSED);
   upper_bound[varid] = UBoundary(r, UBoundary::CLOSED);
   break;
 case Constraint::NONSTRICT_INEQUALITY:
 case Constraint::STRICT_INEQUALITY:
   if (d > 0)


     lower_bound[varid]
       = LBoundary(r, (c_type == Constraint::NONSTRICT_INEQUALITY
         ? LBoundary::CLOSED
         : LBoundary::OPEN));
   else {



     assert(d < 0);
     upper_bound[varid]
       = UBoundary(r, (c_type == Constraint::NONSTRICT_INEQUALITY
         ? UBoundary::CLOSED
         : UBoundary::OPEN));
   }
   break;
 }
      }
    }
  }
  else {






    const GenSys& gs = gen_sys;

    const GenSys::const_iterator gs_begin = gs.begin();
    const GenSys::const_iterator gs_end = gs.end();



    for (GenSys::const_iterator i = gs_begin; i != gs_end; ++i) {
      const Generator& g = *i;
      Generator::Type g_type = g.type();
      switch (g_type) {
      case Generator::LINE:


 for (dimension_type j = space_dim; j-- > 0; )
   if (g.coefficient(Variable(j)) != 0) {
     lower_bound[j] = LBoundary(ERational('-'), LBoundary::OPEN);
     upper_bound[j] = UBoundary(ERational('+'), UBoundary::OPEN);
   }
 break;
      case Generator::RAY:


 for (dimension_type j = space_dim; j-- > 0; ) {
   int sign = sgn(g.coefficient(Variable(j)));
   if (sign < 0)
     lower_bound[j] = LBoundary(ERational('-'), LBoundary::OPEN);
   else if (sign > 0)
     upper_bound[j] = UBoundary(ERational('+'), UBoundary::OPEN);
 }
 break;
      case Generator::POINT:
      case Generator::CLOSURE_POINT:
 {
   const Integer& d = g.divisor();
   for (dimension_type j = space_dim; j-- > 0; ) {
     const Integer& n = g.coefficient(Variable(j));
     ERational r(n, d);
     LBoundary lb(r,(g_type == Generator::CLOSURE_POINT
       ? LBoundary::OPEN
       : LBoundary::CLOSED));
     if (lb < lower_bound[j])
       lower_bound[j] = lb;
     UBoundary ub(r, (g_type == Generator::CLOSURE_POINT
        ? UBoundary::OPEN
        : UBoundary::CLOSED));
     if (ub > upper_bound[j])
       upper_bound[j] = ub;
   }
 }
 break;
      }
    }
  }


  for (dimension_type j = space_dim; j-- > 0; ) {

    const LBoundary& lb = lower_bound[j];
    const ERational& lr = lb.bound();
    if (lr.direction_of_infinity() == 0)
      box.raise_lower_bound(j, lb.is_closed(),
       lr.numerator(), lr.denominator());


    const UBoundary& ub = upper_bound[j];
    const ERational& ur = ub.bound();
    if (ur.direction_of_infinity() == 0)
      box.lower_upper_bound(j, ub.is_closed(),
       ur.numerator(), ur.denominator());
  }
}

template <typename PartialFunction>
void
Polyhedron::map_dimensions(const PartialFunction& pfunc) {
  if (space_dim == 0)
    return;

  if (pfunc.has_empty_codomain()) {

    if (marked_empty()
 || (has_pending_constraints()
     && !remove_pending_to_obtain_generators())
 || (!generators_are_up_to_date() && !update_generators())) {

      space_dim = 0;
      con_sys.clear();
    }
    else

      set_zero_dim_univ();

    assert(OK());
    return;
  }

  dimension_type new_space_dimension = pfunc.max_in_codomain() + 1;

  if (new_space_dimension == space_dim) {
    std::vector<dimension_type> cycles;
    cycles.reserve(space_dim + space_dim/2);


    std::deque<bool> visited(space_dim);

    for (dimension_type i = space_dim; i-- > 0; ) {
      if (!visited[i]) {
 dimension_type j = i;
 do {
   visited[j] = true;
   dimension_type k;
   (void) pfunc.maps(j, k);
   if (k == j)

     goto skip;

   cycles.push_back(j+1);

   j = k;
 } while (!visited[j]);

 cycles.push_back(0);
      skip:
 ;
      }
    }


    if (cycles.empty())
      return;




    if (constraints_are_up_to_date())
      con_sys.permute_columns(cycles);

    if (generators_are_up_to_date())
      gen_sys.permute_columns(cycles);

    assert(OK());
    return;
  }





  const GenSys& old_gensys = generators();

  if (old_gensys.num_rows() == 0) {

    Polyhedron new_polyhedron(topology(), new_space_dimension, EMPTY);
    std::swap(*this, new_polyhedron);
    assert(OK());
    return;
  }


  std::vector<dimension_type> pfunc_maps(space_dim, not_a_dimension());
  for (dimension_type j = space_dim; j-- > 0; ) {
    dimension_type pfunc_j;
    if (pfunc.maps(j, pfunc_j))
      pfunc_maps[j] = pfunc_j;
  }

  GenSys new_gensys;
  for (GenSys::const_iterator i = old_gensys.begin(),
  old_gensys_end = old_gensys.end(); i != old_gensys_end; ++i) {
    const Generator& old_g = *i;
    LinExpression e(0 * Variable(new_space_dimension-1));
    bool all_zeroes = true;
    for (dimension_type j = space_dim; j-- > 0; ) {
      if (old_g.coefficient(Variable(j)) != 0
   && pfunc_maps[j] != not_a_dimension()) {
 e += Variable(pfunc_maps[j]) * old_g.coefficient(Variable(j));
 all_zeroes = false;
      }
    }
    switch (old_g.type()) {
    case Generator::LINE:
      if (!all_zeroes)
 new_gensys.insert(line(e));
      break;
    case Generator::RAY:
      if (!all_zeroes)
 new_gensys.insert(ray(e));
      break;
    case Generator::POINT:

      new_gensys.insert(point(e, old_g.divisor()));
      break;
    case Generator::CLOSURE_POINT:

      new_gensys.insert(closure_point(e, old_g.divisor()));
      break;
    }
  }
  Polyhedron new_polyhedron(topology(), new_gensys);
  std::swap(*this, new_polyhedron);
  assert(OK(true));
}

}
class Parma_Polyhedra_Library::C_Polyhedron : public Polyhedron {
public:

  static dimension_type max_space_dimension();
  explicit C_Polyhedron(dimension_type num_dimensions = 0,
   Degenerate_Kind kind = UNIVERSE);
  C_Polyhedron(const ConSys& cs);
  C_Polyhedron(ConSys& cs);
  C_Polyhedron(const GenSys& gs);
  C_Polyhedron(GenSys& gs);






  explicit C_Polyhedron(const NNC_Polyhedron& y);
  template <typename Box>
  C_Polyhedron(const Box& box, From_Bounding_Box dummy);


  C_Polyhedron(const C_Polyhedron& y);




  C_Polyhedron& operator=(const C_Polyhedron& y);


  ~C_Polyhedron();
};
namespace Parma_Polyhedra_Library {

inline dimension_type
C_Polyhedron::max_space_dimension() {
  using std::min;


  return min(std::numeric_limits<dimension_type>::max() - 1,
      min(ConSys::max_space_dimension(),
   GenSys::max_space_dimension()
   )
      );
}

inline
C_Polyhedron::C_Polyhedron(dimension_type num_dimensions,
      Degenerate_Kind kind)
  : Polyhedron(NECESSARILY_CLOSED,
        num_dimensions > max_space_dimension() ? 0 : num_dimensions,
        kind) {
  if (num_dimensions > max_space_dimension())
    throw std::length_error("PPL::C_Polyhedron::C_Polyhedron(n, k):\n"
       "n exceeds the maximum allowed space dimension.");
}

inline
C_Polyhedron::C_Polyhedron(const ConSys& cs)
  : Polyhedron(NECESSARILY_CLOSED, cs) {
}

inline
C_Polyhedron::C_Polyhedron(ConSys& cs)
  : Polyhedron(NECESSARILY_CLOSED, cs) {
}

inline
C_Polyhedron::C_Polyhedron(const GenSys& gs)
  : Polyhedron(NECESSARILY_CLOSED, gs) {
}

inline
C_Polyhedron::C_Polyhedron(GenSys& gs)
  : Polyhedron(NECESSARILY_CLOSED, gs) {
}

template <typename Box>
C_Polyhedron::C_Polyhedron(const Box& box, From_Bounding_Box)
  : Polyhedron(NECESSARILY_CLOSED, box) {
}

inline
C_Polyhedron::C_Polyhedron(const C_Polyhedron& y)
  : Polyhedron(y) {
}

inline C_Polyhedron&
C_Polyhedron::operator=(const C_Polyhedron& y) {
  Polyhedron::operator=(y);
  return *this;
}

inline
C_Polyhedron::~C_Polyhedron() {
}

}
class Parma_Polyhedra_Library::NNC_Polyhedron : public Polyhedron {
public:

  static dimension_type max_space_dimension();
  explicit NNC_Polyhedron(dimension_type num_dimensions = 0,
     Degenerate_Kind kind = UNIVERSE);
  NNC_Polyhedron(const ConSys& cs);
  NNC_Polyhedron(ConSys& cs);
  NNC_Polyhedron(const GenSys& gs);
  NNC_Polyhedron(GenSys& gs);


  explicit NNC_Polyhedron(const C_Polyhedron& y);
  template <typename Box>
  NNC_Polyhedron(const Box& box, From_Bounding_Box dummy);


  NNC_Polyhedron(const NNC_Polyhedron& y);




  NNC_Polyhedron& operator=(const NNC_Polyhedron& y);


  ~NNC_Polyhedron();

};
namespace Parma_Polyhedra_Library {

inline dimension_type
NNC_Polyhedron::max_space_dimension() {

  return C_Polyhedron::max_space_dimension() - 1;
}

inline
NNC_Polyhedron::NNC_Polyhedron(dimension_type num_dimensions,
          Degenerate_Kind kind)
  : Polyhedron(NOT_NECESSARILY_CLOSED,
        num_dimensions > max_space_dimension() ? 0 : num_dimensions,
        kind) {
  if (num_dimensions > max_space_dimension())
    throw std::length_error("PPL::NNC_Polyhedron::NNC_Polyhedron(n, k):\n"
       "n exceeds the maximum allowed space dimension.");
}

inline
NNC_Polyhedron::NNC_Polyhedron(const ConSys& cs)
  : Polyhedron(NOT_NECESSARILY_CLOSED, cs) {
}

inline
NNC_Polyhedron::NNC_Polyhedron(ConSys& cs)
  : Polyhedron(NOT_NECESSARILY_CLOSED, cs) {
}

inline
NNC_Polyhedron::NNC_Polyhedron(const GenSys& gs)
  : Polyhedron(NOT_NECESSARILY_CLOSED, gs) {
}

inline
NNC_Polyhedron::NNC_Polyhedron(GenSys& gs)
  : Polyhedron(NOT_NECESSARILY_CLOSED, gs) {
}

template <typename Box>
NNC_Polyhedron::NNC_Polyhedron(const Box& box, From_Bounding_Box)
  : Polyhedron(NOT_NECESSARILY_CLOSED, box) {
}

inline
NNC_Polyhedron::NNC_Polyhedron(const NNC_Polyhedron& y)
  : Polyhedron(y) {
}

inline NNC_Polyhedron&
NNC_Polyhedron::operator=(const NNC_Polyhedron& y) {
  Polyhedron::operator=(y);
  return *this;
}

inline
NNC_Polyhedron::~NNC_Polyhedron() {
}

}
namespace Parma_Polyhedra_Library {

template <typename PH>
class Widening_Function;

template <typename PH>
class Limited_Widening_Function;

}
template <typename PH>
class Parma_Polyhedra_Library::Widening_Function {
public:

  typedef void (PH::* Widening_Method)(const PH&, unsigned*);


  Widening_Function(Widening_Method wm);






  void operator()(PH& x, const PH& y, unsigned* tp = 0) const;

private:

  Widening_Method w_method;
};





template <typename PH>
class Parma_Polyhedra_Library::Limited_Widening_Function {
public:

  typedef void (PH::* Limited_Widening_Method)(const PH&,
            const ConSys&,
            unsigned*);
  Limited_Widening_Function(Limited_Widening_Method lwm, const ConSys& cs);







  void operator()(PH& x, const PH& y, unsigned* tp = 0) const;

private:

  Limited_Widening_Method lw_method;

  const ConSys& limiting_cs;
};

namespace Parma_Polyhedra_Library {
template <typename PH>
Widening_Function<PH>
widen_fun(void (PH::* wm)(const PH&, unsigned*));
template <typename PH>
Limited_Widening_Function<PH>
widen_fun(void (PH::* lwm)(const PH&, const ConSys&, unsigned*),
   const ConSys& cs);

}
namespace Parma_Polyhedra_Library {

template <typename PH>
Widening_Function<PH>::Widening_Function(Widening_Method wm)
  : w_method(wm) {
}

template <typename PH>
inline void
Widening_Function<PH>::
operator()(PH& x, const PH& y, unsigned* tp) const {
  (x.*w_method)(y, tp);
}

template <typename PH>
Limited_Widening_Function<PH>::
Limited_Widening_Function(Limited_Widening_Method lwm, const ConSys& cs)
  : lw_method(lwm), limiting_cs(cs) {
}

template <typename PH>
inline void
Limited_Widening_Function<PH>::
operator()(PH& x, const PH& y, unsigned* tp) const {
  (x.*lw_method)(y, limiting_cs, tp);
}

template <typename PH>
inline Widening_Function<PH>
widen_fun(void (PH::* wm)(const PH&, unsigned*)) {
  return Widening_Function<PH>(wm);
}

template <typename PH>
inline Limited_Widening_Function<PH>
widen_fun(void (PH::* lwm)(const PH&, const ConSys&, unsigned*),
   const ConSys& cs) {
  return Limited_Widening_Function<PH>(lwm, cs);
}

}
namespace Parma_Polyhedra_Library {

template <typename PH>
class Determinate;

}
namespace Parma_Polyhedra_Library {
template <typename PH>
bool operator==(const Determinate<PH>& x, const Determinate<PH>& y);
template <typename PH>
bool operator!=(const Determinate<PH>& x, const Determinate<PH>& y);

namespace IO_Operators {



template <typename PH>
std::ostream&
operator<<(std::ostream&, const Determinate<PH>&);

}

}


template <typename PH>
class Parma_Polyhedra_Library::Determinate {
public:
  explicit
  Determinate(dimension_type num_dimensions = 0, bool universe = true);




  Determinate(const PH& p);




  Determinate(const ConSys& cs);


  Determinate(const Determinate& y);


  ~Determinate();







  dimension_type space_dimension() const;


  const ConSys& constraints() const;


  const ConSys& minimized_constraints() const;


  const PH& element() const;


  PH& element();






  void mutate();




  bool is_top() const;




  bool is_bottom() const;




  bool definitely_entails(const Determinate& y) const;




  bool is_definitely_equivalent_to(const Determinate& y) const;


  bool OK() const;




  Determinate& operator <<= (dimension_type n);
  Determinate& hide_assign(dimension_type n);

  friend bool
  operator==<PH>(const Determinate<PH>& x, const Determinate<PH>& y);
  friend bool
  operator!=<PH>(const Determinate<PH>& x, const Determinate<PH>& y);
  void upper_bound_assign(const Determinate& y);


  void meet_assign(const Determinate& y);
  void add_constraint(const Constraint& c);
  void add_constraints(ConSys& cs);







  Determinate& operator=(const Determinate& y);


  void swap(Determinate& y);




  void add_dimensions_and_embed(dimension_type m);




  void add_dimensions_and_project(dimension_type m);




  void concatenate_assign(const Determinate& y);
  void remove_dimensions(const Variables_Set& to_be_removed);
  void remove_higher_dimensions(dimension_type new_dimension);







  template <typename PartialFunction>
  void map_dimensions(const PartialFunction& pfunc);



private:






  class Rep {
  private:





    mutable unsigned long references;


    Rep& operator=(const Rep& y);


    Rep(const Rep& y);


    Rep();

  public:

    PH ph;




    Rep(dimension_type num_dimensions, Polyhedron::Degenerate_Kind kind);


    Rep(const PH& p);


    Rep(const ConSys& cs);


    ~Rep();


    void new_reference() const;




    bool del_reference() const;


    bool is_shared() const;
  };




  Rep* prep;
};


namespace std {



template <typename PH>
void swap(Parma_Polyhedra_Library::Determinate<PH>& x,
   Parma_Polyhedra_Library::Determinate<PH>& y);

}
namespace Parma_Polyhedra_Library {

template <typename PH>
Determinate<PH>::Rep::Rep(dimension_type num_dimensions,
     Polyhedron::Degenerate_Kind kind)
  : references(0), ph(num_dimensions, kind) {
}

template <typename PH>
Determinate<PH>::Rep::Rep(const PH& p)
  : references(0), ph(p) {
}

template <typename PH>
Determinate<PH>::Rep::Rep(const ConSys& cs)
  : references(0), ph(cs) {
}

template <typename PH>
Determinate<PH>::Rep::~Rep() {
  assert(references == 0);
}

template <typename PH>
void
Determinate<PH>::Rep::new_reference() const {
  ++references;
}

template <typename PH>
bool Determinate<PH>::Rep::del_reference() const {
  return --references == 0;
}

template <typename PH>
bool
Determinate<PH>::Rep::is_shared() const {
  return references > 1;
}

template <typename PH>
Determinate<PH>::Determinate(dimension_type num_dimensions, bool universe)
  : prep(new Rep(num_dimensions,
   universe ? Polyhedron::UNIVERSE : Polyhedron::EMPTY)) {
  prep->new_reference();
}

template <typename PH>
Determinate<PH>::Determinate(const PH& ph)
  : prep(new Rep(ph)) {
  prep->new_reference();
}

template <typename PH>
Determinate<PH>::Determinate(const ConSys& cs)
  : prep(new Rep(cs)) {
  prep->new_reference();
}

template <typename PH>
Determinate<PH>::Determinate(const Determinate& y)
  : prep(y.prep) {
  prep->new_reference();
}

template <typename PH>
Determinate<PH>::~Determinate() {
  if (prep->del_reference())
    delete prep;
}

template <typename PH>
Determinate<PH>&
Determinate<PH>::operator=(const Determinate& y) {
  y.prep->new_reference();
  if (prep->del_reference())
    delete prep;
  prep = y.prep;
  return *this;
}

template <typename PH>
inline void
Determinate<PH>::swap(Determinate& y) {
  std::swap(prep, y.prep);
}

template <typename PH>
void
Determinate<PH>::mutate() {
  if (prep->is_shared()) {
    Rep* new_prep = new Rep(prep->ph);
    (void) prep->del_reference();
    new_prep->new_reference();
    prep = new_prep;
  }
}

template <typename PH>
const PH&
Determinate<PH>::element() const {
  return prep->ph;
}

template <typename PH>
PH&
Determinate<PH>::element() {
  mutate();
  return prep->ph;
}

template <typename PH>
void
Determinate<PH>::upper_bound_assign(const Determinate& y) {
  mutate();
  prep->ph.poly_hull_assign(y.prep->ph);
}

template <typename PH>
void
Determinate<PH>::meet_assign(const Determinate& y) {
  mutate();
  prep->ph.intersection_assign(y.prep->ph);
}

template <typename PH>
void
Determinate<PH>::concatenate_assign(const Determinate& y) {
  mutate();
  prep->ph.concatenate_assign(y.prep->ph);
}

template <typename PH>
bool
Determinate<PH>::definitely_entails(const Determinate& y) const {
  return prep == y.prep || y.prep->ph.contains(prep->ph);
}

template <typename PH>
bool
Determinate<PH>::is_definitely_equivalent_to(const Determinate& y) const {
  return prep == y.prep || prep->ph == y.prep->ph;
}

template <typename PH>
inline bool
Determinate<PH>::is_top() const {
  return prep->ph.is_universe();
}

template <typename PH>
inline bool
Determinate<PH>::is_bottom() const {
  return prep->ph.is_empty();
}

template <typename PH>
bool
Determinate<PH>::OK() const {
  return prep->ph.OK();
}

namespace IO_Operators {


template <typename PH>
std::ostream&
operator<<(std::ostream& s, const Determinate<PH>& x) {
  if (x.is_top())
    s << "true";
  else if (x.is_bottom())
    s << "false";
  else {
    const ConSys& cs = x.constraints();
    ConSys::const_iterator i = cs.begin();
    ConSys::const_iterator cs_end = cs.end();
    s << "{ ";
    while (i != cs_end) {
      s << *i++;
      if (i != cs_end)
 s << ", ";
    }
    s << " }";
  }
  return s;
}

}


template <typename PH>
bool
operator==(const Determinate<PH>& x, const Determinate<PH>& y) {
  return x.prep == y.prep || x.prep->ph == y.prep->ph;
}


template <typename PH>
bool
operator!=(const Determinate<PH>& x, const Determinate<PH>& y) {
  return x.prep != y.prep && x.prep->ph != y.prep->ph;
}

template <typename PH>
dimension_type
Determinate<PH>::space_dimension() const {
  return prep->ph.space_dimension();
}

template <typename PH>
const ConSys&
Determinate<PH>::constraints() const {
  return prep->ph.constraints();
}

template <typename PH>
const ConSys&
Determinate<PH>::minimized_constraints() const {
  return prep->ph.minimized_constraints();
}

template <typename PH>
void
Determinate<PH>::add_constraint(const Constraint& c) {
  mutate();
  prep->ph.add_constraint(c);
}

template <typename PH>
void
Determinate<PH>::add_constraints(ConSys& cs) {
  mutate();
  prep->ph.add_constraints(cs);
}

template <typename PH>
void
Determinate<PH>::add_dimensions_and_embed(dimension_type m) {
  mutate();
  prep->ph.add_dimensions_and_embed(m);
}

template <typename PH>
void
Determinate<PH>::add_dimensions_and_project(dimension_type m) {
  mutate();
  prep->ph.add_dimensions_and_project(m);
}

template <typename PH>
void
Determinate<PH>::remove_dimensions(const Variables_Set& to_be_removed) {
  mutate();
  prep->ph.remove_dimensions(to_be_removed);
}

template <typename PH>
void
Determinate<PH>::remove_higher_dimensions(dimension_type new_dimension) {
  mutate();
  prep->ph.remove_higher_dimensions(new_dimension);
}

template <typename PH>
template <typename PartialFunction>
void
Determinate<PH>::map_dimensions(const PartialFunction& pfunc) {
  mutate();
  prep->ph.map_dimensions(pfunc);
}

}


namespace std {


template <typename PH>
inline void
swap(Parma_Polyhedra_Library::Determinate<PH>& x,
     Parma_Polyhedra_Library::Determinate<PH>& y) {
  x.swap(y);
}

}
namespace Parma_Polyhedra_Library {

template <typename CS>
class PowerSet;

}
namespace Parma_Polyhedra_Library {



template <typename CS>
bool
operator==(const PowerSet<CS>& x, const PowerSet<CS>& y);



template <typename CS>
bool
operator!=(const PowerSet<CS>& x, const PowerSet<CS>& y);

namespace IO_Operators {



template <typename CS>
std::ostream&
operator<<(std::ostream& s, const PowerSet<CS>& x);

}

}







template <typename CS>
class Parma_Polyhedra_Library::PowerSet {
public:

  PowerSet();


  PowerSet(const PowerSet& y);


  PowerSet& operator=(const PowerSet& y);


  void swap(PowerSet& y);


  void add_disjunct(const CS& d);


  void upper_bound_assign(const PowerSet& y);


  void meet_assign(const PowerSet& y);






  bool definitely_entails(const PowerSet& y) const;





  bool is_top() const;





  bool is_bottom() const;



  bool OK(bool disallow_bottom = false) const;





  void collapse();

protected:





  typedef std::list<CS> Sequence;


  Sequence sequence;


  mutable bool reduced;


  void omega_reduce() const;




  bool is_omega_reduced() const;

  void collapse(unsigned max_disjuncts);

private:



  bool check_omega_reduced() const;

public:
  typedef typename Sequence::size_type size_type;

  size_type size() const;
  bool empty() const;

  typedef typename Sequence::iterator iterator;
  typedef typename Sequence::const_iterator const_iterator;
  typedef typename Sequence::reverse_iterator reverse_iterator;
  typedef typename Sequence::const_reverse_iterator const_reverse_iterator;
  typedef typename Sequence::value_type value_type;

  iterator begin();
  const_iterator begin() const;

  iterator end();
  const_iterator end() const;

  reverse_iterator rbegin();
  const_reverse_iterator rbegin() const;

  reverse_iterator rend();
  const_reverse_iterator rend() const;

  void push_back(const CS& y);
  void pop_back();
  iterator erase(iterator first, iterator last);
  iterator erase(iterator position);

protected:
  static void add_non_bottom_disjunct(Sequence& s,
          const CS& d,
          iterator& first,
          iterator last);




  static void add_non_bottom_disjunct(Sequence& s, const CS& d);

private:
  void collapse(iterator sink);
};


namespace std {



template <typename CS>
void swap(Parma_Polyhedra_Library::PowerSet<CS>& x,
   Parma_Polyhedra_Library::PowerSet<CS>& y);

}
namespace Parma_Polyhedra_Library {

template <typename CS>
typename PowerSet<CS>::iterator
PowerSet<CS>::begin() {
  return sequence.begin();
}

template <typename CS>
typename PowerSet<CS>::const_iterator
PowerSet<CS>::begin() const {
  return sequence.begin();
}

template <typename CS>
typename PowerSet<CS>::iterator
PowerSet<CS>::end() {
  return sequence.end();
}

template <typename CS>
typename PowerSet<CS>::const_iterator
PowerSet<CS>::end() const {
  return sequence.end();
}

template <typename CS>
typename PowerSet<CS>::reverse_iterator
PowerSet<CS>::rbegin() {
  return sequence.rbegin();
}

template <typename CS>
typename PowerSet<CS>::const_reverse_iterator
PowerSet<CS>::rbegin() const {
  return sequence.rbegin();
}

template <typename CS>
typename PowerSet<CS>::reverse_iterator
PowerSet<CS>::rend() {
  return sequence.rend();
}

template <typename CS>
typename PowerSet<CS>::const_reverse_iterator
PowerSet<CS>::rend() const {
  return sequence.rend();
}

template <typename CS>
typename PowerSet<CS>::size_type
PowerSet<CS>::size() const {
  return sequence.size();
}

template <typename CS>
bool
PowerSet<CS>::empty() const {
  return sequence.empty();
}

template <typename CS>
void
PowerSet<CS>::push_back(const CS& y) {
  sequence.push_back(y);
}

template <typename CS>
void
PowerSet<CS>::pop_back() {
  sequence.pop_back();
}

template <typename CS>
typename PowerSet<CS>::iterator
PowerSet<CS>::erase(iterator first, iterator last) {
  return sequence.erase(first, last);
}

template <typename CS>
typename PowerSet<CS>::iterator
PowerSet<CS>::erase(iterator position) {
  return sequence.erase(position);
}

template <typename CS>
PowerSet<CS>::PowerSet(const PowerSet<CS>& y)
  : sequence(y.sequence), reduced(y.reduced) {
}

template <typename CS>
PowerSet<CS>&
PowerSet<CS>::operator=(const PowerSet<CS>& y) {
  sequence = y.sequence;
  reduced = y.reduced;
  return *this;
}

template <typename CS>
inline void
PowerSet<CS>::swap(PowerSet& y) {
  std::swap(sequence, y.sequence);
  std::swap(reduced, y.reduced);
}

template <typename CS>
PowerSet<CS>::PowerSet()
  : sequence(), reduced(true) {
}

template <typename CS>
void
PowerSet<CS>::collapse(const iterator sink) {
  assert(sink != end());


  CS& d = *sink;
  iterator j = sink;
  iterator send = end();
  for (++j; j != send; ++j)
    d.upper_bound_assign(*j);


  j = sink;
  erase(++j, send);


  for (iterator k = begin(), kn = k; k != sink; k = kn) {
    ++kn;
    if (k->definitely_entails(d))
      erase(k);
  }
  assert(OK());
}

template <typename CS>
void
PowerSet<CS>::omega_reduce() const {
  if (reduced)
    return;

  PowerSet& ps = const_cast<PowerSet&>(*this);

  for (iterator xi = ps.begin(), xin = xi; xi != ps.end(); xi = xin) {
    ++xin;
    if (xi->is_bottom())
      ps.erase(xi);
  }

  for (iterator xi = ps.begin(), xin = xi; xi != ps.end(); xi = xin) {
    ++xin;
    const CS& xv = *xi;
    for (iterator yi = ps.begin(), yin = yi; yi != ps.end(); yi = yin) {
      ++yin;
      if (xi == yi)
 continue;
      const CS& yv = *yi;
      if (yv.definitely_entails(xv)) {
 if (yi == xin)
   ++xin;
 ps.erase(yi);
      }
      else if (xv.definitely_entails(yv)) {
 ps.erase(xi);
 break;
      }
    }
    if (abandon_expensive_computations && xin != ps.end()) {

      ps.collapse(xin);
      break;
    }
  }
  reduced = true;
  assert(OK());
}

template <typename CS>
void
PowerSet<CS>::collapse(const unsigned max_disjuncts) {
  assert(max_disjuncts > 0);

  omega_reduce();
  size_type n = size();
  if (n > max_disjuncts) {
    iterator i = begin();

    for (unsigned m = max_disjuncts-1; m-- > 0; )
      ++i;



    collapse(i);
  }
  assert(OK());
  assert(is_omega_reduced());
}

template <typename CS>
bool
PowerSet<CS>::check_omega_reduced() const {
  for (const_iterator sbegin = begin(), send = end(),
  xi = sbegin; xi != send; ++xi) {
    const CS& xv = *xi;
    if (xv.is_bottom())
      return false;
    for (const_iterator yi = sbegin; yi != send; ++yi) {
      if (xi == yi)
 continue;
      const CS& yv = *yi;
      if (xv.definitely_entails(yv) || yv.definitely_entails(xv))
 return false;
    }
  }
  return true;
}

template <typename CS>
bool
PowerSet<CS>::is_omega_reduced() const {
  if (!reduced && check_omega_reduced())
    reduced = true;
  return reduced;
}

template <typename CS>
void
PowerSet<CS>::add_non_bottom_disjunct(Sequence& s,
          const CS& d,
          iterator& first,
          iterator last) {
  for (iterator xi = first, xin = xi; xi != last; xi = xin) {
    ++xin;
    const CS& xv = *xi;
    if (d.definitely_entails(xv))
      return;
    else if (xv.definitely_entails(d)) {
      if (xi == first)
 first = xin;
      s.erase(xi);
    }
  }
  s.push_back(d);
}

template <typename CS>
void
PowerSet<CS>::add_non_bottom_disjunct(Sequence& s, const CS& d) {
  assert(!d.is_bottom());
  iterator s_begin = s.begin();
  iterator s_end = s.end();
  add_non_bottom_disjunct(s, d, s_begin, s_end);
}

template <typename CS>
void
PowerSet<CS>::add_disjunct(const CS& d) {
  if (!d.is_bottom())
    add_non_bottom_disjunct(sequence, d);
}

template <typename CS>
bool
PowerSet<CS>::definitely_entails(const PowerSet<CS>& y) const {
  const PowerSet<CS>& x = *this;
  bool found = true;
  for (const_iterator xi = x.begin(),
  x_end = x.end(); found && xi != x_end; ++xi) {
    found = false;
    for (const_iterator yi = y.begin(),
    y_end = y.end(); !found && yi != y_end; ++yi)
      found = (*xi).definitely_entails(*yi);
  }
  return found;
}


template <typename CS>
inline
bool operator==(const PowerSet<CS>& x, const PowerSet<CS>& y) {
  return (x.size() == y.size() && equal(x.begin(), x.end(), y.begin()));
}


template <typename CS>
inline
bool operator!=(const PowerSet<CS>& x, const PowerSet<CS>& y) {
  return !(x == y);
}

template <typename CS>
inline bool
PowerSet<CS>::is_top() const {

  omega_reduce();
  const_iterator i = begin();
  const_iterator send = end();
  return i != send && i->is_top() && ++i == send;
}

template <typename CS>
inline bool
PowerSet<CS>::is_bottom() const {

  omega_reduce();
  return sequence.empty();
}

template <typename CS>
inline void
PowerSet<CS>::collapse() {
  if (!empty())
    collapse(begin());
}

template <typename CS>
void
PowerSet<CS>::meet_assign(const PowerSet<CS>& y) {
  const PowerSet<CS>& x = *this;
  Sequence new_sequence;
  for (const_iterator xi = x.begin(), x_end = x.end(); xi != x_end; ++xi)
    for (const_iterator yi = y.begin(), y_end = y.end(); yi != y_end; ++yi) {
      CS zi = *xi;
      zi.meet_assign(*yi);
      if (!zi.is_bottom())
 new_sequence.push_back(zi);
    }
  std::swap(sequence, new_sequence);
  omega_reduce();
}

template <typename CS>
void
PowerSet<CS>::upper_bound_assign(const PowerSet<CS>& y) {

  omega_reduce();
  y.omega_reduce();
  iterator sbegin = begin();
  iterator send = end();
  for (const_iterator i = y.begin(), y_end = y.end(); i != y_end; ++i)
    add_non_bottom_disjunct(sequence, *i, sbegin, send);
}

namespace IO_Operators {


template <typename CS>
std::ostream&
operator<<(std::ostream& s, const PowerSet<CS>& x) {
  if (x.is_bottom())
    s << "false";
  else if (x.is_top())
    s << "true";
  else {
    s << "{ ";
    typename PowerSet<CS>::const_iterator i = x.begin();
    typename PowerSet<CS>::const_iterator x_end = x.end();
    while (i != x_end) {
      s << *i++;
      if (i != x_end)
 s << ", ";
    }
    s << " }";
  }
  return s;
}

}

template <typename CS>
bool
PowerSet<CS>::OK(const bool disallow_bottom) const {
  for (const_iterator i = begin(), send = end(); i != send; ++i) {
    if (!i->OK())
      return false;
    if (disallow_bottom && i->is_bottom()) {




      return false;
    }
  }
  if (reduced && !check_omega_reduced()) {




    return false;
  }
  return true;
}

}


namespace std {


template <typename CS>
inline void
swap(Parma_Polyhedra_Library::PowerSet<CS>& x,
     Parma_Polyhedra_Library::PowerSet<CS>& y) {
  x.swap(y);
}

}
namespace Parma_Polyhedra_Library {

template <typename PH>
class Polyhedra_PowerSet;

}
template <typename PH>
class Parma_Polyhedra_Library::Polyhedra_PowerSet
  : public Parma_Polyhedra_Library::PowerSet
<Parma_Polyhedra_Library::Determinate<PH> > {
public:

  static dimension_type max_space_dimension();
  explicit
  Polyhedra_PowerSet(dimension_type num_dimensions = 0,
       Polyhedron::Degenerate_Kind kind = Polyhedron::UNIVERSE);


  Polyhedra_PowerSet(const Polyhedra_PowerSet& y);


  Polyhedra_PowerSet(const ConSys& cs);







  dimension_type space_dimension() const;
  bool geometrically_covers(const Polyhedra_PowerSet& y) const;
  bool geometrically_equals(const Polyhedra_PowerSet& y) const;


  bool OK() const;
  void add_constraint(const Constraint& c);
  bool add_constraint_and_minimize(const Constraint& c);
  void add_constraints(const ConSys& cs);
  bool add_constraints_and_minimize(const ConSys& cs);
  void pairwise_reduce();
  template <typename Widening>
  void BGP99_extrapolation_assign(const Polyhedra_PowerSet& y,
      Widening wf,
      unsigned max_disjuncts);
  template <typename Cert, typename Widening>
  void BHZ03_widening_assign(const Polyhedra_PowerSet& y, Widening wf);




  template <typename Widening>
  void BHZ03_widening_assign(const Polyhedra_PowerSet& y, Widening wf);
  Polyhedra_PowerSet& operator=(const Polyhedra_PowerSet& y);


  void swap(Polyhedra_PowerSet& y);




  void add_dimensions_and_embed(dimension_type m);




  void add_dimensions_and_project(dimension_type m);







  void concatenate_assign(const Polyhedra_PowerSet& y);
  void remove_dimensions(const Variables_Set& to_be_removed);
  void remove_higher_dimensions(dimension_type new_dimension);







  template <typename PartialFunction>
  void map_dimensions(const PartialFunction& pfunc);



private:
  typedef Determinate<PH> CS;
  typedef PowerSet<CS> Base;

public:
  typedef typename Base::Sequence Sequence;

  typedef typename Sequence::size_type size_type;

  typedef typename Base::iterator iterator;
  typedef typename Base::const_iterator const_iterator;
  typedef typename Base::reverse_iterator reverse_iterator;
  typedef typename Base::const_reverse_iterator const_reverse_iterator;
  typedef typename Base::value_type value_type;






  void ascii_dump(std::ostream& s) const;







  bool ascii_load(std::istream& s);

private:

  dimension_type space_dim;




  template <typename Widening>
  void BGP99_heuristics_assign(const Polyhedra_PowerSet& y, Widening wf);


  template <typename Cert>
  void collect_certificates(std::map<Cert, size_type,
                       typename Cert::Compare>& cert_ms) const;




  template <typename Cert>
  bool is_cert_multiset_stabilizing(const std::map<Cert, size_type,
                                                   typename Cert::Compare>&
        y_cert_ms) const;
};


namespace std {



template <typename PH>
void swap(Parma_Polyhedra_Library::Polyhedra_PowerSet<PH>& x,
   Parma_Polyhedra_Library::Polyhedra_PowerSet<PH>& y);

}
namespace Parma_Polyhedra_Library {

template <typename PH>
dimension_type
Polyhedra_PowerSet<PH>::max_space_dimension() {
  return PH::max_space_dimension();
}

template <typename PH>
Polyhedra_PowerSet<PH>::Polyhedra_PowerSet(dimension_type num_dimensions,
        Polyhedron::Degenerate_Kind kind)
  : space_dim(num_dimensions) {
  if (kind == Polyhedron::UNIVERSE)
    push_back(Determinate<PH>(num_dimensions, true));
}

template <typename PH>
Polyhedra_PowerSet<PH>::Polyhedra_PowerSet(const Polyhedra_PowerSet& y)
  : Base(y), space_dim(y.space_dim) {
}

template <typename PH>
Polyhedra_PowerSet<PH>::Polyhedra_PowerSet(const ConSys& cs)
  : space_dim(cs.space_dimension()) {
  push_back(Determinate<PH>(cs));
}

template <typename PH>
Polyhedra_PowerSet<PH>&
Polyhedra_PowerSet<PH>::operator=(const Polyhedra_PowerSet& y) {
  Base::operator=(y);
  space_dim = y.space_dim;
  return *this;
}

template <typename PH>
inline void
Polyhedra_PowerSet<PH>::swap(Polyhedra_PowerSet& y) {
  Base::swap(y);
  std::swap(space_dim, y.space_dim);
}

template <typename PH>
dimension_type
Polyhedra_PowerSet<PH>::space_dimension() const {
  return space_dim;
}

template <typename PH>
void
Polyhedra_PowerSet<PH>::concatenate_assign(const Polyhedra_PowerSet& y) {

  Base::omega_reduce();
  y.Base::omega_reduce();
  Sequence new_sequence;
  const Polyhedra_PowerSet<PH>& x = *this;
  for (const_iterator xi = x.begin(), x_end = x.end(),
  y_begin = y.begin(), y_end = y.end(); xi != x_end; ) {
    for (const_iterator yi = y_begin; yi != y_end; ++yi) {
      CS zi = *xi;
      zi.concatenate_assign(*yi);
      assert(!zi.is_bottom());
      new_sequence.push_back(zi);
    }
    ++xi;
    if (abandon_expensive_computations && xi != x_end && y_begin != y_end) {

      PH xph = xi->element();
      for (++xi; xi != x_end; ++xi)
 xph.poly_hull_assign(xi->element());
      const_iterator yi = y_begin;
      PH yph = yi->element();
      for (++yi; yi != y_end; ++yi)
 yph.poly_hull_assign(yi->element());
      xph.concatenate_assign(yph);
      std::swap(Base::sequence, new_sequence);
      add_disjunct(xph);
      goto done;
    }
  }
  std::swap(Base::sequence, new_sequence);
 done:
  space_dim += y.space_dim;
  assert(OK());
}

template <typename PH>
void
Polyhedra_PowerSet<PH>::add_constraint(const Constraint& c) {
  for (iterator xi = Base::begin(), x_end = Base::end(); xi != x_end; ++xi)
    xi->element().add_constraint(c);
  Base::reduced = false;
}

template <typename PH>
bool
Polyhedra_PowerSet<PH>::add_constraint_and_minimize(const Constraint& c) {
  for (iterator xi = Base::begin(), xin = xi, x_end = Base::end(); xi != x_end; xi = xin) {
    ++xin;
    if (!xi->element().add_constraint_and_minimize(c)) {
      erase(xi);
      x_end = Base::end();
    }
    else
      Base::reduced = false;
  }
  return !Base::empty();
}

template <typename PH>
void
Polyhedra_PowerSet<PH>::add_constraints(const ConSys& cs) {
  for (iterator xi = Base::begin(), x_end = Base::end(); xi != x_end; ++xi)
    xi->element().add_constraints(cs);
  Base::reduced = false;
}

template <typename PH>
bool
Polyhedra_PowerSet<PH>::add_constraints_and_minimize(const ConSys& cs) {
  for (iterator xi = Base::begin(),
  xin = xi, x_end = Base::end(); xi != x_end; xi = xin) {
    ++xin;
    if (!xi->element().add_constraints_and_minimize(cs)) {
      erase(xi);
      x_end = Base::end();
    }
    else
      Base::reduced = false;
  }
  return !Base::empty();
}

template <typename PH>
void
Polyhedra_PowerSet<PH>::add_dimensions_and_embed(dimension_type m) {
  for (iterator i = Base::begin(), send = Base::end(); i != send; ++i)
    i->add_dimensions_and_embed(m);
  space_dim += m;
  assert(OK());
}

template <typename PH>
void
Polyhedra_PowerSet<PH>::add_dimensions_and_project(dimension_type m) {
  for (iterator i = Base::begin(), send = Base::end(); i != send; ++i)
    i->add_dimensions_and_project(m);
  space_dim += m;
  assert(OK());
}

template <typename PH>
void
Polyhedra_PowerSet<PH>::remove_dimensions(const Variables_Set& to_be_removed) {
  Variables_Set::size_type num_removed = to_be_removed.size();
  if (num_removed > 0) {
    for (iterator i = Base::begin(), send = Base::end(); i != send; ++i) {
      i->remove_dimensions(to_be_removed);
      Base::reduced = false;
    }
    space_dim -= num_removed;
    assert(OK());
  }
}

template <typename PH>
void
Polyhedra_PowerSet<PH>::remove_higher_dimensions(dimension_type
       new_dimension) {
  if (new_dimension < space_dim) {
    for (iterator i = Base::begin(), send = Base::end(); i != send; ++i) {
      i->remove_higher_dimensions(new_dimension);
      Base::reduced = false;
    }
    space_dim = new_dimension;
    assert(OK());
  }
}

template <typename PH>
bool
Polyhedra_PowerSet<PH>::
geometrically_covers(const Polyhedra_PowerSet& y) const {
  for (const_iterator yi = y.Base::begin(),
  y_end = y.Base::end(); yi != y_end; ++yi)
    if (!check_containment(yi->element(), *this))
      return false;
  return true;
}

template <typename PH>
bool
Polyhedra_PowerSet<PH>::
geometrically_equals(const Polyhedra_PowerSet& y) const {
  const Polyhedra_PowerSet& x = * this;
  return x.geometrically_covers(y) && y.geometrically_covers(x);
}

template <typename PH>
template <typename PartialFunction>
void
Polyhedra_PowerSet<PH>::map_dimensions(const PartialFunction& pfunc) {
  if (Base::is_bottom()) {
    dimension_type n = 0;
    for (dimension_type i = space_dim; i-- > 0; ) {
      dimension_type new_i;
      if (pfunc.maps(i, new_i))
 ++n;
    }
    space_dim = n;
  }
  else {
    iterator sbegin = Base::begin();
    for (iterator i = sbegin, send = Base::end(); i != send; ++i)
      i->map_dimensions(pfunc);
    space_dim = sbegin->space_dimension();
    Base::reduced = false;
  }
  assert(OK());
}

template <typename PH>
void
Polyhedra_PowerSet<PH>::pairwise_reduce() {

  Base::omega_reduce();

  size_type n = Base::size();
  size_type deleted;
  do {
    Sequence new_sequence;
    std::deque<bool> marked(n, false);
    deleted = 0;
    iterator sbegin = Base::begin();
    iterator send = Base::end();
    unsigned i_index = 0;
    for (iterator i = sbegin; i != send; ++i, ++i_index) {
      if (marked[i_index])
 continue;
      PH& pi = i->element();
      const_iterator j = i;
      int j_index = i_index;
      for (++j, ++j_index; j != send; ++j, ++j_index) {
 if (marked[j_index])
   continue;
 const PH& pj = j->element();
 if (poly_hull_assign_if_exact(pi, pj)) {
   marked[i_index] = marked[j_index] = true;
   add_non_bottom_disjunct(new_sequence, pi);
   ++deleted;
   goto next;
 }
      }
    next:
      ;
    }
    iterator nsbegin = new_sequence.begin();
    iterator nsend = new_sequence.end();
    i_index = 0;
    for (const_iterator i = sbegin; i != send; ++i, ++i_index)
      if (!marked[i_index])
 add_non_bottom_disjunct(new_sequence, *i, nsbegin, nsend);
    std::swap(Base::sequence, new_sequence);
    n -= deleted;
  } while (deleted > 0);
  assert(OK());
}

template <typename PH>
template <typename Widening>
void
Polyhedra_PowerSet<PH>::
BGP99_heuristics_assign(const Polyhedra_PowerSet& y, Widening wf) {
  size_type n = Base::size();
  Sequence new_sequence;
  std::deque<bool> marked(n, false);
  const_iterator sbegin = Base::begin();
  const_iterator send = Base::end();
  unsigned i_index = 0;
  for (const_iterator i = sbegin; i != send; ++i, ++i_index)
    for (const_iterator j = y.Base::begin(),
    y_end = y.Base::end(); j != y_end; ++j) {
      const PH& pi = i->element();
      const PH& pj = j->element();
      if (pi.contains(pj)) {
 PH pi_copy = pi;
 wf(pi_copy, pj);
 add_non_bottom_disjunct(new_sequence, pi_copy);
 marked[i_index] = true;
      }
    }
  iterator nsbegin = new_sequence.begin();
  iterator nsend = new_sequence.end();
  i_index = 0;
  for (const_iterator i = sbegin; i != send; ++i, ++i_index)
    if (!marked[i_index])
      add_non_bottom_disjunct(new_sequence, *i, nsbegin, nsend);
  std::swap(Base::sequence, new_sequence);
  assert(OK());
  assert(Base::is_omega_reduced());
}

template <typename PH>
template <typename Widening>
void
Polyhedra_PowerSet<PH>::
BGP99_extrapolation_assign(const Polyhedra_PowerSet& y,
      Widening wf,
      unsigned max_disjuncts) {

  Polyhedra_PowerSet<PH>& x = *this;
  x.pairwise_reduce();
  if (max_disjuncts != 0)
    x.collapse(max_disjuncts);
  x.BGP99_heuristics_assign(y, wf);
}

template <typename PH>
template <typename Cert>
void
Polyhedra_PowerSet<PH>::
collect_certificates(std::map<Cert, size_type,
                typename Cert::Compare>& cert_ms) const {
  assert(Base::is_omega_reduced());
  assert(cert_ms.size() == 0);
  for (const_iterator i = Base::begin(),
  iend = Base::end(); i != iend; i++) {
    Cert ph_cert(i->element());
    ++cert_ms[ph_cert];
  }
}

template <typename PH>
template <typename Cert>
bool
Polyhedra_PowerSet<PH>::
is_cert_multiset_stabilizing(const std::map<Cert, size_type,
                       typename Cert::Compare>& y_cert_ms
        ) const {
  typedef std::map<Cert, size_type, typename Cert::Compare> Cert_Multiset;
  Cert_Multiset x_cert_ms;
  collect_certificates(x_cert_ms);
  typename Cert_Multiset::const_iterator
    xi = x_cert_ms.begin(),
    xend = x_cert_ms.end(),
    yi = y_cert_ms.begin(),
    yend = y_cert_ms.end();
  while (xi != xend && yi != yend) {
    const Cert& xi_cert = xi->first;
    const Cert& yi_cert = yi->first;
    switch (xi_cert.compare(yi_cert)) {
    case 0:

      {
 const size_type& xi_count = xi->second;
 const size_type& yi_count = yi->second;
 if (xi_count == yi_count) {

   ++xi;
   ++yi;
 }
 else

   return xi_count < yi_count;
 break;
      }
    case 1:

      return false;

    case -1:

      return true;
    }
  }


  return (yi != yend);
}

template <typename PH>
template <typename Cert, typename Widening>
void
Polyhedra_PowerSet<PH>::BHZ03_widening_assign(const Polyhedra_PowerSet& y,
           Widening wf) {

  Polyhedra_PowerSet<PH>& x = *this;
  assert(x.size() > 0);
  if (y.size() == 0)
    return;


  PH x_hull(x.space_dim, PH::EMPTY);
  for (const_iterator i = x.begin(), iend = x.end(); i != iend; ++i)
    x_hull.poly_hull_assign(i->element());


  PH y_hull(y.space_dim, PH::EMPTY);
  for (const_iterator i = y.begin(), iend = y.end(); i != iend; ++i)
    y_hull.poly_hull_assign(i->element());

  const Cert y_hull_cert(y_hull);


  int hull_stabilization = y_hull_cert.compare(x_hull);
  if (hull_stabilization == 1)
    return;


  const bool y_is_not_a_singleton = y.size() > 1;



  typedef std::map<Cert, size_type, typename Cert::Compare> Cert_Multiset;
  Cert_Multiset y_cert_ms;
  bool y_cert_ms_computed = false;

  if (hull_stabilization == 0 && y_is_not_a_singleton) {

    y.collect_certificates(y_cert_ms);
    y_cert_ms_computed = true;

    if (x.is_cert_multiset_stabilizing(y_cert_ms))
      return;
  }


  Polyhedra_PowerSet<PH> bgp99_heuristics = x;
  bgp99_heuristics.BGP99_heuristics_assign(y, wf);


  PH bgp99_heuristics_hull(x.space_dim, PH::EMPTY);
  for (const_iterator i = bgp99_heuristics.begin(),
  iend = bgp99_heuristics.end(); i != iend; ++i)
    bgp99_heuristics_hull.poly_hull_assign(i->element());



  hull_stabilization = y_hull_cert.compare(bgp99_heuristics_hull);
  if (hull_stabilization == 1) {

    std::swap(x, bgp99_heuristics);
    return;
  }
  else if (hull_stabilization == 0 && y_is_not_a_singleton) {

    if (!y_cert_ms_computed) {
      y.collect_certificates(y_cert_ms);
      y_cert_ms_computed = true;
    }
    if (bgp99_heuristics.is_cert_multiset_stabilizing(y_cert_ms)) {
      std::swap(x, bgp99_heuristics);
      return;
    }




    Polyhedra_PowerSet<PH> reduced_bgp99_heuristics(bgp99_heuristics);
    reduced_bgp99_heuristics.pairwise_reduce();
    if (reduced_bgp99_heuristics.is_cert_multiset_stabilizing(y_cert_ms)) {
      std::swap(x, reduced_bgp99_heuristics);
      return;
    }
  }



  if (bgp99_heuristics_hull.strictly_contains(y_hull)) {

    PH ph = bgp99_heuristics_hull;
    wf(ph, y_hull);

    ph.poly_difference_assign(bgp99_heuristics_hull);
    x.add_disjunct(ph);
    return;
  }


  Polyhedra_PowerSet<PH> x_hull_singleton(x.space_dim, PH::EMPTY);
  x_hull_singleton.add_disjunct(x_hull);
  std::swap(x, x_hull_singleton);
}

template <typename PH>
template <typename Widening>
inline void
Polyhedra_PowerSet<PH>::
BHZ03_widening_assign(const Polyhedra_PowerSet& y, Widening wf) {
  BHZ03_widening_assign<BHRZ03_Certificate>(y, wf);
}

template <typename PH>
void
Polyhedra_PowerSet<PH>::ascii_dump(std::ostream& s) const {
  s << "size " << Base::size()
    << "\nspace_dim " << space_dim
    << std::endl;
  for (const_iterator i = Base::begin(), send = Base::end(); i != send; ++i)
    i->element().ascii_dump(s);
}

template <typename PH>
bool
Polyhedra_PowerSet<PH>::ascii_load(std::istream& s) {
  std::string str;

  if (!(s >> str) || str != "size")
    return false;

  size_type sz;

  if (!(s >> sz))
    return false;

  if (!(s >> str) || str != "space_dim")
    return false;

  if (!(s >> space_dim))
    return false;

  Polyhedra_PowerSet new_pps(space_dim, Polyhedron::EMPTY);
  while (sz-- > 0) {
    PH ph;
    if (!ph.ascii_load(s))
      return false;
    new_pps.add_disjunct(ph);
  }
  swap(new_pps);


  assert(OK());
  return true;
}

template <typename PH>
bool
Polyhedra_PowerSet<PH>::OK() const {
  for (const_iterator i = Base::begin(), send = Base::end(); i != send; ++i)
    if (i->space_dimension() != space_dim) {






      return false;
    }
  return Base::OK();
}

}


namespace std {


template <typename PH>
inline void
swap(Parma_Polyhedra_Library::Polyhedra_PowerSet<PH>& x,
     Parma_Polyhedra_Library::Polyhedra_PowerSet<PH>& y) {
  x.swap(y);
}

}
namespace Parma_Polyhedra_Library {

inline dimension_type
Polyhedron::max_space_dimension() {
  return std::min(C_Polyhedron::max_space_dimension(),
    NNC_Polyhedron::max_space_dimension());
}


inline dimension_type
max_space_dimension() {
  using std::min;
  return
    min(Polyhedron::max_space_dimension(),
 min(Polyhedra_PowerSet<C_Polyhedron>::max_space_dimension(),
     Polyhedra_PowerSet<NNC_Polyhedron>::max_space_dimension()
     )
 );
}

}
namespace Parma_Polyhedra_Library {
template <typename PH>
std::pair<PH, Polyhedra_PowerSet<NNC_Polyhedron> >
linear_partition(const PH& p, const PH& q);

namespace {

template <typename PH>
void
linear_partition_aux(const Constraint& c,
       PH& qq,
       Polyhedra_PowerSet<NNC_Polyhedron>& r) {
  LinExpression le(c);
  Constraint neg_c = c.is_strict_inequality() ? (le <= 0) : (le < 0);
  NNC_Polyhedron qqq(qq);
  if (qqq.add_constraint_and_minimize(neg_c))
    r.add_disjunct(qqq);
  qq.add_constraint(c);
}

}


template <typename PH>
std::pair<PH, Polyhedra_PowerSet<NNC_Polyhedron> >
linear_partition(const PH& p, const PH& q) {
  Polyhedra_PowerSet<NNC_Polyhedron> r(p.space_dimension(),
           Polyhedron::EMPTY);
  PH qq = q;
  const ConSys& pcs = p.constraints();
  for (ConSys::const_iterator i = pcs.begin(),
  pcs_end = pcs.end(); i != pcs_end; ++i) {
    const Constraint c = *i;
    if (c.is_equality()) {
      LinExpression le(c);
      linear_partition_aux(le <= 0, qq, r);
      linear_partition_aux(le >= 0, qq, r);
    }
    else
      linear_partition_aux(c, qq, r);
  }
  return std::pair<PH, Polyhedra_PowerSet<NNC_Polyhedron> >(qq, r);
}



template <typename PH>
bool
poly_hull_assign_if_exact(PH& p, const PH& q);


template <typename PH>
bool
poly_hull_assign_if_exact(PH& p, const PH& q) {
  PH phull = p;
  NNC_Polyhedron nnc_p(p);
  phull.poly_hull_assign(q);
  std::pair<PH, Polyhedra_PowerSet<NNC_Polyhedron> >
    partition = linear_partition(q, phull);
  const Polyhedra_PowerSet<NNC_Polyhedron>& s = partition.second;
  typedef Polyhedra_PowerSet<NNC_Polyhedron>::const_iterator iter;
  for (iter i = s.begin(), s_end = s.end(); i != s_end; ++i)



    if (!nnc_p.contains(i->element()))
      return false;
  p = phull;
  return true;
}

template <typename PH>
bool
check_containment(const PH& ph, const Polyhedra_PowerSet<PH>& ps);

template <typename PH>
bool
check_containment(const PH& ph, const Polyhedra_PowerSet<PH>& ps) {
  Polyhedra_PowerSet<NNC_Polyhedron> tmp(ph.space_dimension(),
      Polyhedron::EMPTY);
  tmp.add_disjunct(NNC_Polyhedron(ph));
  for (typename Polyhedra_PowerSet<PH>::const_iterator i = ps.begin(),
  ps_end = ps.end(); i != ps_end; ++i) {
    const NNC_Polyhedron pi(i->element());
    for (typename Polyhedra_PowerSet<NNC_Polyhedron>::iterator j = tmp.begin(),
    jn = j; j != tmp.end(); j = jn) {
      ++jn;
      const NNC_Polyhedron& pj = j->element();
      if (pi.contains(pj))
 tmp.erase(j);
    }
    if (tmp.empty())
      return true;
    else {
      Polyhedra_PowerSet<NNC_Polyhedron> new_disjuncts(ph.space_dimension(),
             Polyhedron::EMPTY);
      for (Polyhedra_PowerSet<NNC_Polyhedron>::iterator j = tmp.begin(),
      jn = j; j != tmp.end(); j = jn) {
 ++jn;
 const NNC_Polyhedron& pj = j->element();
 if (!pj.is_disjoint_from(pi)) {
   std::pair<NNC_Polyhedron, Polyhedra_PowerSet<NNC_Polyhedron> >
     partition = linear_partition(pi, pj);
   tmp.erase(j);
   new_disjuncts.add_disjunct(partition.first);
   new_disjuncts.upper_bound_assign(partition.second);
 }
      }
      tmp.upper_bound_assign(new_disjuncts);
    }
  }
  return false;
}

}

#ifdef PPL_SAVE_NDEBUG
# define NDEBUG 1
# undef PPL_SAVE_NDEBUG
#else
# undef NDEBUG
#endif
#include <cassert>

#undef PPL_FLEXIBLE_ARRAY
#undef PACKAGE
#undef PACKAGE_BUGREPORT
#undef PACKAGE_NAME
#undef PACKAGE_STRING
#undef PACKAGE_TARNAME
#undef PACKAGE_VERSION
#undef VERSION

#endif
