source: lemon-project-template-glpk/deps/glpk/src/glpnpp06.c @ 10:5545663ca997

/* glpnpp06.c (translate feasibility problem to CNF-SAT) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
*  E-mail: <mao@gnu.org>.
*
*  GLPK 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.
*
*  GLPK 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 GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glpnpp.h"

/***********************************************************************
*  npp_sat_free_row - process free (unbounded) row
*
*  This routine processes row p, which is free (i.e. has no finite
*  bounds):
*
*     -inf < sum a[p,j] x[j] < +inf.                                 (1)
*
*  The constraint (1) cannot be active and therefore it is redundant,
*  so the routine simply removes it from the original problem. */

void npp_sat_free_row(NPP *npp, NPPROW *p)
{     /* the row should be free */
      xassert(p->lb == -DBL_MAX && p->ub == +DBL_MAX);
      /* remove the row from the problem */
      npp_del_row(npp, p);
      return;
}

/***********************************************************************
*  npp_sat_fixed_col - process fixed column
*
*  This routine processes column q, which is fixed:
*
*     x[q] = s[q],                                                   (1)
*
*  where s[q] is a fixed column value.
*
*  The routine substitutes fixed value s[q] into constraint rows and
*  then removes column x[q] from the original problem.
*
*  Substitution of x[q] = s[q] into row i gives:
*
*     L[i] <= sum a[i,j] x[j] <= U[i]   ==>
*              j
*
*     L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i]   ==>
*            j!=q
*
*     L[i] <= sum a[i,j] x[j] + a[i,q] s[q] <= U[i]   ==>
*            j!=q
*
*     L~[i] <= sum a[i,j] x[j] <= U~[i],
*             j!=q
*
*  where
*
*     L~[i] = L[i] - a[i,q] s[q],                                    (2)
*
*     U~[i] = U[i] - a[i,q] s[q]                                     (3)
*
*  are, respectively, lower and upper bound of row i in the transformed
*  problem.
*
*  On recovering solution x[q] is assigned the value of s[q]. */

struct sat_fixed_col
{     /* fixed column */
      int q;
      /* column reference number for variable x[q] */
      int s;
      /* value, at which x[q] is fixed */
};

static int rcv_sat_fixed_col(NPP *, void *);

int npp_sat_fixed_col(NPP *npp, NPPCOL *q)
{     struct sat_fixed_col *info;
      NPPROW *i;
      NPPAIJ *aij;
      int temp;
      /* the column should be fixed */
      xassert(q->lb == q->ub);
      /* create transformation stack entry */
      info = npp_push_tse(npp,
         rcv_sat_fixed_col, sizeof(struct sat_fixed_col));
      info->q = q->j;
      info->s = (int)q->lb;
      xassert((double)info->s == q->lb);
      /* substitute x[q] = s[q] into constraint rows */
      if (info->s == 0)
         goto skip;
      for (aij = q->ptr; aij != NULL; aij = aij->c_next)
      {  i = aij->row;
         if (i->lb != -DBL_MAX)
         {  i->lb -= aij->val * (double)info->s;
            temp = (int)i->lb;
            if ((double)temp != i->lb)
               return 1; /* integer arithmetic error */
         }
         if (i->ub != +DBL_MAX)
         {  i->ub -= aij->val * (double)info->s;
            temp = (int)i->ub;
            if ((double)temp != i->ub)
               return 2; /* integer arithmetic error */
         }
      }
skip: /* remove the column from the problem */
      npp_del_col(npp, q);
      return 0;
}

static int rcv_sat_fixed_col(NPP *npp, void *info_)
{     struct sat_fixed_col *info = info_;
      npp->c_value[info->q] = (double)info->s;
      return 0;
}

/***********************************************************************
*  npp_sat_is_bin_comb - test if row is binary combination
*
*  This routine tests if the specified row is a binary combination,
*  i.e. all its constraint coefficients are +1 and -1 and all variables
*  are binary. If the test was passed, the routine returns non-zero,
*  otherwise zero. */

int npp_sat_is_bin_comb(NPP *npp, NPPROW *row)
{     NPPCOL *col;
      NPPAIJ *aij;
      xassert(npp == npp);
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  if (!(aij->val == +1.0 || aij->val == -1.0))
            return 0; /* non-unity coefficient */
         col = aij->col;
         if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
            return 0; /* non-binary column */
      }
      return 1; /* test was passed */
}

/***********************************************************************
*  npp_sat_num_pos_coef - determine number of positive coefficients
*
*  This routine returns the number of positive coefficients in the
*  specified row. */

int npp_sat_num_pos_coef(NPP *npp, NPPROW *row)
{     NPPAIJ *aij;
      int num = 0;
      xassert(npp == npp);
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  if (aij->val > 0.0)
            num++;
      }
      return num;
}

/***********************************************************************
*  npp_sat_num_neg_coef - determine number of negative coefficients
*
*  This routine returns the number of negative coefficients in the
*  specified row. */

int npp_sat_num_neg_coef(NPP *npp, NPPROW *row)
{     NPPAIJ *aij;
      int num = 0;
      xassert(npp == npp);
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  if (aij->val < 0.0)
            num++;
      }
      return num;
}

/***********************************************************************
*  npp_sat_is_cover_ineq - test if row is covering inequality
*
*  The canonical form of a covering inequality is the following:
*
*     sum x[j] >= 1,                                                 (1)
*   j in J
*
*  where all x[j] are binary variables.
*
*  In general case a covering inequality may have one of the following
*  two forms:
*
*     sum  x[j] -  sum  x[j] >= 1 - |J-|,                            (2)
*   j in J+      j in J-
*
*
*     sum  x[j] -  sum  x[j] <= |J+| - 1.                            (3)
*   j in J+      j in J-
*
*  Obviously, the inequality (2) can be transformed to the form (1) by
*  substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the
*  negation of variable x[j]. And the inequality (3) can be transformed
*  to (2) by multiplying both left- and right-hand sides by -1.
*
*  This routine returns one of the following codes:
*
*  0, if the specified row is not a covering inequality;
*
*  1, if the specified row has the form (2);
*
*  2, if the specified row has the form (3). */

