The base class for convex polyhedra. More...

#include <Polyhedron.defs.hh>

Inheritance diagram for Parma_Polyhedra_Library::Polyhedron:

Collaboration diagram for Parma_Polyhedra_Library::Polyhedron:

Classes
class	Status
	A conjunctive assertion about a polyhedron. More...
Public Types
typedef Coefficient	coefficient_type
	The numeric type of coefficients.
Public Member Functions
template<typename Input >
Input &	check_obj_space_dimension_overflow (Input &input, const Topology topol, const char method, const char reason)
Member Functions that Do Not Modify the Polyhedron
dimension_type	space_dimension () const
	Returns the dimension of the vector space enclosing `*this`.
dimension_type	affine_dimension () const
	Returns , if `this` is empty; otherwise, returns the affine dimension of `this`.
const Constraint_System &	constraints () const
	Returns the system of constraints.
const Constraint_System &	minimized_constraints () const
	Returns the system of constraints, with no redundant constraint.
const Generator_System &	generators () const
	Returns the system of generators.
const Generator_System &	minimized_generators () const
	Returns the system of generators, with no redundant generator.
Congruence_System	congruences () const
	Returns a system of (equality) congruences satisfied by `*this`.
Congruence_System	minimized_congruences () const
	Returns a system of (equality) congruences satisfied by `this`, with no redundant congruences and having the same affine dimension as `this`.
Poly_Con_Relation	relation_with (const Constraint &c) const
	Returns the relations holding between the polyhedron `*this` and the constraint `c`.
Poly_Gen_Relation	relation_with (const Generator &g) const
	Returns the relations holding between the polyhedron `*this` and the generator `g`.
Poly_Con_Relation	relation_with (const Congruence &cg) const
	Returns the relations holding between the polyhedron `*this` and the congruence `c`.
bool	is_empty () const
	Returns `true` if and only if `*this` is an empty polyhedron.
bool	is_universe () const
	Returns `true` if and only if `*this` is a universe polyhedron.
bool	is_topologically_closed () const
	Returns `true` if and only if `*this` is a topologically closed subset of the vector space.
bool	is_disjoint_from (const Polyhedron &y) const
	Returns `true` if and only if `*this` and `y` are disjoint.
bool	is_discrete () const
	Returns `true` if and only if `*this` is discrete.
bool	is_bounded () const
	Returns `true` if and only if `*this` is a bounded polyhedron.
bool	contains_integer_point () const
	Returns `true` if and only if `*this` contains at least one integer point.
bool	constrains (Variable var) const
	Returns `true` if and only if `var` is constrained in `*this`.
bool	bounds_from_above (const Linear_Expression &expr) const
	Returns `true` if and only if `expr` is bounded from above in `*this`.
bool	bounds_from_below (const Linear_Expression &expr) const
	Returns `true` if and only if `expr` is bounded from below in `*this`.
bool	maximize (const Linear_Expression &expr, Coefficient &sup_n, Coefficient &sup_d, bool &maximum) const
	Returns `true` if and only if `this` is not empty and `expr` is bounded from above in `this`, in which case the supremum value is computed.
bool	maximize (const Linear_Expression &expr, Coefficient &sup_n, Coefficient &sup_d, bool &maximum, Generator &g) const
	Returns `true` if and only if `this` is not empty and `expr` is bounded from above in `this`, in which case the supremum value and a point where `expr` reaches it are computed.
bool	minimize (const Linear_Expression &expr, Coefficient &inf_n, Coefficient &inf_d, bool &minimum) const
	Returns `true` if and only if `this` is not empty and `expr` is bounded from below in `this`, in which case the infimum value is computed.
bool	minimize (const Linear_Expression &expr, Coefficient &inf_n, Coefficient &inf_d, bool &minimum, Generator &g) const
	Returns `true` if and only if `this` is not empty and `expr` is bounded from below in `this`, in which case the infimum value and a point where `expr` reaches it are computed.
bool	frequency (const Linear_Expression &expr, Coefficient &freq_n, Coefficient &freq_d, Coefficient &val_n, Coefficient &val_d) const
	Returns `true` if and only if there exist a unique value `val` such that `*this` saturates the equality `expr = val`.
bool	contains (const Polyhedron &y) const
	Returns `true` if and only if `*this` contains `y`.
bool	strictly_contains (const Polyhedron &y) const
	Returns `true` if and only if `*this` strictly contains `y`.
bool	OK (bool check_not_empty=false) const
	Checks if all the invariants are satisfied.
Space Dimension Preserving Member Functions that May Modify the Polyhedron
void	add_constraint (const Constraint &c)
	Adds a copy of constraint `c` to the system of constraints of `*this` (without minimizing the result).
void	add_generator (const Generator &g)
	Adds a copy of generator `g` to the system of generators of `*this` (without minimizing the result).
void	add_congruence (const Congruence &cg)
	Adds a copy of congruence `cg` to `*this`, if `cg` can be exactly represented by a polyhedron.
void	add_constraints (const Constraint_System &cs)
	Adds a copy of the constraints in `cs` to the system of constraints of `*this` (without minimizing the result).
void	add_recycled_constraints (Constraint_System &cs)
	Adds the constraints in `cs` to the system of constraints of `*this` (without minimizing the result).
void	add_generators (const Generator_System &gs)
	Adds a copy of the generators in `gs` to the system of generators of `*this` (without minimizing the result).
void	add_recycled_generators (Generator_System &gs)
	Adds the generators in `gs` to the system of generators of `*this` (without minimizing the result).
void	add_congruences (const Congruence_System &cgs)
	Adds a copy of the congruences in `cgs` to `*this`, if all the congruences can be exactly represented by a polyhedron.
void	add_recycled_congruences (Congruence_System &cgs)
	Adds the congruences in `cgs` to `*this`, if all the congruences can be exactly represented by a polyhedron.
void	refine_with_constraint (const Constraint &c)
	Uses a copy of constraint `c` to refine `*this`.
void	refine_with_congruence (const Congruence &cg)
	Uses a copy of congruence `cg` to refine `*this`.
void	refine_with_constraints (const Constraint_System &cs)
	Uses a copy of the constraints in `cs` to refine `*this`.
void	refine_with_congruences (const Congruence_System &cgs)
	Uses a copy of the congruences in `cgs` to refine `*this`.
template<typename FP_Format , typename Interval_Info >
void	refine_with_linear_form_inequality (const Linear_Form< Interval< FP_Format, Interval_Info > > &left, const Linear_Form< Interval< FP_Format, Interval_Info > > &right, bool is_strict=false)
	Refines `*this` with the constraint expressed by `left` `right` if `is_strict` is set, with the constraint `left` $\leq$ `right` otherwise.
template<typename FP_Format , typename Interval_Info >
void	generalized_refine_with_linear_form_inequality (const Linear_Form< Interval< FP_Format, Interval_Info > > &left, const Linear_Form< Interval< FP_Format, Interval_Info > > &right, Relation_Symbol relsym)
	Refines `*this` with the constraint expressed by `left` $\relsym$ `right`, where $\relsym$ is the relation symbol specified by `relsym`..
template<typename FP_Format , typename Interval_Info >
void	refine_fp_interval_abstract_store (Box< Interval< FP_Format, Interval_Info > > &store) const
	Refines `store` with the constraints defining `*this`.
void	unconstrain (Variable var)
	Computes the cylindrification of `this` with respect to space dimension `var`, assigning the result to `this`.
void	unconstrain (const Variables_Set &vars)
	Computes the cylindrification of `this` with respect to the set of space dimensions `vars`, assigning the result to `this`.
void	intersection_assign (const Polyhedron &y)
	Assigns to `this` the intersection of `this` and `y`.
void	poly_hull_assign (const Polyhedron &y)
	Assigns to `this` the poly-hull of `this` and `y`.
void	upper_bound_assign (const Polyhedron &y)
	Same as poly_hull_assign(y).
void	poly_difference_assign (const Polyhedron &y)
	Assigns to `this` the poly-difference of `this` and `y`.
void	difference_assign (const Polyhedron &y)
	Same as poly_difference_assign(y).
bool	simplify_using_context_assign (const Polyhedron &y)
	Assigns to `this` a meet-preserving simplification of `this` with respect to `y`. If `false` is returned, then the intersection is empty.
void	affine_image (Variable var, const Linear_Expression &expr, Coefficient_traits::const_reference denominator=Coefficient_one())
	Assigns to `this` the affine image of `this` under the function mapping variable `var` to the affine expression specified by `expr` and `denominator`.
template<typename FP_Format , typename Interval_Info >
void	affine_form_image (Variable var, const Linear_Form< Interval< FP_Format, Interval_Info > > &lf)
void	affine_preimage (Variable var, const Linear_Expression &expr, Coefficient_traits::const_reference denominator=Coefficient_one())
	Assigns to `this` the affine preimage of `this` under the function mapping variable `var` to the affine expression specified by `expr` and `denominator`.
void	generalized_affine_image (Variable var, Relation_Symbol relsym, const Linear_Expression &expr, Coefficient_traits::const_reference denominator=Coefficient_one())
	Assigns to `this` the image of `this` with respect to the generalized affine relation $\mathrm{var}' \relsym \frac{\mathrm{expr}}{\mathrm{denominator}}$ , where $\mathord{\relsym}$ is the relation symbol encoded by `relsym`.
void	generalized_affine_preimage (Variable var, Relation_Symbol relsym, const Linear_Expression &expr, Coefficient_traits::const_reference denominator=Coefficient_one())
	Assigns to `this` the preimage of `this` with respect to the generalized affine relation $\mathrm{var}' \relsym \frac{\mathrm{expr}}{\mathrm{denominator}}$ , where $\mathord{\relsym}$ is the relation symbol encoded by `relsym`.
void	generalized_affine_image (const Linear_Expression &lhs, Relation_Symbol relsym, const Linear_Expression &rhs)
	Assigns to `this` the image of `this` with respect to the generalized affine relation $\mathrm{lhs}' \relsym \mathrm{rhs}$ , where $\mathord{\relsym}$ is the relation symbol encoded by `relsym`.
void	generalized_affine_preimage (const Linear_Expression &lhs, Relation_Symbol relsym, const Linear_Expression &rhs)
	Assigns to `this` the preimage of `this` with respect to the generalized affine relation $\mathrm{lhs}' \relsym \mathrm{rhs}$ , where $\mathord{\relsym}$ is the relation symbol encoded by `relsym`.
void	bounded_affine_image (Variable var, const Linear_Expression &lb_expr, const Linear_Expression &ub_expr, Coefficient_traits::const_reference denominator=Coefficient_one())
	Assigns to `this` the image of `this` with respect to the bounded affine relation $\frac{\mathrm{lb\_expr}}{\mathrm{denominator}} \leq \mathrm{var}' \leq \frac{\mathrm{ub\_expr}}{\mathrm{denominator}}$ .
void	bounded_affine_preimage (Variable var, const Linear_Expression &lb_expr, const Linear_Expression &ub_expr, Coefficient_traits::const_reference denominator=Coefficient_one())
	Assigns to `this` the preimage of `this` with respect to the bounded affine relation $\frac{\mathrm{lb\_expr}}{\mathrm{denominator}} \leq \mathrm{var}' \leq \frac{\mathrm{ub\_expr}}{\mathrm{denominator}}$ .
void	time_elapse_assign (const Polyhedron &y)
	Assigns to `this` the result of computing the time-elapse between `this` and `y`.
void	wrap_assign (const Variables_Set &vars, Bounded_Integer_Type_Width w, Bounded_Integer_Type_Representation r, Bounded_Integer_Type_Overflow o, const Constraint_System *cs_p=0, unsigned complexity_threshold=16, bool wrap_individually=true)
	Wraps the specified dimensions of the vector space.
void	drop_some_non_integer_points (Complexity_Class complexity=ANY_COMPLEXITY)
	Possibly tightens `*this` by dropping some points with non-integer coordinates.
void	drop_some_non_integer_points (const Variables_Set &vars, Complexity_Class complexity=ANY_COMPLEXITY)
	Possibly tightens `*this` by dropping some points with non-integer coordinates for the space dimensions corresponding to `vars`.
void	topological_closure_assign ()
	Assigns to `*this` its topological closure.
void	BHRZ03_widening_assign (const Polyhedron &y, unsigned *tp=0)
	Assigns to `this` the result of computing the BHRZ03-widening between `this` and `y`.
void	limited_BHRZ03_extrapolation_assign (const Polyhedron &y, const Constraint_System &cs, unsigned *tp=0)
	Assigns to `this` the result of computing the limited extrapolation between `this` and `y` using the BHRZ03-widening operator.
void	bounded_BHRZ03_extrapolation_assign (const Polyhedron &y, const Constraint_System &cs, unsigned *tp=0)
	Assigns to `this` the result of computing the bounded extrapolation between `this` and `y` using the BHRZ03-widening operator.
void	H79_widening_assign (const Polyhedron &y, unsigned *tp=0)
	Assigns to `this` the result of computing the H79_widening between `this` and `y`.
void	widening_assign (const Polyhedron &y, unsigned *tp=0)
	Same as H79_widening_assign(y, tp).
void	limited_H79_extrapolation_assign (const Polyhedron &y, const Constraint_System &cs, unsigned *tp=0)
	Assigns to `this` the result of computing the limited extrapolation between `this` and `y` using the H79-widening operator.
void	bounded_H79_extrapolation_assign (const Polyhedron &y, const Constraint_System &cs, unsigned *tp=0)
	Assigns to `this` the result of computing the bounded extrapolation between `this` and `y` using the H79-widening operator.
Member Functions that May Modify the Dimension of the Vector Space
void	add_space_dimensions_and_embed (dimension_type m)
	Adds `m` new space dimensions and embeds the old polyhedron in the new vector space.
void	add_space_dimensions_and_project (dimension_type m)
	Adds `m` new space dimensions to the polyhedron and does not embed it in the new vector space.
void	concatenate_assign (const Polyhedron &y)
	Assigns to `this` the concatenation of `this` and `y`, taken in this order.
void	remove_space_dimensions (const Variables_Set &vars)
	Removes all the specified dimensions from the vector space.
void	remove_higher_space_dimensions (dimension_type new_dimension)
	Removes the higher dimensions of the vector space so that the resulting space will have dimension `new_dimension`.
template<typename Partial_Function >
void	map_space_dimensions (const Partial_Function &pfunc)
	Remaps the dimensions of the vector space according to a partial function.
void	expand_space_dimension (Variable var, dimension_type m)
	Creates `m` copies of the space dimension corresponding to `var`.
void	fold_space_dimensions (const Variables_Set &vars, Variable dest)
	Folds the space dimensions in `vars` into `dest`.
Miscellaneous Member Functions
	~Polyhedron ()
	Destructor.
void	m_swap (Polyhedron &y)
	Swaps `this` with polyhedron `y`. (`this` and `y` can be dimension-incompatible.)
void	ascii_dump () const
	Writes to `std::cerr` an ASCII representation of `*this`.
void	ascii_dump (std::ostream &s) const
	Writes to `s` an ASCII representation of `*this`.
void	print () const
	Prints `*this` to `std::cerr` using `operator<<`.
bool	ascii_load (std::istream &s)
	Loads from `s` an ASCII representation (as produced by ascii_dump(std::ostream&) const) and sets `*this` accordingly. Returns `true` if successful, `false` otherwise.
memory_size_type	total_memory_in_bytes () const
	Returns the total size in bytes of the memory occupied by `*this`.
memory_size_type	external_memory_in_bytes () const
	Returns the size in bytes of the memory managed by `*this`.
int32_t	hash_code () const
	Returns a 32-bit hash code for `*this`.
Static Public Member Functions
static dimension_type	max_space_dimension ()
	Returns the maximum space dimension all kinds of Polyhedron can handle.
static bool	can_recycle_constraint_systems ()
	Returns `true` indicating that this domain has methods that can recycle constraints.
static void	initialize ()
	Initializes the class.
static void	finalize ()
	Finalizes the class.
static bool	can_recycle_congruence_systems ()
	Returns `false` indicating that this domain cannot recycle congruences.
Protected Member Functions
	Polyhedron (Topology topol, dimension_type num_dimensions, Degenerate_Element kind)
	Builds a polyhedron having the specified properties.
	Polyhedron (const Polyhedron &y, Complexity_Class complexity=ANY_COMPLEXITY)
	Ordinary copy constructor.
	Polyhedron (Topology topol, const Constraint_System &cs)
	Builds a polyhedron from a system of constraints.
	Polyhedron (Topology topol, Constraint_System &cs, Recycle_Input dummy)
	Builds a polyhedron recycling a system of constraints.
	Polyhedron (Topology topol, const Generator_System &gs)
	Builds a polyhedron from a system of generators.
	Polyhedron (Topology topol, Generator_System &gs, Recycle_Input dummy)
	Builds a polyhedron recycling a system of generators.
template<typename Interval >
	Polyhedron (Topology topol, const Box< Interval > &box, Complexity_Class complexity=ANY_COMPLEXITY)
	Builds a polyhedron from a box.
Polyhedron &	operator= (const Polyhedron &y)
	The assignment operator. (`*this` and `y` can be dimension-incompatible.)
bool	BFT00_poly_hull_assign_if_exact (const Polyhedron &y)
	If the poly-hull of `this` and `y` is exact it is assigned to `this` and `true` is returned, otherwise `false` is returned.
bool	BHZ09_poly_hull_assign_if_exact (const Polyhedron &y)
bool	BHZ09_C_poly_hull_assign_if_exact (const Polyhedron &y)
bool	BHZ09_NNC_poly_hull_assign_if_exact (const Polyhedron &y)
void	drop_some_non_integer_points (const Variables_Set *vars_p, Complexity_Class complexity)
	Possibly tightens `this` by dropping some points with non-integer coordinates for the space dimensions corresponding to `vars_p`.
template<typename FP_Format , typename Interval_Info >
void	overapproximate_linear_form (const Linear_Form< Interval< FP_Format, Interval_Info > > &lf, const dimension_type lf_dimension, Linear_Form< Interval< FP_Format, Interval_Info > > &result)
	Helper function that overapproximates an interval linear form.
Static Protected Member Functions
template<typename FP_Format , typename Interval_Info >
static void	convert_to_integer_expression (const Linear_Form< Interval< FP_Format, Interval_Info > > &lf, const dimension_type lf_dimension, Linear_Expression &result)
	Helper function that makes `result` become a Linear_Expression obtained by normalizing the denominators in `lf`.
template<typename FP_Format , typename Interval_Info >
static void	convert_to_integer_expressions (const Linear_Form< Interval< FP_Format, Interval_Info > > &lf, const dimension_type lf_dimension, Linear_Expression &res, Coefficient &res_low_coeff, Coefficient &res_hi_coeff, Coefficient &denominator)
	Normalization helper function.
Private Types
enum	Three_Valued_Boolean { TVB_TRUE, TVB_FALSE, TVB_DONT_KNOW }
Private Member Functions
Topology	topology () const
	Returns the topological kind of the polyhedron.
bool	is_necessarily_closed () const
	Returns `true` if and only if the polyhedron is necessarily closed.
void	refine_no_check (const Constraint &c)
	Uses a copy of constraint `c` to refine the system of constraints of `*this`.
Three_Valued_Boolean	quick_equivalence_test (const Polyhedron &y) const
	Polynomial but incomplete equivalence test between polyhedra.
bool	is_included_in (const Polyhedron &y) const
	Returns `true` if and only if `*this` is included in `y`.
bool	bounds (const Linear_Expression &expr, bool from_above) const
	Checks if and how `expr` is bounded in `*this`.
bool	max_min (const Linear_Expression &expr, bool maximize, Coefficient &ext_n, Coefficient &ext_d, bool &included, Generator &g) const
	Maximizes or minimizes `expr` subject to `*this`.
Private Verifiers: Verify if Individual Flags are Set
bool	marked_empty () const
	Returns `true` if the polyhedron is known to be empty.
bool	constraints_are_up_to_date () const
	Returns `true` if the system of constraints is up-to-date.
bool	generators_are_up_to_date () const
	Returns `true` if the system of generators is up-to-date.
bool	constraints_are_minimized () const
	Returns `true` if the system of constraints is minimized.
bool	generators_are_minimized () const
	Returns `true` if the system of generators is minimized.
bool	has_pending_constraints () const
	Returns `true` if there are pending constraints.
bool	has_pending_generators () const
	Returns `true` if there are pending generators.
bool	has_something_pending () const
	Returns `true` if there are either pending constraints or pending generators.
bool	can_have_something_pending () const
	Returns `true` if the polyhedron can have something pending.
bool	sat_c_is_up_to_date () const
	Returns `true` if the saturation matrix `sat_c` is up-to-date.
bool	sat_g_is_up_to_date () const
	Returns `true` if the saturation matrix `sat_g` is up-to-date.
State Flag Setters: Set Only the Specified Flags
void	set_zero_dim_univ ()
	Sets `status` to express that the polyhedron is the universe 0-dimension vector space, clearing all corresponding matrices.
void	set_empty ()
	Sets `status` to express that the polyhedron is empty, clearing all corresponding matrices.
void	set_constraints_up_to_date ()
	Sets `status` to express that constraints are up-to-date.
void	set_generators_up_to_date ()
	Sets `status` to express that generators are up-to-date.
void	set_constraints_minimized ()
	Sets `status` to express that constraints are minimized.
void	set_generators_minimized ()
	Sets `status` to express that generators are minimized.
void	set_constraints_pending ()
	Sets `status` to express that constraints are pending.
void	set_generators_pending ()
	Sets `status` to express that generators are pending.
void	set_sat_c_up_to_date ()
	Sets `status` to express that `sat_c` is up-to-date.
void	set_sat_g_up_to_date ()
	Sets `status` to express that `sat_g` is up-to-date.
State Flag Cleaners: Clear Only the Specified Flag
void	clear_empty ()
	Clears the `status` flag indicating that the polyhedron is empty.
void	clear_constraints_up_to_date ()
	Sets `status` to express that constraints are no longer up-to-date.
void	clear_generators_up_to_date ()
	Sets `status` to express that generators are no longer up-to-date.
void	clear_constraints_minimized ()
	Sets `status` to express that constraints are no longer minimized.
void	clear_generators_minimized ()
	Sets `status` to express that generators are no longer minimized.
void	clear_pending_constraints ()
	Sets `status` to express that there are no longer pending constraints.
void	clear_pending_generators ()
	Sets `status` to express that there are no longer pending generators.
void	clear_sat_c_up_to_date ()
	Sets `status` to express that `sat_c` is no longer up-to-date.
void	clear_sat_g_up_to_date ()
	Sets `status` to express that `sat_g` is no longer up-to-date.
The Handling of Pending Rows
bool	process_pending () const
	Processes the pending rows of either description of the polyhedron and obtains a minimized polyhedron.
bool	process_pending_constraints () const
	Processes the pending constraints and obtains a minimized polyhedron.
void	process_pending_generators () const
	Processes the pending generators and obtains a minimized polyhedron.
void	remove_pending_to_obtain_constraints () const
	Lazily integrates the pending descriptions of the polyhedron to obtain a constraint system without pending rows.
bool	remove_pending_to_obtain_generators () const
	Lazily integrates the pending descriptions of the polyhedron to obtain a generator system without pending rows.
Updating and Sorting Matrices
void	update_constraints () const
	Updates constraints starting from generators and minimizes them.
bool	update_generators () const
	Updates generators starting from constraints and minimizes them.
void	update_sat_c () const
	Updates `sat_c` using the updated constraints and generators.
void	update_sat_g () const
	Updates `sat_g` using the updated constraints and generators.
void	obtain_sorted_constraints () const
	Sorts the matrix of constraints keeping status consistency.
void	obtain_sorted_generators () const
	Sorts the matrix of generators keeping status consistency.
void	obtain_sorted_constraints_with_sat_c () const
	Sorts the matrix of constraints and updates `sat_c`.
void	obtain_sorted_generators_with_sat_g () const
	Sorts the matrix of generators and updates `sat_g`.
Weak and Strong Minimization of Descriptions
bool	minimize () const
	Applies (weak) minimization to both the constraints and generators.
bool	strongly_minimize_constraints () const
	Applies strong minimization to the constraints of an NNC polyhedron.
bool	strongly_minimize_generators () const
	Applies strong minimization to the generators of an NNC polyhedron.
Constraint_System	simplified_constraints () const
	If constraints are up-to-date, obtain a simplified copy of them.
Widening- and Extrapolation-Related Functions
void	select_CH78_constraints (const Polyhedron &y, Constraint_System &cs_selection) const
	Copies to `cs_selection` the constraints of `y` corresponding to the definition of the CH78-widening of `*this` and `y`.
void	select_H79_constraints (const Polyhedron &y, Constraint_System &cs_selected, Constraint_System &cs_not_selected) const
	Splits the constraints of `x' into two subsets, depending on whether or not they are selected to compute the H79-widening of `*this` and `y`.
bool	BHRZ03_combining_constraints (const Polyhedron &y, const BHRZ03_Certificate &y_cert, const Polyhedron &H79, const Constraint_System &x_minus_H79_cs)
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 Private Member Functions
static void	add_space_dimensions (Linear_System &sys1, Linear_System &sys2, Bit_Matrix &sat1, Bit_Matrix &sat2, dimension_type add_dim)
	Adds new space dimensions to the given linear systems.
Minimization-Related Static Member Functions
static bool	minimize (bool con_to_gen, Linear_System &source, Linear_System &dest, Bit_Matrix &sat)
	Builds and simplifies constraints from generators (or vice versa).
static bool	add_and_minimize (bool con_to_gen, Linear_System &source1, Linear_System &dest, Bit_Matrix &sat, const Linear_System &source2)
	Adds given constraints and builds minimized corresponding generators or vice versa.
static bool	add_and_minimize (bool con_to_gen, Linear_System &source, Linear_System &dest, Bit_Matrix &sat)
	Adds given constraints and builds minimized corresponding generators or vice versa. The given constraints are in `source`.
static dimension_type	conversion (Linear_System &source, dimension_type start, Linear_System &dest, Bit_Matrix &sat, dimension_type num_lines_or_equalities)
	Performs the conversion from constraints to generators and vice versa.
static dimension_type	simplify (Linear_System &sys, Bit_Matrix &sat)
	Uses Gauss' elimination method to simplify the result of `conversion()`.
Private Attributes
Constraint_System	con_sys
	The system of constraints.
Generator_System	gen_sys
	The system of generators.
Bit_Matrix	sat_c
	The saturation matrix having constraints on its columns.
Bit_Matrix	sat_g
	The saturation matrix having generators on its columns.
Status	status
	The status flags to keep track of the polyhedron's internal state.
dimension_type	space_dim
	The number of dimensions of the enclosing vector space.
Static Private Attributes
static dimension_type *	simplify_num_saturators_p = 0
	Pointer to an array used by simplify().
static size_t	simplify_num_saturators_size = 0
	Dimension of an array used by simplify().
Friends
class	Parma_Polyhedra_Library::Box
class	Parma_Polyhedra_Library::BD_Shape
class	Parma_Polyhedra_Library::Octagonal_Shape
class	Parma_Polyhedra_Library::Grid
class	Parma_Polyhedra_Library::BHRZ03_Certificate
class	Parma_Polyhedra_Library::H79_Certificate
bool	operator== (const Polyhedron &x, const Polyhedron &y)
bool	Parma_Polyhedra_Library::Interfaces::is_necessarily_closed_for_interfaces (const Polyhedron &)
Related Functions
(Note that these are not member functions.)
template<typename PH >
bool	poly_hull_assign_if_exact (PH &p, const PH &q)
	If the poly-hull of `p` and `q` is exact it is assigned to `p` and `true` is returned, otherwise `false` is returned.
template<typename PH >
bool	poly_hull_assign_if_exact (PH &p, const PH &q)
bool	operator== (const Polyhedron &x, const Polyhedron &y)
bool	operator== (const Polyhedron &x, const Polyhedron &y)
std::ostream &	operator<< (std::ostream &s, const Polyhedron &ph)
	Output operator.
void	swap (Polyhedron &x, Polyhedron &y)
	Swaps `x` with `y`.
bool	operator== (const Polyhedron &x, const Polyhedron &y)
	Returns `true` if and only if `x` and `y` are the same polyhedron.
bool	operator!= (const Polyhedron &x, const Polyhedron &y)
	Returns `true` if and only if `x` and `y` are different polyhedra.
void	swap (Polyhedron &x, Polyhedron &y)
bool	operator!= (const Polyhedron &x, const Polyhedron &y)
bool	operator== (const Polyhedron &x, const Polyhedron &y)
std::ostream &	operator<< (std::ostream &s, const Polyhedron &ph)
Exception Throwers
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 Constraint_System &cs) const
void	throw_topology_incompatible (const char method, const char gs_name, const Generator_System &gs) const
void	throw_dimension_incompatible (const char method, const char other_name, dimension_type other_dim) 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 le_name, const Linear_Expression &le) const
void	throw_dimension_incompatible (const char method, const char c_name, const Constraint &c) const
void	throw_dimension_incompatible (const char method, const char g_name, const Generator &g) const
void	throw_dimension_incompatible (const char method, const char cg_name, const Congruence &cg) const
void	throw_dimension_incompatible (const char method, const char cs_name, const Constraint_System &cs) const
void	throw_dimension_incompatible (const char method, const char gs_name, const Generator_System &gs) const
void	throw_dimension_incompatible (const char method, const char cgs_name, const Congruence_System &cgs) const
template<typename C >
void	throw_dimension_incompatible (const char method, const char lf_name, const Linear_Form< C > &lf) const
void	throw_dimension_incompatible (const char method, const char var_name, Variable var) 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
static dimension_type	check_space_dimension_overflow (dimension_type dim, dimension_type max, const Topology topol, const char method, const char reason)
static dimension_type	check_space_dimension_overflow (dimension_type dim, const Topology topol, const char method, const char reason)
template<typename Object >
static Object &	check_obj_space_dimension_overflow (Object &input, Topology topol, const char method, const char reason)

Detailed Description

The base class for convex polyhedra.

An object of the class Polyhedron represents a convex polyhedron in the vector space $\Rset^n$ .

A polyhedron can be specified as either a finite system of constraints or a finite system of generators (see Section Representations of Convex Polyhedra) and it is always possible to obtain either representation. That is, if we know the system of constraints, we can obtain from this the system of generators that define the same polyhedron and vice versa. These systems can contain redundant members: in this case we say that they are not in the minimal form.

Two key attributes of any polyhedron are its topological kind (recording whether it is a C_Polyhedron or an NNC_Polyhedron object) and its space dimension (the dimension $n \in \Nset$ of the enclosing vector space):

all polyhedra, the empty ones included, are endowed with a specific topology and space dimension;
most operations working on a polyhedron and another object (i.e., another polyhedron, a constraint or generator, a set of variables, etc.) will throw an exception if the polyhedron and the object are not both topology-compatible and dimension-compatible (see Section Representations of Convex Polyhedra);
the topology of a polyhedron cannot be changed; rather, there are constructors for each of the two derived classes that will build a new polyhedron with the topology of that class from another polyhedron from either class and any topology;
the only ways in which the space dimension of a polyhedron can be changed are:
- explicit calls to operators provided for that purpose;
- standard copy, assignment and swap operators.

Note that four different polyhedra can be defined on the zero-dimension space: the empty polyhedron, either closed or NNC, and the universe polyhedron $R^0$ , again either closed or NNC.

In all the examples it is assumed that variables x and y are defined (where they are used) as follows:

  Variable x(0);
  Variable y(1);

Example 1

The following code builds a polyhedron corresponding to a square in $\Rset^2$ , given as a system of constraints:

  Constraint_System cs;
  cs.insert(x >= 0);
  cs.insert(x <= 3);
  cs.insert(y >= 0);
  cs.insert(y <= 3);
  C_Polyhedron ph(cs);

The following code builds the same polyhedron as above, but starting from a system of generators specifying the four vertices of the square:

  Generator_System gs;
  gs.insert(point(0*x + 0*y));
  gs.insert(point(0*x + 3*y));
  gs.insert(point(3*x + 0*y));
  gs.insert(point(3*x + 3*y));
  C_Polyhedron ph(gs);

Example 2

The following code builds an unbounded polyhedron corresponding to a half-strip in $\Rset^2$ , given as a system of constraints:

  Constraint_System cs;
  cs.insert(x >= 0);
  cs.insert(x - y <= 0);
  cs.insert(x - y + 1 >= 0);
  C_Polyhedron ph(cs);

The following code builds the same polyhedron as above, but starting from the system of generators specifying the two vertices of the polyhedron and one ray:

  Generator_System gs;
  gs.insert(point(0*x + 0*y));
  gs.insert(point(0*x + y));
  gs.insert(ray(x - y));
  C_Polyhedron ph(gs);

Example 3

The following code builds the polyhedron corresponding to a half-plane by adding a single constraint to the universe polyhedron in $\Rset^2$ :

  C_Polyhedron ph(2);
  ph.add_constraint(y >= 0);

The following code builds the same polyhedron as above, but starting from the empty polyhedron in the space $\Rset^2$ and inserting the appropriate generators (a point, a ray and a line).

  C_Polyhedron ph(2, EMPTY);
  ph.add_generator(point(0*x + 0*y));
  ph.add_generator(ray(y));
  ph.add_generator(line(x));

Note that, although the above polyhedron has no vertices, we must add one point, because otherwise the result of the Minkowski's sum would be an empty polyhedron. To avoid subtle errors related to the minimization process, it is required that the first generator inserted in an empty polyhedron is a point (otherwise, an exception is thrown).

Example 4

The following code shows the use of the function add_space_dimensions_and_embed:

  C_Polyhedron ph(1);
  ph.add_constraint(x == 2);
  ph.add_space_dimensions_and_embed(1);

We build the universe polyhedron in the 1-dimension space $\Rset$ . Then we add a single equality constraint, thus obtaining the polyhedron corresponding to the singleton set $\{ 2 \} \sseq \Rset$ . After the last line of code, the resulting polyhedron is

$\bigl\{\, (2, y)^\transpose \in \Rset^2 \bigm| y \in \Rset \,\bigr\}.$

Example 5

The following code shows the use of the function add_space_dimensions_and_project:

  C_Polyhedron ph(1);
  ph.add_constraint(x == 2);
  ph.add_space_dimensions_and_project(1);

The first two lines of code are the same as in Example 4 for add_space_dimensions_and_embed. After the last line of code, the resulting polyhedron is the singleton set $\bigl\{ (2, 0)^\transpose \bigr\} \sseq \Rset^2$ .

Example 6

The following code shows the use of the function affine_image:

  C_Polyhedron ph(2, EMPTY);
  ph.add_generator(point(0*x + 0*y));
  ph.add_generator(point(0*x + 3*y));
  ph.add_generator(point(3*x + 0*y));
  ph.add_generator(point(3*x + 3*y));
  Linear_Expression expr = x + 4;
  ph.affine_image(x, expr);

In this example the starting polyhedron is a square in $\Rset^2$ , the considered variable is $x$

and the affine expression is $x+4$

. The resulting polyhedron is the same square translated to the right. Moreover, if the affine transformation for the same variable x is $x+y$

:

  Linear_Expression expr = x + y;

the resulting polyhedron is a parallelogram with the height equal to the side of the square and the oblique sides parallel to the line $x-y$

. Instead, if we do not use an invertible transformation for the same variable; for example, the affine expression $y$

:

  Linear_Expression expr = y;

the resulting polyhedron is a diagonal of the square.

Example 7

The following code shows the use of the function affine_preimage:

  C_Polyhedron ph(2);
  ph.add_constraint(x >= 0);
  ph.add_constraint(x <= 3);
  ph.add_constraint(y >= 0);
  ph.add_constraint(y <= 3);
  Linear_Expression expr = x + 4;
  ph.affine_preimage(x, expr);

In this example the starting polyhedron, var and the affine expression and the denominator are the same as in Example 6, while the resulting polyhedron is again the same square, but translated to the left. Moreover, if the affine transformation for x is $x+y$

  Linear_Expression expr = x + y;

the resulting polyhedron is a parallelogram with the height equal to the side of the square and the oblique sides parallel to the line $x+y$

. Instead, if we do not use an invertible transformation for the same variable x, for example, the affine expression $y$

:

  Linear_Expression expr = y;

the resulting polyhedron is a line that corresponds to the $y$

axis.

Example 8

For this example we use also the variables:

  Variable z(2);
  Variable w(3);

The following code shows the use of the function remove_space_dimensions:

  Generator_System gs;
  gs.insert(point(3*x + y +0*z + 2*w));
  C_Polyhedron ph(gs);
  Variables_Set vars;
  vars.insert(y);
  vars.insert(z);
  ph.remove_space_dimensions(vars);

The starting polyhedron is the singleton set $\bigl\{ (3, 1, 0, 2)^\transpose \bigr\} \sseq \Rset^4$ , while the resulting polyhedron is $\bigl\{ (3, 2)^\transpose \bigr\} \sseq \Rset^2$ . Be careful when removing space dimensions incrementally: since dimensions are automatically renamed after each application of the remove_space_dimensions operator, unexpected results can be obtained. For instance, by using the following code we would obtain a different result:

  set<Variable> vars1;
  vars1.insert(y);
  ph.remove_space_dimensions(vars1);
  set<Variable> vars2;
  vars2.insert(z);
  ph.remove_space_dimensions(vars2);

In this case, the result is the polyhedron $\bigl\{(3, 0)^\transpose \bigr\} \sseq \Rset^2$ : when removing the set of dimensions vars2 we are actually removing variable $w$

of the original polyhedron. For the same reason, the operator remove_space_dimensions is not idempotent: removing twice the same non-empty set of dimensions is never the same as removing them just once.

Definition at line 369 of file Polyhedron.defs.hh.

Member Typedef Documentation

typedef Coefficient Parma_Polyhedra_Library::Polyhedron::coefficient_type

The numeric type of coefficients.

Definition at line 372 of file Polyhedron.defs.hh.

Member Enumeration Documentation

enum Parma_Polyhedra_Library::Polyhedron::Three_Valued_Boolean

private

Enumerator:

TVB_TRUE
TVB_FALSE
TVB_DONT_KNOW

Definition at line 2395 of file Polyhedron.defs.hh.

                            {
    TVB_TRUE,
    TVB_FALSE,
    TVB_DONT_KNOW
  };

Constructor & Destructor Documentation

Parma_Polyhedra_Library::Polyhedron::Polyhedron	(	Topology	topol,
		dimension_type	num_dimensions,
		Degenerate_Element	kind
	)

protected

Builds a polyhedron having the specified properties.

Parameters:

topol	The topology of the polyhedron;
num_dimensions	The number of dimensions of the vector space enclosing the polyhedron;
kind	Specifies whether the universe or the empty polyhedron has to be built.

Definition at line 51 of file Polyhedron_nonpublic.cc.

References Parma_Polyhedra_Library::Constraint_System::add_low_level_constraints(), Parma_Polyhedra_Library::Constraint_System::adjust_topology_and_space_dimension(), con_sys, Parma_Polyhedra_Library::EMPTY, max_space_dimension(), OK(), PPL_ASSERT, PPL_ASSERT_HEAVY, set_constraints_minimized(), Parma_Polyhedra_Library::Polyhedron::Status::set_empty(), space_dim, and status.

Referenced by affine_form_image().

  : con_sys(topol),
    gen_sys(topol),
    sat_c(),
    sat_g() {
  // Protecting against space dimension overflow is up to the caller.
  PPL_ASSERT(num_dimensions <= max_space_dimension());

  if (kind == EMPTY)
    status.set_empty();
  else if (num_dimensions > 0) {
    con_sys.add_low_level_constraints();
    con_sys.adjust_topology_and_space_dimension(topol, num_dimensions);
    set_constraints_minimized();
  }
  space_dim = num_dimensions;
  PPL_ASSERT_HEAVY(OK());
}

Parma_Polyhedra_Library::Polyhedron::Polyhedron	(	const Polyhedron &	y,
		Complexity_Class	complexity = `ANY_COMPLEXITY`
	)

protected

Ordinary copy constructor.

The complexity argument is ignored.

Definition at line 72 of file Polyhedron_nonpublic.cc.

References Parma_Polyhedra_Library::Linear_System::assign_with_pending(), con_sys, constraints_are_up_to_date(), gen_sys, generators_are_up_to_date(), PPL_ASSERT, sat_c, sat_c_is_up_to_date(), sat_g, sat_g_is_up_to_date(), and topology().

  : con_sys(y.topology()),
    gen_sys(y.topology()),
    status(y.status),
    space_dim(y.space_dim) {
  // Being a protected method, we simply assert that topologies do match.
  PPL_ASSERT(topology() == y.topology());
  if (y.constraints_are_up_to_date())
    con_sys.assign_with_pending(y.con_sys);
  if (y.generators_are_up_to_date())
    gen_sys.assign_with_pending(y.gen_sys);
  if (y.sat_c_is_up_to_date())
    sat_c = y.sat_c;
  if (y.sat_g_is_up_to_date())
    sat_g = y.sat_g;
}

Parma_Polyhedra_Library::Polyhedron::Polyhedron	(	Topology	topol,
		const Constraint_System &	cs
	)

protected

Builds a polyhedron from a system of constraints.

The polyhedron inherits the space dimension of the constraint system.

Parameters:

topol	The topology of the polyhedron;
cs	The system of constraints defining the polyhedron.

Exceptions:

std::invalid_argument Thrown if the topology of cs is incompatible with topol.

Definition at line 89 of file Polyhedron_nonpublic.cc.

References Parma_Polyhedra_Library::Constraint_System::add_low_level_constraints(), con_sys, max_space_dimension(), Parma_Polyhedra_Library::NECESSARILY_CLOSED, Parma_Polyhedra_Library::Linear_System::num_pending_rows(), OK(), PPL_ASSERT, PPL_ASSERT_HEAVY, set_constraints_up_to_date(), set_empty(), Parma_Polyhedra_Library::Linear_System::set_sorted(), space_dim, Parma_Polyhedra_Library::Constraint_System::space_dimension(), swap(), throw_topology_incompatible(), and Parma_Polyhedra_Library::Linear_System::unset_pending_rows().

  : con_sys(topol),
    gen_sys(topol),
    sat_c(),
    sat_g() {
  // Protecting against space dimension overflow is up to the caller.
  PPL_ASSERT(cs.space_dimension() <= max_space_dimension());

  // TODO: this implementation is just an executable specification.
  Constraint_System cs_copy = cs;

  // Try to adapt `cs_copy' to the required topology.
  const dimension_type cs_copy_space_dim = cs_copy.space_dimension();
  if (!cs_copy.adjust_topology_and_space_dimension(topol, cs_copy_space_dim))
    throw_topology_incompatible((topol == NECESSARILY_CLOSED)
                                ? "C_Polyhedron(cs)"
                                : "NNC_Polyhedron(cs)", "cs", cs_copy);

  // Set the space dimension.
  space_dim = cs_copy_space_dim;

  if (space_dim > 0) {
    // Stealing the rows from `cs_copy'.
    using std::swap;
    swap(con_sys, cs_copy);
    if (con_sys.num_pending_rows() > 0) {
      // Even though `cs_copy' has pending constraints, since the
      // generators of the polyhedron are not up-to-date, the
      // polyhedron cannot have pending constraints. By integrating
      // the pending part of `con_sys' we may loose sortedness.
      con_sys.unset_pending_rows();
      con_sys.set_sorted(false);
    }
    con_sys.add_low_level_constraints();
    set_constraints_up_to_date();
  }
  else {
    // Here `space_dim == 0'.
    if (cs_copy.num_columns() > 0)
      // See if an inconsistent constraint has been passed.
      for (dimension_type i = cs_copy.num_rows(); i-- > 0; )
        if (cs_copy[i].is_inconsistent()) {
          // Inconsistent constraint found: the polyhedron is empty.
          set_empty();
          break;
        }
  }
  PPL_ASSERT_HEAVY(OK());
}

