/* Polyhedron class implementation: conversion().
Copyright (C) 2001-2004 Roberto Bagnara <bagnara@cs.unipr.it>
This file is part of the Parma Polyhedra Library (PPL).
The PPL is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
The PPL is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
USA.
For the most up-to-date information see the Parma Polyhedra Library
site: http://www.cs.unipr.it/ppl/ . */
#include <config.h>
#include "Matrix.defs.hh"
#include "SatMatrix.defs.hh"
#include "Polyhedron.defs.hh"
#include "globals.hh"
#include <cstddef>
namespace PPL = Parma_Polyhedra_Library;
// True if abandon_expensive_computations should be checked often,
// where the meaning of "often" is as stated in the documentation
// of that variable.
#define REACTIVE_ABANDONING 1
/*!
\return
The number of lines of the polyhedron or the number of equality
constraints in the result of conversion.
\param source
The matrix to use to convert \p dest: it may be modified;
\param start
The index of \p source row from which conversion begin;
\param dest
The result of the conversion;
\param sat
The saturation matrix telling us, for each row in \p source, which
are the rows of \p dest that satisfy but do not saturate it;
\param num_lines_or_equalities
The number of rows in the matrix \p dest that are either lines of
the polyhedron (when \p dest is a system of generators) or equality
constraints (when \p dest is a system of constraints).
For simplicity, all the following comments assume we are converting a
constraint system \p source to a generator system \p dest;
the comments for the symmetric case can be obtained by duality.
If some of the constraints in \p source are redundant, they will be removed.
This is why the \p source is not declared to be a constant parameter.
If \p start is 0, then \p source is a sorted matrix; also, \p dest is
a generator system corresponding to an empty constraint system.
If otherwise \p start is greater than 0, then the two sub-matrices of
\p source made by the non-pending rows and the pending rows, respectively,
are both sorted; also, \p dest is the generator system corresponding to
the non-pending constraints of \p source.
Independently from the value of \p start, \p dest has lines from index 0
to index \p num_lines_or_equalities - 1 and rays/points from index
\p num_lines_or_equalities to the last of its rows.
Note that here the rows of \p sat are indexed by rows of \p dest
and its columns are indexed by rows of \p 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 \ref representation).
When we have both descriptions for a polyhedron \f$P\f$
we have what is called a <EM>double description</EM>
(or <EM>DD pair</EM>) for \f$P\f$.
Here, the <EM>representation matrix</EM> refers to the matrix \f$C\f$
whose rows represent the constraints that characterize \f$P\f$
and the <EM>generating matrix</EM>, the matrix \f$G\f$ whose rows
represent the generators of \f$P\f$.
We say that a pair \f$(C, G)\f$ of (real) matrices is
a <EM>double description pair</EM> if
\f[
C\vect{x} \geq \vect{0}
\quad\iff\quad
\exists \vect{\lambda} \geq \vect{0} \mathrel{.}
\vect{x} = G\vect{\lambda}.
\f]
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 \f$P\f$
and we call \f$C\f$ and \f$G\f$ the matrices of constraints and
rays respectively, we have
\f[
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}\, \}.
\f]
Because of the theorem of Minkowski (see Section \ref prelims),
we can say that, given a \f$m \times n\f$ representation matrix
\f$C\f$ such that \f$\mathop{\mathrm{rank}}(C) = n = dimension of
the whole space \f$ for a non-empty polyhedron \f$P\f$,
it is always possible to find a generating matrix \f$G\f$ for \f$P\f$
such that \f$(C, G)\f$ is a DD pair.
Conversely, Weil's theorem ensures that, for each generating matrix
\f$G\f$, it is possible to find a representation matrix \f$C\f$
such that \f$(C, G)\f$ 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 \f$C\f$ find \f$G\f$ such that \f$(C, G)\f$ is a DD pair;
-# given \f$G\f$ find \f$C\f$ such that \f$(C, G)\f$ 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
\f$\Rset^n\f$ that, in a sense, captures the essence of duality.
Consider a matrix \f$A \in \Rset^{m \times n}\f$ and let
\f$\{ \vect{a}_1, \ldots, \vect{a}_m \}\f$ be its set of row vectors.
Consider also another vector \f$\vect{c} \in \Rset^n\f$ such that,
whenever a vector \f$\vect{y} \in \Rset^n\f$ has a non-negative projection
on the \f$\vect{a}_i\f$'s, it also has a non-negative projection
on \f$\vect{c}\f$.
The lemma states that \f$\vect{c}\f$ has this property if and only if
it is in the cone generated by the \f$\vect{a}_i\f$'s.
Formally, the lemma states the equivalence of the two following
assertions:
-# \f$
\forall \vect{y}
\mathrel{:} (A\vect{y} \geq 0 \implies
\langle \vect{y},\vect{c} \rangle \geq 0)
\f$;
-# \f$
\exists \vect{\lambda} \geq \vect{0}
\mathrel{.} \vect{c}^\mathrm{T} = \vect{\lambda}^\mathrm{T}A
\f$.
With this result we can prove that \f$(C, G)\f$ is a DD pair
if and only if \f$(G^\mathrm{T}, C^\mathrm{T})\f$ is a DD pair.
Suppose \f$(C, G)\f$ is a DD pair.
Thus, for each \f$x\f$ of the appropriate dimension,
\f$C\vect{x} \geq \vect{0}\f$ if and only if
\f$\exists \lambda \geq 0 \mathrel{.} \vect{x} = G\vect{\lambda}\f$,
which is of course equivalent to
\f$
\exists \vect{\lambda} \geq \vect{0}
\mathrel{.} \vect{x}^\mathrm{T} = \vect{\lambda}^\mathrm{T}G^\mathrm{T}
\f$.
First, we assume that \f$\vect{z}\f$ is such that
\f$G^\mathrm{T}\vect{z} \geq \vect{0}\f$
and we will show that
\f$\exists \vect{\mu} \geq \vect{0} \mathrel{.}
\vect{z} = C^\mathrm{T}\vect{\mu}\f$.
Let \f$\vect{x}\f$ be such that \f$C\vect{x} \geq \vect{0}\f$.
Since \f$(C, G)\f$ is a DD pair, this is equivalent to
\f$
\exists \vect{\lambda} \geq \vect{0}
\mathrel{.} \vect{x}^\mathrm{T} = \vect{\lambda}^\mathrm{T}G^\mathrm{T}
\f$,
which, by Farkas' lemma is equivalent to
\f$
\forall \vect{y} \mathrel{:} (G^\mathrm{T}\vect{y} \geq \vect{0} \implies
\langle \vect{y}, \vect{x} \rangle \geq 0)
\f$.
Taking \f$\vect{y} = \vect{z}\f$ and recalling our assumption that
\f$G^\mathrm{T}\vect{z} \geq \vect{0}\f$
we can conclude that \f$\langle \vect{z}, \vect{x} \rangle \geq 0\f$,
that is equivalent to \f$\langle \vect{x}, \vect{z} \rangle \geq 0\f$.
We have thus established that
\f$
\forall \vect{x} \mathrel{:} (C\vect{x} \geq \vect{0} \implies
\langle \vect{x}, \vect{z} \rangle \geq 0)
\f$.
By Farkas' lemma, this is equivalent to
\f$\exists \vect{\mu} \geq \vect{0} \mathrel{.}
\vect{z}^\mathrm{T} = \vect{\mu}^\mathrm{T} C\f$,
which is equivalent to what we wanted to prove, that is,
\f$\exists \vect{\mu} \geq \vect{0} \mathrel{.}
\vect{z} = C^\mathrm{T}\vect{\mu}\f$.
In order to prove the reverse implication, the following observation
turns out to be useful:
when \f$(C, G)\f$ is a DD pair, \f$CG \geq 0\f$.
In fact,
let \f$\vect{e}_j\f$ be the vector whose components are all \f$0\f$
apart from the \f$j\f$-th one, which is \f$1\f$.
Clearly \f$\vect{e}_j \geq \vect{0}\f$ and, taking
\f$\vect{\lambda} = \vect{e}_j\f$ and
\f$\vect{x} = G\vect{\lambda} = G \vect{e}_j\f$, we have
\f$C\vect{x} = C(G \vect{e}_j) = (CG)\vect{e}_j \geq \vect{0}\f$,
since \f$(C, G)\f$ is a DD pair.
Thus, as \f$(CG)\vect{e}_j\f$ is the \f$j\f$-th column of \f$CG\f$
and since the choice of \f$j\f$ was arbitrary, \f$CG \geq \vect{0}\f$.
We now assume that \f$\vect{z}\f$ is such that
\f$\exists \vect{\mu} \geq \vect{0} \mathrel{.}
\vect{z} = C^\mathrm{T}\vect{\mu}\f$
and we will prove that \f$G^\mathrm{T}\vect{z} \geq \vect{0}\f$.
By Farkas' lemma, the assumption
\f$\exists \vect{\mu} \geq \vect{0} \mathrel{.}
\vect{z}^\mathrm{T} = \vect{\mu}^\mathrm{T}C\f$,
is equivalent to
\f$\forall \vect{y} \mathrel{:} (C\vect{y} \geq \vect{0}
\implies \langle \vect{y}, \vect{z} \rangle \geq 0)\f$.
If we take \f$\vect{y} = G\vect{e}_j\f$ then \f$C\vect{y}
= CG\vect{e}_j \geq 0\f$,
since \f$CG \geq \vect{0}\f$.
So
\f$
\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
\f$,
that is, the \f$j\f$-th component of \f$G^\mathrm{T}\vect{z}\f$
is non-negative. The arbitrary choice of \f$j\f$ allows us to conclude
that \f$G^\mathrm{T}\vect{z} \geq \vect{0}\f$, 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 \p source and \p dest can be interchanged,
in the sequel we assume the \p source matrix will contain the constraints
that represent the polyhedron and the \p dest matrix will contain
the generator that generates it.
There are some observations that are useful to understand this function:
Observation 1: Let \f$A\f$ be a matrix of constraints that generate
the polyhedron \f$P\f$ and \f$\vect{c}\f$ a new constraint that must
be added. Suppose that there is a line \f$\vect{z}\f$ that does not
saturate the constraint \f$\vect{c}\f$. If we combine the old lines
and rays that do not saturate \f$\vect{c}\f$ (except \f$\vect{z}\f$)
with \f$\vect{z}\f$ such that the new ones saturate \f$\vect{c}\f$,
the new lines and rays also saturate the constraints saturated by
the old constraints.
In fact, if \f$\vect{y}_1\f$ is the old generator that does not saturate
\f$\vect{c}\f$, \f$\vect{y}_2\f$ is the new one such that
\f[
\vect{y}_2 = \lambda \vect{y}_1 + \mu \vect{z}
\f]
and \f$\vect{c}_1\f$ is a previous constraint that \f$\vect{y}_1\f$
and \f$\vect{z}\f$ saturates, we can see
\f[
\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
\f]
and
\f[
\mu \langle \vect{c}_1, \vect{z} \rangle = 0.
\f]
Proposition 1: Let \f$\vect{r}_1\f$ and \f$\vect{r}_2\f$ be distinct
rays of \f$P\f$.
Then the following statements are equivalent:
a) \f$\vect{r}_1\f$ and \f$\vect{r}_2\f$ are adjacent extreme rays
(see Section \ref prelims);
b) \f$\vect{r}_1\f$ and \f$\vect{r}_2\f$ are extreme rays and the
rank of the matrix composed by the constraints saturated by both
\f$\vect{r}_1\f$ and \f$\vect{r}_2\f$ is equal to
\f$d - 2\f$, where \f$d\f$ is the rank of the matrix of constraints.
In fact, let \f$F\f$ be the system of generators that saturate the
constraints saturated by both \f$\vect{r}_1\f$ and \f$\vect{r}_2\f$.
If b) holds, the set \f$F\f$ is 2-dimensional and \f$\vect{r}_1\f$ and
\f$\vect{r}_2\f$ generate this set. So, every generator
\f$\vect{x}\f$ of \f$F\f$ can be built as a combination of
\f$\vect{r}_1\f$ and \f$\vect{r}_2\f$, i.e.
\f[
\vect{x} = \lambda \vect{r}_1 + \mu \vect{r}_2.
\f]
This combination is non-negative because there exists at least a
constraint \f$c\f$ saturated by \f$\vect{r}_1\f$ and not
\f$\vect{r}_2\f$ (or vice versa) (because they are distinct) for which
\f[
\langle \vect{c}, \vect{x} \rangle \geq 0
\f]
and
\f[
\langle \vect{c}, \vect{x} \rangle
= \lambda \langle \vect{c}, \vect{r}_1 \rangle
(or = \mu \langle \vect{c}, \vect{r}_2 \rangle).
\f]
So, there is no other extreme ray in \f$F\f$ and a) holds.
Otherwise, if b) does not hold, the rank of the matrix generated by
the constraints saturated by both \f$\vect{r}_1\f$ and \f$\vect{r}_2\f$
is equal to \f$d - k\f$, with \p k \>= 3, the set \f$F\f$ is
\p k -dimensional and at least \p k extreme rays are necessary
to generate \f$F\f$.
So, \f$\vect{r}_1\f$ and \f$\vect{r}_2\f$ are not adjacent and
a) does not hold.
Proposition 2: When we build the new system of generators starting a
matrix \f$A\f$ of constraints of \f$P\f$, if \f$\vect{c}\f$ is the
constraint to add to \f$A\f$ and all lines of \f$P\f$ saturate
\f$\vect{c}\f$, the new set of rays is the union of those rays that
saturate, of those that satisfy and of a set \f$\overline Q\f$ 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 \f$\vect{r}_1\f$ and
\f$\vect{r}_2\f$ 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 \f$\vect{r}_1\f$ and \f$\vect{r}_2\f$ are not adjacent,
the rank of the matrix composed by the constraints saturated by both
\f$\vect{r}_1\f$ and \f$\vect{r}_2\f$ is different from \f$d - 2\f$
(see the previous proposition) or neither \f$\vect{r}_1\f$ nor
\f$\vect{r}_2\f$ are extreme rays. Since the new ray \f$\vect{r}\f$
is a combination of \f$\vect{r}_1\f$ and \f$\vect{r}_2\f$,
it saturates the same constraints saturated by both \f$\vect{r}_1\f$ and
\f$\vect{r}_2\f$.
If the rank is less than \f$d - 2\f$, the rank of
the matrix composed by \f$\vect{c}\f$ (that is saturated by \f$\vect{r}\f$)
and by the constraints of \f$A\f$ saturated by \f$\vect{r}\f$ is less
than \f$d - 1\f$. It means that \f$r\f$ is redundant (see
Section \ref prelims).
If neither \f$\vect{r}_1\f$ nor \f$\vect{r}_2\f$ are extreme rays,
they belong to a 2-dimensional face containing exactly two extreme rays
of \f$P\f$.
These two adjacent rays built a ray equal to \f$\vect{r}\f$ and so
\f$\vect{r}\f$ is redundant.
*/
PPL::dimension_type
PPL::Polyhedron::conversion(Matrix& source,
const dimension_type start,
Matrix& dest,
SatMatrix& sat,
dimension_type num_lines_or_equalities) {
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();
// 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'.
assert(source_num_rows == sat.num_columns());
assert(dest_num_rows == sat.num_rows());
// If `start > 0', then we are converting the pending constraints.
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;
// Converting the sub-matrix 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).
std::swap(source[k], source[k+source_num_redundant]);
Row& source_k = source[k];
// Constraints and generators must have the same dimension,
// otherwise the scalar product below will bomb.
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 iff the generator saturates the constraint.
static std::vector<Integer> scalar_prod;
const int needed_space = dest_num_rows - scalar_prod.size();
if (needed_space > 0)
scalar_prod.insert(scalar_prod.end(), needed_space, Integer_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) {
scalar_prod[index_non_zero] = source_k * dest[index_non_zero];
if (scalar_prod[index_non_zero] != 0)
// The generator does not saturate the constraint.
break;
#if REACTIVE_ABANDONING
// Check if the client has requested abandoning all exponential
// computations. If so, the exception specified by the client
// is thrown now.
maybe_abandon();
#endif
}
for (dimension_type i = index_non_zero + 1; i < dest_num_rows; ++i) {
scalar_prod[i] = source_k * dest[i];
#if REACTIVE_ABANDONING
maybe_abandon();
#endif
}
// 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.
negate(scalar_prod[index_non_zero]);
for (dimension_type j = dest_num_columns; j-- > 0; )
negate(dest[index_non_zero][j]);
}
// Having changed a line to a ray, we set `dest' to be a
// non-sorted matrix, 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) {
std::swap(dest[index_non_zero],
dest[num_lines_or_equalities]);
std::swap(scalar_prod[index_non_zero],
scalar_prod[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[num_lines_or_equalities]'
// 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.
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
//
// Integer scale = scalar_prod[i];
// scale.gcd_assign(scalar_prod[num_lines_or_equalities]);
// Integer scaled_sp_i = scalar_prod[i] / scale;
// Integer scaled_sp_n
// = scalar_prod[num_lines_or_equalities] / scale;
// for (dimension_type c = dest_num_columns; c-- > 0; ) {
// dest[i][c] *= scaled_sp_n;
// dest[i][c] -= scaled_sp_i * dest[num_lines_or_equalities][c];
// }
gcd_assign(tmp_Integer[1],
scalar_prod[i],
scalar_prod[num_lines_or_equalities]);
exact_div_assign(tmp_Integer[2],
scalar_prod[i],
tmp_Integer[1]);
exact_div_assign(tmp_Integer[3],
scalar_prod[num_lines_or_equalities],
tmp_Integer[1]);
for (dimension_type c = dest_num_columns; c-- > 0; ) {
tmp_Integer[4] = tmp_Integer[3] * dest[i][c];
tmp_Integer[5] = tmp_Integer[2]
* dest[num_lines_or_equalities][c];
dest[i][c] = tmp_Integer[4] - tmp_Integer[5];
}
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[num_lines_or_equalities]'
// 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
//
// Integer scale = scalar_prod[i];
// scale.gcd_assign(scalar_prod[num_lines_or_equalities]);
// Integer scaled_sp_i = scalar_prod[i] / scale;
// Integer scaled_sp_n
// = scalar_prod[num_lines_or_equalities] / scale;
// for (dimension_type c = dest_num_columns; c-- > 0; ) {
// dest[i][c] *= scaled_sp_n;
// dest[i][c] -= scaled_sp_i * dest[num_lines_or_equalities][c];
// }
gcd_assign(tmp_Integer[1],
scalar_prod[i],
scalar_prod[num_lines_or_equalities]);
exact_div_assign(tmp_Integer[2],
scalar_prod[i],
tmp_Integer[1]);
exact_div_assign(tmp_Integer[3],
scalar_prod[num_lines_or_equalities],
tmp_Integer[1]);
for (dimension_type c = dest_num_columns; c-- > 0; ) {
tmp_Integer[4] = tmp_Integer[3] * dest[i][c];
tmp_Integer[5] = tmp_Integer[2]
* dest[num_lines_or_equalities][c];
dest[i][c] = tmp_Integer[4] - tmp_Integer[5];
}
dest[i].strong_normalize();
scalar_prod[i] = 0;
// `dest' has already been set as non-sorted.
}
#if REACTIVE_ABANDONING
maybe_abandon();
#endif
}
// Since the `scalar_prod[num_lines_or_equalities]' 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' ...
if (source_k.is_ray_or_point_or_inequality())
sat[num_lines_or_equalities].set(k);
// ... otherwise, the constraint is an equality which is
// violated by the generator `dest[num_lines_or_equalities]':
// the generator has to be removed from `dest'.
else {
--dest_num_rows;
std::swap(dest[num_lines_or_equalities],
dest[dest_num_rows]);
std::swap(scalar_prod[dest_num_rows],
scalar_prod[num_lines_or_equalities]);
std::swap(sat[num_lines_or_equalities],
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=.
std::swap(dest[sup_bound], dest[lines_or_equal_bound]);
std::swap(scalar_prod[sup_bound], scalar_prod[lines_or_equal_bound]);
std::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;
std::swap(dest[sup_bound], dest[inf_bound]);
std::swap(scalar_prod[sup_bound], scalar_prod[inf_bound]);
std::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.
SatRow new_satrow;
assert(sat[i].last() < 0 || unsigned(sat[i].last()) < k);
assert(sat[j].last() < 0 || unsigned(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]'.
set_union(sat[i], sat[j], new_satrow);
// Computing 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;
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;
}
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(Row::Type(dest.topology(),
Row::RAY_OR_POINT_OR_INEQUALITY));
sat.add_row(new_satrow);
}
else
sat[dest_num_rows] = new_satrow;
Row& new_row = dest[dest_num_rows];
// The following fragment optimizes the computation of
//
// Integer scale = scalar_prod[i];
// scale.gcd_assign(scalar_prod[j]);
// Integer scaled_sp_i = scalar_prod[i] / scale;
// Integer scaled_sp_j = scalar_prod[j] / scale;
// for (dimension_type c = dest_num_columns; c-- > 0; ) {
// new_row[c] = scaled_sp_i * dest[j][c];
// new_row[c] -= scaled_sp_j * dest[i][c];
// }
gcd_assign(tmp_Integer[1],
scalar_prod[i],
scalar_prod[j]);
exact_div_assign(tmp_Integer[2],
scalar_prod[i],
tmp_Integer[1]);
exact_div_assign(tmp_Integer[3],
scalar_prod[j],
tmp_Integer[1]);
for (dimension_type c = dest_num_columns; c-- > 0; ) {
tmp_Integer[4] = tmp_Integer[2] * dest[j][c];
tmp_Integer[5] = tmp_Integer[3] * dest[i][c];
new_row[c] = tmp_Integer[4] - tmp_Integer[5];
}
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.
assert(scalar_prod.size() >= dest_num_rows);
if (scalar_prod.size() <= dest_num_rows)
scalar_prod.push_back(Integer_zero());
else
scalar_prod[dest_num_rows] = Integer_zero();
// Increment the number of generators.
++dest_num_rows;
}
}
}
#if REACTIVE_ABANDONING
maybe_abandon();
#endif
}
// 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;
std::swap(dest[i], dest[j]);
std::swap(scalar_prod[i], scalar_prod[j]);
std::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;
}
}
#if !REACTIVE_ABANDONING
// Check if the client has requested abandoning all exponential
// computations. If so, the exception specified by the client
// is thrown now.
maybe_abandon();
#endif
}
// We may have identified some redundant constraints in `source',
// which have been swapped at the end of the matrix.
if (source_num_redundant > 0) {
assert(source_num_redundant == source.num_rows() - source_num_rows);
source.erase_to_end(source_num_rows);
sat.columns_erase_to_end(source_num_rows);
}
// If `start == 0', then `source' was sorted and remained so.
// If otherwise `start > 0', then the two sub-matrix made by the
// non-pending rows and the pending rows, respectively, were both sorted.
// Thus, the overall matrix is sorted if and only if either
// `start == source_num_rows' (i.e., the second sub-matrix is empty)
// or the row ordering holds for the two rows at the boundary between
// the two sub-matrices.
if (start > 0 && start < source_num_rows)
source.set_sorted(source[start - 1] <= source[start]);
// 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 matrix.
if (dest_num_rows < dest.num_rows()) {
dest.erase_to_end(dest_num_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.rows_erase_to_end(dest_num_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 (dest[i - 1] > dest[i]) {
dest.set_sorted(false);
break;
}
// There are no pending generators in `dest'.
dest.unset_pending_rows();
return num_lines_or_equalities;
}
syntax highlighted by Code2HTML, v. 0.9.1