int npp_sat_is_cover_ineq(NPP *npp, NPPROW *row)
{     xassert(npp == npp);
      if (row->lb != -DBL_MAX && row->ub == +DBL_MAX)
      {  /* row is inequality of '>=' type */
         if (npp_sat_is_bin_comb(npp, row))
         {  /* row is a binary combination */
            if (row->lb == 1.0 - npp_sat_num_neg_coef(npp, row))
            {  /* row has the form (2) */
               return 1;
            }
         }
      }
      else if (row->lb == -DBL_MAX && row->ub != +DBL_MAX)
      {  /* row is inequality of '<=' type */
         if (npp_sat_is_bin_comb(npp, row))
         {  /* row is a binary combination */
            if (row->ub == npp_sat_num_pos_coef(npp, row) - 1.0)
            {  /* row has the form (3) */
               return 2;
            }
         }
      }
      /* row is not a covering inequality */
      return 0;
}

/***********************************************************************
*  npp_sat_is_pack_ineq - test if row is packing inequality
*
*  The canonical form of a packing inequality is the following:
*
*     sum x[j] <= 1,                                                 (1)
*   j in J
*
*  where all x[j] are binary variables.
*
*  In general case a packing inequality may have one of the following
*  two forms:
*
*     sum  x[j] -  sum  x[j] <= 1 - |J-|,                            (2)
*   j in J+      j in J-
*
*
*     sum  x[j] -  sum  x[j] >= |J+| - 1.                            (3)
*   j in J+      j in J-
*
*  Obviously, the inequality (2) can be transformed to the form (1) by
*  substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the
*  negation of variable x[j]. And the inequality (3) can be transformed
*  to (2) by multiplying both left- and right-hand sides by -1.
*
*  This routine returns one of the following codes:
*
*  0, if the specified row is not a packing inequality;
*
*  1, if the specified row has the form (2);
*
*  2, if the specified row has the form (3). */

int npp_sat_is_pack_ineq(NPP *npp, NPPROW *row)
{     xassert(npp == npp);
      if (row->lb == -DBL_MAX && row->ub != +DBL_MAX)
      {  /* row is inequality of '<=' type */
         if (npp_sat_is_bin_comb(npp, row))
         {  /* row is a binary combination */
            if (row->ub == 1.0 - npp_sat_num_neg_coef(npp, row))
            {  /* row has the form (2) */
               return 1;
            }
         }
      }
      else if (row->lb != -DBL_MAX && row->ub == +DBL_MAX)
      {  /* row is inequality of '>=' type */
         if (npp_sat_is_bin_comb(npp, row))
         {  /* row is a binary combination */
            if (row->lb == npp_sat_num_pos_coef(npp, row) - 1.0)
            {  /* row has the form (3) */
               return 2;
            }
         }
      }
      /* row is not a packing inequality */
      return 0;
}

/***********************************************************************
*  npp_sat_is_partn_eq - test if row is partitioning equality
*
*  The canonical form of a partitioning equality is the following:
*
*     sum x[j] = 1,                                                  (1)
*   j in J
*
*  where all x[j] are binary variables.
*
*  In general case a partitioning equality may have one of the following
*  two forms:
*
*     sum  x[j] -  sum  x[j] = 1 - |J-|,                             (2)
*   j in J+      j in J-
*
*
*     sum  x[j] -  sum  x[j] = |J+| - 1.                             (3)
*   j in J+      j in J-
*
*  Obviously, the equality (2) can be transformed to the form (1) by
*  substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the
*  negation of variable x[j]. And the equality (3) can be transformed
*  to (2) by multiplying both left- and right-hand sides by -1.
*
*  This routine returns one of the following codes:
*
*  0, if the specified row is not a partitioning equality;
*
*  1, if the specified row has the form (2);
*
*  2, if the specified row has the form (3). */

int npp_sat_is_partn_eq(NPP *npp, NPPROW *row)
{     xassert(npp == npp);
      if (row->lb == row->ub)
      {  /* row is equality constraint */
         if (npp_sat_is_bin_comb(npp, row))
         {  /* row is a binary combination */
            if (row->lb == 1.0 - npp_sat_num_neg_coef(npp, row))
            {  /* row has the form (2) */
               return 1;
            }
            if (row->ub == npp_sat_num_pos_coef(npp, row) - 1.0)
            {  /* row has the form (3) */
               return 2;
            }
         }
      }
      /* row is not a partitioning equality */
      return 0;
}

/***********************************************************************
*  npp_sat_reverse_row - multiply both sides of row by -1
*
*  This routines multiplies by -1 both left- and right-hand sides of
*  the specified row:
*
*     L <= sum x[j] <= U,
*
*  that results in the following row:
*
*     -U <= sum (-x[j]) <= -L.
*
*  If no integer overflow occured, the routine returns zero, otherwise
*  non-zero. */

int npp_sat_reverse_row(NPP *npp, NPPROW *row)
{     NPPAIJ *aij;
      int temp, ret = 0;
      double old_lb, old_ub;
      xassert(npp == npp);
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  aij->val = -aij->val;
         temp = (int)aij->val;
         if ((double)temp != aij->val)
            ret = 1;
      }
      old_lb = row->lb, old_ub = row->ub;
      if (old_ub == +DBL_MAX)
         row->lb = -DBL_MAX;
      else
      {  row->lb = -old_ub;
         temp = (int)row->lb;
         if ((double)temp != row->lb)
            ret = 2;
      }
      if (old_lb == -DBL_MAX)
         row->ub = +DBL_MAX;
      else
      {  row->ub = -old_lb;
         temp = (int)row->ub;
         if ((double)temp != row->ub)
            ret = 3;
      }
      return ret;
}

/***********************************************************************
*  npp_sat_split_pack - split packing inequality
*
*  Let there be given a packing inequality in canonical form:
*
*     sum  t[j] <= 1,                                                (1)
*   j in J
*
*  where t[j] = x[j] or t[j] = 1 - x[j], x[j] is a binary variable.
*  And let J = J1 U J2 is a partition of the set of literals. Then the
*  inequality (1) is obviously equivalent to the following two packing
*  inequalities:
*
*     sum  t[j] <= y       <-->   sum  t[j] + (1 - y) <= 1,          (2)
*   j in J1                     j in J1
*
*     sum  t[j] <= 1 - y   <-->   sum  t[j] + y       <= 1,          (3)
*   j in J2                     j in J2
*
*  where y is a new binary variable added to the transformed problem.
*
*  Assuming that the specified row is a packing inequality (1), this
*  routine constructs the set J1 by including there first nlit literals
*  (terms) from the specified row, and the set J2 = J \ J1. Then the
*  routine creates a new row, which corresponds to inequality (2), and
*  replaces the specified row with inequality (3). */

