Polyhedron_conversion_templates.hh

Go to the documentation of this file.

/* Polyhedron class implementation: conversion().
   Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
   Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)

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 3 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., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.

For the most up-to-date information see the Parma Polyhedra Library
site: http://bugseng.com/products/ppl/ . */

#ifndef PPL_Polyhedron_conversion_templates_hh
#define PPL_Polyhedron_conversion_templates_hh 1

#include "Bit_Row_defs.hh"
#include "Bit_Matrix_defs.hh"
#include "Polyhedron_defs.hh"
#include "Scalar_Products_defs.hh"
#include "Scalar_Products_inlines.hh"
#include "Temp_defs.hh"
#include "math_utilities_defs.hh"

#include <cstddef>
#include <climits>

// This flag turns on the quick non-adjacency test;
// the performance impact of this test was evaluated by H. Le Verge (1994);
// see also Corollary 4.2 in B. Genov's PhD thesis (2014).
#ifndef PPL_QUICK_NON_ADJ_TEST
#define PPL_QUICK_NON_ADJ_TEST 1
#endif

// This flag turns on the quick adjacency test;
// for a justification, see Corollary 4.3 in B. Genov's PhD thesis (2014),
// where it is also said that the test was implemented in cddlib.
#ifndef PPL_QUICK_ADJ_TEST
#define PPL_QUICK_ADJ_TEST 1
#endif

namespace Parma_Polyhedra_Library {

template <typename Source_Linear_System, typename Dest_Linear_System>
dimension_type
Polyhedron::conversion(Source_Linear_System& source,
                       const dimension_type start,
                       Dest_Linear_System& dest,
                       Bit_Matrix& sat,
                       dimension_type num_lines_or_equalities) {
  typedef typename Dest_Linear_System::row_type dest_row_type;
  typedef typename Source_Linear_System::row_type source_row_type;

  // Constraints and generators must have the same dimension,
  // otherwise the scalar products below will bomb.
  PPL_ASSERT(source.space_dimension() == dest.space_dimension());
  const dimension_type source_space_dim = source.space_dimension();
  const dimension_type source_num_rows = source.num_rows();
  const dimension_type source_num_columns = source_space_dim
    + (source.is_necessarily_closed() ? 1U : 2U);


  dimension_type dest_num_rows = dest.num_rows();
  // The rows removed from `dest' will be placed in this vector, so they
  // can be recycled if needed.
  std::vector<dest_row_type> recyclable_dest_rows;

  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());

  PPL_DIRTY_TEMP_COEFFICIENT(normalized_sp_i);
  PPL_DIRTY_TEMP_COEFFICIENT(normalized_sp_o);

  bool dest_sorted = dest.is_sorted();
  const dimension_type dest_first_pending_row = dest.first_pending_row();

  // This will contain the row indexes of the redundant rows of `source'.
  std::vector<dimension_type> redundant_source_rows;

#if PPL_QUICK_ADJ_TEST
  // This will contain the number of ones in each row of `sat'.
  PPL_DIRTY_TEMP(std::vector<dimension_type>, sat_num_ones);
  sat_num_ones.resize(dest_num_rows, 0);
  for (dimension_type i = dest_num_rows; i-- > 0; ) {
    sat_num_ones[i] = sat[i].count_ones();
  }
#endif // PPL_QUICK_ADJ_TEST

  // 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; ++k) {
    const source_row_type& source_k = source[k];

#ifndef NDEBUG
#if PPL_QUICK_ADJ_TEST
    for (dimension_type i = dest_num_rows; i-- > 0; ) {
      PPL_ASSERT(sat_num_ones[i] == sat[i].count_ones());
    }
#endif // PPL_QUICK_ADJ_TEST
#endif // NDEBUG