Parma_Polyhedra_Library::Polyhedron::Polyhedron	(	Topology	topol,
		Constraint_System &	cs,
		Recycle_Input	dummy
	)

protected

Builds a polyhedron recycling a system of constraints.

The polyhedron inherits the space dimension of the constraint system.

Parameters:

topol	The topology of the polyhedron;
cs	The system of constraints defining the polyhedron. It is not declared `const` because its data-structures may be recycled to build the polyhedron.
dummy	A dummy tag to syntactically differentiate this one from the other constructors.

Exceptions:

std::invalid_argument Thrown if the topology of cs is incompatible with topol.

Definition at line 139 of file Polyhedron_nonpublic.cc.

  : con_sys(topol),
    gen_sys(topol),
    sat_c(),
    sat_g() {
  // Protecting against space dimension overflow is up to the caller.
  PPL_ASSERT(cs.space_dimension() <= max_space_dimension());

  // Try to adapt `cs' to the required topology.
  const dimension_type cs_space_dim = cs.space_dimension();
  if (!cs.adjust_topology_and_space_dimension(topol, cs_space_dim))
    throw_topology_incompatible((topol == NECESSARILY_CLOSED)
                                ? "C_Polyhedron(cs, recycle)"
                                : "NNC_Polyhedron(cs, recycle)", "cs", cs);

  // Set the space dimension.
  space_dim = cs_space_dim;

  if (space_dim > 0) {
    // Stealing the rows from `cs'.
    using std::swap;
    swap(con_sys, cs);
    if (con_sys.num_pending_rows() > 0) {
      // Even though `cs' has pending constraints, since the generators
      // of the polyhedron are not up-to-date, the polyhedron cannot
      // have pending constraints. By integrating the pending part
      // of `con_sys' we may loose sortedness.
      con_sys.unset_pending_rows();
      con_sys.set_sorted(false);
    }
    con_sys.add_low_level_constraints();
    set_constraints_up_to_date();
  }
  else {
    // Here `space_dim == 0'.
    if (cs.num_columns() > 0)
      // See if an inconsistent constraint has been passed.
      for (dimension_type i = cs.num_rows(); i-- > 0; )
        if (cs[i].is_inconsistent()) {
          // Inconsistent constraint found: the polyhedron is empty.
          set_empty();
          break;
        }
  }
  PPL_ASSERT_HEAVY(OK());
}

Parma_Polyhedra_Library::Polyhedron::Polyhedron	(	Topology	topol,
		const Generator_System &	gs
	)

protected

Builds a polyhedron from a system of generators.

The polyhedron inherits the space dimension of the generator system.

Parameters:

topol	The topology of the polyhedron;
gs	The system of generators defining the polyhedron.

Exceptions:

std::invalid_argument Thrown if the topology of gs is incompatible with topol, or if the system of generators is not empty but has no points.

Definition at line 188 of file Polyhedron_nonpublic.cc.

  : con_sys(topol),
    gen_sys(topol),
    sat_c(),
    sat_g() {
  // Protecting against space dimension overflow is up to the caller.
  PPL_ASSERT(gs.space_dimension() <= max_space_dimension());

  // An empty set of generators defines the empty polyhedron.
  if (gs.has_no_rows()) {
    space_dim = gs.space_dimension();
    status.set_empty();
    PPL_ASSERT_HEAVY(OK());
    return;
  }

  // Non-empty valid generator systems have a supporting point, at least.
  if (!gs.has_points())
    throw_invalid_generators((topol == NECESSARILY_CLOSED)
                             ? "C_Polyhedron(gs)"
                             : "NNC_Polyhedron(gs)", "gs");

  // TODO: this implementation is just an executable specification.
  Generator_System gs_copy = gs;

  const dimension_type gs_copy_space_dim = gs_copy.space_dimension();
  // Try to adapt `gs_copy' to the required topology.
  if (!gs_copy.adjust_topology_and_space_dimension(topol, gs_copy_space_dim))
    throw_topology_incompatible((topol == NECESSARILY_CLOSED)
                                ? "C_Polyhedron(gs)"
                                : "NNC_Polyhedron(gs)", "gs", gs_copy);

  if (gs_copy_space_dim > 0) {
    // Stealing the rows from `gs_copy'.
    using std::swap;
    swap(gen_sys, gs_copy);
    // In a generator system describing a NNC polyhedron,
    // for each point we must also have the corresponding closure point.
    if (topol == NOT_NECESSARILY_CLOSED)
      gen_sys.add_corresponding_closure_points();
    if (gen_sys.num_pending_rows() > 0) {
      // Even though `gs_copy' has pending generators, since the
      // constraints of the polyhedron are not up-to-date, the
      // polyhedron cannot have pending generators. By integrating the
      // pending part of `gen_sys' we may loose sortedness.
      gen_sys.unset_pending_rows();
      gen_sys.set_sorted(false);
    }
    // Generators are now up-to-date.
    set_generators_up_to_date();

    // Set the space dimension.
    space_dim = gs_copy_space_dim;
    PPL_ASSERT_HEAVY(OK());
    return;
  }

  // Here `gs_copy.num_rows > 0' and `gs_copy_space_dim == 0':
  // we already checked for both the topology-compatibility
  // and the supporting point.
  space_dim = 0;
  PPL_ASSERT_HEAVY(OK());
}

Parma_Polyhedra_Library::Polyhedron::Polyhedron	(	Topology	topol,
		Generator_System &	gs,
		Recycle_Input	dummy
	)

protected

Builds a polyhedron recycling a system of generators.

The polyhedron inherits the space dimension of the generator system.

Parameters:

topol	The topology of the polyhedron;
gs	The system of generators defining the polyhedron. It is not declared `const` because its data-structures may be recycled to build the polyhedron.
dummy	A dummy tag to syntactically differentiate this one from the other constructors.

Exceptions:

std::invalid_argument Thrown if the topology of gs is incompatible with topol, or if the system of generators is not empty but has no points.

Definition at line 252 of file Polyhedron_nonpublic.cc.

  : con_sys(topol),
    gen_sys(topol),
    sat_c(),
    sat_g() {
  // Protecting against space dimension overflow is up to the caller.
  PPL_ASSERT(gs.space_dimension() <= max_space_dimension());

  // An empty set of generators defines the empty polyhedron.
  if (gs.has_no_rows()) {
    space_dim = gs.space_dimension();
    status.set_empty();
    PPL_ASSERT_HEAVY(OK());
    return;
  }

  // Non-empty valid generator systems have a supporting point, at least.
  if (!gs.has_points())
    throw_invalid_generators((topol == NECESSARILY_CLOSED)
                             ? "C_Polyhedron(gs, recycle)"
                             : "NNC_Polyhedron(gs, recycle)", "gs");

  const dimension_type gs_space_dim = gs.space_dimension();
  // Try to adapt `gs' to the required topology.
  if (!gs.adjust_topology_and_space_dimension(topol, gs_space_dim))
    throw_topology_incompatible((topol == NECESSARILY_CLOSED)
                                ? "C_Polyhedron(gs, recycle)"
                                : "NNC_Polyhedron(gs, recycle)", "gs", gs);

  if (gs_space_dim > 0) {
    // Stealing the rows from `gs'.
    using std::swap;
    swap(gen_sys, gs);
    // In a generator system describing a NNC polyhedron,
    // for each point we must also have the corresponding closure point.
    if (topol == NOT_NECESSARILY_CLOSED)
      gen_sys.add_corresponding_closure_points();
    if (gen_sys.num_pending_rows() > 0) {
      // Even though `gs' has pending generators, since the constraints
      // of the polyhedron are not up-to-date, the polyhedron cannot
      // have pending generators. By integrating the pending part
      // of `gen_sys' we may loose sortedness.
      gen_sys.unset_pending_rows();
      gen_sys.set_sorted(false);
    }
    // Generators are now up-to-date.
    set_generators_up_to_date();

    // Set the space dimension.
    space_dim = gs_space_dim;
    PPL_ASSERT_HEAVY(OK());
    return;
  }

  // Here `gs.num_rows > 0' and `gs_space_dim == 0':
  // we already checked for both the topology-compatibility
  // and the supporting point.
  space_dim = 0;
  PPL_ASSERT_HEAVY(OK());
}

template<typename Interval >

Parma_Polyhedra_Library::Polyhedron::Polyhedron	(	Topology	topol,
		const Box< Interval > &	box,
		Complexity_Class	complexity = `ANY_COMPLEXITY`
	)

protected

Builds a polyhedron from a box.

This will use an algorithm whose complexity is polynomial and build the smallest polyhedron with topology topol containing box.

Parameters:

topol	The topology of the polyhedron;
box	The box representing the polyhedron to be built;
complexity	This argument is ignored.

Definition at line 39 of file Polyhedron.templates.hh.

  : con_sys(topol),
    gen_sys(topol),
    sat_c(),
    sat_g() {
  // Initialize the space dimension as indicated by the box.
  space_dim = box.space_dimension();

  // Check for emptiness.
  if (box.is_empty()) {
    set_empty();
    return;
  }

  // Zero-dim universe polyhedron.
  if (space_dim == 0) {
    set_zero_dim_univ();
    return;
  }

  // Insert a dummy constraint of the highest dimension to avoid the
  // need of resizing the matrix of constraints later;
  // this constraint will be removed at the end.
  con_sys.insert(Variable(space_dim - 1) >= 0);

  PPL_DIRTY_TEMP_COEFFICIENT(l_n);
  PPL_DIRTY_TEMP_COEFFICIENT(l_d);
  PPL_DIRTY_TEMP_COEFFICIENT(u_n);
  PPL_DIRTY_TEMP_COEFFICIENT(u_d);

  if (topol == NECESSARILY_CLOSED) {
    for (dimension_type k = space_dim; k-- > 0; ) {
      const Variable v_k = Variable(k);
      // See if we have a valid lower bound.
      bool l_closed = false;
      bool l_bounded = box.has_lower_bound(v_k, l_n, l_d, l_closed);
      // See if we have a valid upper bound.
      bool u_closed = false;
      bool u_bounded = box.has_upper_bound(v_k, u_n, u_d, u_closed);

      // See if we have an implicit equality constraint.
      if (l_bounded && u_bounded
          && l_closed && u_closed
          && l_n == u_n && l_d == u_d) {
        // Add the constraint `l_d*v_k == l_n'.
        con_sys.insert(l_d * v_k == l_n);
      }
      else {
        if (l_bounded)
          // Add the constraint `l_d*v_k >= l_n'.
          con_sys.insert(l_d * v_k >= l_n);
        if (u_bounded)
          // Add the constraint `u_d*v_k <= u_n'.
          con_sys.insert(u_d * v_k <= u_n);
      }
    }
  }
  else {
    // topol == NOT_NECESSARILY_CLOSED
    for (dimension_type k = space_dim; k-- > 0; ) {
      const Variable v_k = Variable(k);
      // See if we have a valid lower bound.
      bool l_closed = false;
      bool l_bounded = box.has_lower_bound(v_k, l_n, l_d, l_closed);
      // See if we have a valid upper bound.
      bool u_closed = false;
      bool u_bounded = box.has_upper_bound(v_k, u_n, u_d, u_closed);

      // See if we have an implicit equality constraint.
      if (l_bounded && u_bounded
          && l_closed && u_closed
          && l_n == u_n && l_d == u_d) {
        // Add the constraint `l_d*v_k == l_n'.
        con_sys.insert(l_d * v_k == l_n);
      }
      else {
        // Check if a lower bound constraint is required.
        if (l_bounded) {
          if (l_closed)
            // Add the constraint `l_d*v_k >= l_n'.
            con_sys.insert(l_d * v_k >= l_n);
          else
            // Add the constraint `l_d*v_k > l_n'.
            con_sys.insert(l_d * v_k > l_n);
        }
        // Check if an upper bound constraint is required.
        if (u_bounded) {
          if (u_closed)
            // Add the constraint `u_d*v_k <= u_n'.
            con_sys.insert(u_d * v_k <= u_n);
          else
            // Add the constraint `u_d*v_k < u_n'.
            con_sys.insert(u_d * v_k < u_n);
        }
      }
    }
  }

  // Adding the low-level constraints.
  con_sys.add_low_level_constraints();
  // Now removing the dummy constraint inserted before.
  dimension_type n_rows = con_sys.num_rows() - 1;
  con_sys[0].m_swap(con_sys[n_rows]);
  con_sys.set_sorted(false);
  // NOTE: here there are no pending constraints.
  con_sys.set_index_first_pending_row(n_rows);
  con_sys.remove_trailing_rows(1);

  // Constraints are up-to-date.
  set_constraints_up_to_date();
  PPL_ASSERT_HEAVY(OK());
}

Parma_Polyhedra_Library::Polyhedron::~Polyhedron ( )

inline

Destructor.

Definition at line 96 of file Polyhedron.inlines.hh.

{
}

Member Function Documentation

bool Parma_Polyhedra_Library::Polyhedron::add_and_minimize	(	bool	con_to_gen,
		Linear_System &	source1,
		Linear_System &	dest,
		Bit_Matrix &	sat,
		const Linear_System &	source2
	)

staticprivate

Adds given constraints and builds minimized corresponding generators or vice versa.

Returns:: true if the obtained polyhedron is empty, false otherwise.

Parameters:

con_to_gen	`true` if `source1` and `source2` are system of constraints, `false` otherwise;
source1	The first element of the given DD pair;
dest	The second element of the given DD pair;
sat	The saturation matrix that bind `source1` to `dest`;
source2	The new system of generators or constraints.

It is assumed that source1 and source2 are sorted and have no pending rows. It is also assumed that dest has no pending rows. On entry, the rows of sat are indexed by the rows of dest and its columns are indexed by the rows of source1. On exit, the rows of sat are indexed by the rows of dest and its columns are indexed by the rows of the system obtained by merging source1 and source2.

Let us suppose we want to add some constraints to a given system of constraints source1. This method, given a minimized double description pair (source1, dest) and a system of new constraints source2, modifies source1 by adding to it the constraints of source2 that are not in source1. Then, by invoking add_and_minimize(bool, Linear_System&, Linear_System&, Bit_Matrix&), processes the added constraints obtaining a new DD pair.

This method treats also the dual case, i.e., adding new generators to a previous system of generators. In this case source1 contains the old generators, source2 the new ones and dest is the system of constraints in the given minimized DD pair.

Since source2 contains the constraints (or the generators) that will be added to source1, it is constant: it will not be modified.

Definition at line 236 of file minimize.cc.

References Parma_Polyhedra_Library::Linear_System::add_pending_row(), Parma_Polyhedra_Library::cmp(), Parma_Polyhedra_Library::Dense_Matrix::has_no_rows(), Parma_Polyhedra_Library::Linear_System::is_sorted(), Parma_Polyhedra_Library::Dense_Matrix::num_columns(), Parma_Polyhedra_Library::Linear_System::num_pending_rows(), Parma_Polyhedra_Library::Dense_Matrix::num_rows(), and PPL_ASSERT.

                                                                {
  // `source1' and `source2' cannot be empty.
  PPL_ASSERT(!source1.has_no_rows() && !source2.has_no_rows());
  // `source1' and `source2' must have the same number of columns
  // to be merged.
  PPL_ASSERT(source1.num_columns() == source2.num_columns());
  // `source1' and `source2' are fully sorted.
  PPL_ASSERT(source1.is_sorted() && source1.num_pending_rows() == 0);
  PPL_ASSERT(source2.is_sorted() && source2.num_pending_rows() == 0);
  PPL_ASSERT(dest.num_pending_rows() == 0);

  const dimension_type old_source1_num_rows = source1.num_rows();
  // `k1' and `k2' run through the rows of `source1' and `source2', resp.
  dimension_type k1 = 0;
  dimension_type k2 = 0;
  dimension_type source2_num_rows = source2.num_rows();
  while (k1 < old_source1_num_rows && k2 < source2_num_rows) {
    // Add to `source1' the constraints from `source2', as pending rows.
    // We exploit the property that initially both `source1' and `source2'
    // are sorted and index `k1' only scans the non-pending rows of `source1',
    // so that it is not influenced by the pending rows appended to it.
    // This way no duplicate (i.e., trivially redundant) constraint
    // is introduced in `source1'.
    const int cmp = compare(source1[k1], source2[k2]);
    if (cmp == 0) {
      // We found the same row: there is no need to add `source2[k2]'.
      ++k2;
      // By sortedness, since `k1 < old_source1_num_rows',
      // we can increment index `k1' too.
      ++k1;
    }
    else if (cmp < 0)
      // By sortedness, we can increment `k1'.
      ++k1;
    else {
      // Here `cmp > 0'.
      // By sortedness, `source2[k2]' cannot be in `source1'.
      // We add it as a pending row of `source1' (sortedness unaffected).
      source1.add_pending_row(source2[k2]);
      // We can increment `k2'.
      ++k2;
    }
  }
  // Have we scanned all the rows in `source2'?
  if (k2 < source2_num_rows)
    // By sortedness, all the rows in `source2' having indexes
    // greater than or equal to `k2' were not in `source1'.
    // We add them as pending rows of 'source1' (sortedness not affected).
    for ( ; k2 < source2_num_rows; ++k2)
      source1.add_pending_row(source2[k2]);

  if (source1.num_pending_rows() == 0)
    // No row was appended to `source1', because all the constraints
    // in `source2' were already in `source1'.
    // There is nothing left to do ...
    return false;

  return add_and_minimize(con_to_gen, source1, dest, sat);
}

bool Parma_Polyhedra_Library::Polyhedron::add_and_minimize	(	bool	con_to_gen,
		Linear_System &	source,
		Linear_System &	dest,
		Bit_Matrix &	sat
	)

staticprivate

Adds given constraints and builds minimized corresponding generators or vice versa. The given constraints are in source.

Returns:: true if the obtained polyhedron is empty, false otherwise.

Parameters:

con_to_gen	`true` if `source` is a system of constraints, `false` otherwise;
source	The first element of the given DD pair. It also contains the pending rows to be processed;
dest	The second element of the given DD pair. It cannot have pending rows;
sat	The saturation matrix that bind the upper part of `source` to `dest`.

On entry, the rows of sat are indexed by the rows of dest and its columns are indexed by the non-pending rows of source. On exit, the rows of sat are indexed by the rows of dest and its columns are indexed by the rows of source.

Let us suppose that source is a system of constraints. This method assumes that the non-pending part of source and system dest form a double description pair in minimal form and will build a new DD pair in minimal form by processing the pending constraints in source. To this end, it will call conversion()) and simplify.

This method treats also the dual case, i.e., processing pending generators. In this case source contains generators and dest is the system of constraints corresponding to the non-pending part of source.

Definition at line 337 of file minimize.cc.

                                                   {
  PPL_ASSERT(source.num_pending_rows() > 0);
  PPL_ASSERT(source.num_columns() == dest.num_columns());
  PPL_ASSERT(source.is_sorted());

  // First, pad the saturation matrix with new columns (of zeroes)
  // to accommodate for the pending rows of `source'.
  sat.resize(dest.num_rows(), source.num_rows());

  // Incrementally compute the new system of generators.
  // Parameter `start' is set to the index of the first pending constraint.
  const dimension_type num_lines_or_equalities
    = conversion(source, source.first_pending_row(),
                 dest, sat,
                 dest.num_lines_or_equalities());

  // conversion() may have modified the number of rows in `dest'.
  const dimension_type dest_num_rows = dest.num_rows();

  // Checking if the generators in `dest' represent an empty polyhedron:
  // the polyhedron is empty if there are no points
  // (because rays, lines and closure points need a supporting point).
  // Points can be detected by looking at:
  // - the divisor, for necessarily closed polyhedra;
  // - the epsilon coordinate, for NNC polyhedra.
  const dimension_type checking_index
    = dest.is_necessarily_closed()
    ? 0
    : (dest.num_columns() - 1);
  dimension_type first_point;
  for (first_point = num_lines_or_equalities;
       first_point < dest_num_rows;
       ++first_point)
    if (dest[first_point][checking_index] > 0)
      break;

  if (first_point == dest_num_rows)
    if (con_to_gen)
      // No point has been found: the polyhedron is empty.
      return true;
    else {
      // Here `con_to_gen' is false: `dest' is a system of constraints.
      // In this case the condition `first_point == dest_num_rows'
      // actually means that all the constraints in `dest' have their
      // inhomogeneous term equal to 0.
      // This is an ILLEGAL situation, because it implies that
      // the constraint system `dest' lacks the positivity constraint
      // and no linear combination of the constraints in `dest'
      // can reintroduce the positivity constraint.
      PPL_UNREACHABLE;
      return false;
    }
  else {
    // A point has been found: the polyhedron is not empty.
    // Now invoking `simplify()' to remove all the redundant constraints
    // from the system `source'.
    // Since the saturation matrix `sat' returned by `conversion()'
    // has rows indexed by generators (the rows of `dest') and columns
    // indexed by constraints (the rows of `source'), we have to
    // transpose it to obtain the saturation matrix needed by `simplify()'.
    sat.transpose();
    simplify(source, sat);
    // Transposing back.
    sat.transpose();
    return false;
  }
}

void Parma_Polyhedra_Library::Polyhedron::add_congruence ( const Congruence & cg )

Adds a copy of congruence cg to *this, if cg can be exactly represented by a polyhedron.

Exceptions:

std::invalid_argument Thrown if *this and congruence cg are dimension-incompatible, of if cg is a proper congruence which is neither a tautology, nor a contradiction.

Definition at line 1268 of file Polyhedron_public.cc.

References c, Parma_Polyhedra_Library::Constraint::EQUALITY, Parma_Polyhedra_Library::Congruence::is_equality(), Parma_Polyhedra_Library::Congruence::is_inconsistent(), Parma_Polyhedra_Library::Congruence::is_proper_congruence(), Parma_Polyhedra_Library::Congruence::is_tautological(), Parma_Polyhedra_Library::Checked::le, Parma_Polyhedra_Library::NECESSARILY_CLOSED, PPL_ASSERT, Parma_Polyhedra_Library::Congruence::space_dimension(), and Parma_Polyhedra_Library::Linear_Row::strong_normalize().

                                                  {
  // Dimension-compatibility check:
  // the dimension of `cg' can not be greater than space_dim.
  if (space_dim < cg.space_dimension())
    throw_dimension_incompatible("add_congruence(cg)", "cg", cg);

  // Handle the case of proper congruences first.
  if (cg.is_proper_congruence()) {
    if (cg.is_tautological())
      return;
    if (cg.is_inconsistent()) {
      set_empty();
      return;
    }
    // Non-trivial and proper congruences are not allowed.
    throw_invalid_argument("add_congruence(cg)",
                           "cg is a non-trivial, proper congruence");
  }

  PPL_ASSERT(cg.is_equality());
  // Handle empty and 0-dim cases first.
  if (marked_empty())
    return;
  if (space_dim == 0) {
    if (cg.is_inconsistent())
      set_empty();
    return;
  }

  // Add the equality.
  Linear_Expression le(cg);
  Constraint c(le, Constraint::EQUALITY, NECESSARILY_CLOSED);
  // Enforce normalization.
  c.strong_normalize();
  refine_no_check(c);
}

void Parma_Polyhedra_Library::Polyhedron::add_congruences ( const Congruence_System & cgs )

Adds a copy of the congruences in cgs to *this, if all the congruences can be exactly represented by a polyhedron.

Parameters:

cgs	The congruences to be added.

Exceptions:

std::invalid_argument Thrown if *this and cgs are dimension-incompatible, of if there exists in cgs a proper congruence which is neither a tautology, nor a contradiction.

Definition at line 1645 of file Polyhedron_public.cc.

Referenced by add_recycled_congruences(), Parma_Polyhedra_Library::C_Polyhedron::C_Polyhedron(), and Parma_Polyhedra_Library::NNC_Polyhedron::NNC_Polyhedron().

                                                           {
  // Dimension-compatibility check.
  if (space_dim < cgs.space_dimension())
    throw_dimension_incompatible("add_congruences(cgs)", "cgs", cgs);

  Constraint_System cs;
  bool inserted = false;
  for (Congruence_System::const_iterator i = cgs.begin(),
         cgs_end = cgs.end(); i != cgs_end; ++i) {
    const Congruence& cg = *i;
    if (cg.is_equality()) {
      Linear_Expression le(cg);
      Constraint c(le, Constraint::EQUALITY, NECESSARILY_CLOSED);
      // Enforce normalization.
      c.strong_normalize();
      // TODO: Consider stealing the row in c when adding it to cs.
      cs.insert(c);
      inserted = true;
    }
    else {
      PPL_ASSERT(cg.is_proper_congruence());
      if (cg.is_inconsistent()) {
        set_empty();
        return;
      }
      if (!cg.is_tautological())
        throw_invalid_argument("add_congruences(cgs)",
                               "cgs has a non-trivial, proper congruence");
    }
  }
  // Only add cs if it contains something.
  if (inserted)
    add_recycled_constraints(cs);
}

void Parma_Polyhedra_Library::Polyhedron::add_constraint ( const Constraint & c )

Adds a copy of constraint c to the system of constraints of *this (without minimizing the result).

Parameters:

c	The constraint that will be added to the system of constraints of `*this`.

Exceptions:

std::invalid_argument Thrown if *this and constraint c are topology-incompatible or dimension-incompatible.

Definition at line 1244 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Constraint::is_inconsistent(), Parma_Polyhedra_Library::Constraint::is_strict_inequality(), Parma_Polyhedra_Library::Constraint::is_tautological(), and Parma_Polyhedra_Library::Constraint::space_dimension().

Referenced by Parma_Polyhedra_Library::Pointset_Powerset< PSET >::affine_dimension(), Parma_Polyhedra_Library::Implementation::Termination::all_affine_ranking_functions_PR(), Parma_Polyhedra_Library::Implementation::Termination::all_affine_ranking_functions_PR_original(), BHZ09_NNC_poly_hull_assign_if_exact(), Parma_Polyhedra_Library::C_Polyhedron::C_Polyhedron(), and Parma_Polyhedra_Library::Pointset_Powerset< PSET >::linear_partition_aux().

                                                 {
  // Topology-compatibility check.
  if (c.is_strict_inequality() && is_necessarily_closed()) {
    // Trivially true/false strict inequalities are legal.
    if (c.is_tautological())
      return;
    if (c.is_inconsistent()) {
      set_empty();
      return;
    }
    // Here c is a non-trivial strict inequality.
    throw_topology_incompatible("add_constraint(c)", "c", c);
  }

  // Dimension-compatibility check:
  // the dimension of `c' can not be greater than space_dim.
  if (space_dim < c.space_dimension())
    throw_dimension_incompatible("add_constraint(c)", "c", c);

  if (!marked_empty())
    refine_no_check(c);
}

void Parma_Polyhedra_Library::Polyhedron::add_constraints ( const Constraint_System & cs )

Adds a copy of the constraints in cs to the system of constraints of *this (without minimizing the result).

Parameters:

cs	Contains the constraints that will be added to the system of constraints of `*this`.

Exceptions:

std::invalid_argument Thrown if *this and cs are topology-incompatible or dimension-incompatible.

Definition at line 1536 of file Polyhedron_public.cc.

Referenced by Parma_Polyhedra_Library::C_Polyhedron::C_Polyhedron(), Parma_Polyhedra_Library::NNC_Polyhedron::NNC_Polyhedron(), and simplify_using_context_assign().

                                                          {
  // TODO: this is just an executable specification.
  Constraint_System cs_copy = cs;
  add_recycled_constraints(cs_copy);
}

void Parma_Polyhedra_Library::Polyhedron::add_generator ( const Generator & g )

Adds a copy of generator g to the system of generators of *this (without minimizing the result).

Exceptions:

std::invalid_argument Thrown if *this and generator g are topology-incompatible or dimension-incompatible, or if *this is an empty polyhedron and g is not a point.

Definition at line 1306 of file Polyhedron_public.cc.