NPPROW *npp_sat_split_pack(NPP *npp, NPPROW *row, int nlit)
{     NPPROW *rrr;
      NPPCOL *col;
      NPPAIJ *aij;
      int k;
      /* original row should be packing inequality (1) */
      xassert(npp_sat_is_pack_ineq(npp, row) == 1);
      /* and nlit should be less than the number of literals (terms)
         in the original row */
      xassert(0 < nlit && nlit < npp_row_nnz(npp, row));
      /* create new row corresponding to inequality (2) */
      rrr = npp_add_row(npp);
      rrr->lb = -DBL_MAX, rrr->ub = 1.0;
      /* move first nlit literals (terms) from the original row to the
         new row; the original row becomes inequality (3) */
      for (k = 1; k <= nlit; k++)
      {  aij = row->ptr;
         xassert(aij != NULL);
         /* add literal to the new row */
         npp_add_aij(npp, rrr, aij->col, aij->val);
         /* correct rhs */
         if (aij->val < 0.0)
            rrr->ub -= 1.0, row->ub += 1.0;
         /* remove literal from the original row */
         npp_del_aij(npp, aij);
      }
      /* create new binary variable y */
      col = npp_add_col(npp);
      col->is_int = 1, col->lb = 0.0, col->ub = 1.0;
      /* include literal (1 - y) in the new row */
      npp_add_aij(npp, rrr, col, -1.0);
      rrr->ub -= 1.0;
      /* include literal y in the original row */
      npp_add_aij(npp, row, col, +1.0);
      return rrr;
}

/***********************************************************************
*  npp_sat_encode_pack - encode packing inequality
*
*  Given a packing inequality in canonical form:
*
*     sum  t[j] <= 1,                                                (1)
*   j in J
*
*  where t[j] = x[j] or t[j] = 1 - x[j], x[j] is a binary variable,
*  this routine translates it to CNF by replacing it with the following
*  equivalent set of edge packing inequalities:
*
*     t[j] + t[k] <= 1   for all j, k in J, j != k.                  (2)
*
*  Then the routine transforms each edge packing inequality (2) to
*  corresponding covering inequality (that encodes two-literal clause)
*  by multiplying both its part by -1:
*
*     - t[j] - t[k] >= -1   <-->   (1 - t[j]) + (1 - t[k]) >= 1.     (3)
*
*  On exit the routine removes the original row from the problem. */

void npp_sat_encode_pack(NPP *npp, NPPROW *row)
{     NPPROW *rrr;
      NPPAIJ *aij, *aik;
      /* original row should be packing inequality (1) */
      xassert(npp_sat_is_pack_ineq(npp, row) == 1);
      /* create equivalent system of covering inequalities (3) */
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  /* due to symmetry only one of inequalities t[j] + t[k] <= 1
            and t[k] <= t[j] <= 1 can be considered */
         for (aik = aij->r_next; aik != NULL; aik = aik->r_next)
         {  /* create edge packing inequality (2) */
            rrr = npp_add_row(npp);
            rrr->lb = -DBL_MAX, rrr->ub = 1.0;
            npp_add_aij(npp, rrr, aij->col, aij->val);
            if (aij->val < 0.0)
               rrr->ub -= 1.0;
            npp_add_aij(npp, rrr, aik->col, aik->val);
            if (aik->val < 0.0)
               rrr->ub -= 1.0;
            /* and transform it to covering inequality (3) */
            npp_sat_reverse_row(npp, rrr);
            xassert(npp_sat_is_cover_ineq(npp, rrr) == 1);
         }
      }
      /* remove the original row from the problem */
      npp_del_row(npp, row);
      return;
}

/***********************************************************************
*  npp_sat_encode_sum2 - encode 2-bit summation
*
*  Given a set containing two literals x and y this routine encodes
*  the equality
*
*     x + y = s + 2 * c,                                             (1)
*
*  where
*
*     s = (x + y) % 2                                                (2)
*
*  is a binary variable modeling the low sum bit, and
*
*     c = (x + y) / 2                                                (3)
*
*  is a binary variable modeling the high (carry) sum bit. */

void npp_sat_encode_sum2(NPP *npp, NPPLSE *set, NPPSED *sed)
{     NPPROW *row;
      int x, y, s, c;
      /* the set should contain exactly two literals */
      xassert(set != NULL);
      xassert(set->next != NULL);
      xassert(set->next->next == NULL);
      sed->x = set->lit;
      xassert(sed->x.neg == 0 || sed->x.neg == 1);
      sed->y = set->next->lit;
      xassert(sed->y.neg == 0 || sed->y.neg == 1);
      sed->z.col = NULL, sed->z.neg = 0;
      /* perform encoding s = (x + y) % 2 */
      sed->s = npp_add_col(npp);
      sed->s->is_int = 1, sed->s->lb = 0.0, sed->s->ub = 1.0;
      for (x = 0; x <= 1; x++)
      {  for (y = 0; y <= 1; y++)
         {  for (s = 0; s <= 1; s++)
            {  if ((x + y) % 2 != s)
               {  /* generate CNF clause to disable infeasible
                     combination */
                  row = npp_add_row(npp);
                  row->lb = 1.0, row->ub = +DBL_MAX;
                  if (x == sed->x.neg)
                     npp_add_aij(npp, row, sed->x.col, +1.0);
                  else
                  {  npp_add_aij(npp, row, sed->x.col, -1.0);
                     row->lb -= 1.0;
                  }
                  if (y == sed->y.neg)
                     npp_add_aij(npp, row, sed->y.col, +1.0);
                  else
                  {  npp_add_aij(npp, row, sed->y.col, -1.0);
                     row->lb -= 1.0;
                  }
                  if (s == 0)
                     npp_add_aij(npp, row, sed->s, +1.0);
                  else
                  {  npp_add_aij(npp, row, sed->s, -1.0);
                     row->lb -= 1.0;
                  }
               }
            }
         }
      }
      /* perform encoding c = (x + y) / 2 */
      sed->c = npp_add_col(npp);
      sed->c->is_int = 1, sed->c->lb = 0.0, sed->c->ub = 1.0;
      for (x = 0; x <= 1; x++)
      {  for (y = 0; y <= 1; y++)
         {  for (c = 0; c <= 1; c++)
            {  if ((x + y) / 2 != c)
               {  /* generate CNF clause to disable infeasible
                     combination */
                  row = npp_add_row(npp);
                  row->lb = 1.0, row->ub = +DBL_MAX;
                  if (x == sed->x.neg)
                     npp_add_aij(npp, row, sed->x.col, +1.0);
                  else
                  {  npp_add_aij(npp, row, sed->x.col, -1.0);
                     row->lb -= 1.0;
                  }
                  if (y == sed->y.neg)
                     npp_add_aij(npp, row, sed->y.col, +1.0);
                  else
                  {  npp_add_aij(npp, row, sed->y.col, -1.0);
                     row->lb -= 1.0;
                  }
                  if (c == 0)
                     npp_add_aij(npp, row, sed->c, +1.0);
                  else
                  {  npp_add_aij(npp, row, sed->c, -1.0);
                     row->lb -= 1.0;
                  }
               }
            }
         }
      }
      return;
}