    // `scalar_prod[i]' will contain the scalar product of the
    // constraint `source_k' and the generator `dest_rows[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.insert(scalar_prod.end(),
                         dest_num_rows - scalar_prod.size(),
                         Coefficient_zero());
    }
    // `index_non_zero' will indicate the first generator in `dest_rows'
    // 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.sys.rows[index_non_zero]);
      WEIGHT_ADD_MUL(17, source_space_dim);
      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.sys.rows[i]);
      WEIGHT_ADD_MUL(25, source_space_dim);
      // 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_rows' (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_rows[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.sys.rows[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_rows[index_non_zero]' lies on the wrong half-space:
        // we change it to have the opposite direction.
        neg_assign(scalar_prod[index_non_zero]);
        neg_assign(dest.sys.rows[index_non_zero].expr);
        // The modified row may still not be OK(), so don't assert OK here.
        // They are all checked at the end of this function.
      }
      // Having changed a line to a ray, we set `dest_rows' to be a
      // non-sorted system, we decrement the number of lines of `dest_rows'
      // and, if necessary, we move the new ray below all the remaining lines.
      dest_sorted = false;
      --num_lines_or_equalities;
      if (index_non_zero != num_lines_or_equalities) {
        swap(dest.sys.rows[index_non_zero],
             dest.sys.rows[num_lines_or_equalities]);
        swap(scalar_prod[index_non_zero],
             scalar_prod[num_lines_or_equalities]);
      }
      const dest_row_type& dest_nle = dest.sys.rows[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];
      PPL_ASSERT(scalar_prod_nle != 0);
      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);
          dest_row_type& dest_i = dest.sys.rows[i];
          neg_assign(normalized_sp_i);
          dest_i.expr.linear_combine(dest_nle.expr,
                                     normalized_sp_o, normalized_sp_i);
          dest_i.strong_normalize();
          // The modified row may still not be OK(), so don't assert OK here.
          // They are all checked at the end of this function.
          scalar_prod[i] = 0;
          // dest_sorted has already been set to false.
        }
      }

      // 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);
          dest_row_type& dest_i = dest.sys.rows[i];
          WEIGHT_BEGIN();
          neg_assign(normalized_sp_i);
          dest_i.expr.linear_combine(dest_nle.expr,
                                     normalized_sp_o, normalized_sp_i);
          dest_i.strong_normalize();
          // The modified row may still not be OK(), so don't assert OK here.
          // They are all checked at the end of this function.
          scalar_prod[i] = 0;
          // `dest_sorted' has already been set to false.
          WEIGHT_ADD_MUL(41, source_space_dim);
        }
        // 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 - redundant_source_rows.size());
#if PPL_QUICK_ADJ_TEST
        ++sat_num_ones[num_lines_or_equalities];
#endif // PPL_QUICK_ADJ_TEST
      }
      else {
        // ... otherwise, the constraint is an equality which is
        // violated by the generator `dest_nle': the generator has to be
        // removed from `dest_rows'.
        --dest_num_rows;
        swap(dest.sys.rows[num_lines_or_equalities],
             dest.sys.rows[dest_num_rows]);
        recyclable_dest_rows.resize(recyclable_dest_rows.size() + 1);
        swap(dest.sys.rows.back(), recyclable_dest_rows.back());
        dest.sys.rows.pop_back();
        PPL_ASSERT(dest_num_rows == dest.sys.rows.size());

        swap(scalar_prod_nle, scalar_prod[dest_num_rows]);
        swap(sat_nle, sat[dest_num_rows]);
#if PPL_QUICK_ADJ_TEST
        swap(sat_num_ones[num_lines_or_equalities],
             sat_num_ones[dest_num_rows]);
#endif // PPL_QUICK_ADJ_TEST
        // dest_sorted has already been set to false.
      }
      // Finished handling the line or equality case:
      // continue with next `k'.
      continue;
    }

    // Here all the lines in `dest_rows' saturate the constraint `source_k'.
    PPL_ASSERT(index_non_zero >= num_lines_or_equalities);
    // First, we reorder the generators in `dest_rows' 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.sys.rows[sup_bound], dest.sys.rows[lines_or_equal_bound]);
        swap(scalar_prod[sup_bound], scalar_prod[lines_or_equal_bound]);
        swap(sat[sup_bound], sat[lines_or_equal_bound]);
#if PPL_QUICK_ADJ_TEST
        swap(sat_num_ones[sup_bound], sat_num_ones[lines_or_equal_bound]);