Referenced by BHZ09_NNC_poly_hull_assign_if_exact().

                                               {
  // Topology-compatibility check.
  if (g.is_closure_point() && is_necessarily_closed())
    throw_topology_incompatible("add_generator(g)", "g", g);
  // Dimension-compatibility check:
  // the dimension of `g' can not be greater than space_dim.
  const dimension_type g_space_dim = g.space_dimension();
  if (space_dim < g_space_dim)
    throw_dimension_incompatible("add_generator(g)", "g", g);

  // Dealing with a zero-dimensional space polyhedron first.
  if (space_dim == 0) {
    // It is not possible to create 0-dim rays or lines.
    PPL_ASSERT(g.is_point() || g.is_closure_point());
    // Closure points can only be inserted in non-empty polyhedra.
    if (marked_empty()) {
      if (g.type() != Generator::POINT)
        throw_invalid_generator("add_generator(g)", "g");
      else
        set_zero_dim_univ();
    }
    PPL_ASSERT_HEAVY(OK());
    return;
  }

  if (marked_empty()
      || (has_pending_constraints() && !process_pending_constraints())
      || (!generators_are_up_to_date() && !update_generators())) {
    // Here the polyhedron is empty:
    // the specification says we can only insert a point.
    if (!g.is_point())
      throw_invalid_generator("add_generator(g)", "g");
    if (g.is_necessarily_closed() || !is_necessarily_closed()) {
      gen_sys.insert(g);
      // Since `gen_sys' was empty, after inserting `g' we have to resize
      // the system of generators to have the right dimension.
      gen_sys.adjust_topology_and_space_dimension(topology(), space_dim);
      if (!is_necessarily_closed()) {
        // In the NNC topology, each point has to be matched by
        // a corresponding closure point:
        // turn the just inserted point into the corresponding
        // (normalized) closure point.
        Generator& cp = gen_sys[gen_sys.num_rows() - 1];
        cp[space_dim + 1] = 0;
        cp.normalize();
        // Re-insert the point (which is already normalized).
        gen_sys.insert(g);
      }
    }
    else {
      // Note: here we have a _legal_ topology mismatch,
      // because `g' is NOT a closure point (it is a point!)
      // However, by barely invoking `gen_sys.insert(g)' we would
      // cause a change in the topology of `gen_sys', which is wrong.
      // Thus, we insert a "topology corrected" copy of `g'.
      const Linear_Expression nc_expr = Linear_Expression(g);
      gen_sys.insert(Generator::point(nc_expr, g.divisor()));
      // Since `gen_sys' was empty, after inserting `g' we have to resize
      // the system of generators to have the right dimension.
      gen_sys.adjust_topology_and_space_dimension(topology(), space_dim);
    }
    // No longer empty, generators up-to-date and minimized.
    clear_empty();
    set_generators_minimized();
  }
  else {
    PPL_ASSERT(generators_are_up_to_date());
    const bool has_pending = can_have_something_pending();
    if (g.is_necessarily_closed() || !is_necessarily_closed()) {
      // Since `gen_sys' is not empty, the topology and space dimension
      // of the inserted generator are automatically adjusted.
      if (has_pending)
        gen_sys.insert_pending(g);
      else
        gen_sys.insert(g);
      if (!is_necessarily_closed() && g.is_point()) {
        // In the NNC topology, each point has to be matched by
        // a corresponding closure point:
        // turn the just inserted point into the corresponding
        // (normalized) closure point.
        Generator& cp = gen_sys[gen_sys.num_rows() - 1];
        cp[space_dim + 1] = 0;
        cp.normalize();
        // Re-insert the point (which is already normalized).
        if (has_pending)
          gen_sys.insert_pending(g);
        else
          gen_sys.insert(g);
      }
    }
    else {
      PPL_ASSERT(!g.is_closure_point());
      // Note: here we have a _legal_ topology mismatch, because
      // `g' is NOT a closure point.
      // However, by barely invoking `gen_sys.insert(g)' we would
      // cause a change in the topology of `gen_sys', which is wrong.
      // Thus, we insert a "topology corrected" copy of `g'.
      const Linear_Expression nc_expr = Linear_Expression(g);
      switch (g.type()) {
      case Generator::LINE:
        if (has_pending)
          gen_sys.insert_pending(Generator::line(nc_expr));
        else
          gen_sys.insert(Generator::line(nc_expr));
        break;
      case Generator::RAY:
        if (has_pending)
          gen_sys.insert_pending(Generator::ray(nc_expr));
        else
          gen_sys.insert(Generator::ray(nc_expr));
        break;
      case Generator::POINT:
        if (has_pending)
          gen_sys.insert_pending(Generator::point(nc_expr, g.divisor()));
        else
          gen_sys.insert(Generator::point(nc_expr, g.divisor()));
        break;
      case Generator::CLOSURE_POINT:
        PPL_UNREACHABLE;
        break;
      }
    }

    if (has_pending)
      set_generators_pending();
    else {
      // After adding the new generator,
      // constraints are no longer up-to-date.
      clear_generators_minimized();
      clear_constraints_up_to_date();
    }
  }
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::add_generators ( const Generator_System & gs )

Adds a copy of the generators in gs to the system of generators of *this (without minimizing the result).

Parameters:

gs	Contains the generators that will be added to the system of generators of `*this`.

Exceptions:

std::invalid_argument Thrown if *this and gs are topology-incompatible or dimension-incompatible, or if *this is empty and the system of generators gs is not empty, but has no points.

Definition at line 1638 of file Polyhedron_public.cc.

                                                        {
  // TODO: this is just an executable specification.
  Generator_System gs_copy = gs;
  add_recycled_generators(gs_copy);
}

void Parma_Polyhedra_Library::Polyhedron::add_recycled_congruences ( Congruence_System & cgs )

inline

Adds the congruences in cgs to *this, if all the congruences can be exactly represented by a polyhedron.

Parameters:

cgs	The congruences to be added. Its elements may be recycled.

Exceptions:

std::invalid_argument Thrown if *this and cgs are dimension-incompatible, of if there exists in cgs a proper congruence which is neither a tautology, nor a contradiction

Warning:: The only assumption that can be made on cgs upon successful or exceptional return is that it can be safely destroyed.

Definition at line 373 of file Polyhedron.inlines.hh.

References add_congruences().

                                                           {
  add_congruences(cgs);
}

void Parma_Polyhedra_Library::Polyhedron::add_recycled_constraints ( Constraint_System & cs )

Adds the constraints in cs to the system of constraints of *this (without minimizing the result).

Parameters:

cs	The constraint system to be added to `*this`. The constraints in `cs` may be recycled.

Exceptions:

std::invalid_argument Thrown if *this and cs are topology-incompatible or dimension-incompatible.

Warning:: The only assumption that can be made on cs upon successful or exceptional return is that it can be safely destroyed.

Definition at line 1442 of file Polyhedron_public.cc.

Referenced by BHRZ03_combining_constraints(), H79_widening_assign(), limited_BHRZ03_extrapolation_assign(), limited_H79_extrapolation_assign(), and simplify_using_context_assign().

                                                             {
  // Topology compatibility check.
  if (is_necessarily_closed() && cs.has_strict_inequalities()) {
    // We check if _all_ strict inequalities in cs are trivially false.
    // (The iterators already filter away trivially true constraints.)
    for (Constraint_System::const_iterator i = cs.begin(),
           i_end = cs.end(); i != i_end; ++i) {
      if (i->is_strict_inequality()
          && !i->is_inconsistent())
        throw_topology_incompatible("add_recycled_constraints(cs)",
                                    "cs", cs);
    }
    // If we reach this point, all strict inequalities were inconsistent.
    set_empty();
    return;
  }

  // Dimension-compatibility check:
  // the dimension of `cs' can not be greater than space_dim.
  const dimension_type cs_space_dim = cs.space_dimension();
  if (space_dim < cs_space_dim)
    throw_dimension_incompatible("add_recycled_constraints(cs)", "cs", cs);

  // Adding no constraints is a no-op.
  if (cs.has_no_rows())
    return;

  if (space_dim == 0) {
    // In a 0-dimensional space the constraints are
    // tautologies (e.g., 0 == 0 or 1 >= 0 or 1 > 0) or
    // inconsistent (e.g., 1 == 0 or -1 >= 0 or 0 > 0).
    // In a system of constraints `begin()' and `end()' are equal
    // if and only if the system only contains tautologies.
    if (cs.begin() != cs.end())
      // There is a constraint, it must be inconsistent,
      // the polyhedron is empty.
      status.set_empty();
    return;
  }

  if (marked_empty())
    return;

  // The constraints (possibly with pending rows) are required.
  if (has_pending_generators())
    process_pending_generators();
  else if (!constraints_are_up_to_date())
    update_constraints();

  // Adjust `cs' to the right topology and space dimension.
  // NOTE: we already checked for topology compatibility.
  cs.adjust_topology_and_space_dimension(topology(), space_dim);

  const bool adding_pending = can_have_something_pending();

  // Here we do not require `con_sys' to be sorted.
  // also, we _swap_ (instead of copying) the coefficients of `cs'
  // (which is not a const).
  const dimension_type old_num_rows = con_sys.num_rows();
  const dimension_type cs_num_rows = cs.num_rows();
  const dimension_type cs_num_columns = cs.num_columns();
  con_sys.add_zero_rows(cs_num_rows,
                        Linear_Row::Flags(topology(),
                                          Linear_Row::RAY_OR_POINT_OR_INEQUALITY));
  for (dimension_type i = cs_num_rows; i-- > 0; ) {
    // NOTE: we cannot directly swap the rows, since they might have
    // different capacities (besides possibly having different sizes):
    // thus, we steal one coefficient at a time.
    Constraint& new_c = con_sys[old_num_rows + i];
    Constraint& old_c = cs[i];
    if (old_c.is_equality())
      new_c.set_is_equality();
    using std::swap;
    for (dimension_type j = cs_num_columns; j-- > 0; )
      swap(new_c[j], old_c[j]);
  }

  if (adding_pending)
    set_constraints_pending();
  else {
    // The newly added ones are not pending constraints.
    con_sys.unset_pending_rows();
    // They have been simply appended.
    con_sys.set_sorted(false);
    // Constraints are not minimized and generators are not up-to-date.
    clear_constraints_minimized();
    clear_generators_up_to_date();
  }
  // Note: the constraint system may have become unsatisfiable, thus
  // we do not check for satisfiability.
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::add_recycled_generators ( Generator_System & gs )

Adds the generators in gs to the system of generators of *this (without minimizing the result).

Parameters:

gs	The generator system to be added to `*this`. The generators in `gs` may be recycled.

Exceptions:

std::invalid_argument Thrown if *this and gs are topology-incompatible or dimension-incompatible, or if *this is empty and the system of generators gs is not empty, but has no points.

Warning:: The only assumption that can be made on gs upon successful or exceptional return is that it can be safely destroyed.

Definition at line 1543 of file Polyhedron_public.cc.

Referenced by BHRZ03_evolving_points(), and BHRZ03_evolving_rays().

                                                           {
  // Topology compatibility check.
  if (is_necessarily_closed() && gs.has_closure_points())
    throw_topology_incompatible("add_recycled_generators(gs)", "gs", gs);
  // Dimension-compatibility check:
  // the dimension of `gs' can not be greater than space_dim.
  const dimension_type gs_space_dim = gs.space_dimension();
  if (space_dim < gs_space_dim)
    throw_dimension_incompatible("add_recycled_generators(gs)", "gs", gs);

  // Adding no generators is a no-op.
  if (gs.has_no_rows())
    return;

  // Adding valid generators to a zero-dimensional polyhedron
  // transform it in the zero-dimensional universe polyhedron.
  if (space_dim == 0) {
    if (marked_empty() && !gs.has_points())
      throw_invalid_generators("add_recycled_generators(gs)", "gs");
    set_zero_dim_univ();
    PPL_ASSERT_HEAVY(OK(true));
    return;
  }

  // Adjust `gs' to the right topology and dimensions.
  // NOTE: we already checked for topology compatibility.
  gs.adjust_topology_and_space_dimension(topology(), space_dim);
  // For NNC polyhedra, each point must be matched by
  // the corresponding closure point.
  if (!is_necessarily_closed())
    gs.add_corresponding_closure_points();

  // The generators (possibly with pending rows) are required.
  if ((has_pending_constraints() && !process_pending_constraints())
      || (!generators_are_up_to_date() && !minimize())) {
    // We have just discovered that `*this' is empty.
    // So `gs' must contain at least one point.
    if (!gs.has_points())
      throw_invalid_generators("add_recycled_generators(gs)", "gs");
    // The polyhedron is no longer empty and generators are up-to-date.
    using std::swap;
    swap(gen_sys, gs);
    if (gen_sys.num_pending_rows() > 0) {
      // Even though `gs' has pending generators, since the constraints
      // of the polyhedron are not up-to-date, the polyhedron cannot
      // have pending generators. By integrating the pending part
      // of `gen_sys' we may loose sortedness.
      gen_sys.unset_pending_rows();
      gen_sys.set_sorted(false);
    }
    set_generators_up_to_date();
    clear_empty();
    PPL_ASSERT_HEAVY(OK());
    return;
  }

  const bool adding_pending = can_have_something_pending();

  // Here we do not require `gen_sys' to be sorted.
  // also, we _swap_ (instead of copying) the coefficients of `gs'
  // (which is not a const).
  const dimension_type old_num_rows = gen_sys.num_rows();
  const dimension_type gs_num_rows = gs.num_rows();
  const dimension_type gs_num_columns = gs.num_columns();
  gen_sys.add_zero_rows(gs_num_rows,
                        Linear_Row::Flags(topology(),
                                          Linear_Row::RAY_OR_POINT_OR_INEQUALITY));
  for (dimension_type i = gs_num_rows; i-- > 0; ) {
    // NOTE: we cannot directly swap the rows, since they might have
    // different capacities (besides possibly having different sizes):
    // thus, we steal one coefficient at a time.
    Generator& new_g = gen_sys[old_num_rows + i];
    Generator& old_g = gs[i];
    if (old_g.is_line())
      new_g.set_is_line();
    using std::swap;
    for (dimension_type j = gs_num_columns; j-- > 0; )
      swap(new_g[j], old_g[j]);
  }

  if (adding_pending)
    set_generators_pending();
  else {
    // The newly added ones are not pending generators.
    gen_sys.unset_pending_rows();
    // They have been simply appended.
    gen_sys.set_sorted(false);
    // Constraints are not up-to-date and generators are not minimized.
    clear_constraints_up_to_date();
    clear_generators_minimized();
  }
  PPL_ASSERT_HEAVY(OK(true));
}

void Parma_Polyhedra_Library::Polyhedron::add_space_dimensions	(	Linear_System &	sys1,
		Linear_System &	sys2,
		Bit_Matrix &	sat1,
		Bit_Matrix &	sat2,
		dimension_type	add_dim
	)

staticprivate

Adds new space dimensions to the given linear systems.

Parameters:

sys1	The linear system to which columns are added;
sys2	The linear system to which rows and columns are added;
sat1	The saturation matrix whose columns are indexed by the rows of `sys1`. On entry it is up-to-date;
sat2	The saturation matrix whose columns are indexed by the rows of `sys2`;
add_dim	The number of space dimensions to add.

Adds new space dimensions to the vector space modifying the linear systems and saturation matrices. This function is invoked only by add_space_dimensions_and_embed() and add_space_dimensions_and_project(), passing the linear system of constraints and that of generators (and the corresponding saturation matrices) in different order (see those methods for details).

Definition at line 35 of file Polyhedron_chdims.cc.

                                                              {
  PPL_ASSERT(sys1.topology() == sys2.topology());
  PPL_ASSERT(sys1.num_columns() == sys2.num_columns());
  PPL_ASSERT(add_dim != 0);
  using std::swap;

  sys1.add_zero_columns(add_dim);
  dimension_type old_index = sys2.first_pending_row();
  sys2.add_rows_and_columns(add_dim);
  // The added rows are in the non-pending part.
  sys2.set_index_first_pending_row(old_index + add_dim);

  // The resulting saturation matrix will be as follows:
  // from row    0    to      add_dim-1       : only zeroes
  //          add_dim     add_dim+num_rows-1  : old saturation matrix

  // In fact all the old generators saturate all the new constraints
  // because the polyhedron has not been embedded in the new space.
  sat1.resize(sat1.num_rows() + add_dim, sat1.num_columns());
  // The old matrix is moved to the end of the new matrix.
  for (dimension_type i = sat1.num_rows() - add_dim; i-- > 0; )
    swap(sat1[i], sat1[i+add_dim]);
  // Computes the "sat_c", too.
  sat2.transpose_assign(sat1);

  if (!sys1.is_necessarily_closed()) {
    // Moving the epsilon coefficients to the new last column.
    dimension_type new_eps_index = sys1.num_columns() - 1;
    dimension_type old_eps_index = new_eps_index - add_dim;
    // This swap preserves sortedness of `sys1'.
    sys1.swap_columns(old_eps_index, new_eps_index);

    // Try to preserve sortedness of `sys2'.
    if (!sys2.is_sorted())
      sys2.swap_columns(old_eps_index, new_eps_index);
    else {
      for (dimension_type i = sys2.num_rows(); i-- > add_dim; ) {
        Linear_Row& r = sys2[i];
        swap(r[old_eps_index], r[new_eps_index]);
      }
      // The upper-right corner of `sys2' contains the J matrix:
      // swap coefficients to preserve sortedness.
      for (dimension_type i = add_dim; i-- > 0; ++old_eps_index) {
        Linear_Row& r = sys2[i];
        swap(r[old_eps_index], r[old_eps_index + 1]);
      }
    }
    // NOTE: since we swapped columns in both `sys1' and `sys2',
    // no swapping is required for `sat1' and `sat2'.
  }
}

void Parma_Polyhedra_Library::Polyhedron::add_space_dimensions_and_embed ( dimension_type m )

Adds m new space dimensions and embeds the old polyhedron in the new vector space.

Parameters:

m	The number of dimensions to add.

Exceptions:

std::length_error Thrown if adding m new space dimensions would cause the vector space to exceed dimension max_space_dimension().

The new space dimensions will be those having the highest indexes in the new polyhedron, which is characterized by a system of constraints in which the variables running through the new dimensions are not constrained. For instance, when starting from the polyhedron $\cP \sseq \Rset^2$ and adding a third space dimension, the result will be the polyhedron

$\bigl\{\, (x, y, z)^\transpose \in \Rset^3 \bigm| (x, y)^\transpose \in \cP \,\bigr\}.$

Definition at line 92 of file Polyhedron_chdims.cc.

References Parma_Polyhedra_Library::check_space_dimension_overflow(), Parma_Polyhedra_Library::max_space_dimension(), PPL_ASSERT, PPL_ASSERT_HEAVY, Parma_Polyhedra_Library::swap(), and Parma_Polyhedra_Library::UNIVERSE.

Referenced by Parma_Polyhedra_Library::Implementation::Termination::all_affine_quasi_ranking_functions_MS(), Parma_Polyhedra_Library::Implementation::Termination::all_affine_ranking_functions_MS(), Parma_Polyhedra_Library::Implementation::Termination::all_affine_ranking_functions_PR(), and Parma_Polyhedra_Library::Implementation::Termination::all_affine_ranking_functions_PR_original().

                                                              {
  // The space dimension of the resulting polyhedron should not
  // overflow the maximum allowed space dimension.
  check_space_dimension_overflow(m, max_space_dimension() - space_dimension(),
                                 topology(),
                                 "add_space_dimensions_and_embed(m)",
                                 "adding m new space dimensions exceeds "
                                 "the maximum allowed space dimension");

  // Adding no dimensions to any polyhedron is a no-op.
  if (m == 0)
    return;

  // Adding dimensions to an empty polyhedron is obtained by adjusting
  // `space_dim' and clearing `con_sys' (since it can contain the
  // unsatisfiable constraint system of the wrong dimension).
  if (marked_empty()) {
    space_dim += m;
    con_sys.clear();
    return;
  }

  // The case of a zero-dimensional space polyhedron.
  if (space_dim == 0) {
    // Since it is not empty, it has to be the universe polyhedron.
    PPL_ASSERT(status.test_zero_dim_univ());
    // We swap `*this' with a newly created
    // universe polyhedron of dimension `m'.
    Polyhedron ph(topology(), m, UNIVERSE);
    m_swap(ph);
    return;
  }

  // To embed an n-dimension space polyhedron in a (n+m)-dimension space,
  // we just add `m' zero-columns to the rows in the system of constraints;
  // in contrast, the system of generators needs additional rows,
  // corresponding to the vectors of the canonical basis
  // for the added dimensions. That is, for each new dimension `x[k]'
  // we add the line having that direction. This is done by invoking
  // the function add_space_dimensions() giving the system of generators
  // as the second argument.
  if (constraints_are_up_to_date())
    if (generators_are_up_to_date()) {
      // `sat_c' must be up to date for add_space_dimensions().
      if (!sat_c_is_up_to_date())
        update_sat_c();
      // Adds rows and/or columns to both matrices.
      // `add_space_dimensions' correctly handles pending constraints
      // or generators.
      add_space_dimensions(con_sys, gen_sys, sat_c, sat_g, m);
    }
    else {
      // Only constraints are up-to-date: no need to modify the generators.
      con_sys.add_zero_columns(m);
      // If the polyhedron is not necessarily closed,
      // move the epsilon coefficients to the last column.
      if (!is_necessarily_closed())
        con_sys.swap_columns(space_dim + 1, space_dim + 1 + m);
    }
  else {
    // Only generators are up-to-date: no need to modify the constraints.
    PPL_ASSERT(generators_are_up_to_date());
    gen_sys.add_rows_and_columns(m);
    // The polyhedron does not support pending generators.
    gen_sys.unset_pending_rows();
    // If the polyhedron is not necessarily closed,
    // move the epsilon coefficients to the last column.
    if (!is_necessarily_closed()) {
      // Try to preserve sortedness of `gen_sys'.
      if (!gen_sys.is_sorted())
        gen_sys.swap_columns(space_dim + 1, space_dim + 1 + m);
      else {
        using std::swap;
        dimension_type old_eps_index = space_dim + 1;
        dimension_type new_eps_index = old_eps_index + m;
        for (dimension_type i = gen_sys.num_rows(); i-- > m; ) {
          Generator& r = gen_sys[i];
          swap(r[old_eps_index], r[new_eps_index]);
        }
        // The upper-right corner of `gen_sys' contains the J matrix:
        // swap coefficients to preserve sortedness.
        for (dimension_type i = m; i-- > 0; ++old_eps_index) {
          Generator& r = gen_sys[i];
          swap(r[old_eps_index], r[old_eps_index + 1]);
        }
      }
    }
  }
  // Update the space dimension.
  space_dim += m;

  // Note: we do not check for satisfiability, because the system of
  // constraints may be unsatisfiable.
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::add_space_dimensions_and_project ( dimension_type m )

Adds m new space dimensions to the polyhedron and does not embed it in the new vector space.

Parameters:

m	The number of space dimensions to add.

Exceptions:

std::length_error Thrown if adding m new space dimensions would cause the vector space to exceed dimension max_space_dimension().

The new space dimensions will be those having the highest indexes in the new polyhedron, which is characterized by a system of constraints in which the variables running through the new dimensions are all constrained to be equal to 0. For instance, when starting from the polyhedron $\cP \sseq \Rset^2$ and adding a third space dimension, the result will be the polyhedron

$\bigl\{\, (x, y, 0)^\transpose \in \Rset^3 \bigm| (x, y)^\transpose \in \cP \,\bigr\}.$

Definition at line 189 of file Polyhedron_chdims.cc.

References Parma_Polyhedra_Library::check_space_dimension_overflow(), Parma_Polyhedra_Library::max_space_dimension(), PPL_ASSERT, PPL_ASSERT_HEAVY, Parma_Polyhedra_Library::swap(), Parma_Polyhedra_Library::Generator::zero_dim_closure_point(), and Parma_Polyhedra_Library::Generator::zero_dim_point().

                                                                {
  // The space dimension of the resulting polyhedron should not
  // overflow the maximum allowed space dimension.
  check_space_dimension_overflow(m, max_space_dimension() - space_dimension(),
                                 topology(),
                                 "add_space_dimensions_and_project(m)",
                                 "adding m new space dimensions exceeds "
                                 "the maximum allowed space dimension");

  // Adding no dimensions to any polyhedron is a no-op.
  if (m == 0)
    return;

  // Adding dimensions to an empty polyhedron is obtained
  // by merely adjusting `space_dim'.
  if (marked_empty()) {
    space_dim += m;
    con_sys.clear();
    return;
  }

  if (space_dim == 0) {
    PPL_ASSERT(status.test_zero_dim_univ() && gen_sys.has_no_rows());
    // The system of generators for this polyhedron has only
    // the origin as a point.
    // In an NNC polyhedron, all points have to be accompanied
    // by the corresponding closure points
    // (this time, dimensions are automatically adjusted).
    if (!is_necessarily_closed())
      gen_sys.insert(Generator::zero_dim_closure_point());
    gen_sys.insert(Generator::zero_dim_point());
    gen_sys.adjust_topology_and_space_dimension(topology(), m);
    set_generators_minimized();
    space_dim = m;
    PPL_ASSERT_HEAVY(OK());
    return;
  }

  // To project an n-dimension space polyhedron in a (n+m)-dimension space,
  // we just add to the system of generators `m' zero-columns;
  // In contrast, in the system of constraints, new rows are needed
  // in order to avoid embedding the old polyhedron in the new space.
  // Thus, for each new dimensions `x[k]', we add the constraint
  // x[k] = 0: this is done by invoking the function add_space_dimensions()
  // giving the system of constraints as the second argument.
  if (constraints_are_up_to_date())
    if (generators_are_up_to_date()) {
      // `sat_g' must be up to date for add_space_dimensions().
      if (!sat_g_is_up_to_date())
        update_sat_g();
      // Adds rows and/or columns to both matrices.
      // `add_space_dimensions()' correctly handles pending constraints
      // or generators.
      add_space_dimensions(gen_sys, con_sys, sat_g, sat_c, m);
    }
    else {
      // Only constraints are up-to-date: no need to modify the generators.
      con_sys.add_rows_and_columns(m);
      // The polyhedron does not support pending constraints.
      con_sys.unset_pending_rows();
      // If the polyhedron is not necessarily closed,
      // move the epsilon coefficients to the last column.
      if (!is_necessarily_closed()) {
        // Try to preserve sortedness of `con_sys'.
        if (!con_sys.is_sorted())
          con_sys.swap_columns(space_dim + 1, space_dim + 1 + m);
        else {
          using std::swap;
          dimension_type old_eps_index = space_dim + 1;
          dimension_type new_eps_index = old_eps_index + m;
          for (dimension_type i = con_sys.num_rows(); i-- > m; ) {
            Constraint& r = con_sys[i];
            swap(r[old_eps_index], r[new_eps_index]);
          }
          // The upper-right corner of `con_sys' contains the J matrix:
          // swap coefficients to preserve sortedness.
          for (dimension_type i = m; i-- > 0; ++old_eps_index) {
            Constraint& r = con_sys[i];
            swap(r[old_eps_index], r[old_eps_index + 1]);
          }
        }
      }
    }
  else {
    // Only generators are up-to-date: no need to modify the constraints.
    PPL_ASSERT(generators_are_up_to_date());
    gen_sys.add_zero_columns(m);
    // If the polyhedron is not necessarily closed,
    // move the epsilon coefficients to the last column.
    if (!is_necessarily_closed())
      gen_sys.swap_columns(space_dim + 1, space_dim + 1 + m);
  }
  // Now we update the space dimension.
  space_dim += m;

  // Note: we do not check for satisfiability, because the system of
  // constraints may be unsatisfiable.
  PPL_ASSERT_HEAVY(OK());
}

PPL::dimension_type Parma_Polyhedra_Library::Polyhedron::affine_dimension ( ) const

Returns $0$ , if *this is empty; otherwise, returns the affine dimension of *this.

Definition at line 60 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Constraint_System::begin(), and Parma_Polyhedra_Library::Constraint_System::end().

Referenced by Parma_Polyhedra_Library::Pointset_Powerset< PSET >::affine_dimension(), BHZ09_C_poly_hull_assign_if_exact(), and is_discrete().

                                      {
  if (is_empty())
    return 0;

  const Constraint_System& cs = minimized_constraints();
  dimension_type d = space_dim;
  for (Constraint_System::const_iterator i = cs.begin(),
         cs_end = cs.end(); i != cs_end; ++i)
    if (i->is_equality())
      --d;
  return d;
}

template<typename FP_Format , typename Interval_Info >

void Parma_Polyhedra_Library::Polyhedron::affine_form_image	(	Variable	var,
		const Linear_Form< Interval< FP_Format, Interval_Info > > &	lf
	)

Assigns to *this the affine form image of *this under the function mapping variable var into the affine expression(s) specified by lf.

Parameters:

var	The variable to which the affine expression is assigned.
lf	The linear form on intervals with floating point boundaries that defines the affine expression(s). ALL of its coefficients MUST be bounded.

Exceptions:

std::invalid_argument Thrown if lf and *this are dimension-incompatible or if var is not a space dimension of *this.

This function is used in abstract interpretation to model an assignment of a value that is correctly overapproximated by lf to the floating point variable represented by var.

Definition at line 369 of file Polyhedron.templates.hh.

References bounded_affine_image(), convert_to_integer_expressions(), Parma_Polyhedra_Library::Variable::id(), marked_empty(), overapproximate_linear_form(), Polyhedron(), PPL_ASSERT, PPL_COMPILE_TIME_CHECK, PPL_DIRTY_TEMP_COEFFICIENT, space_dim, throw_dimension_incompatible(), topology(), and Parma_Polyhedra_Library::UNIVERSE.

                                                             {

  // Check that FP_Format is indeed a floating point type.
  PPL_COMPILE_TIME_CHECK(!std::numeric_limits<FP_Format>::is_exact,
                         "Polyhedron::affine_form_image:"
                         " FP_Format not a floating point type.");

  // Dimension compatibility checks.
  // The dimension of lf should not be greater than the dimension of *this.
  const dimension_type lf_space_dim = lf.space_dimension();
  if (space_dim < lf_space_dim)
    throw_dimension_incompatible("affine_form_image(v, l, s)", "l", lf);

  // `var' should be one of the dimensions of the polyhedron.
  const dimension_type var_id = var.id();
  if (space_dim < var_id + 1)
    throw_dimension_incompatible("affine_form_image(v, l, s)", "v", var);

  // We assume that the analyzer will not perform an unreachable assignment.
  PPL_ASSERT(!marked_empty());

  typedef Interval<FP_Format, Interval_Info> FP_Interval_Type;
  typedef Linear_Form<FP_Interval_Type> FP_Linear_Form;

  if (Floating_Point_Expression<FP_Interval_Type, float_ieee754_single>::
      overflows(lf)) {
    *this = Polyhedron(topology(), space_dim, UNIVERSE);
    return;
  }

  // Overapproximate lf.
  FP_Linear_Form lf_approx;
  overapproximate_linear_form(lf, lf_space_dim, lf_approx);

  if (Floating_Point_Expression<FP_Interval_Type, float_ieee754_single>::
      overflows(lf_approx)) {
    *this = Polyhedron(topology(), space_dim, UNIVERSE);
    return;
  }

  // Normalize lf.
  Linear_Expression lf_approx_le;
  PPL_DIRTY_TEMP_COEFFICIENT(lo_coeff);
  PPL_DIRTY_TEMP_COEFFICIENT(hi_coeff);
  PPL_DIRTY_TEMP_COEFFICIENT(denominator);
  convert_to_integer_expressions(lf_approx, lf_space_dim, lf_approx_le,
                                 lo_coeff, hi_coeff, denominator);

  // Finally, do the assignment.
  bounded_affine_image(var, lf_approx_le + lo_coeff, lf_approx_le + hi_coeff,
                       denominator);
}

void Parma_Polyhedra_Library::Polyhedron::affine_image	(	Variable	var,
		const Linear_Expression &	expr,
		Coefficient_traits::const_reference	denominator = `Coefficient_one()`
	)

Assigns to *this the affine image of *this under the function mapping variable var to the affine expression specified by expr and denominator.

Parameters:

var	The variable to which the affine expression is assigned;
expr	The numerator of the affine expression;
denominator	The denominator of the affine expression (optional argument with default value 1).

Exceptions:

std::invalid_argument Thrown if denominator is zero or if expr and *this are dimension-incompatible or if var is not a space dimension of *this.

When considering the generators of a polyhedron, the affine transformation

$\frac{\sum_{i=0}^{n-1} a_i x_i + b}{\mathrm{denominator}}$

is assigned to var where expr is $\sum_{i=0}^{n-1} a_i x_i + b$ ( $b$ is the inhomogeneous term).

If constraints are up-to-date, it uses the specialized function affine_preimage() (for the system of constraints) and inverse transformation to reach the same result. To obtain the inverse transformation we use the following observation.

Observation:

The affine transformation is invertible if the coefficient of var in this transformation (i.e., $a_\mathrm{var}$ ) is different from zero.
If the transformation is invertible, then we can write
$\mathrm{denominator} * {x'}_\mathrm{var} = \sum_{i = 0}^{n - 1} a_i x_i + b = a_\mathrm{var} x_\mathrm{var} + \sum_{i \neq var} a_i x_i + b,$
so that the inverse transformation is
$a_\mathrm{var} x_\mathrm{var} = \mathrm{denominator} * {x'}_\mathrm{var} - \sum_{i \neq j} a_i x_i - b.$

Then, if the transformation is invertible, all the entities that were up-to-date remain up-to-date. Otherwise only generators remain up-to-date.

In other words, if $R$ is a $m_1 \times n$ matrix representing the rays of the polyhedron, $V$ is a $m_2 \times n$ matrix representing the points of the polyhedron and

$P = \bigl\{\, \vect{x} = (x_0, \ldots, x_{n-1})^\mathrm{T} \bigm| \vect{x} = \vect{\lambda} R + \vect{\mu} V, \vect{\lambda} \in \Rset^{m_1}_+, \vect{\mu} \in \Rset^{m_2}_+, \sum_{i = 0}^{m_2 - 1} \mu_i = 1 \,\bigr\}$

and $T$ is the affine transformation to apply to $P$ , then the resulting polyhedron is

$P' = \bigl\{\, (x_0, \ldots, T(x_0, \ldots, x_{n-1}), \ldots, x_{n-1})^\mathrm{T} \bigm| (x_0, \ldots, x_{n-1})^\mathrm{T} \in P \,\bigr\}.$

Affine transformations are, for example:

translations
rotations
symmetries.

Definition at line 2560 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Linear_Expression::coefficient(), Parma_Polyhedra_Library::inverse(), Parma_Polyhedra_Library::neg_assign(), PPL_ASSERT_HEAVY, Parma_Polyhedra_Library::Variable::space_dimension(), and Parma_Polyhedra_Library::Linear_Expression::space_dimension().

Referenced by fold_space_dimensions().

                                                            {
  // The denominator cannot be zero.
  if (denominator == 0)
    throw_invalid_argument("affine_image(v, e, d)", "d == 0");

  // Dimension-compatibility checks.
  // The dimension of `expr' should not be greater than the dimension
  // of `*this'.
  if (space_dim < expr.space_dimension())
    throw_dimension_incompatible("affine_image(v, e, d)", "e", expr);
  // `var' should be one of the dimensions of the polyhedron.
  const dimension_type var_space_dim = var.space_dimension();
  if (space_dim < var_space_dim)
    throw_dimension_incompatible("affine_image(v, e, d)", "v", var);

  if (marked_empty())
    return;

  if (expr.coefficient(var) != 0) {
    // The transformation is invertible:
    // minimality and saturators are preserved, so that
    // pending rows, if present, are correctly handled.
    if (generators_are_up_to_date()) {
      // Generator_System::affine_image() requires the third argument
      // to be a positive Coefficient.
      if (denominator > 0)
        gen_sys.affine_image(var_space_dim, expr, denominator);
      else
        gen_sys.affine_image(var_space_dim, -expr, -denominator);
    }
    if (constraints_are_up_to_date()) {
      // To build the inverse transformation,
      // after copying and negating `expr',
      // we exchange the roles of `expr[var_space_dim]' and `denominator'.
      Linear_Expression inverse;
      if (expr[var_space_dim] > 0) {
        inverse = -expr;
        inverse[var_space_dim] = denominator;
        con_sys.affine_preimage(var_space_dim, inverse, expr[var_space_dim]);
      }
      else {
        // The new denominator is negative: we negate everything once
        // more, as Constraint_System::affine_preimage() requires the
        // third argument to be positive.
        inverse = expr;
        inverse[var_space_dim] = denominator;
        neg_assign(inverse[var_space_dim]);
        con_sys.affine_preimage(var_space_dim, inverse, -expr[var_space_dim]);
      }
    }
  }
  else {
    // The transformation is not invertible.
    // We need an up-to-date system of generators.
    if (has_something_pending())
      remove_pending_to_obtain_generators();
    else if (!generators_are_up_to_date())
      minimize();
    if (!marked_empty()) {
      // Generator_System::affine_image() requires the third argument
      // to be a positive Coefficient.
      if (denominator > 0)
        gen_sys.affine_image(var_space_dim, expr, denominator);
      else
        gen_sys.affine_image(var_space_dim, -expr, -denominator);

      clear_constraints_up_to_date();
      clear_generators_minimized();
      clear_sat_c_up_to_date();
      clear_sat_g_up_to_date();
    }
  }
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::affine_preimage	(	Variable	var,
		const Linear_Expression &	expr,
		Coefficient_traits::const_reference	denominator = `Coefficient_one()`
	)

Assigns to *this the affine preimage of *this under the function mapping variable var to the affine expression specified by expr and denominator.

Parameters:

var	The variable to which the affine expression is substituted;
expr	The numerator of the affine expression;
denominator	The denominator of the affine expression (optional argument with default value 1).

Exceptions:

std::invalid_argument Thrown if denominator is zero or if expr and *this are dimension-incompatible or if var is not a space dimension of *this.

When considering constraints of a polyhedron, the affine transformation

$\frac{\sum_{i=0}^{n-1} a_i x_i + b}{denominator},$

is assigned to var where expr is $\sum_{i=0}^{n-1} a_i x_i + b$ ( $b$ is the inhomogeneous term).

If generators are up-to-date, then the specialized function affine_image() is used (for the system of generators) and inverse transformation to reach the same result. To obtain the inverse transformation, we use the following observation.

Observation:

The affine transformation is invertible if the coefficient of var in this transformation (i.e. $a_\mathrm{var}$ ) is different from zero.
If the transformation is invertible, then we can write
$\mathrm{denominator} * {x'}_\mathrm{var} = \sum_{i = 0}^{n - 1} a_i x_i + b = a_\mathrm{var} x_\mathrm{var} + \sum_{i \neq \mathrm{var}} a_i x_i + b,$
, the inverse transformation is
$a_\mathrm{var} x_\mathrm{var} = \mathrm{denominator} * {x'}_\mathrm{var} - \sum_{i \neq j} a_i x_i - b.$
.

Then, if the transformation is invertible, all the entities that were up-to-date remain up-to-date. Otherwise only constraints remain up-to-date.

In other words, if $A$ is a $m \times n$ matrix representing the constraints of the polyhedron, $T$ is the affine transformation to apply to $P$ and

$P = \bigl\{\, \vect{x} = (x_0, \ldots, x_{n-1})^\mathrm{T} \bigm| A\vect{x} \geq \vect{0} \,\bigr\}.$

The resulting polyhedron is

$P' = \bigl\{\, \vect{x} = (x_0, \ldots, x_{n-1}))^\mathrm{T} \bigm| A'\vect{x} \geq \vect{0} \,\bigr\},$

where $A'$ is defined as follows:

${a'}_{ij} = \begin{cases} a_{ij} * \mathrm{denominator} + a_{i\mathrm{var}}*\mathrm{expr}[j] \quad \mathrm{for } j \neq \mathrm{var}; \\ \mathrm{expr}[\mathrm{var}] * a_{i\mathrm{var}}, \quad \text{for } j = \mathrm{var}. \end{cases}$

Definition at line 2640 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Linear_Expression::coefficient(), Parma_Polyhedra_Library::inverse(), Parma_Polyhedra_Library::neg_assign(), PPL_ASSERT_HEAVY, Parma_Polyhedra_Library::Variable::space_dimension(), and Parma_Polyhedra_Library::Linear_Expression::space_dimension().

                                                               {
  // The denominator cannot be zero.
  if (denominator == 0)
    throw_invalid_argument("affine_preimage(v, e, d)", "d == 0");

  // Dimension-compatibility checks.
  // The dimension of `expr' should not be greater than the dimension
  // of `*this'.
  if (space_dim < expr.space_dimension())
    throw_dimension_incompatible("affine_preimage(v, e, d)", "e", expr);
  // `var' should be one of the dimensions of the polyhedron.
  const dimension_type var_space_dim = var.space_dimension();
  if (space_dim < var_space_dim)
    throw_dimension_incompatible("affine_preimage(v, e, d)", "v", var);

  if (marked_empty())
    return;

  if (expr.coefficient(var) != 0) {
    // The transformation is invertible:
    // minimality and saturators are preserved.
    if (constraints_are_up_to_date()) {
      // Constraint_System::affine_preimage() requires the third argument
      // to be a positive Coefficient.
      if (denominator > 0)
        con_sys.affine_preimage(var_space_dim, expr, denominator);
      else
        con_sys.affine_preimage(var_space_dim, -expr, -denominator);
    }
    if (generators_are_up_to_date()) {
      // To build the inverse transformation,
      // after copying and negating `expr',
      // we exchange the roles of `expr[var_space_dim]' and `denominator'.
      Linear_Expression inverse;
      if (expr[var_space_dim] > 0) {
        inverse = -expr;
        inverse[var_space_dim] = denominator;
        gen_sys.affine_image(var_space_dim, inverse, expr[var_space_dim]);
      }
      else {
        // The new denominator is negative:
        // we negate everything once more, as Generator_System::affine_image()
        // requires the third argument to be positive.
        inverse = expr;
        inverse[var_space_dim] = denominator;
        neg_assign(inverse[var_space_dim]);
        gen_sys.affine_image(var_space_dim, inverse, -expr[var_space_dim]);
      }
    }
  }
  else {
    // The transformation is not invertible.
    // We need an up-to-date system of constraints.
    if (has_something_pending())
      remove_pending_to_obtain_constraints();
    else if (!constraints_are_up_to_date())
      minimize();
    // Constraint_System::affine_preimage() requires the third argument
    // to be a positive Coefficient.
    if (denominator > 0)
      con_sys.affine_preimage(var_space_dim, expr, denominator);
    else
      con_sys.affine_preimage(var_space_dim, -expr, -denominator);
    // Generators, minimality and saturators are no longer valid.
    clear_generators_up_to_date();
    clear_constraints_minimized();
    clear_sat_c_up_to_date();
    clear_sat_g_up_to_date();
  }
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::ascii_dump ( ) const

Writes to std::cerr an ASCII representation of *this.

void Parma_Polyhedra_Library::Polyhedron::ascii_dump ( std::ostream & s ) const

Writes to s an ASCII representation of *this.

Definition at line 3689 of file Polyhedron_public.cc.

                                             {
  s << "space_dim " << space_dim << "\n";
  status.ascii_dump(s);
  s << "\ncon_sys ("
    << (constraints_are_up_to_date() ? "" : "not_")
    << "up-to-date)"
    << "\n";
  con_sys.ascii_dump(s);
  s << "\ngen_sys ("
    << (generators_are_up_to_date() ? "" : "not_")
    << "up-to-date)"
    << "\n";
  gen_sys.ascii_dump(s);
  s << "\nsat_c\n";
  sat_c.ascii_dump(s);
  s << "\nsat_g\n";
  sat_g.ascii_dump(s);
  s << "\n";
}

bool Parma_Polyhedra_Library::Polyhedron::ascii_load ( std::istream & s )

Loads from s an ASCII representation (as produced by ascii_dump(std::ostream&) const) and sets *this accordingly. Returns true if successful, false otherwise.

Definition at line 3712 of file Polyhedron_public.cc.

References PPL_ASSERT_HEAVY.

                                       {
  std::string str;

  if (!(s >> str) || str != "space_dim")
    return false;

  if (!(s >> space_dim))
    return false;

  if (!status.ascii_load(s))
    return false;

  if (!(s >> str) || str != "con_sys")
    return false;

  if (!(s >> str) || (str != "(not_up-to-date)" && str != "(up-to-date)"))
    return false;

  if (!con_sys.ascii_load(s))
    return false;

  if (!(s >> str) || str != "gen_sys")
    return false;

  if (!(s >> str) || (str != "(not_up-to-date)" && str != "(up-to-date)"))
    return false;

  if (!gen_sys.ascii_load(s))
    return false;

  if (!(s >> str) || str != "sat_c")
    return false;

  if (!sat_c.ascii_load(s))
    return false;

  if (!(s >> str) || str != "sat_g")
    return false;

  if (!sat_g.ascii_load(s))
    return false;

  // Check invariants.
  PPL_ASSERT_HEAVY(OK());
  return true;
}

bool Parma_Polyhedra_Library::Polyhedron::BFT00_poly_hull_assign_if_exact ( const Polyhedron & y )

protected

If the poly-hull of *this and y is exact it is assigned to *this and true is returned, otherwise false is returned.

Current implementation is based on (a variant of) Algorithm 8.1 in A. Bemporad, K. Fukuda, and F. D. Torrisi Convexity Recognition of the Union of Polyhedra Technical Report AUT00-13, ETH Zurich, 2000

Note:: It is assumed that *this and y are topologically closed and dimension-compatible; if the assumption does not hold, the behavior is undefined.

Definition at line 1902 of file Polyhedron_nonpublic.cc.

                                                                  {
  // Declare a const reference to *this (to avoid accidental modifications).
  const Polyhedron& x = *this;
  // Private method: the caller must ensure the following.
  PPL_ASSERT(x.is_necessarily_closed());
  PPL_ASSERT(x.topology() == y.topology());
  PPL_ASSERT(x.space_dim == y.space_dim);

  // The zero-dim case is trivial.
  if (x.space_dim == 0) {
    upper_bound_assign(y);
    return true;
  }
  // If `x' or `y' is (known to be) empty, the convex union is exact.
  if (x.marked_empty()) {
    *this = y;
    return true;
  }
  else if (y.is_empty())
    return true;
  else if (x.is_empty()) {
    *this = y;
    return true;
  }

  // Here both `x' and `y' are known to be non-empty.

  // Implementation based on Algorithm 8.1 (page 15) in [BemporadFT00TR],
  // generalized so as to also allow for unbounded polyhedra.
  // The extension to unbounded polyhedra is obtained by mimicking
  // what done in Algorithm 8.2 (page 19) with respect to
  // Algorithm 6.2 (page 13).
  // We also apply a couple of improvements (see steps 2.1, 3.1, 6.1, 7.1)
  // so as to quickly handle special cases and avoid the splitting
  // of equalities/lines into pairs of inequalities/rays.

  (void) x.minimize();
  (void) y.minimize();
  const Constraint_System& x_cs = x.con_sys;
  const Constraint_System& y_cs = y.con_sys;
  const Generator_System& x_gs = x.gen_sys;
  const Generator_System& y_gs = y.gen_sys;
  const dimension_type x_gs_num_rows = x_gs.num_rows();
  const dimension_type y_gs_num_rows = y_gs.num_rows();

  // Step 1: generators of `x' that are redundant in `y', and vice versa.
  std::vector<bool> x_gs_red_in_y(x_gs_num_rows, false);
  dimension_type num_x_gs_red_in_y = 0;
  for (dimension_type i = x_gs_num_rows; i-- > 0; )
    if (y.relation_with(x_gs[i]).implies(Poly_Gen_Relation::subsumes())) {
      x_gs_red_in_y[i] = true;
      ++num_x_gs_red_in_y;
    }
  std::vector<bool> y_gs_red_in_x(y_gs_num_rows, false);
  dimension_type num_y_gs_red_in_x = 0;
  for (dimension_type i = y_gs_num_rows; i-- > 0; )
    if (x.relation_with(y_gs[i]).implies(Poly_Gen_Relation::subsumes())) {
      y_gs_red_in_x[i] = true;
      ++num_y_gs_red_in_x;
    }

  // Step 2: if no redundant generator has been identified,
  // then the union is not convex. CHECKME: why?
  if (num_x_gs_red_in_y == 0 && num_y_gs_red_in_x == 0)
    return false;

  // Step 2.1: while at it, also perform quick inclusion tests.
  if (num_y_gs_red_in_x == y_gs_num_rows)
    // `y' is included into `x': union is convex.
    return true;
  if (num_x_gs_red_in_y == x_gs_num_rows) {
    // `x' is included into `y': union is convex.
    *this = y;
    return true;
  }

  // Here we know that `x' is not included in `y', and vice versa.

  // Step 3: constraints of `x' that are satisfied by `y', and vice versa.
  const dimension_type x_cs_num_rows = x_cs.num_rows();
  std::vector<bool> x_cs_red_in_y(x_cs_num_rows, false);
  for (dimension_type i = x_cs_num_rows; i-- > 0; ) {
    const Constraint& x_cs_i = x_cs[i];
    if (y.relation_with(x_cs_i).implies(Poly_Con_Relation::is_included()))
      x_cs_red_in_y[i] = true;
    else if (x_cs_i.is_equality())
      // Step 3.1: `x' has an equality not satisfied by `y':
      // union is not convex (recall that `y' does not contain `x').
      // NOTE: this would be false for NNC polyhedra.
      // Example: x = { A == 0 }, y = { 0 < A <= 1 }.
      return false;
  }
  const dimension_type y_cs_num_rows = y_cs.num_rows();
  std::vector<bool> y_cs_red_in_x(y_cs_num_rows, false);
  for (dimension_type i = y_cs_num_rows; i-- > 0; ) {
    const Constraint& y_cs_i = y_cs[i];
    if (x.relation_with(y_cs_i).implies(Poly_Con_Relation::is_included()))
      y_cs_red_in_x[i] = true;
    else if (y_cs_i.is_equality())
      // Step 3.1: `y' has an equality not satisfied by `x':
      // union is not convex (see explanation above).
      return false;
  }

  // Loop in steps 5-9: for each pair of non-redundant generators,
  // compute their "mid-point" and check if it is both in `x' and `y'.

  // Note: reasoning at the polyhedral cone level.
  // CHECKME, FIXME: Polyhedron is a (deprecated) friend of Generator.
  // Here below we systematically exploit such a friendship, so as to
  // freely reinterpret a Generator as a Linear_Row and vice versa.
  Linear_Row mid_row;
  const Generator& mid_g = static_cast<const Generator&>(mid_row);

  for (dimension_type i = x_gs_num_rows; i-- > 0; ) {
    if (x_gs_red_in_y[i])
      continue;
    const Linear_Row& x_row = static_cast<const Linear_Row&>(x_gs[i]);
    const dimension_type row_sz = x_row.size();
    const bool x_row_is_line = x_row.is_line_or_equality();
    for (dimension_type j = y_gs_num_rows; j-- > 0; ) {
      if (y_gs_red_in_x[j])
        continue;
      const Linear_Row& y_row = static_cast<const Linear_Row&>(y_gs[j]);
      const bool y_row_is_line = y_row.is_line_or_equality();

      // Step 6: compute mid_row = x_row + y_row.
      // NOTE: no need to actually compute the "mid-point",
      // since any strictly positive combination would do.
      mid_row = x_row;
      for (dimension_type k = row_sz; k-- > 0; )
        mid_row[k] += y_row[k];
      // A zero ray is not a well formed generator.
      const bool illegal_ray
        = (mid_row[0] == 0 && mid_row.all_homogeneous_terms_are_zero());
      // A zero ray cannot be generated from a line: this holds
      // because x_row (resp., y_row) is not subsumed by y (resp., x).
      PPL_ASSERT(!(illegal_ray && (x_row_is_line || y_row_is_line)));
      if (illegal_ray)
        continue;
      // Normalize mid_row (strongly, if needed).
      mid_row.normalize();
      if (x_row_is_line) {
        if (y_row_is_line)
          // mid_row is a line too: sign normalization is needed.
          mid_row.sign_normalize();
        else
          // mid_row is a ray/point.
          mid_row.set_is_ray_or_point_or_inequality();
      }
      PPL_ASSERT(mid_g.OK());

      // Step 7: check if mid_g is in the union of x and y.
      if (x.relation_with(mid_g) == Poly_Gen_Relation::nothing()
          && y.relation_with(mid_g) == Poly_Gen_Relation::nothing())
        return false;

      // If either x_row or y_row is a line, we should use its
      // negation to produce another generator to be tested too.
      // NOTE: exclusive-or is meant.
      if (!x_row_is_line && y_row_is_line) {
        // Step 6.1: (re-)compute mid_row = x_row - y_row.
        mid_row = x_row;
        for (dimension_type k = row_sz; k-- > 0; )
          mid_row[k] -= y_row[k];
        mid_row.normalize();
        // Step 7.1: check if mid_g is in the union of x and y.
        if (x.relation_with(mid_g) == Poly_Gen_Relation::nothing()
            && y.relation_with(mid_g) == Poly_Gen_Relation::nothing())
          return false;
      }
      else if (x_row_is_line && !y_row_is_line) {
        // Step 6.1: (re-)compute mid_row = - x_row + y_row.
        mid_row = y_row;
        for (dimension_type k = row_sz; k-- > 0; )
          mid_row[k] -= x_row[k];
        mid_row.normalize();
        // Step 7.1: check if mid_g is in the union of x and y.
        if (x.relation_with(mid_g) == Poly_Gen_Relation::nothing()
            && y.relation_with(mid_g) == Poly_Gen_Relation::nothing())
          return false;
      }
    }
  }

  // Here we know that the union of x and y is convex.
  // TODO: exploit knowledge on the cardinality of non-redundant
  // constraints/generators to improve the convex-hull computation.
  // Using generators allows for exploiting incrementality.
  for (dimension_type j = 0; j < y_gs_num_rows; ++j) {
    if (!y_gs_red_in_x[j])
      add_generator(y_gs[j]);
  }
  PPL_ASSERT_HEAVY(OK());
  return true;
}

bool Parma_Polyhedra_Library::Polyhedron::BHRZ03_combining_constraints	(	const Polyhedron &	y,
		const BHRZ03_Certificate &	y_cert,
		const Polyhedron &	H79,
		const Constraint_System &	x_minus_H79_cs
	)

private

Definition at line 374 of file Polyhedron_widenings.cc.

Referenced by BHRZ03_widening_assign().

                                                                        {
  Polyhedron& x = *this;
  // It is assumed that `y <= x <= H79'.
  PPL_ASSERT(x.topology() == y.topology()
         && x.topology() == H79.topology()
         && x.topology() == x_minus_H79_cs.topology());
  PPL_ASSERT(x.space_dim == y.space_dim
         && x.space_dim == H79.space_dim
         && x.space_dim == x_minus_H79_cs.space_dimension());
  PPL_ASSERT(!x.marked_empty() && !x.has_something_pending()
         && x.constraints_are_minimized() && x.generators_are_minimized());
  PPL_ASSERT(!y.marked_empty() && !y.has_something_pending()
         && y.constraints_are_minimized() && y.generators_are_minimized());
  PPL_ASSERT(!H79.marked_empty() && !H79.has_something_pending()
         && H79.constraints_are_minimized() && H79.generators_are_minimized());

  // We will choose from `x_minus_H79_cs' many subsets of constraints,
  // that will be collected (one at a time) in `combining_cs'.
  // For each group collected, we compute an average constraint,
  // that will be stored in `new_cs'.

  // There is no point in applying this technique when `x_minus_H79_cs'
  // has one constraint at most (no ``new'' constraint can be computed).
  const dimension_type x_minus_H79_cs_num_rows = x_minus_H79_cs.num_rows();
  if (x_minus_H79_cs_num_rows <= 1)
    return false;

  const Topology topol = x.topology();
  Constraint_System combining_cs(topol);
  Constraint_System new_cs(topol);

  // Consider the points that belong to both `x.gen_sys' and `y.gen_sys'.
  // For NNC polyhedra, the role of points is played by closure points.
  const bool closed = x.is_necessarily_closed();
  for (dimension_type i = y.gen_sys.num_rows(); i-- > 0; ) {
    const Generator& g = y.gen_sys[i];
    if ((g.is_point() && closed) || (g.is_closure_point() && !closed)) {
      // If in `H79.con_sys' there is already an inequality constraint
      // saturating this point, then there is no need to produce another
      // constraint.
      bool lies_on_the_boundary_of_H79 = false;
      const Constraint_System& H79_cs = H79.con_sys;
      for (dimension_type j = H79_cs.num_rows(); j-- > 0; ) {
        const Constraint& c = H79_cs[j];
        if (c.is_inequality() && Scalar_Products::sign(c, g) == 0) {
          lies_on_the_boundary_of_H79 = true;
          break;
        }
      }
      if (lies_on_the_boundary_of_H79)
        continue;

      // Consider all the constraints in `x_minus_H79_cs'
      // that are saturated by the point `g'.
      combining_cs.clear();
      for (dimension_type j = x_minus_H79_cs_num_rows; j-- > 0; ) {
        const Constraint& c = x_minus_H79_cs[j];
        if (Scalar_Products::sign(c, g) == 0)
          combining_cs.insert(c);
      }
      // Build a new constraint by combining all the chosen constraints.
      const dimension_type combining_cs_num_rows = combining_cs.num_rows();
      if (combining_cs_num_rows > 0) {
        if (combining_cs_num_rows == 1)
          // No combination is needed.
          new_cs.insert(combining_cs[0]);
        else {
          Linear_Expression e(0);
          bool strict_inequality = false;
          for (dimension_type h = combining_cs_num_rows; h-- > 0; ) {
            if (combining_cs[h].is_strict_inequality())
              strict_inequality = true;
            e += Linear_Expression(combining_cs[h]);
          }

          if (!e.all_homogeneous_terms_are_zero()) {
            if (strict_inequality)
              new_cs.insert(e > 0);
            else
              new_cs.insert(e >= 0);
          }
        }
      }
    }
  }

  // If none of the collected constraints strictly intersects `H79',
  // then the technique was unsuccessful.
  bool improves_upon_H79 = false;
  const Poly_Con_Relation si = Poly_Con_Relation::strictly_intersects();
  for (dimension_type i = new_cs.num_rows(); i-- > 0; )
    if (H79.relation_with(new_cs[i]) == si) {
      improves_upon_H79 = true;
      break;
    }
  if (!improves_upon_H79)
    return false;

  // The resulting polyhedron is obtained by adding the constraints
  // in `new_cs' to polyhedron `H79'.
  Polyhedron result = H79;
  result.add_recycled_constraints(new_cs);
  // Force minimization.
  result.minimize();

  // Check for stabilization with respect to `y_cert' and improvement
  // over `H79'.
  if (y_cert.is_stabilizing(result) && !result.contains(H79)) {
    // The technique was successful.
    x.m_swap(result);
    PPL_ASSERT_HEAVY(x.OK(true));
    return true;
  }
  else
    // The technique was unsuccessful.
    return false;
}

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

private

Definition at line 496 of file Polyhedron_widenings.cc.

References add_recycled_generators(), constraints_are_minimized(), contains(), gen_sys, generators_are_minimized(), has_something_pending(), intersection_assign(), Parma_Polyhedra_Library::Generator::is_closure_point(), is_necessarily_closed(), Parma_Polyhedra_Library::Generator::is_point(), Parma_Polyhedra_Library::BHRZ03_Certificate::is_stabilizing(), Parma_Polyhedra_Library::Linear_Row::linear_combine(), m_swap(), marked_empty(), minimize(), Parma_Polyhedra_Library::Poly_Gen_Relation::nothing(), Parma_Polyhedra_Library::Dense_Matrix::num_rows(), OK(), PPL_ASSERT, PPL_ASSERT_HEAVY, relation_with(), space_dim, and topology().

Referenced by BHRZ03_widening_assign().

                                                               {
  Polyhedron& x = *this;
  // It is assumed that `y <= x <= H79'.
  PPL_ASSERT(x.topology() == y.topology()
         && x.topology() == H79.topology());
  PPL_ASSERT(x.space_dim == y.space_dim
         && x.space_dim == H79.space_dim);
  PPL_ASSERT(!x.marked_empty() && !x.has_something_pending()
         && x.constraints_are_minimized() && x.generators_are_minimized());
  PPL_ASSERT(!y.marked_empty() && !y.has_something_pending()
         && y.constraints_are_minimized() && y.generators_are_minimized());
  PPL_ASSERT(!H79.marked_empty() && !H79.has_something_pending()
         && H79.constraints_are_minimized() && H79.generators_are_minimized());

  // For each point in `x.gen_sys' that is not in `y',
  // this technique tries to identify a set of rays that:
  //  - are included in polyhedron `H79';
  //  - when added to `y' will subsume the point.
  Generator_System candidate_rays;

  const dimension_type x_gen_sys_num_rows = x.gen_sys.num_rows();
  const dimension_type y_gen_sys_num_rows = y.gen_sys.num_rows();
  const bool closed = x.is_necessarily_closed();
  for (dimension_type i = x_gen_sys_num_rows; i-- > 0; ) {
    Generator& g1 = x.gen_sys[i];
    // For C polyhedra, we choose a point of `x.gen_sys'
    // that is not included in `y'.
    // In the case of NNC polyhedra, we can restrict attention to
    // closure points (considering also points will only add redundancy).
    if (((g1.is_point() && closed) || (g1.is_closure_point() && !closed))
        && y.relation_with(g1) == Poly_Gen_Relation::nothing()) {
      // For each point (resp., closure point) `g2' in `y.gen_sys',
      // where `g1' and `g2' are different,
      // build the candidate ray `g1 - g2'.
      for (dimension_type j = y_gen_sys_num_rows; j-- > 0; ) {
        const Generator& g2 = y.gen_sys[j];
        if ((g2.is_point() && closed)
            || (g2.is_closure_point() && !closed)) {
          PPL_ASSERT(compare(g1, g2) != 0);
          Generator ray_from_g2_to_g1 = g1;
          ray_from_g2_to_g1.linear_combine(g2, 0);
          candidate_rays.insert(ray_from_g2_to_g1);
        }
      }
    }
  }

  // Be non-intrusive.
  Polyhedron result = x;
  result.add_recycled_generators(candidate_rays);
  result.intersection_assign(H79);
  // Force minimization.
  result.minimize();

  // Check for stabilization with respect to `y_cert' and improvement
  // over `H79'.
  if (y_cert.is_stabilizing(result) && !result.contains(H79)) {
    // The technique was successful.
    x.m_swap(result);
    PPL_ASSERT_HEAVY(x.OK(true));
    return true;
  }
  else
    // The technique was unsuccessful.
    return false;
}

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

private

Definition at line 566 of file Polyhedron_widenings.cc.

Referenced by BHRZ03_widening_assign().

                                                             {
  Polyhedron& x = *this;
  // It is assumed that `y <= x <= H79'.
  PPL_ASSERT(x.topology() == y.topology()
         && x.topology() == H79.topology());
  PPL_ASSERT(x.space_dim == y.space_dim
         && x.space_dim == H79.space_dim);
  PPL_ASSERT(!x.marked_empty() && !x.has_something_pending()
         && x.constraints_are_minimized() && x.generators_are_minimized());
  PPL_ASSERT(!y.marked_empty() && !y.has_something_pending()
         && y.constraints_are_minimized() && y.generators_are_minimized());
  PPL_ASSERT(!H79.marked_empty() && !H79.has_something_pending()
         && H79.constraints_are_minimized() && H79.generators_are_minimized());

  const dimension_type x_gen_sys_num_rows = x.gen_sys.num_rows();
  const dimension_type y_gen_sys_num_rows = y.gen_sys.num_rows();

  // Candidate rays are kept in a temporary generator system.
  Generator_System candidate_rays;
  PPL_DIRTY_TEMP_COEFFICIENT(tmp);
  for (dimension_type i = x_gen_sys_num_rows; i-- > 0; ) {
    const Generator& x_g = x.gen_sys[i];
    // We choose a ray of `x' that does not belong to `y'.
    if (x_g.is_ray() && y.relation_with(x_g) == Poly_Gen_Relation::nothing()) {
      for (dimension_type j = y_gen_sys_num_rows; j-- > 0; ) {
        const Generator& y_g = y.gen_sys[j];
        if (y_g.is_ray()) {
          Generator new_ray(x_g);
          // Modify `new_ray' according to the evolution of `x_g' with
          // respect to `y_g'.
          std::deque<bool> considered(x.space_dim + 1);
          for (dimension_type k = 1; k < x.space_dim; ++k)
            if (!considered[k])
              for (dimension_type h = k + 1; h <= x.space_dim; ++h)
                if (!considered[h]) {
                  tmp = x_g[k] * y_g[h];
                  // The following line optimizes the computation of
                  // <CODE> tmp -= x_g[h] * y_g[k]; </CODE>
                  sub_mul_assign(tmp, x_g[h], y_g[k]);
                  const int clockwise
                    = sgn(tmp);
                  const int first_or_third_quadrant
                    = sgn(x_g[k]) * sgn(x_g[h]);
                  switch (clockwise * first_or_third_quadrant) {
                  case -1:
                    new_ray[k] = 0;
                    considered[k] = true;
                    break;
                  case 1:
                    new_ray[h] = 0;
                    considered[h] = true;
                    break;
                  default:
                    break;
                  }
                }
          new_ray.normalize();
          candidate_rays.insert(new_ray);
        }
      }
    }
  }

  // If there are no candidate rays, we cannot obtain stabilization.
  if (candidate_rays.has_no_rows())
    return false;

  // Be non-intrusive.
  Polyhedron result = x;
  result.add_recycled_generators(candidate_rays);
  result.intersection_assign(H79);
  // Force minimization.
  result.minimize();

  // Check for stabilization with respect to `y' and improvement over `H79'.
  if (y_cert.is_stabilizing(result) && !result.contains(H79)) {
    // The technique was successful.
    x.m_swap(result);
    PPL_ASSERT_HEAVY(x.OK(true));
    return true;
  }
  else
    // The technique was unsuccessful.
    return false;
}

void Parma_Polyhedra_Library::Polyhedron::BHRZ03_widening_assign	(	const Polyhedron &	y,
		unsigned *	tp = `0`
	)

Assigns to *this the result of computing the BHRZ03-widening between *this and y.

Parameters:

y	A polyhedron that must be contained in `*this`;
tp	An optional pointer to an unsigned variable storing the number of available tokens (to be used when applying the widening with tokens delay technique).

Exceptions:

std::invalid_argument Thrown if *this and y are topology-incompatible or dimension-incompatible.

Definition at line 655 of file Polyhedron_widenings.cc.

References BHRZ03_combining_constraints(), BHRZ03_evolving_points(), BHRZ03_evolving_rays(), contains(), Parma_Polyhedra_Library::copy_contains(), Parma_Polyhedra_Library::Dense_Matrix::has_no_rows(), Parma_Polyhedra_Library::BHRZ03_Certificate::is_stabilizing(), m_swap(), marked_empty(), minimize(), OK(), PPL_ASSERT, PPL_ASSERT_HEAVY, PPL_EXPECT_HEAVY, select_H79_constraints(), space_dim, topology(), and Parma_Polyhedra_Library::UNIVERSE.

Referenced by limited_BHRZ03_extrapolation_assign().

                                                                       {
  Polyhedron& x = *this;
  // Topology compatibility check.
  if (x.topology() != y.topology())
    throw_topology_incompatible("BHRZ03_widening_assign(y)", "y", y);
  // Dimension-compatibility check.
  if (x.space_dim != y.space_dim)
    throw_dimension_incompatible("BHRZ03_widening_assign(y)", "y", y);

  // Assume `y' is contained in or equal to `x'.
  PPL_EXPECT_HEAVY(copy_contains(x, y));

  // If any argument is zero-dimensional or empty,
  // the BHRZ03-widening behaves as the identity function.
  if (x.space_dim == 0 || x.marked_empty() || y.marked_empty())
    return;

  // `y.con_sys' and `y.gen_sys' should be in minimal form.
  if (!y.minimize())
    // `y' is empty: the result is `x'.
    return;
  // `x.con_sys' and `x.gen_sys' should be in minimal form.
  x.minimize();

  // Compute certificate info for polyhedron `y'.
  BHRZ03_Certificate y_cert(y);

  // If the iteration is stabilizing, the resulting polyhedron is `x'.
  // At this point, also check if the two polyhedra are the same
  // (exploiting the knowledge that `y <= x').
  if (y_cert.is_stabilizing(x) || y.contains(x)) {
    PPL_ASSERT_HEAVY(OK());
    return;
  }

  // Here the iteration is not immediately stabilizing.
  // If we are using the widening-with-tokens technique and
  // there are tokens available, use one of them and return `x'.
  if (tp != 0 && *tp > 0) {
    --(*tp);
    PPL_ASSERT_HEAVY(OK());
    return;
  }

  // Copy into `H79_cs' the constraints that are common to `x' and `y',
  // according to the definition of the H79 widening.
  // The other ones are copied into `x_minus_H79_cs'.
  const Topology topol = x.topology();
  Constraint_System H79_cs(topol);
  Constraint_System x_minus_H79_cs(topol);
  x.select_H79_constraints(y, H79_cs, x_minus_H79_cs);

  // We cannot have selected all of the rows, since otherwise
  // the iteration should have been immediately stabilizing.
  PPL_ASSERT(!x_minus_H79_cs.has_no_rows());
  // Be careful to obtain the right space dimension
  // (because `H79_cs' may be empty).
  Polyhedron H79(topol, x.space_dim, UNIVERSE);
  H79.add_recycled_constraints(H79_cs);
  // Force minimization.
  H79.minimize();

  // NOTE: none of the following widening heuristics is intrusive:
  // they will modify `x' only when returning successfully.
  if (x.BHRZ03_combining_constraints(y, y_cert, H79, x_minus_H79_cs))
    return;

  PPL_ASSERT_HEAVY(H79.OK() && x.OK() && y.OK());

  if (x.BHRZ03_evolving_points(y, y_cert, H79))
    return;

  PPL_ASSERT_HEAVY(H79.OK() && x.OK() && y.OK());

  if (x.BHRZ03_evolving_rays(y, y_cert, H79))
    return;

  PPL_ASSERT_HEAVY(H79.OK() && x.OK() && y.OK());

  // No previous technique was successful: fall back to the H79 widening.
  x.m_swap(H79);
  PPL_ASSERT_HEAVY(x.OK(true));

  // The H79 widening is always stabilizing.
  PPL_ASSERT(y_cert.is_stabilizing(x));
}

bool Parma_Polyhedra_Library::Polyhedron::BHZ09_C_poly_hull_assign_if_exact ( const Polyhedron & y )

protected

Definition at line 1454 of file Polyhedron_nonpublic.cc.

References affine_dimension(), con_sys, gen_sys, Parma_Polyhedra_Library::Poly_Con_Relation::implies(), is_empty(), Parma_Polyhedra_Library::Poly_Con_Relation::is_included(), is_included_in(), is_necessarily_closed(), minimize(), Parma_Polyhedra_Library::Dense_Matrix::num_rows(), PPL_ASSERT, PPL_ASSERT_HEAVY, relation_with(), sat_g, sat_g_is_up_to_date(), Parma_Polyhedra_Library::Bit_Row::set(), Parma_Polyhedra_Library::Bit_Row::set_until(), space_dim, Parma_Polyhedra_Library::Poly_Gen_Relation::subsumes(), Parma_Polyhedra_Library::Bit_Row::union_assign(), and update_sat_g().

Referenced by BHZ09_poly_hull_assign_if_exact().

                                                                    {
  Polyhedron& x = *this;
  // Private method: the caller must ensure the following.
  PPL_ASSERT(x.is_necessarily_closed() && y.is_necessarily_closed());
  PPL_ASSERT(x.space_dim > 0 && x.space_dim == y.space_dim);
  PPL_ASSERT(!x.is_empty() && !y.is_empty());

  // Minimization is not really required, but it is probably the best
  // way of getting constraints, generators and saturation matrices
  // up-to-date.  Minimization it also removes redundant
  // constraints/generators.
  (void) x.minimize();
  (void) y.minimize();

  // Handle a special case: for topologically closed polyhedra P and Q,
  // if the affine dimension of P is greater than that of Q, then
  // their upper bound is exact if and only if P includes Q.
  const dimension_type x_affine_dim = x.affine_dimension();
  const dimension_type y_affine_dim = y.affine_dimension();
  if (x_affine_dim > y_affine_dim)
    return y.is_included_in(x);
  else if (x_affine_dim < y_affine_dim) {
    if (x.is_included_in(y)) {
      x = y;
      return true;
    }
    else
      return false;
  }

  const Constraint_System& x_cs = x.con_sys;
  const Generator_System& x_gs = x.gen_sys;
  const Generator_System& y_gs = y.gen_sys;
  const dimension_type x_gs_num_rows = x_gs.num_rows();
  const dimension_type y_gs_num_rows = y_gs.num_rows();

  // Step 1: generators of `x' that are redundant in `y', and vice versa.
  Bit_Row x_gs_red_in_y;
  dimension_type num_x_gs_red_in_y = 0;
  for (dimension_type i = x_gs_num_rows; i-- > 0; )
    if (y.relation_with(x_gs[i]).implies(Poly_Gen_Relation::subsumes())) {
      x_gs_red_in_y.set(i);
      ++num_x_gs_red_in_y;
    }
  Bit_Row y_gs_red_in_x;
  dimension_type num_y_gs_red_in_x = 0;
  for (dimension_type i = y_gs_num_rows; i-- > 0; )
    if (x.relation_with(y_gs[i]).implies(Poly_Gen_Relation::subsumes())) {
      y_gs_red_in_x.set(i);
      ++num_y_gs_red_in_x;
    }

  // Step 2: filter away special cases.

  // Step 2.1: inclusion tests.
  if (num_y_gs_red_in_x == y_gs_num_rows)
    // `y' is included into `x': upper bound `x' is exact.
    return true;
  if (num_x_gs_red_in_y == x_gs_num_rows) {
    // `x' is included into `y': upper bound `y' is exact.
    x = y;
    return true;
  }

  // Step 2.2: if no generator of `x' is redundant for `y', then
  // (as by 2.1 there exists a constraint of `x' non-redundant for `y')
  // the upper bound is not exact; the same if exchanging `x' and `y'.
  if (num_x_gs_red_in_y == 0 || num_y_gs_red_in_x == 0)
    return false;

  // Step 3: see if `x' has a non-redundant constraint `c_x' that is not
  // satisfied by `y' and a non-redundant generator in `y' (see Step 1)
  // saturating `c_x'. If so, the upper bound is not exact.

  // Make sure the saturation matrix for `x' is up to date.
  // Any sat matrix would do: we choose `sat_g' because it matches
  // the two nested loops (constraints on rows and generators on columns).
  if (!x.sat_g_is_up_to_date())
    x.update_sat_g();
  const Bit_Matrix& x_sat = x.sat_g;

  Bit_Row all_ones;
  all_ones.set_until(x_gs_num_rows);
  Bit_Row row_union;
  for (dimension_type i = x_cs.num_rows(); i-- > 0; ) {
    const bool included
      = y.relation_with(x_cs[i]).implies(Poly_Con_Relation::is_included());
    if (!included) {
      row_union.union_assign(x_gs_red_in_y, x_sat[i]);
      if (row_union != all_ones)
        return false;
    }
  }

  // Here we know that the upper bound is exact: compute it.
  for (dimension_type j = y_gs_num_rows; j-- > 0; )
    if (!y_gs_red_in_x[j])
      add_generator(y_gs[j]);

  PPL_ASSERT_HEAVY(OK());
  return true;
}

bool Parma_Polyhedra_Library::Polyhedron::BHZ09_NNC_poly_hull_assign_if_exact ( const Polyhedron & y )

protected

Definition at line 1558 of file Polyhedron_nonpublic.cc.

Referenced by BHZ09_poly_hull_assign_if_exact().

                                                                      {
  const Polyhedron& x = *this;
  // Private method: the caller must ensure the following.
  PPL_ASSERT(!x.is_necessarily_closed() && !y.is_necessarily_closed());
  PPL_ASSERT(x.space_dim > 0 && x.space_dim == y.space_dim);
  PPL_ASSERT(!x.is_empty() && !y.is_empty());

  // Minimization is not really required, but it is probably the best
  // way of getting constraints, generators and saturation matrices
  // up-to-date.  Minimization also removes redundant
  // constraints/generators.
  (void) x.minimize();
  (void) y.minimize();

  const Generator_System& x_gs = x.gen_sys;
  const Generator_System& y_gs = y.gen_sys;
  const dimension_type x_gs_num_rows = x_gs.num_rows();
  const dimension_type y_gs_num_rows = y_gs.num_rows();

  // Compute generators of `x' that are non-redundant in `y' ...
  Bit_Row x_gs_non_redundant_in_y;
  Bit_Row x_points_non_redundant_in_y;
  Bit_Row x_closure_points;
  dimension_type num_x_gs_non_redundant_in_y = 0;
  for (dimension_type i = x_gs_num_rows; i-- > 0; ) {
    const Generator& x_gs_i = x_gs[i];
    if (x_gs_i.is_closure_point())
      x_closure_points.set(i);
    if (y.relation_with(x_gs[i]).implies(Poly_Gen_Relation::subsumes()))
      continue;
    x_gs_non_redundant_in_y.set(i);
    ++num_x_gs_non_redundant_in_y;
    if (x_gs_i.is_point())
      x_points_non_redundant_in_y.set(i);
  }

  // If `x' is included into `y', the upper bound `y' is exact.
  if (num_x_gs_non_redundant_in_y == 0) {
    *this = y;
    return true;
  }

  // ... and vice versa, generators of `y' that are non-redundant in `x'.
  Bit_Row y_gs_non_redundant_in_x;
  Bit_Row y_points_non_redundant_in_x;
  Bit_Row y_closure_points;
  dimension_type num_y_gs_non_redundant_in_x = 0;
  for (dimension_type i = y_gs_num_rows; i-- > 0; ) {
    const Generator& y_gs_i = y_gs[i];
    if (y_gs_i.is_closure_point())
      y_closure_points.set(i);
    if (x.relation_with(y_gs_i).implies(Poly_Gen_Relation::subsumes()))
      continue;
    y_gs_non_redundant_in_x.set(i);
    ++num_y_gs_non_redundant_in_x;
    if (y_gs_i.is_point())
      y_points_non_redundant_in_x.set(i);
  }

  // If `y' is included into `x', the upper bound `x' is exact.
  if (num_y_gs_non_redundant_in_x == 0)
    return true;

  Bit_Row x_non_points_non_redundant_in_y;
  x_non_points_non_redundant_in_y
    .difference_assign(x_gs_non_redundant_in_y,
                       x_points_non_redundant_in_y);

  const Constraint_System& x_cs = x.con_sys;
  const Constraint_System& y_cs = y.con_sys;
  const dimension_type x_cs_num_rows = x_cs.num_rows();
  const dimension_type y_cs_num_rows = y_cs.num_rows();

  // Filter away the points of `x_gs' that would be redundant
  // in the topological closure of `y'.
  Bit_Row x_points_non_redundant_in_y_closure;
  for (dimension_type i = x_points_non_redundant_in_y.first();
       i != C_Integer<unsigned long>::max;
       i = x_points_non_redundant_in_y.next(i)) {
    const Generator& x_p = x_gs[i];
    PPL_ASSERT(x_p.is_point());
    // NOTE: we cannot use Constraint_System::relation_with()
    // as we need to treat strict inequalities as if they were nonstrict.
    for (dimension_type j = y_cs_num_rows; j-- > 0; ) {
      const Constraint& y_c = y_cs[j];
      const int sp_sign = Scalar_Products::reduced_sign(y_c, x_p);
      if (sp_sign < 0 || (y_c.is_equality() && sp_sign > 0)) {
        x_points_non_redundant_in_y_closure.set(i);
        break;
      }
    }
  }

  // Make sure the saturation matrix for `x' is up to date.
  // Any sat matrix would do: we choose `sat_g' because it matches
  // the two nested loops (constraints on rows and generators on columns).
  if (!x.sat_g_is_up_to_date())
    x.update_sat_g();
  const Bit_Matrix& x_sat = x.sat_g;

  Bit_Row x_cs_condition_3;
  Bit_Row x_gs_condition_3;
  Bit_Row all_ones;
  all_ones.set_until(x_gs_num_rows);
  Bit_Row saturators;
  Bit_Row tmp_set;
  for (dimension_type i = x_cs_num_rows; i-- > 0; ) {
    const Constraint& x_c = x_cs[i];
    // Skip constraint if it is not violated by `y'.
    if (y.relation_with(x_c).implies(Poly_Con_Relation::is_included()))
      continue;
    saturators.difference_assign(all_ones, x_sat[i]);
    // Check condition 1.
    tmp_set.intersection_assign(x_non_points_non_redundant_in_y, saturators);
    if (!tmp_set.empty())
      return false;
    if (x_c.is_strict_inequality()) {
      // Postpone check for condition 3.
      x_cs_condition_3.set(i);
      tmp_set.intersection_assign(x_closure_points, saturators);
      x_gs_condition_3.union_assign(x_gs_condition_3, tmp_set);
    }
    else {
      // Check condition 2.
      tmp_set.intersection_assign(x_points_non_redundant_in_y_closure,
                                  saturators);
      if (!tmp_set.empty())
        return false;
    }
  }

  // Now exchange the roles of `x' and `y'
  // (the statement of the NNC theorem in BHZ09 is symmetric).

  Bit_Row y_non_points_non_redundant_in_x;
  y_non_points_non_redundant_in_x
    .difference_assign(y_gs_non_redundant_in_x,
                       y_points_non_redundant_in_x);

  // Filter away the points of `y_gs' that would be redundant
  // in the topological closure of `x'.
  Bit_Row y_points_non_redundant_in_x_closure;
  for (dimension_type i = y_points_non_redundant_in_x.first();
       i != C_Integer<unsigned long>::max;
       i = y_points_non_redundant_in_x.next(i)) {
    const Generator& y_p = y_gs[i];
    PPL_ASSERT(y_p.is_point());
    // NOTE: we cannot use Constraint_System::relation_with()
    // as we need to treat strict inequalities as if they were nonstrict.
    for (dimension_type j = x_cs_num_rows; j-- > 0; ) {
      const Constraint& x_c = x_cs[j];
      const int sp_sign = Scalar_Products::reduced_sign(x_c, y_p);
      if (sp_sign < 0 || (x_c.is_equality() && sp_sign > 0)) {
        y_points_non_redundant_in_x_closure.set(i);
        break;
      }
    }
  }

  // Make sure the saturation matrix `sat_g' for `y' is up to date.
  if (!y.sat_g_is_up_to_date())
    y.update_sat_g();
  const Bit_Matrix& y_sat = y.sat_g;

  Bit_Row y_cs_condition_3;
  Bit_Row y_gs_condition_3;
  all_ones.clear();
  all_ones.set_until(y_gs_num_rows);
  for (dimension_type i = y_cs_num_rows; i-- > 0; ) {
    const Constraint& y_c = y_cs[i];
    // Skip constraint if it is not violated by `x'.
    if (x.relation_with(y_c).implies(Poly_Con_Relation::is_included()))
      continue;
    saturators.difference_assign(all_ones, y_sat[i]);
    // Check condition 1.
    tmp_set.intersection_assign(y_non_points_non_redundant_in_x, saturators);
    if (!tmp_set.empty())
      return false;
    if (y_c.is_strict_inequality()) {
      // Postpone check for condition 3.
      y_cs_condition_3.set(i);
      tmp_set.intersection_assign(y_closure_points, saturators);
      y_gs_condition_3.union_assign(y_gs_condition_3, tmp_set);
    }
    else {
      // Check condition 2.
      tmp_set.intersection_assign(y_points_non_redundant_in_x_closure,
                                  saturators);
      if (!tmp_set.empty())
        return false;
    }
  }

  // Now considering condition 3.

  if (x_cs_condition_3.empty() && y_cs_condition_3.empty()) {
    // No test for condition 3 is needed.
    // The hull is exact: compute it.
    for (dimension_type j = y_gs_num_rows; j-- > 0; )
      if (y_gs_non_redundant_in_x[j])
        add_generator(y_gs[j]);
    return true;
  }

  // We have anyway to compute the upper bound and its constraints too.
  Polyhedron ub(x);
  for (dimension_type j = y_gs_num_rows; j-- > 0; )
    if (y_gs_non_redundant_in_x[j])
      ub.add_generator(y_gs[j]);
  (void) ub.minimize();
  PPL_ASSERT(!ub.is_empty());

  // NOTE: the following computation of x_gs_condition_3_not_in_y
  // (resp., y_gs_condition_3_not_in_x) is not required for correctness.
  // It is done so as to later apply a speculative test
  // (i.e., a non-conclusive but computationally lighter test).

  // Filter away from `x_gs_condition_3' those closure points
  // that, when considered as points, would belong to `y',
  // i.e., those that violate no strict constraint in `y_cs'.
  Bit_Row x_gs_condition_3_not_in_y;
  for (dimension_type i = y_cs_num_rows; i-- > 0; ) {
    const Constraint& y_c = y_cs[i];
    if (y_c.is_strict_inequality()) {
      for (dimension_type j = x_gs_condition_3.first();
           j != C_Integer<unsigned long>::max; j = x_gs_condition_3.next(j)) {
        const Generator& x_cp = x_gs[j];
        PPL_ASSERT(x_cp.is_closure_point());
        const int sp_sign = Scalar_Products::reduced_sign(y_c, x_cp);
        PPL_ASSERT(sp_sign >= 0);
        if (sp_sign == 0) {
          x_gs_condition_3.clear(j);
          x_gs_condition_3_not_in_y.set(j);
        }
      }
      if (x_gs_condition_3.empty())
        break;
    }
  }
  // Symmetrically, filter away from `y_gs_condition_3' those
  // closure points that, when considered as points, would belong to `x',
  // i.e., those that violate no strict constraint in `x_cs'.
  Bit_Row y_gs_condition_3_not_in_x;
  for (dimension_type i = x_cs_num_rows; i-- > 0; ) {
    if (x_cs[i].is_strict_inequality()) {
      const Constraint& x_c = x_cs[i];
      for (dimension_type j = y_gs_condition_3.first();
           j != C_Integer<unsigned long>::max; j = y_gs_condition_3.next(j)) {
        const Generator& y_cp = y_gs[j];
        PPL_ASSERT(y_cp.is_closure_point());
        const int sp_sign = Scalar_Products::reduced_sign(x_c, y_cp);
        PPL_ASSERT(sp_sign >= 0);
        if (sp_sign == 0) {
          y_gs_condition_3.clear(j);
          y_gs_condition_3_not_in_x.set(j);
        }
      }
      if (y_gs_condition_3.empty())
        break;
    }
  }

  // NOTE: here we apply the speculative test.
  // Check if there exists a closure point in `x_gs_condition_3_not_in_y'
  // or `y_gs_condition_3_not_in_x' that belongs (as point) to the hull.
  // If so, the hull is not exact.
  const Constraint_System& ub_cs = ub.constraints();
  for (dimension_type i = ub_cs.num_rows(); i-- > 0; ) {
    if (ub_cs[i].is_strict_inequality()) {
      const Constraint& ub_c = ub_cs[i];
      for (dimension_type j = x_gs_condition_3_not_in_y.first();
           j != C_Integer<unsigned long>::max;
           j = x_gs_condition_3_not_in_y.next(j)) {
        const Generator& x_cp = x_gs[j];
        PPL_ASSERT(x_cp.is_closure_point());
        const int sp_sign = Scalar_Products::reduced_sign(ub_c, x_cp);
        PPL_ASSERT(sp_sign >= 0);
        if (sp_sign == 0)
          x_gs_condition_3_not_in_y.clear(j);
      }
      for (dimension_type j = y_gs_condition_3_not_in_x.first();
           j != C_Integer<unsigned long>::max;
           j = y_gs_condition_3_not_in_x.next(j)) {
        const Generator& y_cp = y_gs[j];
        PPL_ASSERT(y_cp.is_closure_point());
        const int sp_sign = Scalar_Products::reduced_sign(ub_c, y_cp);
        PPL_ASSERT(sp_sign >= 0);
        if (sp_sign == 0)
          y_gs_condition_3_not_in_x.clear(j);
      }
    }
  }

  if (!(x_gs_condition_3_not_in_y.empty()
        && y_gs_condition_3_not_in_x.empty()))
    // There exist a closure point satisfying condition 3,
    // hence the hull is not exact.
    return false;

  // The speculative test was not successful:
  // apply the expensive (but conclusive) test for condition 3.

  // Consider strict inequalities in `x' violated by `y'.
  for (dimension_type i = x_cs_condition_3.first();
       i != C_Integer<unsigned long>::max; i = x_cs_condition_3.next(i)) {
    const Constraint& x_cs_i = x_cs[i];
    PPL_ASSERT(x_cs_i.is_strict_inequality());
    // Build the equality constraint induced by x_cs_i.
    Constraint eq_i(Linear_Expression(x_cs_i) == 0);
    PPL_ASSERT(!(ub.relation_with(eq_i)
                 .implies(Poly_Con_Relation::is_disjoint())));
    Polyhedron ub_inters_hyperplane(ub);
    ub_inters_hyperplane.add_constraint(eq_i);
    Polyhedron y_inters_hyperplane(y);
    y_inters_hyperplane.add_constraint(eq_i);
    if (!y_inters_hyperplane.contains(ub_inters_hyperplane))
      // The hull is not exact.
      return false;
  }

  // Consider strict inequalities in `y' violated by `x'.
  for (dimension_type i = y_cs_condition_3.first();
       i != C_Integer<unsigned long>::max; i = y_cs_condition_3.next(i)) {
    const Constraint& y_cs_i = y_cs[i];
    PPL_ASSERT(y_cs_i.is_strict_inequality());
    // Build the equality constraint induced by y_cs_i.
    Constraint eq_i(Linear_Expression(y_cs_i) == 0);
    PPL_ASSERT(!(ub.relation_with(eq_i)
                 .implies(Poly_Con_Relation::is_disjoint())));
    Polyhedron ub_inters_hyperplane(ub);
    ub_inters_hyperplane.add_constraint(eq_i);
    Polyhedron x_inters_hyperplane(x);
    x_inters_hyperplane.add_constraint(eq_i);
    if (!x_inters_hyperplane.contains(ub_inters_hyperplane))
      // The hull is not exact.
      return false;
  }

  // The hull is exact.
  m_swap(ub);
  return true;
}

bool Parma_Polyhedra_Library::Polyhedron::BHZ09_poly_hull_assign_if_exact ( const Polyhedron & y )

protected

Definition at line 1422 of file Polyhedron_nonpublic.cc.

References BHZ09_C_poly_hull_assign_if_exact(), BHZ09_NNC_poly_hull_assign_if_exact(), is_empty(), is_necessarily_closed(), marked_empty(), PPL_ASSERT, space_dim, topology(), and upper_bound_assign().

                                                                  {
  Polyhedron& x = *this;

  // Private method: the caller must ensure the following.
  PPL_ASSERT(x.topology() == y.topology());
  PPL_ASSERT(x.space_dim == y.space_dim);

  // The zero-dim case is trivial.
  if (x.space_dim == 0) {
    x.upper_bound_assign(y);
    return true;
  }

  // If `x' or `y' are (known to be) empty, the upper bound is exact.
  if (x.marked_empty()) {
    x = y;
    return true;
  }
  else if (y.is_empty())
    return true;
  else if (x.is_empty()) {
    x = y;
    return true;
  }

  if (x.is_necessarily_closed())
    return x.BHZ09_C_poly_hull_assign_if_exact(y);
  else
    return x.BHZ09_NNC_poly_hull_assign_if_exact(y);
}

void Parma_Polyhedra_Library::Polyhedron::bounded_affine_image	(	Variable	var,
		const Linear_Expression &	lb_expr,
		const Linear_Expression &	ub_expr,
		Coefficient_traits::const_reference	denominator = `Coefficient_one()`
	)

Assigns to *this the image of *this with respect to the bounded affine relation $\frac{\mathrm{lb\_expr}}{\mathrm{denominator}} \leq \mathrm{var}' \leq \frac{\mathrm{ub\_expr}}{\mathrm{denominator}}$ .

Parameters:

var	The variable updated by the affine relation;
lb_expr	The numerator of the lower bounding affine expression;
ub_expr	The numerator of the upper bounding affine expression;
denominator	The (common) denominator for the lower and upper bounding affine expressions (optional argument with default value 1).

Exceptions:

std::invalid_argument Thrown if denominator is zero or if lb_expr (resp., ub_expr) and *this are dimension-incompatible or if var is not a space dimension of *this.

Definition at line 2716 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Linear_Expression::coefficient(), Parma_Polyhedra_Library::GREATER_OR_EQUAL, Parma_Polyhedra_Library::LESS_OR_EQUAL, PPL_ASSERT_HEAVY, Parma_Polyhedra_Library::Variable::space_dimension(), and Parma_Polyhedra_Library::Linear_Expression::space_dimension().

Referenced by affine_form_image().

                                                                    {
  // The denominator cannot be zero.
  if (denominator == 0)
    throw_invalid_argument("bounded_affine_image(v, lb, ub, d)", "d == 0");

  // Dimension-compatibility checks.
  // `var' should be one of the dimensions of the polyhedron.
  const dimension_type var_space_dim = var.space_dimension();
  if (space_dim < var_space_dim)
    throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
                                 "v", var);
  // The dimension of `lb_expr' and `ub_expr' should not be
  // greater than the dimension of `*this'.
  const dimension_type lb_space_dim = lb_expr.space_dimension();
  if (space_dim < lb_space_dim)
    throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
                                 "lb", lb_expr);
  const dimension_type ub_space_dim = ub_expr.space_dimension();
  if (space_dim < ub_space_dim)
    throw_dimension_incompatible("bounded_affine_image(v, lb, ub, d)",
                                 "ub", ub_expr);

  // Any image of an empty polyhedron is empty.
  if (marked_empty())
    return;

  // Check whether `var' occurs in `lb_expr' and/or `ub_expr'.
  if (lb_expr.coefficient(var) == 0) {
    // Here `var' may only occur in `ub_expr'.
    generalized_affine_image(var,
                             LESS_OR_EQUAL,
                             ub_expr,
                             denominator);
    if (denominator > 0)
      refine_no_check(lb_expr <= denominator*var);
    else
      refine_no_check(denominator*var <= lb_expr);
  }
  else if (ub_expr.coefficient(var) == 0) {
    // Here `var' only occurs in `lb_expr'.
    generalized_affine_image(var,
                             GREATER_OR_EQUAL,
                             lb_expr,
                             denominator);
    if (denominator > 0)
      refine_no_check(denominator*var <= ub_expr);
    else
      refine_no_check(ub_expr <= denominator*var);
  }
  else {
    // Here `var' occurs in both `lb_expr' and `ub_expr'.
    // To ease the computation, we add an additional dimension.
    const Variable new_var(space_dim);
    add_space_dimensions_and_embed(1);
    // Constrain the new dimension to be equal to `ub_expr'.
    refine_no_check(denominator*new_var == ub_expr);
    // Apply the affine lower bound.
    generalized_affine_image(var,
                             GREATER_OR_EQUAL,
                             lb_expr,
                             denominator);
    if (!marked_empty())
      // Now apply the affine upper bound, as recorded in `new_var'.
      refine_no_check(new_var >= var);
    // Remove the temporarily added dimension.
    remove_higher_space_dimensions(space_dim-1);
  }
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::bounded_affine_preimage	(	Variable	var,
		const Linear_Expression &	lb_expr,
		const Linear_Expression &	ub_expr,
		Coefficient_traits::const_reference	denominator = `Coefficient_one()`
	)

Assigns to *this the preimage of *this with respect to the bounded affine relation $\frac{\mathrm{lb\_expr}}{\mathrm{denominator}} \leq \mathrm{var}' \leq \frac{\mathrm{ub\_expr}}{\mathrm{denominator}}$ .

Parameters:

var	The variable updated by the affine relation;
lb_expr	The numerator of the lower bounding affine expression;
ub_expr	The numerator of the upper bounding affine expression;
denominator	The (common) denominator for the lower and upper bounding affine expressions (optional argument with default value 1).

Exceptions:

std::invalid_argument Thrown if denominator is zero or if lb_expr (resp., ub_expr) and *this are dimension-incompatible or if var is not a space dimension of *this.

Definition at line 2791 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Linear_Expression::coefficient(), PPL_ASSERT_HEAVY, Parma_Polyhedra_Library::Variable::space_dimension(), and Parma_Polyhedra_Library::Linear_Expression::space_dimension().

                                                                       {
  // The denominator cannot be zero.
  if (denominator == 0)
    throw_invalid_argument("bounded_affine_preimage(v, lb, ub, d)", "d == 0");

  // Dimension-compatibility checks.
  // `var' should be one of the dimensions of the polyhedron.
  const dimension_type var_space_dim = var.space_dimension();
  if (space_dim < var_space_dim)
    throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
                                 "v", var);
  // The dimension of `lb_expr' and `ub_expr' should not be
  // greater than the dimension of `*this'.
  const dimension_type lb_space_dim = lb_expr.space_dimension();
  if (space_dim < lb_space_dim)
    throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
                                 "lb", lb_expr);
  const dimension_type ub_space_dim = ub_expr.space_dimension();
  if (space_dim < ub_space_dim)
    throw_dimension_incompatible("bounded_affine_preimage(v, lb, ub, d)",
                                 "ub", ub_expr);

  // Any preimage of an empty polyhedron is empty.
  if (marked_empty())
    return;

  // Check whether `var' occurs in neither `lb_expr' nor `ub_expr'.
  if (lb_expr.coefficient(var) == 0 && ub_expr.coefficient(var) == 0) {
    if (denominator > 0) {
      refine_no_check(lb_expr <= denominator*var);
      refine_no_check(denominator*var <= ub_expr);
    }
    else {
      refine_no_check(ub_expr <= denominator*var);
      refine_no_check(denominator*var <= lb_expr);
    }
    unconstrain(var);
  }
  else {
    // Here `var' occurs in `lb_expr' or `ub_expr'.
    // To ease the computation, add an additional dimension.
    const Variable new_var(space_dim);
    add_space_dimensions_and_embed(1);
    // Swap dimensions `var' and `new_var'.
    std::vector<dimension_type> swapping_cycle;
    swapping_cycle.push_back(var_space_dim);
    swapping_cycle.push_back(space_dim);
    swapping_cycle.push_back(0);
    if (constraints_are_up_to_date())
      con_sys.permute_columns(swapping_cycle);
    if (generators_are_up_to_date())
      gen_sys.permute_columns(swapping_cycle);
    // Constrain the new dimension as dictated by `lb_expr' and `ub_expr'.
    // (we force minimization because we will need the generators).
    if (denominator > 0) {
      refine_no_check(lb_expr <= denominator*new_var);
      refine_no_check(denominator*new_var <= ub_expr);
    }
    else {
      refine_no_check(ub_expr <= denominator*new_var);
      refine_no_check(denominator*new_var <= lb_expr);
    }
    // Remove the temporarily added dimension.
    remove_higher_space_dimensions(space_dim-1);
  }
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::bounded_BHRZ03_extrapolation_assign	(	const Polyhedron &	y,
		const Constraint_System &	cs,
		unsigned *	tp = `0`
	)

Assigns to *this the result of computing the bounded extrapolation between *this and y using the BHRZ03-widening operator.

Parameters:

y	A polyhedron that must be contained in `*this`;
cs	The system of constraints used to improve the widened polyhedron;
tp	An optional pointer to an unsigned variable storing the number of available tokens (to be used when applying the widening with tokens delay technique).

Exceptions:

std::invalid_argument Thrown if *this, y and cs are topology-incompatible or dimension-incompatible.

Definition at line 822 of file Polyhedron_widenings.cc.

References Parma_Polyhedra_Library::ANY_COMPLEXITY, Parma_Polyhedra_Library::Box< ITV >::CC76_widening_assign(), and Parma_Polyhedra_Library::Box< ITV >::constraints().

                                                    {
  Rational_Box x_box(*this, ANY_COMPLEXITY);
  Rational_Box y_box(y, ANY_COMPLEXITY);
  x_box.CC76_widening_assign(y_box);
  limited_BHRZ03_extrapolation_assign(y, cs, tp);
  Constraint_System x_box_cs = x_box.constraints();
  add_recycled_constraints(x_box_cs);
}

void Parma_Polyhedra_Library::Polyhedron::bounded_H79_extrapolation_assign	(	const Polyhedron &	y,
		const Constraint_System &	cs,
		unsigned *	tp = `0`
	)

Assigns to *this the result of computing the bounded extrapolation between *this and y using the H79-widening operator.

Parameters:

y	A polyhedron that must be contained in `*this`;
cs	The system of constraints used to improve the widened polyhedron;
tp	An optional pointer to an unsigned variable storing the number of available tokens (to be used when applying the widening with tokens delay technique).

Exceptions:

std::invalid_argument Thrown if *this, y and cs are topology-incompatible or dimension-incompatible.

Definition at line 361 of file Polyhedron_widenings.cc.

References Parma_Polyhedra_Library::ANY_COMPLEXITY, Parma_Polyhedra_Library::Box< ITV >::CC76_widening_assign(), and Parma_Polyhedra_Library::Box< ITV >::constraints().

                                                                {
  Rational_Box x_box(*this, ANY_COMPLEXITY);
  Rational_Box y_box(y, ANY_COMPLEXITY);
  x_box.CC76_widening_assign(y_box);
  limited_H79_extrapolation_assign(y, cs, tp);
  Constraint_System x_box_cs = x_box.constraints();
  add_recycled_constraints(x_box_cs);
}

bool Parma_Polyhedra_Library::Polyhedron::bounds	(	const Linear_Expression &	expr,
		bool	from_above
	)		const

private

Checks if and how expr is bounded in *this.

Returns true if and only if from_above is true and expr is bounded from above in *this, or from_above is false and expr is bounded from below in *this.

Parameters:

expr	The linear expression to test;
from_above	`true` if and only if the boundedness of interest is "from above".

Exceptions:

std::invalid_argument Thrown if expr and *this are dimension-incompatible.

Definition at line 524 of file Polyhedron_nonpublic.cc.

References Parma_Polyhedra_Library::Scalar_Products::homogeneous_sign(), Parma_Polyhedra_Library::Generator::is_line(), Parma_Polyhedra_Library::Generator::is_line_or_ray(), and Parma_Polyhedra_Library::Linear_Expression::space_dimension().

Referenced by bounds_from_above(), and bounds_from_below().

                                                     {
  // The dimension of `expr' should not be greater than the dimension
  // of `*this'.
  const dimension_type expr_space_dim = expr.space_dimension();
  if (space_dim < expr_space_dim)
    throw_dimension_incompatible((from_above
                                  ? "bounds_from_above(e)"
                                  : "bounds_from_below(e)"), "e", expr);

  // A zero-dimensional or empty polyhedron bounds everything.
  if (space_dim == 0
      || marked_empty()
      || (has_pending_constraints() && !process_pending_constraints())
      || (!generators_are_up_to_date() && !update_generators()))
    return true;

  // The polyhedron has updated, possibly pending generators.
  for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
    const Generator& g = gen_sys[i];
    // Only lines and rays in `*this' can cause `expr' to be unbounded.
    if (g.is_line_or_ray()) {
      const int sp_sign = Scalar_Products::homogeneous_sign(expr, g);
      if (sp_sign != 0
          && (g.is_line()
              || (from_above && sp_sign > 0)
              || (!from_above && sp_sign < 0)))
        // `*this' does not bound `expr'.
        return false;
    }
  }
  // No sources of unboundedness have been found for `expr'
  // in the given direction.
  return true;
}

bool Parma_Polyhedra_Library::Polyhedron::bounds_from_above ( const Linear_Expression & expr ) const

inline

Returns true if and only if expr is bounded from above in *this.

Exceptions:

std::invalid_argument Thrown if expr and *this are dimension-incompatible.

Definition at line 312 of file Polyhedron.inlines.hh.

References bounds().

                                                                 {
  return bounds(expr, true);
}

bool Parma_Polyhedra_Library::Polyhedron::bounds_from_below ( const Linear_Expression & expr ) const

inline

Returns true if and only if expr is bounded from below in *this.

Exceptions:

std::invalid_argument Thrown if expr and *this are dimension-incompatible.

Definition at line 317 of file Polyhedron.inlines.hh.

References bounds().

                                                                 {
  return bounds(expr, false);
}

bool Parma_Polyhedra_Library::Polyhedron::can_have_something_pending ( ) const

inlineprivate

Returns true if the polyhedron can have something pending.

Definition at line 179 of file Polyhedron.inlines.hh.

References constraints_are_minimized(), generators_are_minimized(), sat_c_is_up_to_date(), and sat_g_is_up_to_date().

Referenced by intersection_assign(), poly_hull_assign(), and time_elapse_assign().

                                             {
  return constraints_are_minimized()
    && generators_are_minimized()
    && (sat_c_is_up_to_date() || sat_g_is_up_to_date());
}

bool Parma_Polyhedra_Library::Polyhedron::can_recycle_congruence_systems ( )

inlinestatic

Returns false indicating that this domain cannot recycle congruences.

Definition at line 124 of file Polyhedron.inlines.hh.

                                           {
  return false;
}

bool Parma_Polyhedra_Library::Polyhedron::can_recycle_constraint_systems ( )

inlinestatic

Returns true indicating that this domain has methods that can recycle constraints.

Definition at line 119 of file Polyhedron.inlines.hh.

                                           {
  return true;
}

template<typename Input >

Input& Parma_Polyhedra_Library::Polyhedron::check_obj_space_dimension_overflow	(	Input &	input,
		const Topology	topol,
		const char *	method,
		const char *	reason
	)

Definition at line 587 of file Polyhedron.templates.hh.

References check_space_dimension_overflow(), and max_space_dimension().

                                                                   {
  check_space_dimension_overflow(input.space_dimension(),
                                 max_space_dimension(),
                                 topol, method, reason);
  return input;
}

template<typename Object >

static Object& Parma_Polyhedra_Library::Polyhedron::check_obj_space_dimension_overflow	(	Object &	input,
		Topology	topol,
		const char *	method,
		const char *	reason
	)

staticprotected

PPL::dimension_type Parma_Polyhedra_Library::Polyhedron::check_space_dimension_overflow	(	dimension_type	dim,
		dimension_type	max,
		const Topology	topol,
		const char *	method,
		const char *	reason
	)

staticprotected

Definition at line 2424 of file Polyhedron_nonpublic.cc.

References Parma_Polyhedra_Library::check_space_dimension_overflow(), and Parma_Polyhedra_Library::NECESSARILY_CLOSED.

Referenced by check_obj_space_dimension_overflow().

                                                                    {
  const char* domain = (topol == NECESSARILY_CLOSED)
    ? "PPL::C_Polyhedron::" : "PPL::NNC_Polyhedron::";
  return PPL::check_space_dimension_overflow(dim, max, domain, method, reason);
}

PPL::dimension_type Parma_Polyhedra_Library::Polyhedron::check_space_dimension_overflow	(	dimension_type	dim,
		const Topology	topol,
		const char *	method,
		const char *	reason
	)

staticprotected

Definition at line 2435 of file Polyhedron_nonpublic.cc.

References Parma_Polyhedra_Library::check_space_dimension_overflow(), and max_space_dimension().

                                                                    {
  return check_space_dimension_overflow(dim, Polyhedron::max_space_dimension(),
                                        topol, method, reason);
}

void Parma_Polyhedra_Library::Polyhedron::clear_constraints_minimized ( )

inlineprivate

Sets status to express that constraints are no longer minimized.

Definition at line 245 of file Polyhedron.inlines.hh.

References Parma_Polyhedra_Library::Polyhedron::Status::reset_c_minimized(), and status.

Referenced by clear_constraints_up_to_date(), intersection_assign(), and remove_pending_to_obtain_constraints().

                                        {
  status.reset_c_minimized();
}

void Parma_Polyhedra_Library::Polyhedron::clear_constraints_up_to_date ( )

inlineprivate

Sets status to express that constraints are no longer up-to-date.

This also implies that they are neither minimized and both saturation matrices are no longer meaningful.

Definition at line 277 of file Polyhedron.inlines.hh.

References clear_constraints_minimized(), clear_pending_constraints(), clear_sat_c_up_to_date(), clear_sat_g_up_to_date(), Parma_Polyhedra_Library::Polyhedron::Status::reset_c_up_to_date(), and status.

Referenced by poly_hull_assign(), remove_pending_to_obtain_generators(), strongly_minimize_generators(), and time_elapse_assign().

                                         {
  clear_pending_constraints();
  clear_constraints_minimized();
  clear_sat_c_up_to_date();
  clear_sat_g_up_to_date();
  status.reset_c_up_to_date();
  // Can get rid of con_sys here.
}

void Parma_Polyhedra_Library::Polyhedron::clear_empty ( )

inlineprivate

Clears the status flag indicating that the polyhedron is empty.

Definition at line 240 of file Polyhedron.inlines.hh.

References Parma_Polyhedra_Library::Polyhedron::Status::reset_empty(), and status.

                        {
  status.reset_empty();
}

void Parma_Polyhedra_Library::Polyhedron::clear_generators_minimized ( )

inlineprivate

Sets status to express that generators are no longer minimized.

Definition at line 250 of file Polyhedron.inlines.hh.

References Parma_Polyhedra_Library::Polyhedron::Status::reset_g_minimized(), and status.

Referenced by clear_generators_up_to_date(), poly_hull_assign(), remove_pending_to_obtain_generators(), and time_elapse_assign().

                                       {
  status.reset_g_minimized();
}

void Parma_Polyhedra_Library::Polyhedron::clear_generators_up_to_date ( )

inlineprivate

Sets status to express that generators are no longer up-to-date.

This also implies that they are neither minimized and both saturation matrices are no longer meaningful.

Definition at line 287 of file Polyhedron.inlines.hh.

References clear_generators_minimized(), clear_pending_generators(), clear_sat_c_up_to_date(), clear_sat_g_up_to_date(), Parma_Polyhedra_Library::Polyhedron::Status::reset_g_up_to_date(), and status.

Referenced by intersection_assign(), remove_pending_to_obtain_constraints(), and strongly_minimize_constraints().

                                        {
  clear_pending_generators();
  clear_generators_minimized();
  clear_sat_c_up_to_date();
  clear_sat_g_up_to_date();
  status.reset_g_up_to_date();
  // Can get rid of gen_sys here.
}

void Parma_Polyhedra_Library::Polyhedron::clear_pending_constraints ( )

inlineprivate

Sets status to express that there are no longer pending constraints.

Definition at line 255 of file Polyhedron.inlines.hh.

References Parma_Polyhedra_Library::Polyhedron::Status::reset_c_pending(), and status.

Referenced by clear_constraints_up_to_date(), process_pending_constraints(), and remove_pending_to_obtain_constraints().

                                      {
  status.reset_c_pending();
}

void Parma_Polyhedra_Library::Polyhedron::clear_pending_generators ( )

inlineprivate

Sets status to express that there are no longer pending generators.

Definition at line 260 of file Polyhedron.inlines.hh.

References Parma_Polyhedra_Library::Polyhedron::Status::reset_g_pending(), and status.

Referenced by clear_generators_up_to_date(), process_pending_generators(), and remove_pending_to_obtain_generators().

                                     {
  status.reset_g_pending();
}

void Parma_Polyhedra_Library::Polyhedron::clear_sat_c_up_to_date ( )

inlineprivate

Sets status to express that sat_c is no longer up-to-date.

Definition at line 265 of file Polyhedron.inlines.hh.

References Parma_Polyhedra_Library::Polyhedron::Status::reset_sat_c_up_to_date(), and status.

Referenced by clear_constraints_up_to_date(), clear_generators_up_to_date(), obtain_sorted_constraints(), process_pending_generators(), and update_generators().

                                   {
  status.reset_sat_c_up_to_date();
  // Can get rid of sat_c here.
}

void Parma_Polyhedra_Library::Polyhedron::clear_sat_g_up_to_date ( )

inlineprivate

Sets status to express that sat_g is no longer up-to-date.

Definition at line 271 of file Polyhedron.inlines.hh.

References Parma_Polyhedra_Library::Polyhedron::Status::reset_sat_g_up_to_date(), and status.

Referenced by clear_constraints_up_to_date(), clear_generators_up_to_date(), obtain_sorted_generators(), process_pending_constraints(), and update_constraints().

                                   {
  status.reset_sat_g_up_to_date();
  // Can get rid of sat_g here.
}

void Parma_Polyhedra_Library::Polyhedron::concatenate_assign ( const Polyhedron & y )

Assigns to *this the concatenation of *this and y, taken in this order.

Exceptions:

std::invalid_argument	Thrown if `*this` and `y` are topology-incompatible.
std::length_error	Thrown if the concatenation would cause the vector space to exceed dimension `max_space_dimension()`.

Definition at line 290 of file Polyhedron_chdims.cc.

                                                     {
  if (topology() != y.topology())
    throw_topology_incompatible("concatenate_assign(y)", "y", y);

  // The space dimension of the resulting polyhedron should not
  // overflow the maximum allowed space dimension.
  const dimension_type added_columns = y.space_dim;
  check_space_dimension_overflow(added_columns,
                                 max_space_dimension() - space_dimension(),
                                 topology(),
                                 "concatenate_assign(y)",
                                 "concatenation exceeds the maximum "
                                 "allowed space dimension");

  // If `*this' or `y' are empty polyhedra, it is sufficient to adjust
  // the dimension of the space.
  if (marked_empty() || y.marked_empty()) {
    space_dim += added_columns;
    set_empty();
    return;
  }

  // If `y' is a non-empty 0-dim space polyhedron, the result is `*this'.
  if (added_columns == 0)
    return;

  // If `*this' is a non-empty 0-dim space polyhedron, the result is `y'.
  if (space_dim == 0) {
    *this = y;
    return;
  }

  // TODO: this implementation is just an executable specification.
  Constraint_System cs = y.constraints();

  // The constraints of `x' (possibly with pending rows) are required.
  if (has_pending_generators())
    process_pending_generators();
  else if (!constraints_are_up_to_date())
    update_constraints();

  // The matrix for the new system of constraints is obtained
  // by leaving the old system of constraints in the upper left-hand side
  // and placing the constraints of `cs' in the lower right-hand side.
  // NOTE: here topologies agree, whereas dimensions may not agree.
  dimension_type old_num_rows = con_sys.num_rows();
  dimension_type old_num_columns = con_sys.num_columns();
  dimension_type added_rows = cs.num_rows();

  // We already dealt with the cases of an empty or zero-dim `y' polyhedron;
  // also, `cs' contains the low-level constraints, at least.
  PPL_ASSERT(added_rows > 0 && added_columns > 0);

  con_sys.add_zero_rows_and_columns(added_rows, added_columns,
                                    Linear_Row::Flags(topology(),
                                                      Linear_Row::RAY_OR_POINT_OR_INEQUALITY));
  // Move the epsilon coefficient to the last column, if needed.
  if (!is_necessarily_closed())
    con_sys.swap_columns(old_num_columns - 1,
                         old_num_columns - 1 + added_columns);
  dimension_type cs_num_columns = cs.num_columns();
  // Steal the constraints from `cs' and put them in `con_sys'
  // using the right displacement for coefficients.
  for (dimension_type i = added_rows; i-- > 0; ) {
    Constraint& c_old = cs[i];
    Constraint& c_new = con_sys[old_num_rows + i];
    // Method `add_zero_rows_and_columns', by default, added inequalities.
    if (c_old.is_equality())
      c_new.set_is_equality();
    // The inhomogeneous term is not displaced.
    using std::swap;
    swap(c_new[0], c_old[0]);
    // All homogeneous terms (included the epsilon coefficient,
    // if present) are displaced by `space_dim' columns.
    for (dimension_type j = 1; j < cs_num_columns; ++j)
      swap(c_old[j], c_new[space_dim + j]);
  }

  if (can_have_something_pending()) {
    // If `*this' can support pending constraints, then, since we have
    // resized the system of constraints, we must also add to the generator
    // system those lines corresponding to the newly added dimensions,
    // because the non-pending parts of `con_sys' and `gen_sys' must still
    // be a DD pair in minimal form.
    gen_sys.add_rows_and_columns(added_columns);
    gen_sys.set_sorted(false);
    if (!is_necessarily_closed())
      gen_sys.swap_columns(old_num_columns - 1,
                           old_num_columns - 1 + added_columns);
    // The added lines are not pending.
    gen_sys.unset_pending_rows();
    // Since we added new lines at the beginning of `x.gen_sys',
    // we also have to adjust the saturation matrix `sat_c'.
    // FIXME: if `sat_c' is not up-to-date, could not we directly update
    // `sat_g' by resizing it and shifting its columns?
    if (!sat_c_is_up_to_date()) {
      sat_c.transpose_assign(sat_g);
      set_sat_c_up_to_date();
    }
    clear_sat_g_up_to_date();
    sat_c.resize(sat_c.num_rows() + added_columns, sat_c.num_columns());
    // The old saturation rows are copied at the end of the matrix.
    // The newly introduced lines saturate all the non-pending constraints,
    // thus their saturation rows are made of zeroes.
    using std::swap;
    for (dimension_type i = sat_c.num_rows() - added_columns; i-- > 0; )
      swap(sat_c[i], sat_c[i+added_columns]);
    // Since `added_rows > 0', we now have pending constraints.
    set_constraints_pending();
  }
  else {
    // The polyhedron cannot have pending constraints.
    con_sys.unset_pending_rows();
#if BE_LAZY
    con_sys.set_sorted(false);
#else
    con_sys.sort_rows();
#endif
    clear_constraints_minimized();
    clear_generators_up_to_date();
    clear_sat_g_up_to_date();
    clear_sat_c_up_to_date();
  }
  // Update space dimension.
  space_dim += added_columns;

  // The system of constraints may be unsatisfiable,
  // thus we do not check for satisfiability.
  PPL_ASSERT_HEAVY(OK());
}

Congruence_System Parma_Polyhedra_Library::Polyhedron::congruences ( ) const

inline

Returns a system of (equality) congruences satisfied by *this.

Definition at line 363 of file Polyhedron.inlines.hh.

References minimized_constraints().

                              {
  return Congruence_System(minimized_constraints());
}

bool Parma_Polyhedra_Library::Polyhedron::constrains ( Variable var ) const

Returns true if and only if var is constrained in *this.

Exceptions:

std::invalid_argument Thrown if var is not a space dimension of *this.

Definition at line 678 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Linear_Row::coefficient(), Parma_Polyhedra_Library::Variable::id(), Parma_Polyhedra_Library::Generator::is_line(), Parma_Polyhedra_Library::Generator::is_line_or_ray(), Parma_Polyhedra_Library::Boundary_NS::sgn(), and Parma_Polyhedra_Library::Variable::space_dimension().

                                                  {
  // `var' should be one of the dimensions of the polyhedron.
  const dimension_type var_space_dim = var.space_dimension();
  if (space_dim < var_space_dim)
    throw_dimension_incompatible("constrains(v)", "v", var);

  // An empty polyhedron constrains all variables.
  if (marked_empty())
    return true;

  if (generators_are_up_to_date() && !has_pending_constraints()) {
    // Since generators are up-to-date and there are no pending
    // constraints, the generator system (since it is well formed)
    // contains a point.  Hence the polyhedron is not empty.
    if (constraints_are_up_to_date() && !has_pending_generators())
      // Here a variable is constrained if and only if it is
      // syntactically constrained.
      goto syntactic_check;

    if (generators_are_minimized()) {
      // Try a quick, incomplete check for the universe polyhedron:
      // a universe polyhedron constrains no variable.
      // Count the number of non-pending
      // (hence, linearly independent) lines.
      dimension_type num_lines = 0;
      const dimension_type first_pending = gen_sys.first_pending_row();
      for (dimension_type i = first_pending; i-- > 0; )
        if (gen_sys[i].is_line())
          ++num_lines;

      if (num_lines == space_dim)
        return false;
    }

    // Scan generators: perhaps we will find a generator equivalent to
    // line(var) or a pair of generators equivalent to ray(-var) and
    // ray(var).
    bool have_positive_ray = false;
    bool have_negative_ray = false;
    const dimension_type var_id = var.id();
    for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
      const Generator& gen_sys_i = gen_sys[i];
      if (gen_sys_i.is_line_or_ray()) {
        const Linear_Row& row = gen_sys_i;
        const int sign = sgn(row.coefficient(var_id));
        if (sign != 0) {
          for (dimension_type j = space_dim+1; --j > 0; )
            if (j != var_id && row[j] != 0)
              goto next;
          if (gen_sys_i.is_line())
            return true;
          if (sign > 0)
            if (have_negative_ray)
              return true;
            else
              have_positive_ray = true;
          else if (have_positive_ray)
            return true;
          else
            have_negative_ray = true;
        }
      }
    next:
      ;
    }

    // We are still here: at least we know that, since generators are
    // up-to-date and there are no pending constraints, then the
    // generator system (since it is well formed) contains a point.
    // Hence the polyhedron is not empty.
    if (has_pending_generators())
      process_pending_generators();
    else if (!constraints_are_up_to_date())
      update_constraints();
    goto syntactic_check;
  }

  // We must minimize to detect emptiness and obtain constraints.
  if (!minimize())
    return true;

 syntactic_check:
  for (dimension_type i = con_sys.num_rows(); i-- > 0; )
    if (con_sys[i].coefficient(var) != 0)
      return true;
  return false;
}

const PPL::Constraint_System & Parma_Polyhedra_Library::Polyhedron::constraints ( ) const

Returns the system of constraints.

Definition at line 74 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Constraint_System::adjust_topology_and_space_dimension(), PPL_ASSERT, and Parma_Polyhedra_Library::Constraint_System::zero_dim_empty().

Referenced by Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), BHZ09_NNC_poly_hull_assign_if_exact(), Parma_Polyhedra_Library::Box< ITV >::Box(), Parma_Polyhedra_Library::C_Polyhedron::C_Polyhedron(), concatenate_assign(), Parma_Polyhedra_Library::Grid::Grid(), Parma_Polyhedra_Library::NNC_Polyhedron::NNC_Polyhedron(), Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape(), and poly_difference_assign().

                                 {
  if (marked_empty()) {
    // We want `con_sys' to only contain the unsatisfiable constraint
    // of the appropriate dimension.
    if (con_sys.has_no_rows()) {
      // The 0-dim unsatisfiable constraint is extended to
      // the appropriate dimension and then stored in `con_sys'.
      Constraint_System unsat_cs = Constraint_System::zero_dim_empty();
      unsat_cs.adjust_topology_and_space_dimension(topology(), space_dim);
      const_cast<Constraint_System&>(con_sys).m_swap(unsat_cs);
    }
    else {
      // Checking that `con_sys' contains the right thing.
      PPL_ASSERT(con_sys.space_dimension() == space_dim);
      PPL_ASSERT(con_sys.num_rows() == 1);
      PPL_ASSERT(con_sys[0].is_inconsistent());
    }
    return con_sys;
  }

  if (space_dim == 0) {
    // Zero-dimensional universe.
    PPL_ASSERT(con_sys.num_rows() == 0 && con_sys.num_columns() == 0);
    return con_sys;
  }

  // If the polyhedron has pending generators, we process them to obtain
  // the constraints. No processing is needed if the polyhedron has
  // pending constraints.
  if (has_pending_generators())
    process_pending_generators();
  else if (!constraints_are_up_to_date())
    update_constraints();

  // TODO: reconsider whether to really sort constraints at this stage.
#if ENSURE_SORTEDNESS
  // We insist in returning a sorted system of constraints,
  // but sorting is useless if there are pending constraints.
  if (!has_pending_constraints())
    obtain_sorted_constraints();
#endif
  return con_sys;
}

bool Parma_Polyhedra_Library::Polyhedron::constraints_are_minimized ( ) const

inlineprivate

Returns true if the system of constraints is minimized.

Note that only weak minimization is entailed, so that an NNC polyhedron may still have $\epsilon$ -redundant constraints.

Definition at line 144 of file Polyhedron.inlines.hh.

References status, and Parma_Polyhedra_Library::Polyhedron::Status::test_c_minimized().

Referenced by Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), BHRZ03_combining_constraints(), BHRZ03_evolving_points(), BHRZ03_evolving_rays(), can_have_something_pending(), Parma_Polyhedra_Library::Grid::Grid(), is_included_in(), Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape(), quick_equivalence_test(), select_CH78_constraints(), select_H79_constraints(), simplified_constraints(), and simplify_using_context_assign().

                                            {
  return status.test_c_minimized();
}

bool Parma_Polyhedra_Library::Polyhedron::constraints_are_up_to_date ( ) const

inlineprivate

Returns true if the system of constraints is up-to-date.

Definition at line 134 of file Polyhedron.inlines.hh.

References status, and Parma_Polyhedra_Library::Polyhedron::Status::test_c_up_to_date().

Referenced by Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), Parma_Polyhedra_Library::Box< ITV >::Box(), Parma_Polyhedra_Library::Grid::Grid(), H79_widening_assign(), intersection_assign(), is_included_in(), map_space_dimensions(), Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape(), operator=(), Polyhedron(), simplified_constraints(), and simplify_using_context_assign().

                                             {
  return status.test_c_up_to_date();
}

bool Parma_Polyhedra_Library::Polyhedron::contains ( const Polyhedron & y ) const

Returns true if and only if *this contains y.

Exceptions:

std::invalid_argument Thrown if *this and y are topology-incompatible or dimension-incompatible.

Definition at line 3658 of file Polyhedron_public.cc.

References is_empty(), is_included_in(), marked_empty(), quick_equivalence_test(), space_dim, topology(), and TVB_TRUE.

Referenced by BHRZ03_combining_constraints(), BHRZ03_evolving_points(), BHRZ03_evolving_rays(), BHRZ03_widening_assign(), BHZ09_NNC_poly_hull_assign_if_exact(), Parma_Polyhedra_Library::Pointset_Powerset< PSET >::check_containment(), H79_widening_assign(), poly_difference_assign(), poly_hull_assign_if_exact(), and strictly_contains().

                                                 {
  const Polyhedron& x = *this;

  // Topology compatibility check.
  if (x.topology() != y.topology())
    throw_topology_incompatible("contains(y)", "y", y);

  // Dimension-compatibility check.
  if (x.space_dim != y.space_dim)
    throw_dimension_incompatible("contains(y)", "y", y);

  if (y.marked_empty())
    return true;
  else if (x.marked_empty())
    return y.is_empty();
  else if (y.space_dim == 0)
    return true;
  else if (x.quick_equivalence_test(y) == Polyhedron::TVB_TRUE)
    return true;
  else
    return y.is_included_in(x);
}

bool Parma_Polyhedra_Library::Polyhedron::contains_integer_point ( ) const

Returns true if and only if *this contains at least one integer point.

Definition at line 548 of file Polyhedron_public.cc.

                                            {
  // Any empty polyhedron does not contain integer points.
  if (marked_empty())
    return false;

  // A zero-dimensional, universe polyhedron has, by convention, an
  // integer point.
  if (space_dim == 0)
    return true;

  // CHECKME: do we really want to call conversion to check for emptiness?
  if (has_pending_constraints() && !process_pending())
    // Empty again.
    return false;

  // FIXME: do also exploit info regarding rays and lines, if possible.
  // Is any integer point already available?
  PPL_ASSERT(!has_pending_constraints());
  if (generators_are_up_to_date())
    for (dimension_type i = gen_sys.num_rows(); i-- > 0; )
      if (gen_sys[i].is_point() && gen_sys[i].divisor() == 1)
        return true;

  const Constraint_System& cs = constraints();
#if 0 // TEMPORARILY DISABLED.
  MIP_Problem mip(space_dim,
                  cs.begin(), cs.end(),
                  Variables_Set(Variable(0), Variable(space_dim-1)));
#else
  // FIXME: temporary workaround, to be removed as soon as the MIP
  // problem class will correctly and precisely handle
  // ((strict) in-) equality constraints having all integer variables.
  MIP_Problem mip(space_dim);
  mip.add_to_integer_space_dimensions(Variables_Set(Variable(0),
                                                    Variable(space_dim-1)));
  PPL_DIRTY_TEMP_COEFFICIENT(homogeneous_gcd);
  PPL_DIRTY_TEMP_COEFFICIENT(gcd);
  PPL_DIRTY_TEMP(mpq_class, rational_inhomogeneous);
  PPL_DIRTY_TEMP_COEFFICIENT(tightened_inhomogeneous);
  for (Constraint_System::const_iterator cs_i = cs.begin(),
         cs_end = cs.end(); cs_i != cs_end; ++cs_i) {
    const Constraint& c = *cs_i;
    const Constraint::Type c_type = c.type();
    const Coefficient& inhomogeneous = c.inhomogeneous_term();
    if (c_type == Constraint::STRICT_INEQUALITY) {
      // CHECKME: should we change the behavior of Linear_Expression(c) ?
      // Compute the GCD of the coefficients of c
      // (disregarding the inhomogeneous term and the epsilon dimension).
      homogeneous_gcd = 0;
      for (dimension_type i = space_dim; i-- > 0; )
        gcd_assign(homogeneous_gcd,
                   homogeneous_gcd, c.coefficient(Variable(i)));
      if (homogeneous_gcd == 0) {
        // NOTE: since tautological constraints are already filtered away
        // by iterators, here we must have an inconsistent constraint.
        PPL_ASSERT(c.is_inconsistent());
        return false;
      }
      Linear_Expression le;
      for (dimension_type i = space_dim; i-- > 0; )
        le += (c.coefficient(Variable(i)) / homogeneous_gcd) * Variable(i);
      // Add the integer part of `inhomogeneous'.
      le += (inhomogeneous / homogeneous_gcd);
      // Further tighten the constraint if the inhomogeneous term
      // was integer, i.e., if `homogeneous_gcd' divides `inhomogeneous'.
      gcd_assign(gcd, homogeneous_gcd, inhomogeneous);
      if (gcd == homogeneous_gcd)
        le -= 1;
      mip.add_constraint(le >= 0);
    }
    else {
      // Equality or non-strict inequality.
      // If possible, avoid useless gcd computations.
      if (inhomogeneous == 0)
        // The inhomogeneous term cannot be tightened.
        mip.add_constraint(c);
      else {
        // Compute the GCD of the coefficients of c
        // (disregarding the inhomogeneous term)
        // to see whether or not the inhomogeneous term can be tightened.
        homogeneous_gcd = 0;
        for (dimension_type i = space_dim; i-- > 0; )
          gcd_assign(homogeneous_gcd,
                     homogeneous_gcd, c.coefficient(Variable(i)));
        if (homogeneous_gcd == 0) {
          // NOTE: since tautological constraints are already filtered away
          // by iterators, here we must have an inconsistent constraint.
          PPL_ASSERT(c.is_inconsistent());
          return false;
        }
        else if (homogeneous_gcd == 1)
          // The normalized inhomogeneous term is integer:
          // add the constraint as-is.
          mip.add_constraint(c);
        else {
          PPL_ASSERT(homogeneous_gcd > 1);
          // Here the normalized inhomogeneous term is rational:
          // the constraint has to be tightened.
#ifndef NDEBUG
          // `homogeneous_gcd' does not divide `inhomogeneous'.
          // FIXME: add a divisibility test for Coefficient.
          gcd_assign(gcd, homogeneous_gcd, inhomogeneous);
          PPL_ASSERT(gcd == 1);
#endif
          if (c.type() == Constraint::EQUALITY)
            return false;
          // Extract the homogeneous part of the constraint.
          Linear_Expression le = Linear_Expression(c);
          le -= inhomogeneous;
          // Tighten the inhomogeneous term.
          assign_r(rational_inhomogeneous.get_num(),
                   inhomogeneous, ROUND_NOT_NEEDED);
          assign_r(rational_inhomogeneous.get_den(),
                   homogeneous_gcd, ROUND_NOT_NEEDED);
          // Note: canonicalization is not needed (as gcd == 1).
          PPL_ASSERT(is_canonical(rational_inhomogeneous));
          assign_r(tightened_inhomogeneous,
                   rational_inhomogeneous, ROUND_DOWN);
          tightened_inhomogeneous *= homogeneous_gcd;
          le += tightened_inhomogeneous;
          mip.add_constraint(le >= 0);
        }
      }
    }
  }
#endif // TEMPORARY WORKAROUND.
  return mip.is_satisfiable();
}

PPL::dimension_type Parma_Polyhedra_Library::Polyhedron::conversion	(	Linear_System &	source,
		dimension_type	start,
		Linear_System &	dest,
		Bit_Matrix &	sat,
		dimension_type	num_lines_or_equalities
	)

staticprivate

Performs the conversion from constraints to generators and vice versa.

Returns:: The number of lines of the polyhedron or the number of equality constraints in the result of conversion.

Parameters:

source	The system to use to convert `dest:` it may be modified;
start	The index of `source` row from which conversion begin;
dest	The result of the conversion;
sat	The saturation matrix telling us, for each row in `source`, which are the rows of `dest` that satisfy but do not saturate it;
num_lines_or_equalities	The number of rows in the system `dest` that are either lines of the polyhedron (when `dest` is a system of generators) or equality constraints (when `dest` is a system of constraints).

For simplicity, all the following comments assume we are converting a constraint system source to a generator system dest; the comments for the symmetric case can be obtained by duality.

If some of the constraints in source are redundant, they will be removed. This is why the source is not declared to be a constant parameter.

If start is 0, then source is a sorted system; also, dest is a generator system corresponding to an empty constraint system. If otherwise start is greater than 0, then the two sub-systems of source made by the non-pending rows and the pending rows, respectively, are both sorted; also, dest is the generator system corresponding to the non-pending constraints of source.

Independently from the value of start, dest has lines from index 0 to index num_lines_or_equalities - 1 and rays/points from index num_lines_or_equalities to the last of its rows.

Note that here the rows of sat are indexed by rows of dest and its columns are indexed by rows of source.

We know that polyhedra can be represented by both a system of constraints or a system of generators (points, rays and lines) (see Section Representations of Convex Polyhedra). When we have both descriptions for a polyhedron $P$ we have what is called a double description (or DD pair) for $P$ .

Here, the representation system refers to the system $C$ whose rows represent the constraints that characterize $P$ and the generating system, the system $G$ whose rows represent the generators of $P$ . We say that a pair $(C, G)$ of (real) systems is a double description pair if

$C\vect{x} \geq \vect{0} \quad\iff\quad \exists \vect{\lambda} \geq \vect{0} \mathrel{.} \vect{x} = G\vect{\lambda}.$

The term "double description" is quite natural in the sense that such a pair contains two different description of the same object. In fact, if we refer to the cone representation of a polyhedron $P$ and we call $C$ and $G$ the systems of constraints and rays respectively, we have

$P = \{\, \vect{x} \in \Rset^n \mid C\vect{x} \geq \vect{0}\, \} = \{\, \vect{x} \in \Rset^n \mid \vect{x} = G\vect{\lambda} \text{ for some } \vect{\lambda} \geq \vect{0}\, \}.$

Because of the theorem of Minkowski (see Section Further Notation and Terminology), we can say that, given a $m \times n$ representation system $C$ such that $\mathop{\mathrm{rank}}(C) = n = \mathit{dimension of the whole space}$ for a non-empty polyhedron $P$ , it is always possible to find a generating system $G$ for $P$ such that $(C, G)$ is a DD pair. Conversely, Weyl's theorem ensures that, for each generating system $G$ , it is possible to find a representation system $C$ such that is a DD pair.

For efficiency reasons, our representation of polyhedra makes use of a double description. We are thus left with two problems:

given find such that is a DD pair;
given find such that is a DD pair.

Using Farkas' Lemma we can prove that these two problems are computationally equivalent (i.e., linear-time reducible to each other). Farkas' Lemma establishes a fundamental property of vectors in $\Rset^n$ that, in a sense, captures the essence of duality. Consider a matrix $A \in \Rset^{m \times n}$ and let $\{ \vect{a}_1, \ldots, \vect{a}_m \}$ be its set of row vectors. Consider also another vector $\vect{c} \in \Rset^n$ such that, whenever a vector $\vect{y} \in \Rset^n$ has a non-negative projection on the $\vect{a}_i$ 's, it also has a non-negative projection on $\vect{c}$ . The lemma states that $\vect{c}$ has this property if and only if it is in the cone generated by the $\vect{a}_i$ 's. Formally, the lemma states the equivalence of the two following assertions:

$\forall \vect{y} \mathrel{:} (A\vect{y} \geq 0 \implies \langle \vect{y},\vect{c} \rangle \geq 0)$ ;
$\exists \vect{\lambda} \geq \vect{0} \mathrel{.} \vect{c}^\mathrm{T} = \vect{\lambda}^\mathrm{T}A$ .

With this result we can prove that $(C, G)$ is a DD pair if and only if $(G^\mathrm{T}, C^\mathrm{T})$ is a DD pair.

Suppose $(C, G)$ is a DD pair. Thus, for each $x$ of the appropriate dimension, $C\vect{x} \geq \vect{0}$ if and only if $\exists \lambda \geq 0 \mathrel{.} \vect{x} = G\vect{\lambda}$ , which is of course equivalent to $\exists \vect{\lambda} \geq \vect{0} \mathrel{.} \vect{x}^\mathrm{T} = \vect{\lambda}^\mathrm{T}G^\mathrm{T}$ .

First, we assume that $\vect{z}$ is such that $G^\mathrm{T}\vect{z} \geq \vect{0}$ and we will show that $\exists \vect{\mu} \geq \vect{0} \mathrel{.} \vect{z} = C^\mathrm{T}\vect{\mu}$ . Let $\vect{x}$ be such that $C\vect{x} \geq \vect{0}$ . Since $(C, G)$ is a DD pair, this is equivalent to $\exists \vect{\lambda} \geq \vect{0} \mathrel{.} \vect{x}^\mathrm{T} = \vect{\lambda}^\mathrm{T}G^\mathrm{T}$ , which, by Farkas' Lemma is equivalent to $\forall \vect{y} \mathrel{:} (G^\mathrm{T}\vect{y} \geq \vect{0} \implies \langle \vect{y}, \vect{x} \rangle \geq 0)$ . Taking $\vect{y} = \vect{z}$ and recalling our assumption that $G^\mathrm{T}\vect{z} \geq \vect{0}$ we can conclude that $\langle \vect{z}, \vect{x} \rangle \geq 0$ , that is equivalent to $\langle \vect{x}, \vect{z} \rangle \geq 0$ . We have thus established that $\forall \vect{x} \mathrel{:} (C\vect{x} \geq \vect{0} \implies \langle \vect{x}, \vect{z} \rangle \geq 0)$ . By Farkas' Lemma, this is equivalent to $\exists \vect{\mu} \geq \vect{0} \mathrel{.} \vect{z}^\mathrm{T} = \vect{\mu}^\mathrm{T} C$ , which is equivalent to what we wanted to prove, that is, $\exists \vect{\mu} \geq \vect{0} \mathrel{.} \vect{z} = C^\mathrm{T}\vect{\mu}$ .

In order to prove the reverse implication, the following observation turns out to be useful: when $(C, G)$ is a DD pair, $CG \geq 0$ . In fact, let $\vect{e}_j$ be the vector whose components are all $0$ apart from the $j$ -th one, which is $1$ . Clearly $\vect{e}_j \geq \vect{0}$ and, taking $\vect{\lambda} = \vect{e}_j$ and $\vect{x} = G\vect{\lambda} = G \vect{e}_j$ , we have $C\vect{x} = C(G \vect{e}_j) = (CG)\vect{e}_j \geq \vect{0}$ , since $(C, G)$ is a DD pair. Thus, as $(CG)\vect{e}_j$ is the $j$ -th column of $CG$ and since the choice of $j$ was arbitrary, $CG \geq \vect{0}$ .

We now assume that $\vect{z}$ is such that $\exists \vect{\mu} \geq \vect{0} \mathrel{.} \vect{z} = C^\mathrm{T}\vect{\mu}$ and we will prove that $G^\mathrm{T}\vect{z} \geq \vect{0}$ . By Farkas' Lemma, the assumption $\exists \vect{\mu} \geq \vect{0} \mathrel{.} \vect{z}^\mathrm{T} = \vect{\mu}^\mathrm{T}C$ , is equivalent to $\forall \vect{y} \mathrel{:} (C\vect{y} \geq \vect{0} \implies \langle \vect{y}, \vect{z} \rangle \geq 0)$ . If we take $\vect{y} = G\vect{e}_j$ then $C\vect{y} = CG\vect{e}_j \geq 0$ , since $CG \geq \vect{0}$ . So $\langle \vect{y}, \vect{z} \rangle = (\vect{e}_j^\mathrm{T}G^\mathrm{T}) \vect{z} = \vect{e}_j^\mathrm{T}(G^\mathrm{T} \vect{z}) \geq 0$ , that is, the $j$ -th component of $G^\mathrm{T}\vect{z}$ is non-negative. The arbitrary choice of $j$ allows us to conclude that $G^\mathrm{T}\vect{z} \geq \vect{0}$ , as required.

In view of this result, the following exposition assumes, for clarity, that the conversion being performed is from constraints to generators. Thus, even if the roles of source and dest can be interchanged, in the sequel we assume the source system will contain the constraints that represent the polyhedron and the dest system will contain the generator that generates it.

There are some observations that are useful to understand this function:

Observation 1: Let $A$ be a system of constraints that generate the polyhedron $P$ and $\vect{c}$ a new constraint that must be added. Suppose that there is a line $\vect{z}$ that does not saturate the constraint $\vect{c}$ . If we combine the old lines and rays that do not saturate $\vect{c}$ (except $\vect{z}$ ) with $\vect{z}$ such that the new ones saturate $\vect{c}$ , the new lines and rays also saturate the constraints saturated by the old lines and rays.

In fact, if $\vect{y}_1$ is the old generator that does not saturate $\vect{c}$ , $\vect{y}_2$ is the new one such that

$\vect{y}_2 = \lambda \vect{y}_1 + \mu \vect{z}$

and $\vect{c}_1$ is a previous constraint that $\vect{y}_1$ and $\vect{z}$ saturates, we can see

$\langle \vect{c}_1, \vect{y}_2 \rangle = \langle \vect{c}_1, (\lambda \vect{y}_1 + \mu \vect{z}) \rangle = \lambda \langle \vect{c}_1, \vect{y}_1 \rangle + \mu \langle \vect{c}_1, \vect{z} \rangle = 0 + \mu \langle \vect{c}_1, \vect{z} \rangle = \mu \langle \vect{c}_1, \vect{z} \rangle$

and

$\mu \langle \vect{c}_1, \vect{z} \rangle = 0.$

Proposition 1: Let $\vect{r}_1$ and $\vect{r}_2$ be distinct rays of $P$ . Then the following statements are equivalent: a) $\vect{r}_1$ and $\vect{r}_2$ are adjacent extreme rays (see Section Further Notation and Terminology); b) $\vect{r}_1$ and $\vect{r}_2$ are extreme rays and the rank of the system composed by the constraints saturated by both $\vect{r}_1$ and $\vect{r}_2$ is equal to $d - 2$ , where $d$ is the rank of the system of constraints.