/***********************************************************************
*  npp_sat_encode_sum3 - encode 3-bit summation
*
*  Given a set containing at least three literals this routine chooses
*  some literals x, y, z from that set and encodes the equality
*
*     x + y + z = s + 2 * c,                                         (1)
*
*  where
*
*     s = (x + y + z) % 2                                            (2)
*
*  is a binary variable modeling the low sum bit, and
*
*     c = (x + y + z) / 2                                            (3)
*
*  is a binary variable modeling the high (carry) sum bit. */

void npp_sat_encode_sum3(NPP *npp, NPPLSE *set, NPPSED *sed)
{     NPPROW *row;
      int x, y, z, s, c;
      /* the set should contain at least three literals */
      xassert(set != NULL);
      xassert(set->next != NULL);
      xassert(set->next->next != NULL);
      sed->x = set->lit;
      xassert(sed->x.neg == 0 || sed->x.neg == 1);
      sed->y = set->next->lit;
      xassert(sed->y.neg == 0 || sed->y.neg == 1);
      sed->z = set->next->next->lit;
      xassert(sed->z.neg == 0 || sed->z.neg == 1);
      /* perform encoding s = (x + y + z) % 2 */
      sed->s = npp_add_col(npp);
      sed->s->is_int = 1, sed->s->lb = 0.0, sed->s->ub = 1.0;
      for (x = 0; x <= 1; x++)
      {  for (y = 0; y <= 1; y++)
         {  for (z = 0; z <= 1; z++)
            {  for (s = 0; s <= 1; s++)
               {  if ((x + y + z) % 2 != s)
                  {  /* generate CNF clause to disable infeasible
                        combination */
                     row = npp_add_row(npp);
                     row->lb = 1.0, row->ub = +DBL_MAX;
                     if (x == sed->x.neg)
                        npp_add_aij(npp, row, sed->x.col, +1.0);
                     else
                     {  npp_add_aij(npp, row, sed->x.col, -1.0);
                        row->lb -= 1.0;
                     }
                     if (y == sed->y.neg)
                        npp_add_aij(npp, row, sed->y.col, +1.0);
                     else
                     {  npp_add_aij(npp, row, sed->y.col, -1.0);
                        row->lb -= 1.0;
                     }
                     if (z == sed->z.neg)
                        npp_add_aij(npp, row, sed->z.col, +1.0);
                     else
                     {  npp_add_aij(npp, row, sed->z.col, -1.0);
                        row->lb -= 1.0;
                     }
                     if (s == 0)
                        npp_add_aij(npp, row, sed->s, +1.0);
                     else
                     {  npp_add_aij(npp, row, sed->s, -1.0);
                        row->lb -= 1.0;
                     }
                  }
               }
            }
         }
      }
      /* perform encoding c = (x + y + z) / 2 */
      sed->c = npp_add_col(npp);
      sed->c->is_int = 1, sed->c->lb = 0.0, sed->c->ub = 1.0;
      for (x = 0; x <= 1; x++)
      {  for (y = 0; y <= 1; y++)
         {  for (z = 0; z <= 1; z++)
            {  for (c = 0; c <= 1; c++)
               {  if ((x + y + z) / 2 != c)
                  {  /* generate CNF clause to disable infeasible
                        combination */
                     row = npp_add_row(npp);
                     row->lb = 1.0, row->ub = +DBL_MAX;
                     if (x == sed->x.neg)
                        npp_add_aij(npp, row, sed->x.col, +1.0);
                     else
                     {  npp_add_aij(npp, row, sed->x.col, -1.0);
                        row->lb -= 1.0;
                     }
                     if (y == sed->y.neg)
                        npp_add_aij(npp, row, sed->y.col, +1.0);
                     else
                     {  npp_add_aij(npp, row, sed->y.col, -1.0);
                        row->lb -= 1.0;
                     }
                     if (z == sed->z.neg)
                        npp_add_aij(npp, row, sed->z.col, +1.0);
                     else
                     {  npp_add_aij(npp, row, sed->z.col, -1.0);
                        row->lb -= 1.0;
                     }
                     if (c == 0)
                        npp_add_aij(npp, row, sed->c, +1.0);
                     else
                     {  npp_add_aij(npp, row, sed->c, -1.0);
                        row->lb -= 1.0;
                     }
                  }
               }
            }
         }
      }
      return;
}