#endif // PPL_QUICK_ADJ_TEST
        ++lines_or_equal_bound;
        ++sup_bound;
        dest_sorted = false;
      }
      else if (sp_sign < 0) {
        // This generator has to be moved in Q-.
        --inf_bound;
        swap(dest.sys.rows[sup_bound], dest.sys.rows[inf_bound]);
        swap(sat[sup_bound], sat[inf_bound]);
        swap(scalar_prod[sup_bound], scalar_prod[inf_bound]);
#if PPL_QUICK_ADJ_TEST
        swap(sat_num_ones[sup_bound], sat_num_ones[inf_bound]);
#endif // PPL_QUICK_ADJ_TEST
        dest_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()) {
        redundant_source_rows.push_back(k);
      }
      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_rows' all the generators of Q+.
        PPL_ASSERT(dest_num_rows >= lines_or_equal_bound);
        while (dest_num_rows != lines_or_equal_bound) {
          recyclable_dest_rows.resize(recyclable_dest_rows.size() + 1);
          swap(dest.sys.rows.back(), recyclable_dest_rows.back());
          dest.sys.rows.pop_back();
          --dest_num_rows;
          }
        PPL_ASSERT(dest_num_rows == dest.sys.rows.size());
      }
      // Finished handling the case when Q- is empty:
      // continue with next `k'.
      continue;
    }

    // The set Q- is not empty, i.e., at least one generator
    // violates the constraint `source_k'.
    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_rows' all the generators in Q-.
      PPL_ASSERT(dest_num_rows >= sup_bound);
      while (dest_num_rows != sup_bound) {
        recyclable_dest_rows.resize(recyclable_dest_rows.size() + 1);
        swap(dest.sys.rows.back(), recyclable_dest_rows.back());
        dest.sys.rows.pop_back();
        --dest_num_rows;
      }
      PPL_ASSERT(dest_num_rows == dest.sys.rows.size());
      // Finished handling the case when Q+ is empty:
      // continue with next `k'.
      continue;
    }

    // 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-.

    const dimension_type bound = dest_num_rows;

#if PPL_QUICK_NON_ADJ_TEST
    // For the quick non-adjacency test, 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.
    const dimension_type min_saturators
      = source_num_columns - num_lines_or_equalities - 2;
    // NOTE: we are treating the `k'-th constraint.
    const dimension_type max_saturators = k - redundant_source_rows.size();
#endif // PPL_QUICK_NON_ADJ_TEST

    // 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_rows[i]' and `dest_rows[j]' are
        // adjacent.
        // If there exist another generator that saturates
        // all the constraints saturated by both `dest_rows[i]' and
        // `dest_rows[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_rows[i]' and
        // `dest_rows[j]'.
        Bit_Row new_satrow(sat[i], sat[j]);

        // Even before actually creating the new ray as a
        // positive combination of `dest_rows[i]' and `dest_rows[j]',
        // we exploit saturation information to perform:
        //  - a quick non-adjacency test;
        //  - a quick adjacency test.

#if (PPL_QUICK_NON_ADJ_TEST || PPL_QUICK_ADJ_TEST)
        // Compute the number of common saturators.
        dimension_type new_satrow_ones = new_satrow.count_ones();
#endif // (PPL_QUICK_NON_ADJ_TEST || PPL_QUICK_ADJ_TEST)

#if PPL_QUICK_NON_ADJ_TEST
        const dimension_type num_common_satur
          = max_saturators - new_satrow_ones;
        if (num_common_satur < min_saturators) {
          // Quick non-adjacency test succeded: consider next `j'.
          continue;
        }
#endif // PPL_QUICK_NON_ADJ_TEST

#if PPL_QUICK_ADJ_TEST
        // If either `sat[i]' or `sat[j]' has exactly one more zeroes
        // than `new_satrow', then `dest_rows[i]' and `dest_rows[j]'
        // are adjacent. Equivalently, adjacency holds if `new_satrow_ones'
        // is equal to 1 plus the maximum of `sat_num_ones[i]' and
        // `sat_num_ones[j]'.
        const dimension_type max_ones_i_j
          = std::max(sat_num_ones[i], sat_num_ones[j]);
        if (max_ones_i_j + 1 == new_satrow_ones) {
          // Quick adjacency test succeded: skip the full test.
          goto are_adjacent;
        }
#endif // PPL_QUICK_ADJ_TEST

        // Perform the full (combinatorial) adjacency test.
        {
          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_rows[i]' and `dest_rows[j]'.
              redundant = true;
              break;
            }
          }
          PPL_ASSERT(bound >= num_lines_or_equalities);
          WEIGHT_ADD_MUL(15, bound - num_lines_or_equalities);
          if (redundant) {
            // Full non-adjacency test succeded: consider next `j'.
            continue;
          }
        }

#if PPL_QUICK_ADJ_TEST
      are_adjacent:
#endif // PPL_QUICK_ADJ_TEST
        // Adding the new ray to `dest_rows' and the corresponding
        // saturation row to `sat'.
        dest_row_type new_row;
        if (recyclable_dest_rows.empty()) {
          sat.add_recycled_row(new_satrow);
#if PPL_QUICK_ADJ_TEST
          sat_num_ones.push_back(new_satrow_ones);
#endif // PPL_QUICK_ADJ_TEST
        }
        else {
          swap(new_row, recyclable_dest_rows.back());
          recyclable_dest_rows.pop_back();
          new_row.set_space_dimension_no_ok(source_space_dim);
          swap(sat[dest_num_rows], new_satrow);
#if PPL_QUICK_ADJ_TEST
          swap(sat_num_ones[dest_num_rows], new_satrow_ones);
#endif // PPL_QUICK_ADJ_TEST
        }

        // 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();

        neg_assign(normalized_sp_o);
        new_row = dest.sys.rows[j];
        // TODO: Check if the following assertions hold.
        PPL_ASSERT(normalized_sp_i != 0);
        PPL_ASSERT(normalized_sp_o != 0);
        new_row.expr.linear_combine(dest.sys.rows[i].expr,
                                    normalized_sp_i, normalized_sp_o);

        WEIGHT_ADD_MUL(86, source_space_dim);
        new_row.strong_normalize();
        // Don't assert new_row.OK() here, because it may fail if
        // the parameter `dest' contained a row that wasn't ok.
        // Since we added a new generator to `dest_rows',
        // 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();
        }
        dest.sys.rows.resize(dest.sys.rows.size() + 1);
        swap(dest.sys.rows.back(), new_row);
        // Increment the number of generators.
        ++dest_num_rows;
      } // End of loop on `j'.
      // Check if the client has requested abandoning all expensive
      // computations.  If so, the exception specified by the client
      // is thrown now.
      maybe_abandon();
    } // End of loop on `i'.
    // 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.

      // After the removal of redundant rows in `source', the k-th
      // row will have index `new_k'.
      const dimension_type new_k = k - redundant_source_rows.size();
      for (dimension_type l = lines_or_equal_bound;
           l < sup_bound; ++l) {
        sat[l].set(new_k);
#if PPL_QUICK_ADJ_TEST
        ++sat_num_ones[l];
#endif // PPL_PPL_QUICK_ADJ_TEST
      }
    }
    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.sys.rows[i], dest.sys.rows[j]);
      swap(scalar_prod[i], scalar_prod[j]);
      swap(sat[i], sat[j]);
#if PPL_QUICK_ADJ_TEST
      swap(sat_num_ones[i], sat_num_ones[j]);
#endif // PPL_QUICK_ADJ_TEST
      ++j;
      dest_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.
    const dimension_type new_num_rows = (j == bound) ? i : j;
    PPL_ASSERT(dest_num_rows >= new_num_rows);
    while (dest_num_rows != new_num_rows) {
      recyclable_dest_rows.resize(recyclable_dest_rows.size() + 1);
      swap(dest.sys.rows.back(), recyclable_dest_rows.back());
      dest.sys.rows.pop_back();
      --dest_num_rows;
    }
    PPL_ASSERT(dest_num_rows == dest.sys.rows.size());
  } // End of loop on `k'.

  // We may have identified some redundant constraints in `source',
  // which have been swapped at the end of the system.
  if (redundant_source_rows.size() > 0) {
    source.remove_rows(redundant_source_rows);
    sat.remove_trailing_columns(redundant_source_rows.size());
  }

  // 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_rows',
  // which have been swapped into recyclable_dest_rows.
  if (!recyclable_dest_rows.empty()) {
    const dimension_type num_removed_rows = recyclable_dest_rows.size();
    sat.remove_trailing_rows(num_removed_rows);
  }
  if (dest_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.sys.rows[i - 1], dest.sys.rows[i]) > 0) {
        dest_sorted = false;
        break;
      }
    }
  }
#ifndef NDEBUG
  // The previous code can modify the rows' fields, exploiting the friendness.
  // Check that all rows are OK now.
  for (dimension_type i = dest.num_rows(); i-- > 0; ) {
    PPL_ASSERT(dest.sys.rows[i].OK());
  }
#endif

  dest.sys.index_first_pending = dest.num_rows();
  dest.set_sorted(dest_sorted);
  PPL_ASSERT(dest.sys.OK());

  return num_lines_or_equalities;
}

} // namespace Parma_Polyhedra_Library

#endif // !defined(PPL_Polyhedron_conversion_templates_hh)