In fact, let $F$ be the system of generators that saturate the constraints saturated by both $\vect{r}_1$ and $\vect{r}_2$ . If b) holds, the set $F$ is 2-dimensional and $\vect{r}_1$ and $\vect{r}_2$ generate this set. So, every generator $\vect{x}$ of $F$ can be built as a combination of $\vect{r}_1$ and $\vect{r}_2$ , i.e.

$\vect{x} = \lambda \vect{r}_1 + \mu \vect{r}_2.$

This combination is non-negative because there exists at least a constraint $c$ saturated by $\vect{r}_1$ and not $\vect{r}_2$ (or vice versa) (because they are distinct) for which

$\langle \vect{c}, \vect{x} \rangle \geq 0$

and

$\langle \vect{c}, \vect{x} \rangle = \lambda \langle \vect{c}, \vect{r}_1 \rangle (or = \mu \langle \vect{c}, \vect{r}_2 \rangle).$

So, there is no other extreme ray in $F$ and a) holds. Otherwise, if b) does not hold, the rank of the system generated by the constraints saturated by both $\vect{r}_1$ and $\vect{r}_2$ is equal to $d - k$ , with k >= 3, the set $F$ is k -dimensional and at least k extreme rays are necessary to generate $F$ . So, $\vect{r}_1$ and $\vect{r}_2$ are not adjacent and a) does not hold.