/***********************************************************************
*  npp_sat_encode_sum_ax - encode linear combination of 0-1 variables
*
*  PURPOSE
*
*  Given a linear combination of binary variables:
*
*     sum a[j] x[j],                                                 (1)
*      j
*
*  which is the linear form of the specified row, this routine encodes
*  (i.e. translates to CNF) the following equality:
*
*                        n
*     sum |a[j]| t[j] = sum 2**(k-1) * y[k],                         (2)
*      j                k=1
*
*  where t[j] = x[j] (if a[j] > 0) or t[j] = 1 - x[j] (if a[j] < 0),
*  and y[k] is either t[j] or a new literal created by the routine or
*  a constant zero. Note that the sum in the right-hand side of (2) can
*  be thought as a n-bit representation of the sum in the left-hand
*  side, which is a non-negative integer number.
*
*  ALGORITHM
*
*  First, the number of bits, n, sufficient to represent any value in
*  the left-hand side of (2) is determined. Obviously, n is the number
*  of bits sufficient to represent the sum (sum |a[j]|).
*
*  Let
*
*               n
*     |a[j]| = sum 2**(k-1) b[j,k],                                  (3)
*              k=1
*
*  where b[j,k] is k-th bit in a n-bit representation of |a[j]|. Then
*
*                          m            n
*     sum |a[j]| * t[j] = sum 2**(k-1) sum b[j,k] * t[j].            (4)
*      j                  k=1          j=1
*
*  Introducing the set
*
*     J[k] = { j : b[j,k] = 1 }                                      (5)
*
*  allows rewriting (4) as follows:
*
*                          n
*     sum |a[j]| * t[j] = sum 2**(k-1)  sum    t[j].                 (6)
*      j                  k=1         j in J[k]
*
*  Thus, our goal is to provide |J[k]| <= 1 for all k, in which case
*  we will have the representation (1).
*
*  Let |J[k]| = 2, i.e. J[k] has exactly two literals u and v. In this
*  case we can apply the following transformation:
*
*     u + v = s + 2 * c,                                             (7)
*
*  where s and c are, respectively, low (sum) and high (carry) bits of
*  the sum of two bits. This allows to replace two literals u and v in
*  J[k] by one literal s, and carry out literal c to J[k+1].
*
*  If |J[k]| >= 3, i.e. J[k] has at least three literals u, v, and w,
*  we can apply the following transformation:
*
*     u + v + w = s + 2 * c.                                         (8)
*
*  Again, literal s replaces literals u, v, and w in J[k], and literal
*  c goes into J[k+1].
*
*  On exit the routine stores each literal from J[k] in element y[k],
*  1 <= k <= n. If J[k] is empty, y[k] is set to constant false.
*
*  RETURNS
*
*  The routine returns n, the number of literals in the right-hand side
*  of (2), 0 <= n <= NBIT_MAX. If the sum (sum |a[j]|) is too large, so
*  more than NBIT_MAX (= 31) literals are needed to encode the original
*  linear combination, the routine returns a negative value. */

#define NBIT_MAX 31
/* maximal number of literals in the right hand-side of (2) */

static NPPLSE *remove_lse(NPP *npp, NPPLSE *set, NPPCOL *col)
{     /* remove specified literal from specified literal set */
      NPPLSE *lse, *prev = NULL;
      for (lse = set; lse != NULL; prev = lse, lse = lse->next)
         if (lse->lit.col == col) break;
      xassert(lse != NULL);
      if (prev == NULL)
         set = lse->next;
      else
         prev->next = lse->next;
      dmp_free_atom(npp->pool, lse, sizeof(NPPLSE));
      return set;
}

int npp_sat_encode_sum_ax(NPP *npp, NPPROW *row, NPPLIT y[])
{     NPPAIJ *aij;
      NPPLSE *set[1+NBIT_MAX], *lse;
      NPPSED sed;
      int k, n, temp;
      double sum;
      /* compute the sum (sum |a[j]|) */
      sum = 0.0;
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
         sum += fabs(aij->val);
      /* determine n, the number of bits in the sum */
      temp = (int)sum;
      if ((double)temp != sum)
         return -1; /* integer arithmetic error */
      for (n = 0; temp > 0; n++, temp >>= 1);
      xassert(0 <= n && n <= NBIT_MAX);
      /* build initial sets J[k], 1 <= k <= n; see (5) */
      /* set[k] is a pointer to the list of literals in J[k] */
      for (k = 1; k <= n; k++)
         set[k] = NULL;
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  temp = (int)fabs(aij->val);
         xassert((int)temp == fabs(aij->val));
         for (k = 1; temp > 0; k++, temp >>= 1)
         {  if (temp & 1)
            {  xassert(k <= n);
               lse = dmp_get_atom(npp->pool, sizeof(NPPLSE));
               lse->lit.col = aij->col;
               lse->lit.neg = (aij->val > 0.0 ? 0 : 1);
               lse->next = set[k];
               set[k] = lse;
            }
         }
      }
      /* main transformation loop */
      for (k = 1; k <= n; k++)
      {  /* reduce J[k] and set y[k] */
         for (;;)
         {  if (set[k] == NULL)
            {  /* J[k] is empty */
               /* set y[k] to constant false */
               y[k].col = NULL, y[k].neg = 0;
               break;
            }
            if (set[k]->next == NULL)
            {  /* J[k] contains one literal */
               /* set y[k] to that literal */
               y[k] = set[k]->lit;
               dmp_free_atom(npp->pool, set[k], sizeof(NPPLSE));
               break;
            }
            if (set[k]->next->next == NULL)
            {  /* J[k] contains two literals */
               /* apply transformation (7) */
               npp_sat_encode_sum2(npp, set[k], &sed);
            }
            else
            {  /* J[k] contains at least three literals */
               /* apply transformation (8) */
               npp_sat_encode_sum3(npp, set[k], &sed);
               /* remove third literal from set[k] */
               set[k] = remove_lse(npp, set[k], sed.z.col);
            }
            /* remove second literal from set[k] */
            set[k] = remove_lse(npp, set[k], sed.y.col);
            /* remove first literal from set[k] */
            set[k] = remove_lse(npp, set[k], sed.x.col);
            /* include new literal s to set[k] */
            lse = dmp_get_atom(npp->pool, sizeof(NPPLSE));
            lse->lit.col = sed.s, lse->lit.neg = 0;
            lse->next = set[k];
            set[k] = lse;
            /* include new literal c to set[k+1] */
            xassert(k < n); /* FIXME: can "overflow" happen? */
            lse = dmp_get_atom(npp->pool, sizeof(NPPLSE));
            lse->lit.col = sed.c, lse->lit.neg = 0;
            lse->next = set[k+1];
            set[k+1] = lse;
         }
      }
      return n;
}

/***********************************************************************
*  npp_sat_normalize_clause - normalize clause
*
*  This routine normalizes the specified clause, which is a disjunction
*  of literals, by replacing multiple literals, which refer to the same
*  binary variable, with a single literal.
*
*  On exit the routine returns the number of literals in the resulting
*  clause. However, if the specified clause includes both a literal and
*  its negation, the routine returns a negative value meaning that the
*  clause is equivalent to the value true. */

int npp_sat_normalize_clause(NPP *npp, int size, NPPLIT lit[])
{     int j, k, new_size;
      xassert(npp == npp);
      xassert(size >= 0);
      new_size = 0;
      for (k = 1; k <= size; k++)
      {  for (j = 1; j <= new_size; j++)
         {  if (lit[k].col == lit[j].col)
            {  /* lit[k] refers to the same variable as lit[j], which
                  is already included in the resulting clause */
               if (lit[k].neg == lit[j].neg)
               {  /* ignore lit[k] due to the idempotent law */
                  goto skip;
               }
               else
               {  /* lit[k] is NOT lit[j]; the clause is equivalent to
                     the value true */
                  return -1;
               }
            }
         }
         /* include lit[k] in the resulting clause */
         lit[++new_size] = lit[k];
skip:    ;
      }
      return new_size;
}

/***********************************************************************
*  npp_sat_encode_clause - translate clause to cover inequality
*
*  Given a clause
*
*     OR  t[j],                                                      (1)
*   j in J
*
*  where t[j] is a literal, i.e. t[j] = x[j] or t[j] = NOT x[j], this
*  routine translates it to the following equivalent cover inequality,
*  which is added to the transformed problem:
*
*     sum t[j] >= 1,                                                 (2)
*   j in J
*
*  where t[j] = x[j] or t[j] = 1 - x[j].
*
*  If necessary, the clause should be normalized before a call to this
*  routine. */

NPPROW *npp_sat_encode_clause(NPP *npp, int size, NPPLIT lit[])
{     NPPROW *row;
      int k;
      xassert(size >= 1);
      row = npp_add_row(npp);
      row->lb = 1.0, row->ub = +DBL_MAX;
      for (k = 1; k <= size; k++)
      {  xassert(lit[k].col != NULL);
         if (lit[k].neg == 0)
            npp_add_aij(npp, row, lit[k].col, +1.0);
         else if (lit[k].neg == 1)
         {  npp_add_aij(npp, row, lit[k].col, -1.0);
            row->lb -= 1.0;
         }
         else
            xassert(lit != lit);
      }
      return row;
}

/***********************************************************************
*  npp_sat_encode_geq - encode "not less than" constraint
*
*  PURPOSE
*
*  This routine translates to CNF the following constraint:
*
*      n
*     sum 2**(k-1) * y[k] >= b,                                      (1)
*     k=1
*
*  where y[k] is either a literal (i.e. y[k] = x[k] or y[k] = 1 - x[k])
*  or constant false (zero), b is a given lower bound.
*
*  ALGORITHM
*
*  If b < 0, the constraint is redundant, so assume that b >= 0. Let
*
*          n
*     b = sum 2**(k-1) b[k],                                         (2)
*         k=1
*
*  where b[k] is k-th binary digit of b. (Note that if b >= 2**n and
*  therefore cannot be represented in the form (2), the constraint (1)
*  is infeasible.) In this case the condition (1) is equivalent to the
*  following condition:
*
*     y[n] y[n-1] ... y[2] y[1] >= b[n] b[n-1] ... b[2] b[1],        (3)
*
*  where ">=" is understood lexicographically.
*
*  Algorithmically the condition (3) can be tested as follows:
*
*     for (k = n; k >= 1; k--)
*     {  if (y[k] < b[k])
*           y is less than b;
*        if (y[k] > b[k])
*           y is greater than b;
*     }
*     y is equal to b;
*
*  Thus, y is less than b iff there exists k, 1 <= k <= n, for which
*  the following condition is satisfied:
*
*     y[n] = b[n] AND ... AND y[k+1] = b[k+1] AND y[k] < b[k].       (4)
*
*  Negating the condition (4) we have that y is not less than b iff for
*  all k, 1 <= k <= n, the following condition is satisfied:
*
*     y[n] != b[n] OR ... OR y[k+1] != b[k+1] OR y[k] >= b[k].       (5)
*
*  Note that if b[k] = 0, the literal y[k] >= b[k] is always true, in
*  which case the entire clause (5) is true and can be omitted.
*
*  RETURNS
*
*  Normally the routine returns zero. However, if the constraint (1) is
*  infeasible, the routine returns non-zero. */

int npp_sat_encode_geq(NPP *npp, int n, NPPLIT y[], int rhs)
{     NPPLIT lit[1+NBIT_MAX];
      int j, k, size, temp, b[1+NBIT_MAX];
      xassert(0 <= n && n <= NBIT_MAX);
      /* if the constraint (1) is redundant, do nothing */
      if (rhs < 0)
         return 0;
      /* determine binary digits of b according to (2) */
      for (k = 1, temp = rhs; k <= n; k++, temp >>= 1)
         b[k] = temp & 1;
      if (temp != 0)
      {  /* b >= 2**n; the constraint (1) is infeasible */
         return 1;
      }
      /* main transformation loop */
      for (k = 1; k <= n; k++)
      {  /* build the clause (5) for current k */
         size = 0; /* clause size = number of literals */
         /* add literal y[k] >= b[k] */
         if (b[k] == 0)
         {  /* b[k] = 0 -> the literal is true */
            goto skip;
         }
         else if (y[k].col == NULL)
         {  /* y[k] = 0, b[k] = 1 -> the literal is false */
            xassert(y[k].neg == 0);
         }
         else
         {  /* add literal y[k] = 1 */
            lit[++size] = y[k];
         }
         for (j = k+1; j <= n; j++)
         {  /* add literal y[j] != b[j] */
            if (y[j].col == NULL)
            {  xassert(y[j].neg == 0);
               if (b[j] == 0)
               {  /* y[j] = 0, b[j] = 0 -> the literal is false */
                  continue;
               }
               else
               {  /* y[j] = 0, b[j] = 1 -> the literal is true */
                  goto skip;
               }
            }
            else
            {  lit[++size] = y[j];
               if (b[j] != 0)
                  lit[size].neg = 1 - lit[size].neg;
            }
         }
         /* normalize the clause */
         size = npp_sat_normalize_clause(npp, size, lit);
         if (size < 0)
         {  /* the clause is equivalent to the value true */
            goto skip;
         }
         if (size == 0)
         {  /* the clause is equivalent to the value false; this means
               that the constraint (1) is infeasible */
            return 2;
         }
         /* translate the clause to corresponding cover inequality */
         npp_sat_encode_clause(npp, size, lit);
skip:    ;
      }
      return 0;
}

/***********************************************************************
*  npp_sat_encode_leq - encode "not greater than" constraint
*
*  PURPOSE
*
*  This routine translates to CNF the following constraint:
*
*      n
*     sum 2**(k-1) * y[k] <= b,                                      (1)
*     k=1
*
*  where y[k] is either a literal (i.e. y[k] = x[k] or y[k] = 1 - x[k])
*  or constant false (zero), b is a given upper bound.
*
*  ALGORITHM
*
*  If b < 0, the constraint is infeasible, so assume that b >= 0. Let
*
*          n
*     b = sum 2**(k-1) b[k],                                         (2)
*         k=1
*
*  where b[k] is k-th binary digit of b. (Note that if b >= 2**n and
*  therefore cannot be represented in the form (2), the constraint (1)
*  is redundant.) In this case the condition (1) is equivalent to the
*  following condition:
*
*     y[n] y[n-1] ... y[2] y[1] <= b[n] b[n-1] ... b[2] b[1],        (3)
*
*  where "<=" is understood lexicographically.
*
*  Algorithmically the condition (3) can be tested as follows:
*
*     for (k = n; k >= 1; k--)
*     {  if (y[k] < b[k])
*           y is less than b;
*        if (y[k] > b[k])
*           y is greater than b;
*     }
*     y is equal to b;
*
*  Thus, y is greater than b iff there exists k, 1 <= k <= n, for which
*  the following condition is satisfied:
*
*     y[n] = b[n] AND ... AND y[k+1] = b[k+1] AND y[k] > b[k].       (4)
*
*  Negating the condition (4) we have that y is not greater than b iff
*  for all k, 1 <= k <= n, the following condition is satisfied:
*
*     y[n] != b[n] OR ... OR y[k+1] != b[k+1] OR y[k] <= b[k].       (5)
*
*  Note that if b[k] = 1, the literal y[k] <= b[k] is always true, in
*  which case the entire clause (5) is true and can be omitted.
*
*  RETURNS
*
*  Normally the routine returns zero. However, if the constraint (1) is
*  infeasible, the routine returns non-zero. */

int npp_sat_encode_leq(NPP *npp, int n, NPPLIT y[], int rhs)
{     NPPLIT lit[1+NBIT_MAX];
      int j, k, size, temp, b[1+NBIT_MAX];
      xassert(0 <= n && n <= NBIT_MAX);
      /* check if the constraint (1) is infeasible */
      if (rhs < 0)
         return 1;
      /* determine binary digits of b according to (2) */
      for (k = 1, temp = rhs; k <= n; k++, temp >>= 1)
         b[k] = temp & 1;
      if (temp != 0)
      {  /* b >= 2**n; the constraint (1) is redundant */
         return 0;
      }
      /* main transformation loop */
      for (k = 1; k <= n; k++)
      {  /* build the clause (5) for current k */
         size = 0; /* clause size = number of literals */
         /* add literal y[k] <= b[k] */
         if (b[k] == 1)
         {  /* b[k] = 1 -> the literal is true */
            goto skip;
         }
         else if (y[k].col == NULL)
         {  /* y[k] = 0, b[k] = 0 -> the literal is true */
            xassert(y[k].neg == 0);
            goto skip;
         }
         else
         {  /* add literal y[k] = 0 */
            lit[++size] = y[k];
            lit[size].neg = 1 - lit[size].neg;
         }
         for (j = k+1; j <= n; j++)
         {  /* add literal y[j] != b[j] */
            if (y[j].col == NULL)
            {  xassert(y[j].neg == 0);
               if (b[j] == 0)
               {  /* y[j] = 0, b[j] = 0 -> the literal is false */
                  continue;
               }
               else
               {  /* y[j] = 0, b[j] = 1 -> the literal is true */
                  goto skip;
               }
            }
            else
            {  lit[++size] = y[j];
               if (b[j] != 0)
                  lit[size].neg = 1 - lit[size].neg;
            }
         }
         /* normalize the clause */
         size = npp_sat_normalize_clause(npp, size, lit);
         if (size < 0)
         {  /* the clause is equivalent to the value true */
            goto skip;
         }
         if (size == 0)
         {  /* the clause is equivalent to the value false; this means
               that the constraint (1) is infeasible */
            return 2;
         }
         /* translate the clause to corresponding cover inequality */
         npp_sat_encode_clause(npp, size, lit);
skip:    ;
      }
      return 0;
}

/***********************************************************************
*  npp_sat_encode_row - encode constraint (row) of general type
*
*  PURPOSE
*
*  This routine translates to CNF the following constraint (row):
*
*     L <= sum a[j] x[j] <= U,                                       (1)
*           j
*
*  where all x[j] are binary variables.
*
*  ALGORITHM
*
*  First, the routine performs substitution x[j] = t[j] for j in J+
*  and x[j] = 1 - t[j] for j in J-, where J+ = { j : a[j] > 0 } and
*  J- = { j : a[j] < 0 }. This gives:
*
*     L <=  sum  a[j] t[j] +   sum  a[j] (1 - t[j]) <= U  ==>
*         j in J+            j in J-
*
*     L' <= sum |a[j]| t[j] <= U',                                   (2)
*            j
*
*  where
*
*     L' = L -   sum  a[j],   U' = U -   sum  a[j].                  (3)
*              j in J-                 j in J-
*
*  (Actually only new bounds L' and U' are computed.)
*
*  Then the routine translates to CNF the following equality:
*
*                        n
*     sum |a[j]| t[j] = sum 2**(k-1) * y[k],                         (4)
*      j                k=1
*
*  where y[k] is either some t[j] or a new literal or a constant zero
*  (see the routine npp_sat_encode_sum_ax).
*
*  Finally, the routine translates to CNF the following conditions:
*
*      n
*     sum 2**(k-1) * y[k] >= L'                                      (5)
*     k=1
*
*  and
*
*      n
*     sum 2**(k-1) * y[k] <= U'                                      (6)
*     k=1
*
*  (see the routines npp_sat_encode_geq and npp_sat_encode_leq).
*
*  All resulting clauses are encoded as cover inequalities and included
*  into the transformed problem.
*
*  Note that on exit the routine removes the specified constraint (row)
*  from the original problem.
*
*  RETURNS
*
*  The routine returns one of the following codes:
*
*  0 - translation was successful;
*  1 - constraint (1) was found infeasible;
*  2 - integer arithmetic error occured. */

int npp_sat_encode_row(NPP *npp, NPPROW *row)
{     NPPAIJ *aij;
      NPPLIT y[1+NBIT_MAX];
      int n, rhs;
      double lb, ub;
      /* the row should not be free */
      xassert(!(row->lb == -DBL_MAX && row->ub == +DBL_MAX));
      /* compute new bounds L' and U' (3) */
      lb = row->lb;
      ub = row->ub;
      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
      {  if (aij->val < 0.0)
         {  if (lb != -DBL_MAX)
               lb -= aij->val;
            if (ub != -DBL_MAX)
               ub -= aij->val;
         }
      }
      /* encode the equality (4) */
      n = npp_sat_encode_sum_ax(npp, row, y);
      if (n < 0)
         return 2; /* integer arithmetic error */
      /* encode the condition (5) */
      if (lb != -DBL_MAX)
      {  rhs = (int)lb;
         if ((double)rhs != lb)
            return 2; /* integer arithmetic error */
         if (npp_sat_encode_geq(npp, n, y, rhs) != 0)
            return 1; /* original constraint is infeasible */
      }
      /* encode the condition (6) */
      if (ub != +DBL_MAX)
      {  rhs = (int)ub;
         if ((double)rhs != ub)
            return 2; /* integer arithmetic error */
         if (npp_sat_encode_leq(npp, n, y, rhs) != 0)
            return 1; /* original constraint is infeasible */
      }
      /* remove the specified row from the problem */
      npp_del_row(npp, row);
      return 0;
}

/***********************************************************************
*  npp_sat_encode_prob - encode 0-1 feasibility problem
*
*  This routine translates the specified 0-1 feasibility problem to an
*  equivalent SAT-CNF problem.
*
*  N.B. Currently this is a very crude implementation.
*
*  RETURNS
*
*  0           success;
*
*  GLP_ENOPFS  primal/integer infeasibility detected;
*
*  GLP_ERANGE  integer overflow occured. */

int npp_sat_encode_prob(NPP *npp)
{     NPPROW *row, *next_row, *prev_row;
      NPPCOL *col, *next_col;
      int cover = 0, pack = 0, partn = 0, ret;
      /* process and remove free rows */
      for (row = npp->r_head; row != NULL; row = next_row)
      {  next_row = row->next;
         if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
            npp_sat_free_row(npp, row);
      }
      /* process and remove fixed columns */
      for (col = npp->c_head; col != NULL; col = next_col)
      {  next_col = col->next;
         if (col->lb == col->ub)
            xassert(npp_sat_fixed_col(npp, col) == 0);
      }
      /* only binary variables should remain */
      for (col = npp->c_head; col != NULL; col = col->next)
         xassert(col->is_int && col->lb == 0.0 && col->ub == 1.0);
      /* new rows may be added to the end of the row list, so we walk
         from the end to beginning of the list */
      for (row = npp->r_tail; row != NULL; row = prev_row)
      {  prev_row = row->prev;
         /* process special cases */
         ret = npp_sat_is_cover_ineq(npp, row);
         if (ret != 0)
         {  /* row is covering inequality */
            cover++;
            /* since it already encodes a clause, just transform it to
               canonical form */
            if (ret == 2)
            {  xassert(npp_sat_reverse_row(npp, row) == 0);
               ret = npp_sat_is_cover_ineq(npp, row);
            }
            xassert(ret == 1);
            continue;
         }
         ret = npp_sat_is_partn_eq(npp, row);
         if (ret != 0)
         {  /* row is partitioning equality */
            NPPROW *cov;
            NPPAIJ *aij;
            partn++;
            /* transform it to canonical form */
            if (ret == 2)
            {  xassert(npp_sat_reverse_row(npp, row) == 0);
               ret = npp_sat_is_partn_eq(npp, row);
            }
            xassert(ret == 1);
            /* and split it into covering and packing inequalities,
               both in canonical forms */
            cov = npp_add_row(npp);
            cov->lb = row->lb, cov->ub = +DBL_MAX;
            for (aij = row->ptr; aij != NULL; aij = aij->r_next)
               npp_add_aij(npp, cov, aij->col, aij->val);
            xassert(npp_sat_is_cover_ineq(npp, cov) == 1);
            /* the cover inequality already encodes a clause and do
               not need any further processing */
            row->lb = -DBL_MAX;
            xassert(npp_sat_is_pack_ineq(npp, row) == 1);
            /* the packing inequality will be processed below */
            pack--;
         }
         ret = npp_sat_is_pack_ineq(npp, row);
         if (ret != 0)
         {  /* row is packing inequality */
            NPPROW *rrr;
            int nlit, desired_nlit = 4;
            pack++;
            /* transform it to canonical form */
            if (ret == 2)
            {  xassert(npp_sat_reverse_row(npp, row) == 0);
               ret = npp_sat_is_pack_ineq(npp, row);
            }
            xassert(ret == 1);
            /* process the packing inequality */
            for (;;)
            {  /* determine the number of literals in the remaining
                  inequality */
               nlit = npp_row_nnz(npp, row);
               if (nlit <= desired_nlit)
                  break;
               /* split the current inequality into one having not more
                  than desired_nlit literals and remaining one */
               rrr = npp_sat_split_pack(npp, row, desired_nlit-1);
               /* translate the former inequality to CNF and remove it
                  from the original problem */
               npp_sat_encode_pack(npp, rrr);
            }
            /* translate the remaining inequality to CNF and remove it
               from the original problem */
            npp_sat_encode_pack(npp, row);
            continue;
         }
         /* translate row of general type to CNF and remove it from the
            original problem */
         ret = npp_sat_encode_row(npp, row);
         if (ret == 0)
            ;
         else if (ret == 1)
            ret = GLP_ENOPFS;
         else if (ret == 2)
            ret = GLP_ERANGE;
         else
            xassert(ret != ret);
         if (ret != 0)
            goto done;
      }
      ret = 0;
      if (cover != 0)
         xprintf("%d covering inequalities\n", cover);
      if (pack != 0)
         xprintf("%d packing inequalities\n", pack);
      if (partn != 0)
         xprintf("%d partitioning equalities\n", partn);
done: return ret;
}

/* eof */