Proposition 2: When we build the new system of generators starting from a system $A$ of constraints of $P$ , if $\vect{c}$ is the constraint to add to $A$ and all lines of $P$ saturate $\vect{c}$ , the new set of rays is the union of those rays that saturate, of those that satisfy and of a set $\overline Q$ of rays such that each of them

lies on the hyper-plane represented by the k-th constraint,
is a positive combination of two adjacent rays $\vect{r}_1$ and $\vect{r}_2$ such that the first one satisfies the constraint and the other does not satisfy it. If the adjacency property is not taken in account, the new set of rays is not irredundant, in general.

In fact, if $\vect{r}_1$ and $\vect{r}_2$ are not adjacent, the rank of the system composed by the constraints saturated by both $\vect{r}_1$ and $\vect{r}_2$ is different from $d - 2$ (see the previous proposition) or neither $\vect{r}_1$ nor $\vect{r}_2$ are extreme rays. Since the new ray $\vect{r}$ is a combination of $\vect{r}_1$ and $\vect{r}_2$ , it saturates the same constraints saturated by both $\vect{r}_1$ and $\vect{r}_2$ . If the rank is less than $d - 2$ , the rank of the system composed by $\vect{c}$ (that is saturated by $\vect{r}$ ) and by the constraints of $A$ saturated by $\vect{r}$ is less than $d - 1$ . It means that $r$ is redundant (see Section Further Notation and Terminology). If neither $\vect{r}_1$ nor $\vect{r}_2$ are extreme rays, they belong to a 2-dimensional face containing exactly two extreme rays of $P$ . These two adjacent rays build a ray equal to $\vect{r}$ and so $\vect{r}$ is redundant.

Definition at line 347 of file conversion.cc.

Referenced by minimize().

                                                                    {
  dimension_type source_num_rows = source.num_rows();
  dimension_type dest_num_rows = dest.num_rows();
  const dimension_type source_num_columns = source.num_columns();
  const dimension_type dest_num_columns = dest.num_columns();

  using std::swap;

  // By construction, the number of columns of `sat' is the same as
  // the number of rows of `source'; also, the number of rows of `sat'
  // is the same as the number of rows of `dest'.
  PPL_ASSERT(source_num_rows == sat.num_columns());
  PPL_ASSERT(dest_num_rows == sat.num_rows());

  // If `start > 0', then we are converting the pending constraints.
  PPL_ASSERT(start == 0 || start == source.first_pending_row());

  // During the iteration on the constraints in `source' we may identify
  // constraints that are redundant: these have to be removed by swapping
  // the rows of `source', taking care not to compromise the sortedness
  // of the constraints that still have to be considered.
  // To this end, the following counter keeps the number of redundant
  // constraints seen so far, to be used as a displacement when swapping rows.
  dimension_type source_num_redundant = 0;

  PPL_DIRTY_TEMP_COEFFICIENT(normalized_sp_i);
  PPL_DIRTY_TEMP_COEFFICIENT(normalized_sp_o);

  // Converting the sub-system of `source' having rows with indexes
  // from `start' to the last one (i.e., `source_num_rows' - 1).
  for (dimension_type k = start; k < source_num_rows; ) {

    // All the `source_num_redundant' redundant constraints identified so far
    // have consecutive indices starting from `k'.
    if (source_num_redundant > 0)
      // Let the next constraint have index `k'.
      // There is no need to swap the columns of `sat' (all zeroes).
      swap(source[k], source[k+source_num_redundant]);

    Linear_Row& source_k = source[k];

    // Constraints and generators must have the same dimension,
    // otherwise the scalar product below will bomb.
    PPL_ASSERT(source_num_columns == dest_num_columns);

    // `scalar_prod[i]' will contain the scalar product of the
    // constraint `source_k' and the generator `dest[i]'.  This
    // product is 0 if and only if the generator saturates the
    // constraint.
    PPL_DIRTY_TEMP(std::vector<Coefficient>, scalar_prod);
    if (dest_num_rows > scalar_prod.size())
      scalar_prod.resize(dest_num_rows, Coefficient_zero());
    // `index_non_zero' will indicate the first generator in `dest'
    // that does not saturate the constraint `source_k'.
    dimension_type index_non_zero = 0;
    for ( ; index_non_zero < dest_num_rows; ++index_non_zero) {
      WEIGHT_BEGIN();
      Scalar_Products::assign(scalar_prod[index_non_zero],
                              source_k,
                              dest[index_non_zero]);
      WEIGHT_ADD_MUL(17, source_num_columns);
      if (scalar_prod[index_non_zero] != 0)
        // The generator does not saturate the constraint.
        break;
      // Check if the client has requested abandoning all expensive
      // computations.  If so, the exception specified by the client
      // is thrown now.
      maybe_abandon();
    }
    for (dimension_type i = index_non_zero + 1; i < dest_num_rows; ++i) {
      WEIGHT_BEGIN();
      Scalar_Products::assign(scalar_prod[i], source_k, dest[i]);
      WEIGHT_ADD_MUL(25, source_num_columns);
      // Check if the client has requested abandoning all expensive
      // computations.  If so, the exception specified by the client
      // is thrown now.
      maybe_abandon();
    }

    // We first treat the case when `index_non_zero' is less than
    // `num_lines_or_equalities', i.e., when the generator that
    // does not saturate the constraint `source_k' is a line.
    // The other case (described later) is when all the lines
    // in `dest' (i.e., all the rows having indexes less than
    // `num_lines_or_equalities') do saturate the constraint.

    if (index_non_zero < num_lines_or_equalities) {
      // Since the generator `dest[index_non_zero]' does not saturate
      // the constraint `source_k', it can no longer be a line
      // (see saturation rule in Section \ref prelims).
      // Therefore, we first transform it to a ray.
      dest[index_non_zero].set_is_ray_or_point_or_inequality();
      // Of the two possible choices, we select the ray satisfying
      // the constraint (namely, the ray whose scalar product
      // with the constraint gives a positive result).
      if (scalar_prod[index_non_zero] < 0) {
        // The ray `dest[index_non_zero]' lies on the wrong half-space:
        // we change it to have the opposite direction.
        neg_assign(scalar_prod[index_non_zero]);
        for (dimension_type j = dest_num_columns; j-- > 0; )
          neg_assign(dest[index_non_zero][j]);
      }
      // Having changed a line to a ray, we set `dest' to be a
      // non-sorted system, we decrement the number of lines of `dest' and,
      // if necessary, we move the new ray below all the remaining lines.
      dest.set_sorted(false);
      --num_lines_or_equalities;
      if (index_non_zero != num_lines_or_equalities) {
        swap(dest[index_non_zero],
             dest[num_lines_or_equalities]);
        swap(scalar_prod[index_non_zero],
             scalar_prod[num_lines_or_equalities]);
      }
      Linear_Row& dest_nle = dest[num_lines_or_equalities];

      // Computing the new lineality space.
      // Since each line must lie on the hyper-plane corresponding to
      // the constraint `source_k', the scalar product between
      // the line and the constraint must be 0.
      // This property already holds for the lines having indexes
      // between 0 and `index_non_zero' - 1.
      // We have to consider the remaining lines, having indexes
      // between `index_non_zero' and `num_lines_or_equalities' - 1.
      // Each line that does not saturate the constraint has to be
      // linearly combined with generator `dest_nle' so that the
      // resulting new line saturates the constraint.
      // Note that, by Observation 1 above, the resulting new line
      // will still saturate all the constraints that were saturated by
      // the old line.

      Coefficient& scalar_prod_nle = scalar_prod[num_lines_or_equalities];
      for (dimension_type
             i = index_non_zero; i < num_lines_or_equalities; ++i) {
        if (scalar_prod[i] != 0) {
          // The following fragment optimizes the computation of
          //
          // <CODE>
          //   Coefficient scale = scalar_prod[i];
          //   scale.gcd_assign(scalar_prod_nle);
          //   Coefficient normalized_sp_i = scalar_prod[i] / scale;
          //   Coefficient normalized_sp_n = scalar_prod_nle / scale;
          //   for (dimension_type c = dest_num_columns; c-- > 0; ) {
          //     dest[i][c] *= normalized_sp_n;
          //     dest[i][c] -= normalized_sp_i * dest_nle[c];
          //   }
          // </CODE>
          normalize2(scalar_prod[i],
                     scalar_prod_nle,
                     normalized_sp_i,
                     normalized_sp_o);
          Linear_Row& dest_i = dest[i];
          for (dimension_type c = dest_num_columns; c-- > 0; ) {
            Coefficient& dest_i_c = dest_i[c];
            dest_i_c *= normalized_sp_o;
            sub_mul_assign(dest_i_c, normalized_sp_i, dest_nle[c]);
          }
          dest_i.strong_normalize();
          scalar_prod[i] = 0;
          // `dest' has already been set as non-sorted.
        }
      }

      // Computing the new pointed cone.
      // Similarly to what we have done during the computation of
      // the lineality space, we consider all the remaining rays
      // (having indexes strictly greater than `num_lines_or_equalities')
      // that do not saturate the constraint `source_k'. These rays
      // are positively combined with the ray `dest_nle' so that the
      // resulting new rays saturate the constraint.
      for (dimension_type
             i = num_lines_or_equalities + 1; i < dest_num_rows; ++i) {
        if (scalar_prod[i] != 0) {
          // The following fragment optimizes the computation of
          //
          // <CODE>
          //   Coefficient scale = scalar_prod[i];
          //   scale.gcd_assign(scalar_prod_nle);
          //   Coefficient normalized_sp_i = scalar_prod[i] / scale;
          //   Coefficient normalized_sp_n = scalar_prod_nle / scale;
          //   for (dimension_type c = dest_num_columns; c-- > 0; ) {
          //     dest[i][c] *= normalized_sp_n;
          //     dest[i][c] -= normalized_sp_i * dest_nle[c];
          //   }
          // </CODE>
          normalize2(scalar_prod[i],
                     scalar_prod_nle,
                     normalized_sp_i,
                     normalized_sp_o);
          Linear_Row& dest_i = dest[i];
          WEIGHT_BEGIN();
          for (dimension_type c = dest_num_columns; c-- > 0; ) {
            Coefficient& dest_i_c = dest_i[c];
            dest_i_c *= normalized_sp_o;
            sub_mul_assign(dest_i_c, normalized_sp_i, dest_nle[c]);
          }
          dest_i.strong_normalize();
          scalar_prod[i] = 0;
          // `dest' has already been set as non-sorted.
          WEIGHT_ADD_MUL(41, dest_num_columns);
        }
        // Check if the client has requested abandoning all expensive
        // computations.  If so, the exception specified by the client
        // is thrown now.
        maybe_abandon();
      }
      // Since the `scalar_prod_nle' is positive (by construction), it
      // does not saturate the constraint `source_k'.  Therefore, if
      // the constraint is an inequality, we set to 1 the
      // corresponding element of `sat' ...
      Bit_Row& sat_nle = sat[num_lines_or_equalities];
      if (source_k.is_ray_or_point_or_inequality())
        sat_nle.set(k);
      // ... otherwise, the constraint is an equality which is
      // violated by the generator `dest_nle': the generator has to be
      // removed from `dest'.
      else {
        --dest_num_rows;
        swap(dest_nle, dest[dest_num_rows]);
        swap(scalar_prod_nle, scalar_prod[dest_num_rows]);
        swap(sat_nle, sat[dest_num_rows]);
        // `dest' has already been set as non-sorted.
      }
      // We continue with the next constraint.
      ++k;
    }
    // Here we have `index_non_zero' >= `num_lines_or_equalities',
    // so that all the lines in `dest' saturate the constraint `source_k'.
    else {
      // First, we reorder the generators in `dest' as follows:
      // -# all the lines should have indexes between 0 and
      //    `num_lines_or_equalities' - 1 (this already holds);
      // -# all the rays that saturate the constraint should have
      //    indexes between `num_lines_or_equalities' and
      //    `lines_or_equal_bound' - 1; these rays form the set Q=.
      // -# all the rays that have a positive scalar product with the
      //    constraint should have indexes between `lines_or_equal_bound'
      //    and `sup_bound' - 1; these rays form the set Q+.
      // -# all the rays that have a negative scalar product with the
      //    constraint should have indexes between `sup_bound' and
      //    `dest_num_rows' - 1; these rays form the set Q-.
      dimension_type lines_or_equal_bound = num_lines_or_equalities;
      dimension_type inf_bound = dest_num_rows;
      // While we find saturating generators, we simply increment
      // `lines_or_equal_bound'.
      while (inf_bound > lines_or_equal_bound
             && scalar_prod[lines_or_equal_bound] == 0)
        ++lines_or_equal_bound;
      dimension_type sup_bound = lines_or_equal_bound;
      while (inf_bound > sup_bound) {
        const int sp_sign = sgn(scalar_prod[sup_bound]);
        if (sp_sign == 0) {
          // This generator has to be moved in Q=.
          swap(dest[sup_bound], dest[lines_or_equal_bound]);
          swap(scalar_prod[sup_bound], scalar_prod[lines_or_equal_bound]);
          swap(sat[sup_bound], sat[lines_or_equal_bound]);
          ++lines_or_equal_bound;
          ++sup_bound;
          dest.set_sorted(false);
        }
        else if (sp_sign < 0) {
          // This generator has to be moved in Q-.
          --inf_bound;
          swap(dest[sup_bound], dest[inf_bound]);
          swap(scalar_prod[sup_bound], scalar_prod[inf_bound]);
          swap(sat[sup_bound], sat[inf_bound]);
          dest.set_sorted(false);
        }
        else
          // sp_sign > 0: this generator has to be moved in Q+.
          ++sup_bound;
      }

      if (sup_bound == dest_num_rows) {
        // Here the set Q- is empty.
        // If the constraint is an inequality, then all the generators
        // in Q= and Q+ satisfy the constraint. The constraint is redundant
        // and it can be safely removed from the constraint system.
        // This is why the `source' parameter is not declared `const'.
        if (source_k.is_ray_or_point_or_inequality()) {
          ++source_num_redundant;
          --source_num_rows;
          // NOTE: we continue with the next cycle of the loop
          // without incrementing the index `k', because:
          // -# either `k == source_num_rows', and we will exit the loop;
          // -# or, having increased `source_num_redundant', we will swap
          //    in position `k' a constraint that still has to be examined.
        }
        else {
          // The constraint is an equality, so that all the generators
          // in Q+ violate it. Since the set Q- is empty, we can simply
          // remove from `dest' all the generators of Q+.
          dest_num_rows = lines_or_equal_bound;
          // We continue with the next constraint.
          ++k;
        }
      }
      else {
        // The set Q- is not empty, i.e., at least one generator
        // violates the constraint `source_k'.
        // We have to further distinguish two cases:
        if (sup_bound == num_lines_or_equalities)
          // The set Q+ is empty, so that all generators that satisfy
          // the constraint also saturate it.
          // We can simply remove from `dest' all the generators in Q-.
          dest_num_rows = sup_bound;
        else {
          // The sets Q+ and Q- are both non-empty.
          // The generators of the new pointed cone are all those satisfying
          // the constraint `source_k' plus a set of new rays enjoying
          // the following properties:
          // -# they lie on the hyper-plane represented by the constraint
          // -# they are obtained as a positive combination of two
          //    adjacent rays, the first taken from Q+ and the second
          //    taken from Q-.

          // The adjacency property is necessary to have an irredundant
          // set of new rays (see proposition 2).
          const dimension_type bound = dest_num_rows;

          // In the following loop,
          // `i' runs through the generators in the set Q+ and
          // `j' runs through the generators in the set Q-.
          for (dimension_type i = lines_or_equal_bound; i < sup_bound; ++i) {
            for(dimension_type j = sup_bound; j < bound; ++j) {
              // Checking if generators `dest[i]' and `dest[j]' are adjacent.
              // If there exist another generator that saturates
              // all the constraints saturated by both `dest[i]' and
              // `dest[j]', then they are NOT adjacent.
              PPL_ASSERT(sat[i].last() == C_Integer<unsigned long>::max
                         || sat[i].last() < k);
              PPL_ASSERT(sat[j].last() == C_Integer<unsigned long>::max
                         || sat[j].last() < k);

              // Being the union of `sat[i]' and `sat[j]',
              // `new_satrow' corresponds to a ray that saturates all the
              // constraints saturated by both `dest[i]' and `dest[j]'.
              Bit_Row new_satrow(sat[i], sat[j]);

              // Compute the number of common saturators.
              // NOTE: this number has to be less than `k' because
              // we are treating the `k'-th constraint.
              const dimension_type
                num_common_satur = k - new_satrow.count_ones();

              // Even before actually creating the new ray as a
              // positive combination of `dest[i]' and `dest[j]',
              // we exploit saturation information to check if
              // it can be an extremal ray. To this end, we refer
              // to the definition of a minimal proper face
              // (see comments in Polyhedron.defs.hh):
              // an extremal ray saturates at least `n' - `t' - 1
              // constraints, where `n' is the dimension of the space
              // and `t' is the dimension of the lineality space.
              // Since `n == source_num_columns - 1' and
              // `t == num_lines_or_equalities', we obtain that
              // an extremal ray saturates at least
              // `source_num_columns - num_lines_or_equalities - 2'
              // constraints.
              if (num_common_satur
                  >= source_num_columns - num_lines_or_equalities - 2) {
                // The minimal proper face rule is satisfied.
                // Now we actually check for redundancy by computing
                // adjacency information.
                bool redundant = false;
                WEIGHT_BEGIN();
                for (dimension_type
                       l = num_lines_or_equalities; l < bound; ++l)
                  if (l != i && l != j
                      && subset_or_equal(sat[l], new_satrow)) {
                    // Found another generator saturating all the
                    // constraints saturated by both `dest[i]' and `dest[j]'.
                    redundant = true;
                    break;
                  }
                PPL_ASSERT(bound >= num_lines_or_equalities);
                WEIGHT_ADD_MUL(15, bound - num_lines_or_equalities);
                if (!redundant) {
                  // Adding the new ray to `dest' and the corresponding
                  // saturation row to `sat'.
                  if (dest_num_rows == dest.num_rows()) {
                    // Make room for one more row.
                    dest.add_pending_row(Linear_Row::Flags(dest.topology(),
                                                           Linear_Row::RAY_OR_POINT_OR_INEQUALITY));
                    sat.add_recycled_row(new_satrow);
                  }
                  else
                    sat[dest_num_rows].m_swap(new_satrow);

                  Linear_Row& new_row = dest[dest_num_rows];
                  // The following fragment optimizes the computation of
                  //
                  // <CODE>
                  //   Coefficient scale = scalar_prod[i];
                  //   scale.gcd_assign(scalar_prod[j]);
                  //   Coefficient normalized_sp_i = scalar_prod[i] / scale;
                  //   Coefficient normalized_sp_j = scalar_prod[j] / scale;
                  //   for (dimension_type c = dest_num_columns; c-- > 0; ) {
                  //     new_row[c] = normalized_sp_i * dest[j][c];
                  //     new_row[c] -= normalized_sp_j * dest[i][c];
                  //   }
                  // </CODE>
                  normalize2(scalar_prod[i],
                             scalar_prod[j],
                             normalized_sp_i,
                             normalized_sp_o);
                  WEIGHT_BEGIN();
                  for (dimension_type c = dest_num_columns; c-- > 0; ) {
                    Coefficient& new_row_c = new_row[c];
                    new_row_c = normalized_sp_i * dest[j][c];
                    sub_mul_assign(new_row_c, normalized_sp_o, dest[i][c]);
                  }
                  WEIGHT_ADD_MUL(86, dest_num_columns);
                  new_row.strong_normalize();
                  // Since we added a new generator to `dest',
                  // we also add a new element to `scalar_prod';
                  // by construction, the new ray lies on the hyper-plane
                  // represented by the constraint `source_k'.
                  // Thus, the added scalar product is 0.
                  PPL_ASSERT(scalar_prod.size() >= dest_num_rows);
                  if (scalar_prod.size() <= dest_num_rows)
                    scalar_prod.push_back(Coefficient_zero());
                  else
                    scalar_prod[dest_num_rows] = Coefficient_zero();
                  // Increment the number of generators.
                  ++dest_num_rows;
                } // if (!redundant)
              }
            }
            // Check if the client has requested abandoning all expensive
            // computations.  If so, the exception specified by the client
            // is thrown now.
            maybe_abandon();
          }
          // Now we substitute the rays in Q- (i.e., the rays violating
          // the constraint) with the newly added rays.
          dimension_type j;
          if (source_k.is_ray_or_point_or_inequality()) {
            // The constraint is an inequality:
            // the violating generators are those in Q-.
            j = sup_bound;
            // For all the generators in Q+, set to 1 the corresponding
            // entry for the constraint `source_k' in the saturation matrix.
            for (dimension_type l = lines_or_equal_bound; l < sup_bound; ++l)
              sat[l].set(k);
          }
          else
            // The constraint is an equality:
            // the violating generators are those in the union of Q+ and Q-.
            j = lines_or_equal_bound;

          // Swapping the newly added rays
          // (index `i' running through `dest_num_rows - 1' down-to `bound')
          // with the generators violating the constraint
          // (index `j' running through `j' up-to `bound - 1').
          dimension_type i = dest_num_rows;
          while (j < bound && i > bound) {
            --i;
            swap(dest[i], dest[j]);
            swap(scalar_prod[i], scalar_prod[j]);
            swap(sat[i], sat[j]);
            ++j;
            dest.set_sorted(false);
          }
          // Setting the number of generators in `dest':
          // - if the number of generators violating the constraint
          //   is less than or equal to the number of the newly added
          //   generators, we assign `i' to `dest_num_rows' because
          //   all generators above this index are significant;
          // - otherwise, we assign `j' to `dest_num_rows' because
          //   all generators below index `j-1' violates the constraint.
          dest_num_rows = (j == bound) ? i : j;
        }
        // We continue with the next constraint.
        ++k;
      }
    }
  }

  // We may have identified some redundant constraints in `source',
  // which have been swapped at the end of the system.
  if (source_num_redundant > 0) {
    PPL_ASSERT(source_num_redundant == source.num_rows() - source_num_rows);
    source.remove_trailing_rows(source_num_redundant);
    sat.remove_trailing_columns(source_num_redundant);
  }
  // If `start == 0', then `source' was sorted and remained so.
  // If otherwise `start > 0', then the two sub-system made by the
  // non-pending rows and the pending rows, respectively, were both sorted.
  // Thus, the overall system is sorted if and only if either
  // `start == source_num_rows' (i.e., the second sub-system is empty)
  // or the row ordering holds for the two rows at the boundary between
  // the two sub-systems.
  if (start > 0 && start < source_num_rows)
    source.set_sorted(compare(source[start - 1], source[start]) <= 0);
  // There are no longer pending constraints in `source'.
  source.unset_pending_rows();

  // We may have identified some redundant rays in `dest',
  // which have been swapped at the end of the system.
  if (dest_num_rows < dest.num_rows()) {
    const dimension_type num_removed_rows = dest.num_rows() - dest_num_rows;
    dest.remove_trailing_rows(num_removed_rows);
    // Be careful: we might have erased some of the non-pending rows.
    if (dest.first_pending_row() > dest_num_rows)
      dest.unset_pending_rows();
    sat.remove_trailing_rows(num_removed_rows);
  }
  if (dest.is_sorted())
    // If the non-pending generators in `dest' are still declared to be
    // sorted, then we have to also check for the sortedness of the
    // pending generators.
    for (dimension_type i = dest.first_pending_row(); i < dest_num_rows; ++i)
      if (compare(dest[i - 1], dest[i]) > 0) {
        dest.set_sorted(false);
        break;
      }
  // There are no pending generators in `dest'.
  dest.unset_pending_rows();

  return num_lines_or_equalities;
}

template<typename FP_Format , typename Interval_Info >

void Parma_Polyhedra_Library::Polyhedron::convert_to_integer_expression	(	const Linear_Form< Interval< FP_Format, Interval_Info > > &	lf,
		const dimension_type	lf_dimension,
		Linear_Expression &	result
	)

staticprotected

Helper function that makes result become a Linear_Expression obtained by normalizing the denominators in lf.

Parameters:

lf	The linear form on intervals with floating point boundaries to normalize. It should be the result of an application of static method `overapproximate_linear_form`.
lf_dimension	Must be the space dimension of `lf`.
result	Used to store the result.

This function ignores the upper bound of intervals in lf, so that in fact result can be seen as lf multiplied by a proper normalization constant.

Definition at line 471 of file Polyhedron.templates.hh.

References Parma_Polyhedra_Library::exact_div_assign(), Parma_Polyhedra_Library::lcm_assign(), Parma_Polyhedra_Library::numer_denom(), and PPL_DIRTY_TEMP_COEFFICIENT.

Referenced by refine_with_linear_form_inequality().

                                           {
  result = Linear_Expression();

  typedef Interval<FP_Format, Interval_Info> FP_Interval_Type;
  std::vector<Coefficient> numerators(lf_dimension+1);
  std::vector<Coefficient> denominators(lf_dimension+1);

  // Convert each floating point number to a pair <numerator, denominator>
  // and compute the lcm of all denominators.
  PPL_DIRTY_TEMP_COEFFICIENT(lcm);
  lcm = 1;
  const FP_Interval_Type& b = lf.inhomogeneous_term();
  // FIXME: are these checks numerator[i] != 0 really necessary?
  numer_denom(b.lower(), numerators[lf_dimension],
                         denominators[lf_dimension]);
  if (numerators[lf_dimension] != 0)
      lcm_assign(lcm, lcm, denominators[lf_dimension]);

  for (dimension_type i = 0; i < lf_dimension; ++i) {
    const FP_Interval_Type& curr_int = lf.coefficient(Variable(i));
    numer_denom(curr_int.lower(), numerators[i], denominators[i]);
    if (numerators[i] != 0)
      lcm_assign(lcm, lcm, denominators[i]);
  }

  for (dimension_type i = 0; i < lf_dimension; ++i) {
    if (numerators[i] != 0) {
      exact_div_assign(denominators[i], lcm, denominators[i]);
      numerators[i] *= denominators[i];
      result += numerators[i] * Variable(i);
    }
  }

  if (numerators[lf_dimension] != 0) {
    exact_div_assign(denominators[lf_dimension],
                     lcm, denominators[lf_dimension]);
    numerators[lf_dimension] *= denominators[lf_dimension];
    result += numerators[lf_dimension];
  }
}

template<typename FP_Format , typename Interval_Info >

void Parma_Polyhedra_Library::Polyhedron::convert_to_integer_expressions	(	const Linear_Form< Interval< FP_Format, Interval_Info > > &	lf,
		const dimension_type	lf_dimension,
		Linear_Expression &	res,
		Coefficient &	res_low_coeff,
		Coefficient &	res_hi_coeff,
		Coefficient &	denominator
	)

staticprotected

Normalization helper function.

Parameters:

lf	The linear form on intervals with floating point boundaries to normalize. It should be the result of an application of static method `overapproximate_linear_form`.
lf_dimension	Must be the space dimension of `lf`.
res	Stores the normalized linear form, except its inhomogeneous term.
res_low_coeff	Stores the lower boundary of the inhomogeneous term of the result.
res_hi_coeff	Stores the higher boundary of the inhomogeneous term of the result.
denominator	Becomes the common denominator of `res_low_coeff`, `res_hi_coeff` and all coefficients in `res`.

Results are obtained by normalizing denominators in lf, ignoring the upper bounds of variable coefficients in lf.

Definition at line 517 of file Polyhedron.templates.hh.

References Parma_Polyhedra_Library::exact_div_assign(), Parma_Polyhedra_Library::lcm_assign(), and Parma_Polyhedra_Library::numer_denom().

Referenced by affine_form_image().

                                          {
  res = Linear_Expression();

  typedef Interval<FP_Format, Interval_Info> FP_Interval_Type;
  std::vector<Coefficient> numerators(lf_dimension+2);
  std::vector<Coefficient> denominators(lf_dimension+2);

  // Convert each floating point number to a pair <numerator, denominator>
  // and compute the lcm of all denominators.
  Coefficient& lcm = denominator;
  lcm = 1;
  const FP_Interval_Type& b = lf.inhomogeneous_term();
  numer_denom(b.lower(), numerators[lf_dimension], denominators[lf_dimension]);
  // FIXME: are these checks numerator[i] != 0 really necessary?
  if (numerators[lf_dimension] != 0)
      lcm_assign(lcm, lcm, denominators[lf_dimension]);

  numer_denom(b.upper(), numerators[lf_dimension+1],
                         denominators[lf_dimension+1]);
  if (numerators[lf_dimension+1] != 0)
      lcm_assign(lcm, lcm, denominators[lf_dimension+1]);

  for (dimension_type i = 0; i < lf_dimension; ++i) {
    const FP_Interval_Type& curr_int = lf.coefficient(Variable(i));
    numer_denom(curr_int.lower(), numerators[i], denominators[i]);
    if (numerators[i] != 0)
      lcm_assign(lcm, lcm, denominators[i]);
  }

  for (dimension_type i = 0; i < lf_dimension; ++i) {
    if (numerators[i] != 0) {
      exact_div_assign(denominators[i], lcm, denominators[i]);
      numerators[i] *= denominators[i];
      res += numerators[i] * Variable(i);
    }
  }

  if (numerators[lf_dimension] != 0) {
    exact_div_assign(denominators[lf_dimension],
                     lcm, denominators[lf_dimension]);
    numerators[lf_dimension] *= denominators[lf_dimension];
    res_low_coeff = numerators[lf_dimension];
  }
  else
    res_low_coeff = Coefficient(0);

  if (numerators[lf_dimension+1] != 0) {
    exact_div_assign(denominators[lf_dimension+1],
                     lcm, denominators[lf_dimension+1]);
    numerators[lf_dimension+1] *= denominators[lf_dimension+1];
    res_hi_coeff = numerators[lf_dimension+1];
  }
  else
    res_hi_coeff = Coefficient(0);
}

void Parma_Polyhedra_Library::Polyhedron::difference_assign ( const Polyhedron & y )

inline

Same as poly_difference_assign(y).

Definition at line 86 of file Polyhedron.inlines.hh.

References poly_difference_assign().

                                                 {
  poly_difference_assign(y);
}

void Parma_Polyhedra_Library::Polyhedron::drop_some_non_integer_points ( Complexity_Class complexity = ANY_COMPLEXITY )

inline

Possibly tightens *this by dropping some points with non-integer coordinates.

Parameters:

complexity The maximal complexity of any algorithms used.

Note:: Currently there is no optimality guarantee, not even if complexity is ANY_COMPLEXITY.

Definition at line 438 of file Polyhedron.inlines.hh.

Referenced by drop_some_non_integer_points().

                                                                    {
  const Variables_Set* p_vs = 0;
  drop_some_non_integer_points(p_vs, complexity);
}

void Parma_Polyhedra_Library::Polyhedron::drop_some_non_integer_points	(	const Variables_Set &	vars,
		Complexity_Class	complexity = `ANY_COMPLEXITY`
	)

inline

Possibly tightens *this by dropping some points with non-integer coordinates for the space dimensions corresponding to vars.

Parameters:

vars	Points with non-integer coordinates for these variables/space-dimensions can be discarded.
complexity	The maximal complexity of any algorithms used.

Note:: Currently there is no optimality guarantee, not even if complexity is ANY_COMPLEXITY.

Definition at line 444 of file Polyhedron.inlines.hh.

References drop_some_non_integer_points().

                                                                      {
  drop_some_non_integer_points(&vars, complexity);
}

void Parma_Polyhedra_Library::Polyhedron::drop_some_non_integer_points	(	const Variables_Set *	vars_p,
		Complexity_Class	complexity
	)

protected

Possibly tightens *this by dropping some points with non-integer coordinates for the space dimensions corresponding to *vars_p.

Parameters:

vars_p	When nonzero, points with non-integer coordinates for the variables/space-dimensions contained in `*vars_p` can be discarded.
complexity	The maximal complexity of any algorithms used.

Note:: Currently there is no optimality guarantee, not even if complexity is ANY_COMPLEXITY.

Definition at line 2100 of file Polyhedron_nonpublic.cc.

References Parma_Polyhedra_Library::ANY_COMPLEXITY, c, Parma_Polyhedra_Library::Constraint::epsilon_leq_one(), Parma_Polyhedra_Library::exact_div_assign(), Parma_Polyhedra_Library::gcd_assign(), Parma_Polyhedra_Library::Constraint::is_equality(), Parma_Polyhedra_Library::Constraint::is_tautological(), Parma_Polyhedra_Library::neg_assign(), Parma_Polyhedra_Library::Dense_Row::normalize(), PPL_ASSERT, PPL_ASSERT_HEAVY, PPL_DIRTY_TEMP_COEFFICIENT, and Parma_Polyhedra_Library::Boundary_NS::sgn().

                                                                           {
  // There is nothing to do for an empty set of variables.
  if (vars_p != 0 && vars_p->empty())
    return;

  // Any empty polyhedron does not contain integer points.
  if (marked_empty())
    return;

  // A zero-dimensional, universe polyhedron has, by convention, an
  // integer point.
  if (space_dim == 0) {
    set_empty();
    return;
  }

  // The constraints (possibly with pending rows) are required.
  if (has_pending_generators()) {
    // Processing of pending generators is exponential in the worst case.
    if (complexity != ANY_COMPLEXITY)
      return;
    else
      process_pending_generators();
  }
  if (!constraints_are_up_to_date()) {
    // Constraints update is exponential in the worst case.
    if (complexity != ANY_COMPLEXITY)
      return;
    else
      update_constraints();
  }
  // For NNC polyhedra we need to process any pending constraints.
  if (!is_necessarily_closed() && has_pending_constraints()) {
    if (complexity != ANY_COMPLEXITY)
      return;
    else if (!process_pending_constraints())
      // We just discovered the polyhedron is empty.
      return;
  }

  PPL_ASSERT(!has_pending_generators() && constraints_are_up_to_date());
  PPL_ASSERT(is_necessarily_closed() || !has_pending_constraints());

  bool changed = false;
  const dimension_type eps_index = space_dim + 1;
  PPL_DIRTY_TEMP_COEFFICIENT(gcd);

  for (dimension_type j = con_sys.num_rows(); j-- > 0; ) {
    Constraint& c = con_sys[j];
    if (c.is_tautological())
      goto next_constraint;

    if (vars_p != 0) {
      for (dimension_type i = space_dim; i-- > 0; )
        if (c[i+1] != 0 && vars_p->find(i) == vars_p->end())
          goto next_constraint;
    }

    if (!is_necessarily_closed()) {
      // Transform all strict inequalities into non-strict ones,
      // with the inhomogeneous term incremented by 1.
      if (c[eps_index] < 0) {
        c[eps_index] = 0;
        --c[0];
        // Enforce normalization.
        // FIXME: is this really necessary?
        c.normalize();
        changed = true;
      }
    }

    {
      // Compute the GCD of all the homogeneous terms.
      dimension_type i = space_dim+1;
      while (i > 1) {
        const Coefficient& c_i = c[--i];
        if (const int c_i_sign = sgn(c_i)) {
          gcd = c_i;
          if (c_i_sign < 0)
            neg_assign(gcd);
          goto compute_gcd;
        }
      }
      // We reach this point only if all the coefficients were zero.
      goto next_constraint;

    compute_gcd:
      if (gcd == 1)
        goto next_constraint;
      while (i > 1) {
        const Coefficient& c_i = c[--i];
        if (c_i != 0) {
          // See the comment in Dense_Row::normalize().
          gcd_assign(gcd, c_i, gcd);
          if (gcd == 1)
            goto next_constraint;
        }
      }
      PPL_ASSERT(gcd != 1);
      PPL_ASSERT(c[0] % gcd != 0);

      // If we have an equality, the polyhedron becomes empty.
      if (c.is_equality()) {
        set_empty();
        return;
      }

      // Divide the inhomogeneous coefficients by the GCD.
      for (dimension_type k = space_dim+1; --k > 0; ) {
        Coefficient& c_k = c[k];
        exact_div_assign(c_k, c_k, gcd);
      }
      Coefficient& c_0 = c[0];
      const int c_0_sign = sgn(c_0);
      c_0 /= gcd;
      if (c_0_sign < 0)
        --c_0;
      changed = true;
    }

  next_constraint:
    ;
  }

  if (changed) {
    if (!is_necessarily_closed()) {
      con_sys.insert(Constraint::epsilon_leq_one());
      // FIXME: make sure that the following line really can stay here
      // and should not be moved below the brace.
      con_sys.set_sorted(false);
    }

    // After changing the system of constraints, the generators
    // are no longer up-to-date and the constraints are no longer
    // minimized.
    clear_generators_up_to_date();
    clear_constraints_minimized();
  }
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::expand_space_dimension	(	Variable	var,
		dimension_type	m
	)

Creates m copies of the space dimension corresponding to var.

Parameters:

var	The variable corresponding to the space dimension to be replicated;
m	The number of replicas to be created.

Exceptions:

std::invalid_argument	Thrown if `var` does not correspond to a dimension of the vector space.
std::length_error	Thrown if adding `m` new space dimensions would cause the vector space to exceed dimension `max_space_dimension()`.

If *this has space dimension $n$ , with $n > 0$ , and var has space dimension $k \leq n$ , then the $k$ -th space dimension is expanded to m new space dimensions $n$ , $n+1$ , $\dots$ , $n+m-1$ .

Definition at line 553 of file Polyhedron_chdims.cc.

                                                                    {
  // TODO: this implementation is _really_ an executable specification.

  // `var' should be one of the dimensions of the vector space.
  if (var.space_dimension() > space_dim)
    throw_dimension_incompatible("expand_space_dimension(v, m)", "v", var);

  // The space dimension of the resulting polyhedron should not
  // overflow the maximum allowed space dimension.
  check_space_dimension_overflow(m, max_space_dimension() - space_dimension(),
                                 topology(),
                                 "expand_dimension(v, m)",
                                 "adding m new space dimensions exceeds "
                                 "the maximum allowed space dimension");

  // Nothing to do, if no dimensions must be added.
  if (m == 0)
    return;

  // Keep track of the dimension before adding the new ones.
  dimension_type old_dim = space_dim;

  // Add the required new dimensions.
  add_space_dimensions_and_embed(m);

  const dimension_type src_d = var.id();
  const Constraint_System& cs = constraints();
  Constraint_System new_constraints;
  for (Constraint_System::const_iterator i = cs.begin(),
         cs_end = cs.end(); i != cs_end; ++i) {
    const Constraint& c = *i;

    // If `c' does not constrain `var', skip it.
    if (c.coefficient(var) == 0)
      continue;

    // Each relevant constraint results in `m' new constraints.
    for (dimension_type dst_d = old_dim; dst_d < old_dim+m; ++dst_d) {
      Linear_Expression e;
      for (dimension_type j = old_dim; j-- > 0; )
        e +=
          c.coefficient(Variable(j))
          * ((j == src_d) ? Variable(dst_d) : Variable(j));
      e += c.inhomogeneous_term();
      new_constraints.insert(c.is_equality()
                             ? (e == 0)
                             : (c.is_nonstrict_inequality()
                                ? (e >= 0)
                                : (e > 0)));
    }
  }
  add_recycled_constraints(new_constraints);
  PPL_ASSERT_HEAVY(OK());
}

PPL::memory_size_type Parma_Polyhedra_Library::Polyhedron::external_memory_in_bytes ( ) const

Returns the size in bytes of the memory managed by *this.

Definition at line 3760 of file Polyhedron_public.cc.

Referenced by total_memory_in_bytes().

                                              {
  return
    con_sys.external_memory_in_bytes()
    + gen_sys.external_memory_in_bytes()
    + sat_c.external_memory_in_bytes()
    + sat_g.external_memory_in_bytes();
}

void Parma_Polyhedra_Library::Polyhedron::finalize ( )

static

Finalizes the class.

Definition at line 53 of file Polyhedron_public.cc.

Referenced by Parma_Polyhedra_Library::Init::~Init().

                        {
  delete [] simplify_num_saturators_p;
  simplify_num_saturators_p = 0;
  simplify_num_saturators_size = 0;
}

void Parma_Polyhedra_Library::Polyhedron::fold_space_dimensions	(	const Variables_Set &	vars,
		Variable	dest
	)

Folds the space dimensions in vars into dest.

Parameters:

vars	The set of Variable objects corresponding to the space dimensions to be folded;
dest	The variable corresponding to the space dimension that is the destination of the folding operation.

Exceptions:

std::invalid_argument Thrown if *this is dimension-incompatible with dest or with one of the Variable objects contained in vars. Also thrown if dest is contained in vars.

If *this has space dimension $n$ , with $n > 0$ , dest has space dimension $k \leq n$ , vars is a set of variables whose maximum space dimension is also less than or equal to $n$ , and dest is not a member of vars, then the space dimensions corresponding to variables in vars are folded into the $k$ -th space dimension.

Definition at line 609 of file Polyhedron_chdims.cc.

References affine_image(), Parma_Polyhedra_Library::Variable::id(), PPL_ASSERT_HEAVY, Parma_Polyhedra_Library::Variables_Set::space_dimension(), and Parma_Polyhedra_Library::Variable::space_dimension().

                                                      {
  // TODO: this implementation is _really_ an executable specification.

  // `dest' should be one of the dimensions of the polyhedron.
  if (dest.space_dimension() > space_dim)
    throw_dimension_incompatible("fold_space_dimensions(vs, v)", "v", dest);

  // The folding of no dimensions is a no-op.
  if (vars.empty())
    return;

  // All variables in `vars' should be dimensions of the polyhedron.
  if (vars.space_dimension() > space_dim)
    throw_dimension_incompatible("fold_space_dimensions(vs, v)",
                                 "vs.space_dimension()",
                                 vars.space_dimension());

  // Moreover, `dest.id()' should not occur in `vars'.
  if (vars.find(dest.id()) != vars.end())
    throw_invalid_argument("fold_space_dimensions(vs, v)",
                           "v should not occur in vs");

  // All of the affine images we are going to compute are not invertible,
  // hence we will need to compute the generators of the polyhedron.
  // Since we keep taking copies, make sure that a single conversion
  // from constraints to generators is computed.
  (void) generators();
  // Having generators, we now know if the polyhedron is empty:
  // in that case, folding is equivalent to just removing space dimensions.
  if (!marked_empty()) {
    for (Variables_Set::const_iterator i = vars.begin(),
           vs_end = vars.end(); i != vs_end; ++i) {
      Polyhedron copy = *this;
      copy.affine_image(dest, Linear_Expression(Variable(*i)));
      poly_hull_assign(copy);
    }
  }
  remove_space_dimensions(vars);
  PPL_ASSERT_HEAVY(OK());
}

bool Parma_Polyhedra_Library::Polyhedron::frequency	(	const Linear_Expression &	expr,
		Coefficient &	freq_n,
		Coefficient &	freq_d,
		Coefficient &	val_n,
		Coefficient &	val_d
	)		const

Returns true if and only if there exist a unique value val such that *this saturates the equality expr = val.

Parameters:

expr	The linear expression for which the frequency is needed;
freq_n	If `true` is returned, the value is set to ; Present for interface compatibility with class Grid, where the frequency can have a non-zero value;
freq_d	If `true` is returned, the value is set to ;
val_n	The numerator of `val`;
val_d	The denominator of `val`;

Exceptions:

std::invalid_argument Thrown if expr and *this are dimension-incompatible.

If false is returned, then freq_n, freq_d, val_n and val_d are left untouched.

Definition at line 3475 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::assign_r(), Parma_Polyhedra_Library::Scalar_Products::homogeneous_assign(), Parma_Polyhedra_Library::Linear_Expression::inhomogeneous_term(), Parma_Polyhedra_Library::Generator::is_closure_point(), Parma_Polyhedra_Library::Generator::is_line_or_ray(), Parma_Polyhedra_Library::Generator::is_point(), PPL_ASSERT, PPL_DIRTY_TEMP, PPL_DIRTY_TEMP_COEFFICIENT, Parma_Polyhedra_Library::ROUND_NOT_NEEDED, Parma_Polyhedra_Library::Boundary_NS::sgn(), Parma_Polyhedra_Library::Linear_Expression::space_dimension(), and value.

                                                                         {
  // The dimension of `expr' must be at most the dimension of *this.
  if (space_dim < expr.space_dimension())
    throw_dimension_incompatible("frequency(e, ...)", "e", expr);

  // If the `expr' has a constant value, then the frequency
  // `freq_n' is 0. Otherwise the values for \p expr are not discrete
  // and we return false.

  // Space dimension is 0: if empty, then return false;
  // otherwise the frequency is 1 and the value is 0.
  if (space_dim == 0) {
    if (is_empty())
      return false;
    freq_n = 0;
    freq_d = 1;
    val_n = expr.inhomogeneous_term();
    val_d = 1;
    return true;
  }

  // For an empty polyhedron, we simply return false.
  if (marked_empty()
      || (has_pending_constraints() && !process_pending_constraints())
      || (!generators_are_up_to_date() && !update_generators()))
    return false;

  // The polyhedron has updated, possibly pending generators.
  // The following loop will iterate through the generator
  // to see if `expr' has a constant value.
  PPL_DIRTY_TEMP(mpq_class, value);

  // True if we have no other candidate value to compare with.
  bool first_candidate = true;

  PPL_DIRTY_TEMP_COEFFICIENT(sp);
  PPL_DIRTY_TEMP(mpq_class, candidate);
  for (dimension_type i = gen_sys.num_rows(); i-- > 0; ) {
    const Generator& gen_sys_i = gen_sys[i];
    Scalar_Products::homogeneous_assign(sp, expr, gen_sys_i);
    // Lines and rays in `*this' can cause `expr' to be non-constant.
    if (gen_sys_i.is_line_or_ray()) {
      const int sp_sign = sgn(sp);
      if (sp_sign != 0)
        // `expr' is unbounded in `*this'.
        return false;
    }
    else {
      // We have a point or a closure point.
      PPL_ASSERT(gen_sys_i.is_point() || gen_sys_i.is_closure_point());
      // Notice that we are ignoring the constant term in `expr' here.
      // We will add it to the value if there is a constant value.
      assign_r(candidate.get_num(), sp, ROUND_NOT_NEEDED);
      assign_r(candidate.get_den(), gen_sys_i[0], ROUND_NOT_NEEDED);
      candidate.canonicalize();
      if (first_candidate) {
        // We have a (new) candidate value.
        first_candidate = false;
        value = candidate;
      }
      else if (candidate != value)
        return false;
    }
  }

  // Add in the constant term in `expr'.
  PPL_DIRTY_TEMP(mpz_class, n);
  assign_r(n, expr.inhomogeneous_term(), ROUND_NOT_NEEDED);
  value += n;
  // FIXME: avoid these temporaries, if possible.
  // This can be done adding an `assign' function working on native
  // and checked or an operator= that have on one side a checked and
  // on the other a native or checked.
  // The reason why now we can't use operator= is the fact that we
  // still can have Coefficient defined to mpz_class (and not
  // Checked_Number<mpz_class>).
  val_n = Coefficient(value.get_num());
  val_d = Coefficient(value.get_den());

  freq_n = 0;
  freq_d = 1;
  return true;
}

void Parma_Polyhedra_Library::Polyhedron::generalized_affine_image	(	Variable	var,
		Relation_Symbol	relsym,
		const Linear_Expression &	expr,
		Coefficient_traits::const_reference	denominator = `Coefficient_one()`
	)

Assigns to *this the image of *this with respect to the generalized affine relation $\mathrm{var}' \relsym \frac{\mathrm{expr}}{\mathrm{denominator}}$ , where $\mathord{\relsym}$ is the relation symbol encoded by relsym.

Parameters:

var	The left hand side variable of the generalized affine relation;
relsym	The relation symbol;
expr	The numerator of the right hand side affine expression;
denominator	The denominator of the right hand side affine expression (optional argument with default value 1).

Exceptions:

std::invalid_argument Thrown if denominator is zero or if expr and *this are dimension-incompatible or if var is not a space dimension of *this or if *this is a C_Polyhedron and relsym is a strict relation symbol.

Definition at line 2864 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::EQUAL, Parma_Polyhedra_Library::GREATER_OR_EQUAL, Parma_Polyhedra_Library::GREATER_THAN, Parma_Polyhedra_Library::LESS_OR_EQUAL, Parma_Polyhedra_Library::LESS_THAN, Parma_Polyhedra_Library::NOT_EQUAL, PPL_ASSERT, PPL_ASSERT_HEAVY, PPL_UNREACHABLE, Parma_Polyhedra_Library::Variable::space_dimension(), and Parma_Polyhedra_Library::Linear_Expression::space_dimension().

                                                                        {
  // The denominator cannot be zero.
  if (denominator == 0)
    throw_invalid_argument("generalized_affine_image(v, r, e, d)", "d == 0");

  // Dimension-compatibility checks.
  // The dimension of `expr' should not be greater than the dimension
  // of `*this'.
  if (space_dim < expr.space_dimension())
    throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
                                 "e", expr);
  // `var' should be one of the dimensions of the polyhedron.
  const dimension_type var_space_dim = var.space_dimension();
  if (space_dim < var_space_dim)
    throw_dimension_incompatible("generalized_affine_image(v, r, e, d)",
                                 "v", var);

  // Strict relation symbols are only admitted for NNC polyhedra.
  if (is_necessarily_closed()
      && (relsym == LESS_THAN || relsym == GREATER_THAN))
    throw_invalid_argument("generalized_affine_image(v, r, e, d)",
                           "r is a strict relation symbol");
  // The relation symbol cannot be a disequality.
  if (relsym == NOT_EQUAL)
    throw_invalid_argument("generalized_affine_image(v, r, e, d)",
                           "r is the disequality relation symbol");

  // First compute the affine image.
  affine_image(var, expr, denominator);

  if (relsym == EQUAL)
    // The affine relation is indeed an affine function.
    return;

  // Any image of an empty polyhedron is empty.
  // Note: DO check for emptiness here, as we will later add a ray.
  if (is_empty())
    return;

  switch (relsym) {
  case LESS_OR_EQUAL:
    add_generator(ray(-var));
    break;
  case GREATER_OR_EQUAL:
    add_generator(ray(var));
    break;
  case LESS_THAN:
  // Intentionally fall through.
  case GREATER_THAN:
    {
      // The relation symbol is strict.
      PPL_ASSERT(!is_necessarily_closed());
      // While adding the ray, we minimize the generators
      // in order to avoid adding too many redundant generators later.
      add_generator(ray((relsym == GREATER_THAN) ? var : -var));
      minimize();
      // We split each point of the generator system into two generators:
      // a closure point, having the same coordinates of the given point,
      // and another point, having the same coordinates for all but the
      // `var' dimension, which is displaced along the direction of the
      // newly introduced ray.
      const dimension_type eps_index = space_dim + 1;
      for (dimension_type i =  gen_sys.num_rows(); i-- > 0; )
        if (gen_sys[i].is_point()) {
          // Add a `var'-displaced copy of `gen_sys[i]' to the generator system.
          gen_sys.add_row(gen_sys[i]);
          if (relsym == GREATER_THAN)
            ++gen_sys[gen_sys.num_rows()-1][var_space_dim];
          else
            --gen_sys[gen_sys.num_rows()-1][var_space_dim];
          // Transform `gen_sys[i]' into a closure point.
          gen_sys[i][eps_index] = 0;
        }
      clear_constraints_up_to_date();
      clear_generators_minimized();
      gen_sys.set_sorted(false);
      clear_sat_c_up_to_date();
      clear_sat_g_up_to_date();
    }
    break;
  default:
    // The EQUAL and NOT_EQUAL cases have been already dealt with.
    PPL_UNREACHABLE;
    break;
  }
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::generalized_affine_image	(	const Linear_Expression &	lhs,
		Relation_Symbol	relsym,
		const Linear_Expression &	rhs
	)

Assigns to *this the image of *this with respect to the generalized affine relation $\mathrm{lhs}' \relsym \mathrm{rhs}$ , where $\mathord{\relsym}$ is the relation symbol encoded by relsym.

Parameters:

lhs	The left hand side affine expression;
relsym	The relation symbol;
rhs	The right hand side affine expression.

Exceptions:

std::invalid_argument Thrown if *this is dimension-incompatible with lhs or rhs or if *this is a C_Polyhedron and relsym is a strict relation symbol.

Definition at line 3060 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Linear_Expression::coefficient(), Parma_Polyhedra_Library::EQUAL, Parma_Polyhedra_Library::GREATER_OR_EQUAL, Parma_Polyhedra_Library::GREATER_THAN, Parma_Polyhedra_Library::Generator_System::insert(), Parma_Polyhedra_Library::LESS_OR_EQUAL, Parma_Polyhedra_Library::LESS_THAN, Parma_Polyhedra_Library::NOT_EQUAL, PPL_ASSERT_HEAVY, PPL_UNREACHABLE, and Parma_Polyhedra_Library::Linear_Expression::space_dimension().

                                                                        {
  // Dimension-compatibility checks.
  // The dimension of `lhs' should not be greater than the dimension
  // of `*this'.
  dimension_type lhs_space_dim = lhs.space_dimension();
  if (space_dim < lhs_space_dim)
    throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
                                 "e1", lhs);
  // The dimension of `rhs' should not be greater than the dimension
  // of `*this'.
  const dimension_type rhs_space_dim = rhs.space_dimension();
  if (space_dim < rhs_space_dim)
    throw_dimension_incompatible("generalized_affine_image(e1, r, e2)",
                                 "e2", rhs);

  // Strict relation symbols are only admitted for NNC polyhedra.
  if (is_necessarily_closed()
      && (relsym == LESS_THAN || relsym == GREATER_THAN))
    throw_invalid_argument("generalized_affine_image(e1, r, e2)",
                           "r is a strict relation symbol");
  // The relation symbol cannot be a disequality.
  if (relsym == NOT_EQUAL)
    throw_invalid_argument("generalized_affine_image(e1, r, e2)",
                           "r is the disequality relation symbol");

  // Any image of an empty polyhedron is empty.
  if (marked_empty())
    return;

  // Compute the actual space dimension of `lhs',
  // i.e., the highest dimension having a non-zero coefficient in `lhs'.
  for ( ; lhs_space_dim > 0; lhs_space_dim--)
    if (lhs.coefficient(Variable(lhs_space_dim - 1)) != 0)
      break;
  // If all variables have a zero coefficient, then `lhs' is a constant:
  // we can simply add the constraint `lhs relsym rhs'.
  if (lhs_space_dim == 0) {
    switch (relsym) {
    case LESS_THAN:
      refine_no_check(lhs < rhs);
      break;
    case LESS_OR_EQUAL:
      refine_no_check(lhs <= rhs);
      break;
    case EQUAL:
      refine_no_check(lhs == rhs);
      break;
    case GREATER_OR_EQUAL:
      refine_no_check(lhs >= rhs);
      break;
    case GREATER_THAN:
      refine_no_check(lhs > rhs);
      break;
    case NOT_EQUAL:
      // The NOT_EQUAL case has been already dealt with.
      PPL_UNREACHABLE;
      break;
    }
    return;
  }

  // Gather in `new_lines' the collections of all the lines having
  // the direction of variables occurring in `lhs'.
  // While at it, check whether or not there exists a variable
  // occurring in both `lhs' and `rhs'.
  Generator_System new_lines;
  bool lhs_vars_intersects_rhs_vars = false;
  for (dimension_type i = lhs_space_dim; i-- > 0; )
    if (lhs.coefficient(Variable(i)) != 0) {
      new_lines.insert(line(Variable(i)));
      if (rhs.coefficient(Variable(i)) != 0)
        lhs_vars_intersects_rhs_vars = true;
    }

  if (lhs_vars_intersects_rhs_vars) {
    // Some variables in `lhs' also occur in `rhs'.
    // To ease the computation, we add an additional dimension.
    const Variable new_var(space_dim);
    add_space_dimensions_and_embed(1);

    // Constrain the new dimension to be equal to the right hand side.
    // (check for emptiness because we will add lines).
    refine_no_check(new_var == rhs);
    if (!is_empty()) {
      // Existentially quantify the variables in the left hand side.
      add_recycled_generators(new_lines);

      // Constrain the new dimension so that it is related to
      // the left hand side as dictated by `relsym'
      // (we force minimization because we will need the generators).
      switch (relsym) {
      case LESS_THAN:
        refine_no_check(lhs < new_var);
        break;
      case LESS_OR_EQUAL:
        refine_no_check(lhs <= new_var);
        break;
      case EQUAL:
        refine_no_check(lhs == new_var);
        break;
      case GREATER_OR_EQUAL:
        refine_no_check(lhs >= new_var);
        break;
      case GREATER_THAN:
        refine_no_check(lhs > new_var);
        break;
      case NOT_EQUAL:
        // The NOT_EQUAL case has been already dealt with.
        PPL_UNREACHABLE;
        break;
      }
    }
    // Remove the temporarily added dimension.
    remove_higher_space_dimensions(space_dim-1);
  }
  else {
    // `lhs' and `rhs' variables are disjoint:
    // there is no need to add a further dimension.

    // Any image of an empty polyhedron is empty.
    // Note: DO check for emptiness here, as we will add lines.
    if (is_empty())
      return;

    // Existentially quantify the variables in the left hand side.
    add_recycled_generators(new_lines);

    // Constrain the left hand side expression so that it is related to
    // the right hand side expression as dictated by `relsym'.
    switch (relsym) {
    case LESS_THAN:
      refine_no_check(lhs < rhs);
      break;
    case LESS_OR_EQUAL:
      refine_no_check(lhs <= rhs);
      break;
    case EQUAL:
      refine_no_check(lhs == rhs);
      break;
    case GREATER_OR_EQUAL:
      refine_no_check(lhs >= rhs);
      break;
    case GREATER_THAN:
      refine_no_check(lhs > rhs);
      break;
    case NOT_EQUAL:
      // The NOT_EQUAL case has been already dealt with.
      PPL_UNREACHABLE;
      break;
    }
  }
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::generalized_affine_preimage	(	Variable	var,
		Relation_Symbol	relsym,
		const Linear_Expression &	expr,
		Coefficient_traits::const_reference	denominator = `Coefficient_one()`
	)

Assigns to *this the preimage of *this with respect to the generalized affine relation $\mathrm{var}' \relsym \frac{\mathrm{expr}}{\mathrm{denominator}}$ , where $\mathord{\relsym}$ is the relation symbol encoded by relsym.

Parameters:

var	The left hand side variable of the generalized affine relation;
relsym	The relation symbol;
expr	The numerator of the right hand side affine expression;
denominator	The denominator of the right hand side affine expression (optional argument with default value 1).

Exceptions:

std::invalid_argument Thrown if denominator is zero or if expr and *this are dimension-incompatible or if var is not a space dimension of *this or if *this is a C_Polyhedron and relsym is a strict relation symbol.

Definition at line 2957 of file Polyhedron_public.cc.

                                                                           {
  // The denominator cannot be zero.
  if (denominator == 0)
    throw_invalid_argument("generalized_affine_preimage(v, r, e, d)",
                           "d == 0");

  // Dimension-compatibility checks.
  // The dimension of `expr' should not be greater than the dimension
  // of `*this'.
  if (space_dim < expr.space_dimension())
    throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
                                 "e", expr);
  // `var' should be one of the dimensions of the polyhedron.
  const dimension_type var_space_dim = var.space_dimension();
  if (space_dim < var_space_dim)
    throw_dimension_incompatible("generalized_affine_preimage(v, r, e, d)",
                                 "v", var);

  // Strict relation symbols are only admitted for NNC polyhedra.
  if (is_necessarily_closed()
      && (relsym == LESS_THAN || relsym == GREATER_THAN))
    throw_invalid_argument("generalized_affine_preimage(v, r, e, d)",
                           "r is a strict relation symbol");
  // The relation symbol cannot be a disequality.
  if (relsym == NOT_EQUAL)
    throw_invalid_argument("generalized_affine_preimage(v, r, e, d)",
                           "r is the disequality relation symbol");

  // Check whether the affine relation is indeed an affine function.
  if (relsym == EQUAL) {
    affine_preimage(var, expr, denominator);
    return;
  }

  // Compute the reversed relation symbol to simplify later coding.
  Relation_Symbol reversed_relsym;
  switch (relsym) {
  case LESS_THAN:
    reversed_relsym = GREATER_THAN;
    break;
  case LESS_OR_EQUAL:
    reversed_relsym = GREATER_OR_EQUAL;
    break;
  case GREATER_OR_EQUAL:
    reversed_relsym = LESS_OR_EQUAL;
    break;
  case GREATER_THAN:
    reversed_relsym = LESS_THAN;
    break;
  default:
    // The EQUAL and NOT_EQUAL cases have been already dealt with.
    PPL_UNREACHABLE;
    return;
  }

  // Check whether the preimage of this affine relation can be easily
  // computed as the image of its inverse relation.
  const Coefficient& var_coefficient = expr.coefficient(var);
  if (var_coefficient != 0) {
    Linear_Expression inverse_expr
      = expr - (denominator + var_coefficient) * var;
    PPL_DIRTY_TEMP_COEFFICIENT(inverse_denominator);
    neg_assign(inverse_denominator, var_coefficient);
    Relation_Symbol inverse_relsym
      = (sgn(denominator) == sgn(inverse_denominator))
      ? relsym : reversed_relsym;
    generalized_affine_image(var, inverse_relsym, inverse_expr,
                             inverse_denominator);
    return;
  }

  // Here `var_coefficient == 0', so that the preimage cannot
  // be easily computed by inverting the affine relation.
  // Shrink the polyhedron by adding the constraint induced
  // by the affine relation.
  const Relation_Symbol corrected_relsym
    = (denominator > 0) ? relsym : reversed_relsym;
  switch (corrected_relsym) {
  case LESS_THAN:
    refine_no_check(denominator*var < expr);
    break;
  case LESS_OR_EQUAL:
    refine_no_check(denominator*var <= expr);
    break;
  case GREATER_OR_EQUAL:
    refine_no_check(denominator*var >= expr);
    break;
  case GREATER_THAN:
    refine_no_check(denominator*var > expr);
    break;
  default:
    // The EQUAL and NOT_EQUAL cases have been already dealt with.
    PPL_UNREACHABLE;
    break;
  }
  unconstrain(var);
  PPL_ASSERT_HEAVY(OK());
}

void Parma_Polyhedra_Library::Polyhedron::generalized_affine_preimage	(	const Linear_Expression &	lhs,
		Relation_Symbol	relsym,
		const Linear_Expression &	rhs
	)

Assigns to *this the preimage of *this with respect to the generalized affine relation $\mathrm{lhs}' \relsym \mathrm{rhs}$ , where $\mathord{\relsym}$ is the relation symbol encoded by relsym.

Parameters:

lhs	The left hand side affine expression;
relsym	The relation symbol;
rhs	The right hand side affine expression.

Exceptions:

std::invalid_argument Thrown if *this is dimension-incompatible with lhs or rhs or if *this is a C_Polyhedron and relsym is a strict relation symbol.

Definition at line 3217 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Linear_Expression::coefficient(), Parma_Polyhedra_Library::EQUAL, Parma_Polyhedra_Library::GREATER_OR_EQUAL, Parma_Polyhedra_Library::GREATER_THAN, Parma_Polyhedra_Library::Generator_System::insert(), Parma_Polyhedra_Library::LESS_OR_EQUAL, Parma_Polyhedra_Library::LESS_THAN, Parma_Polyhedra_Library::NOT_EQUAL, PPL_ASSERT_HEAVY, PPL_UNREACHABLE, and Parma_Polyhedra_Library::Linear_Expression::space_dimension().

                                                                           {
  // Dimension-compatibility checks.
  // The dimension of `lhs' should not be greater than the dimension
  // of `*this'.
  dimension_type lhs_space_dim = lhs.space_dimension();
  if (space_dim < lhs_space_dim)
    throw_dimension_incompatible("generalized_affine_preimage(e1, r, e2)",
                                 "e1", lhs);
  // The dimension of `rhs' should not be greater than the dimension
  // of `*this'.
  const dimension_type rhs_space_dim = rhs.space_dimension();
  if (space_dim < rhs_space_dim)
    throw_dimension_incompatible("generalized_affine_preimage(e1, r, e2)",
                                 "e2", rhs);

  // Strict relation symbols are only admitted for NNC polyhedra.
  if (is_necessarily_closed()
      && (relsym == LESS_THAN || relsym == GREATER_THAN))
    throw_invalid_argument("generalized_affine_preimage(e1, r, e2)",
                           "r is a strict relation symbol");
  // The relation symbol cannot be a disequality.
  if (relsym == NOT_EQUAL)
    throw_invalid_argument("generalized_affine_preimage(e1, r, e2)",
                           "r is the disequality relation symbol");

  // Any preimage of an empty polyhedron is empty.
  if (marked_empty())
    return;

  // Compute the actual space dimension of `lhs',
  // i.e., the highest dimension having a non-zero coefficient in `lhs'.
  for ( ; lhs_space_dim > 0; lhs_space_dim--)
    if (lhs.coefficient(Variable(lhs_space_dim - 1)) != 0)
      break;

  // If all variables have a zero coefficient, then `lhs' is a constant:
  // in this case, preimage and image happen to be the same.
  if (lhs_space_dim == 0) {
    generalized_affine_image(lhs, relsym, rhs);
    return;
  }

  // Gather in `new_lines' the collections of all the lines having
  // the direction of variables occurring in `lhs'.
  // While at it, check whether or not there exists a variable
  // occurring in both `lhs' and `rhs'.
  Generator_System new_lines;
  bool lhs_vars_intersects_rhs_vars = false;
  for (dimension_type i = lhs_space_dim; i-- > 0; )
    if (lhs.coefficient(Variable(i)) != 0) {
      new_lines.insert(line(Variable(i)));
      if (rhs.coefficient(Variable(i)) != 0)
        lhs_vars_intersects_rhs_vars = true;
    }

  if (lhs_vars_intersects_rhs_vars) {
    // Some variables in `lhs' also occur in `rhs'.
    // To ease the computation, we add an additional dimension.
    const Variable new_var(space_dim);
    add_space_dimensions_and_embed(1);

    // Constrain the new dimension to be equal to `lhs'
    // (also check for emptiness because we have to add lines).
    refine_no_check(new_var == lhs);
    if (!is_empty()) {
      // Existentially quantify the variables in the left hand side.
      add_recycled_generators(new_lines);

      // Constrain the new dimension so that it is related to
      // the right hand side as dictated by `relsym'.
      switch (relsym) {
      case LESS_THAN:
        refine_no_check(new_var < rhs);
        break;
      case LESS_OR_EQUAL:
        refine_no_check(new_var <= rhs);
        break;
      case EQUAL:
        refine_no_check(new_var == rhs);
        break;
      case GREATER_OR_EQUAL:
        refine_no_check(new_var >= rhs);
        break;
      case GREATER_THAN:
        refine_no_check(new_var > rhs);
        break;
      case NOT_EQUAL:
        // The NOT_EQUAL case has been already dealt with.
        PPL_UNREACHABLE;
        break;
      }
    }
    // Remove the temporarily added dimension.
    remove_higher_space_dimensions(space_dim-1);
  }
  else {
    // `lhs' and `rhs' variables are disjoint:
    // there is no need to add a further dimension.

    // Constrain the left hand side expression so that it is related to
    // the right hand side expression as dictated by `relsym'.
    switch (relsym) {
    case LESS_THAN:
      refine_no_check(lhs < rhs);
      break;
    case LESS_OR_EQUAL:
      refine_no_check(lhs <= rhs);
      break;
    case EQUAL:
      refine_no_check(lhs == rhs);
      break;
    case GREATER_OR_EQUAL:
      refine_no_check(lhs >= rhs);
      break;
    case GREATER_THAN:
      refine_no_check(lhs > rhs);
      break;
    case NOT_EQUAL:
      // The NOT_EQUAL case has been already dealt with.
      PPL_UNREACHABLE;
      break;
    }
    // Any image of an empty polyhedron is empty.
    // Note: DO check for emptiness here, as we will add lines.
    if (is_empty())
      return;
    // Existentially quantify all the variables occurring in `lhs'.
    add_recycled_generators(new_lines);
  }
  PPL_ASSERT_HEAVY(OK());
}

template<typename FP_Format , typename Interval_Info >

void Parma_Polyhedra_Library::Polyhedron::generalized_refine_with_linear_form_inequality	(	const Linear_Form< Interval< FP_Format, Interval_Info > > &	left,
		const Linear_Form< Interval< FP_Format, Interval_Info > > &	right,
		Relation_Symbol	relsym
	)

inline

Refines *this with the constraint expressed by left $\relsym$ right, where $\relsym$ is the relation symbol specified by relsym..

Parameters:

left	The linear form on intervals with floating point boundaries that is on the left of the comparison operator. All of its coefficients MUST be bounded.
right	The linear form on intervals with floating point boundaries that is on the right of the comparison operator. All of its coefficients MUST be bounded.
relsym	The relation symbol.

Exceptions:

std::invalid_argument	Thrown if `left` (or `right`) is dimension-incompatible with `*this`.
std::runtime_error	Thrown if `relsym` is not a valid relation symbol.

This function is used in abstract interpretation to model a filter that is generated by a comparison of two expressions that are correctly approximated by left and right respectively.

Definition at line 379 of file Polyhedron.inlines.hh.

References Parma_Polyhedra_Library::EQUAL, Parma_Polyhedra_Library::GREATER_OR_EQUAL, Parma_Polyhedra_Library::GREATER_THAN, Parma_Polyhedra_Library::LESS_OR_EQUAL, Parma_Polyhedra_Library::LESS_THAN, Parma_Polyhedra_Library::NOT_EQUAL, PPL_UNREACHABLE, and refine_with_linear_form_inequality().

                                          {
  switch (relsym) {
  case EQUAL:
    // TODO: see if we can handle this case more efficiently.
    refine_with_linear_form_inequality(left, right, false);
    refine_with_linear_form_inequality(right, left, false);
    break;
  case LESS_THAN:
    refine_with_linear_form_inequality(left, right, true);
    break;
  case LESS_OR_EQUAL:
    refine_with_linear_form_inequality(left, right, false);
    break;
  case GREATER_THAN:
    refine_with_linear_form_inequality(right, left, true);
    break;
  case GREATER_OR_EQUAL:
    refine_with_linear_form_inequality(right, left, false);
    break;
  case NOT_EQUAL:
    break;
  default:
    PPL_UNREACHABLE;
    break;
  }
}

const PPL::Generator_System & Parma_Polyhedra_Library::Polyhedron::generators ( ) const

Returns the system of generators.

Definition at line 130 of file Polyhedron_public.cc.

References Parma_Polyhedra_Library::Generator_System::adjust_topology_and_space_dimension(), PPL_ASSERT, and Parma_Polyhedra_Library::Generator_System::zero_dim_univ().

Referenced by Parma_Polyhedra_Library::Implementation::Termination::all_affine_ranking_functions_PR(), Parma_Polyhedra_Library::Implementation::Termination::all_affine_ranking_functions_PR_original(), Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), Parma_Polyhedra_Library::Box< ITV >::Box(), Parma_Polyhedra_Library::Grid::Grid(), map_space_dimensions(), and Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape().

                                {
  if (marked_empty()) {
    PPL_ASSERT(gen_sys.has_no_rows());
    // We want `gen_sys' to have the appropriate space dimension,
    // even though it is an empty generator system.
    if (gen_sys.space_dimension() != space_dim) {
      Generator_System gs;
      gs.adjust_topology_and_space_dimension(topology(), space_dim);
      const_cast<Generator_System&>(gen_sys).m_swap(gs);
    }
    return gen_sys;
  }

  if (space_dim == 0) {
    PPL_ASSERT(gen_sys.num_rows() == 0 && gen_sys.num_columns() == 0);
    return Generator_System::zero_dim_univ();
  }

  // If the polyhedron has pending constraints, we process them to obtain
  // the generators (we may discover that the polyhedron is empty).
  // No processing is needed if the polyhedron has pending generators.
  if ((has_pending_constraints() && !process_pending_constraints())
      || (!generators_are_up_to_date() && !update_generators())) {
    // We have just discovered that `*this' is empty.
    PPL_ASSERT(gen_sys.has_no_rows());
    // We want `gen_sys' to have the appropriate space dimension,
    // even though it is an empty generator system.
    if (gen_sys.space_dimension() != space_dim) {
      Generator_System gs;
      gs.adjust_topology_and_space_dimension(topology(), space_dim);
      const_cast<Generator_System&>(gen_sys).m_swap(gs);
    }
    return gen_sys;
  }

  // TODO: reconsider whether to really sort generators at this stage.
#if ENSURE_SORTEDNESS
  // We insist in returning a sorted system of generators,
  // but sorting is useless if there are pending generators.
  if (!has_pending_generators())
    obtain_sorted_generators();
#else
  // In the case of an NNC polyhedron, if the generator system is fully
  // minimized (i.e., minimized and with no pending generator), then
  // return a sorted system of generators: this is needed so that the
  // const_iterator could correctly filter out the matched closure points.
  if (!is_necessarily_closed()
      && generators_are_minimized() && !

Classes

Public Types

Public Member Functions

Static Public Member Functions

Protected Member Functions

Static Protected Member Functions

Private Types

Private Member Functions

Static Private Member Functions

Private Attributes

Static Private Attributes

Friends

Related Functions

Exception Throwers

Detailed Description

Member Typedef Documentation

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation