1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpnpp06.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,1500 @@
     1.4 +/* glpnpp06.c (translate feasibility problem to CNF-SAT) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpnpp.h"
    1.29 +
    1.30 +/***********************************************************************
    1.31 +*  npp_sat_free_row - process free (unbounded) row
    1.32 +*
    1.33 +*  This routine processes row p, which is free (i.e. has no finite
    1.34 +*  bounds):
    1.35 +*
    1.36 +*     -inf < sum a[p,j] x[j] < +inf.                                 (1)
    1.37 +*
    1.38 +*  The constraint (1) cannot be active and therefore it is redundant,
    1.39 +*  so the routine simply removes it from the original problem. */
    1.40 +
    1.41 +void npp_sat_free_row(NPP *npp, NPPROW *p)
    1.42 +{     /* the row should be free */
    1.43 +      xassert(p->lb == -DBL_MAX && p->ub == +DBL_MAX);
    1.44 +      /* remove the row from the problem */
    1.45 +      npp_del_row(npp, p);
    1.46 +      return;
    1.47 +}
    1.48 +
    1.49 +/***********************************************************************
    1.50 +*  npp_sat_fixed_col - process fixed column
    1.51 +*
    1.52 +*  This routine processes column q, which is fixed:
    1.53 +*
    1.54 +*     x[q] = s[q],                                                   (1)
    1.55 +*
    1.56 +*  where s[q] is a fixed column value.
    1.57 +*
    1.58 +*  The routine substitutes fixed value s[q] into constraint rows and
    1.59 +*  then removes column x[q] from the original problem.
    1.60 +*
    1.61 +*  Substitution of x[q] = s[q] into row i gives:
    1.62 +*
    1.63 +*     L[i] <= sum a[i,j] x[j] <= U[i]   ==>
    1.64 +*              j
    1.65 +*
    1.66 +*     L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i]   ==>
    1.67 +*            j!=q
    1.68 +*
    1.69 +*     L[i] <= sum a[i,j] x[j] + a[i,q] s[q] <= U[i]   ==>
    1.70 +*            j!=q
    1.71 +*
    1.72 +*     L~[i] <= sum a[i,j] x[j] <= U~[i],
    1.73 +*             j!=q
    1.74 +*
    1.75 +*  where
    1.76 +*
    1.77 +*     L~[i] = L[i] - a[i,q] s[q],                                    (2)
    1.78 +*
    1.79 +*     U~[i] = U[i] - a[i,q] s[q]                                     (3)
    1.80 +*
    1.81 +*  are, respectively, lower and upper bound of row i in the transformed
    1.82 +*  problem.
    1.83 +*
    1.84 +*  On recovering solution x[q] is assigned the value of s[q]. */
    1.85 +
    1.86 +struct sat_fixed_col
    1.87 +{     /* fixed column */
    1.88 +      int q;
    1.89 +      /* column reference number for variable x[q] */
    1.90 +      int s;
    1.91 +      /* value, at which x[q] is fixed */
    1.92 +};
    1.93 +
    1.94 +static int rcv_sat_fixed_col(NPP *, void *);
    1.95 +
    1.96 +int npp_sat_fixed_col(NPP *npp, NPPCOL *q)
    1.97 +{     struct sat_fixed_col *info;
    1.98 +      NPPROW *i;
    1.99 +      NPPAIJ *aij;
   1.100 +      int temp;
   1.101 +      /* the column should be fixed */
   1.102 +      xassert(q->lb == q->ub);
   1.103 +      /* create transformation stack entry */
   1.104 +      info = npp_push_tse(npp,
   1.105 +         rcv_sat_fixed_col, sizeof(struct sat_fixed_col));
   1.106 +      info->q = q->j;
   1.107 +      info->s = (int)q->lb;
   1.108 +      xassert((double)info->s == q->lb);
   1.109 +      /* substitute x[q] = s[q] into constraint rows */
   1.110 +      if (info->s == 0)
   1.111 +         goto skip;
   1.112 +      for (aij = q->ptr; aij != NULL; aij = aij->c_next)
   1.113 +      {  i = aij->row;
   1.114 +         if (i->lb != -DBL_MAX)
   1.115 +         {  i->lb -= aij->val * (double)info->s;
   1.116 +            temp = (int)i->lb;
   1.117 +            if ((double)temp != i->lb)
   1.118 +               return 1; /* integer arithmetic error */
   1.119 +         }
   1.120 +         if (i->ub != +DBL_MAX)
   1.121 +         {  i->ub -= aij->val * (double)info->s;
   1.122 +            temp = (int)i->ub;
   1.123 +            if ((double)temp != i->ub)
   1.124 +               return 2; /* integer arithmetic error */
   1.125 +         }
   1.126 +      }
   1.127 +skip: /* remove the column from the problem */
   1.128 +      npp_del_col(npp, q);
   1.129 +      return 0;
   1.130 +}
   1.131 +
   1.132 +static int rcv_sat_fixed_col(NPP *npp, void *info_)
   1.133 +{     struct sat_fixed_col *info = info_;
   1.134 +      npp->c_value[info->q] = (double)info->s;
   1.135 +      return 0;
   1.136 +}
   1.137 +
   1.138 +/***********************************************************************
   1.139 +*  npp_sat_is_bin_comb - test if row is binary combination
   1.140 +*
   1.141 +*  This routine tests if the specified row is a binary combination,
   1.142 +*  i.e. all its constraint coefficients are +1 and -1 and all variables
   1.143 +*  are binary. If the test was passed, the routine returns non-zero,
   1.144 +*  otherwise zero. */
   1.145 +
   1.146 +int npp_sat_is_bin_comb(NPP *npp, NPPROW *row)
   1.147 +{     NPPCOL *col;
   1.148 +      NPPAIJ *aij;
   1.149 +      xassert(npp == npp);
   1.150 +      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
   1.151 +      {  if (!(aij->val == +1.0 || aij->val == -1.0))
   1.152 +            return 0; /* non-unity coefficient */
   1.153 +         col = aij->col;
   1.154 +         if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
   1.155 +            return 0; /* non-binary column */
   1.156 +      }
   1.157 +      return 1; /* test was passed */
   1.158 +}
   1.159 +
   1.160 +/***********************************************************************
   1.161 +*  npp_sat_num_pos_coef - determine number of positive coefficients
   1.162 +*
   1.163 +*  This routine returns the number of positive coefficients in the
   1.164 +*  specified row. */
   1.165 +
   1.166 +int npp_sat_num_pos_coef(NPP *npp, NPPROW *row)
   1.167 +{     NPPAIJ *aij;
   1.168 +      int num = 0;
   1.169 +      xassert(npp == npp);
   1.170 +      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
   1.171 +      {  if (aij->val > 0.0)
   1.172 +            num++;
   1.173 +      }
   1.174 +      return num;
   1.175 +}
   1.176 +
   1.177 +/***********************************************************************
   1.178 +*  npp_sat_num_neg_coef - determine number of negative coefficients
   1.179 +*
   1.180 +*  This routine returns the number of negative coefficients in the
   1.181 +*  specified row. */
   1.182 +
   1.183 +int npp_sat_num_neg_coef(NPP *npp, NPPROW *row)
   1.184 +{     NPPAIJ *aij;
   1.185 +      int num = 0;
   1.186 +      xassert(npp == npp);
   1.187 +      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
   1.188 +      {  if (aij->val < 0.0)
   1.189 +            num++;
   1.190 +      }
   1.191 +      return num;
   1.192 +}
   1.193 +
   1.194 +/***********************************************************************
   1.195 +*  npp_sat_is_cover_ineq - test if row is covering inequality
   1.196 +*
   1.197 +*  The canonical form of a covering inequality is the following:
   1.198 +*
   1.199 +*     sum x[j] >= 1,                                                 (1)
   1.200 +*   j in J
   1.201 +*
   1.202 +*  where all x[j] are binary variables.
   1.203 +*
   1.204 +*  In general case a covering inequality may have one of the following
   1.205 +*  two forms:
   1.206 +*
   1.207 +*     sum  x[j] -  sum  x[j] >= 1 - |J-|,                            (2)
   1.208 +*   j in J+      j in J-
   1.209 +*
   1.210 +*
   1.211 +*     sum  x[j] -  sum  x[j] <= |J+| - 1.                            (3)
   1.212 +*   j in J+      j in J-
   1.213 +*
   1.214 +*  Obviously, the inequality (2) can be transformed to the form (1) by
   1.215 +*  substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the
   1.216 +*  negation of variable x[j]. And the inequality (3) can be transformed
   1.217 +*  to (2) by multiplying both left- and right-hand sides by -1.
   1.218 +*
   1.219 +*  This routine returns one of the following codes:
   1.220 +*
   1.221 +*  0, if the specified row is not a covering inequality;
   1.222 +*
   1.223 +*  1, if the specified row has the form (2);
   1.224 +*
   1.225 +*  2, if the specified row has the form (3). */
   1.226 +
   1.227 +int npp_sat_is_cover_ineq(NPP *npp, NPPROW *row)
   1.228 +{     xassert(npp == npp);
   1.229 +      if (row->lb != -DBL_MAX && row->ub == +DBL_MAX)
   1.230 +      {  /* row is inequality of '>=' type */
   1.231 +         if (npp_sat_is_bin_comb(npp, row))
   1.232 +         {  /* row is a binary combination */
   1.233 +            if (row->lb == 1.0 - npp_sat_num_neg_coef(npp, row))
   1.234 +            {  /* row has the form (2) */
   1.235 +               return 1;
   1.236 +            }
   1.237 +         }
   1.238 +      }
   1.239 +      else if (row->lb == -DBL_MAX && row->ub != +DBL_MAX)
   1.240 +      {  /* row is inequality of '<=' type */
   1.241 +         if (npp_sat_is_bin_comb(npp, row))
   1.242 +         {  /* row is a binary combination */
   1.243 +            if (row->ub == npp_sat_num_pos_coef(npp, row) - 1.0)
   1.244 +            {  /* row has the form (3) */
   1.245 +               return 2;
   1.246 +            }
   1.247 +         }
   1.248 +      }
   1.249 +      /* row is not a covering inequality */
   1.250 +      return 0;
   1.251 +}
   1.252 +
   1.253 +/***********************************************************************
   1.254 +*  npp_sat_is_pack_ineq - test if row is packing inequality
   1.255 +*
   1.256 +*  The canonical form of a packing inequality is the following:
   1.257 +*
   1.258 +*     sum x[j] <= 1,                                                 (1)
   1.259 +*   j in J
   1.260 +*
   1.261 +*  where all x[j] are binary variables.
   1.262 +*
   1.263 +*  In general case a packing inequality may have one of the following
   1.264 +*  two forms:
   1.265 +*
   1.266 +*     sum  x[j] -  sum  x[j] <= 1 - |J-|,                            (2)
   1.267 +*   j in J+      j in J-
   1.268 +*
   1.269 +*
   1.270 +*     sum  x[j] -  sum  x[j] >= |J+| - 1.                            (3)
   1.271 +*   j in J+      j in J-
   1.272 +*
   1.273 +*  Obviously, the inequality (2) can be transformed to the form (1) by
   1.274 +*  substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the
   1.275 +*  negation of variable x[j]. And the inequality (3) can be transformed
   1.276 +*  to (2) by multiplying both left- and right-hand sides by -1.
   1.277 +*
   1.278 +*  This routine returns one of the following codes:
   1.279 +*
   1.280 +*  0, if the specified row is not a packing inequality;
   1.281 +*
   1.282 +*  1, if the specified row has the form (2);
   1.283 +*
   1.284 +*  2, if the specified row has the form (3). */
   1.285 +
   1.286 +int npp_sat_is_pack_ineq(NPP *npp, NPPROW *row)
   1.287 +{     xassert(npp == npp);
   1.288 +      if (row->lb == -DBL_MAX && row->ub != +DBL_MAX)
   1.289 +      {  /* row is inequality of '<=' type */
   1.290 +         if (npp_sat_is_bin_comb(npp, row))
   1.291 +         {  /* row is a binary combination */
   1.292 +            if (row->ub == 1.0 - npp_sat_num_neg_coef(npp, row))
   1.293 +            {  /* row has the form (2) */
   1.294 +               return 1;
   1.295 +            }
   1.296 +         }
   1.297 +      }
   1.298 +      else if (row->lb != -DBL_MAX && row->ub == +DBL_MAX)
   1.299 +      {  /* row is inequality of '>=' type */
   1.300 +         if (npp_sat_is_bin_comb(npp, row))
   1.301 +         {  /* row is a binary combination */
   1.302 +            if (row->lb == npp_sat_num_pos_coef(npp, row) - 1.0)
   1.303 +            {  /* row has the form (3) */
   1.304 +               return 2;
   1.305 +            }
   1.306 +         }
   1.307 +      }
   1.308 +      /* row is not a packing inequality */
   1.309 +      return 0;
   1.310 +}
   1.311 +
   1.312 +/***********************************************************************
   1.313 +*  npp_sat_is_partn_eq - test if row is partitioning equality
   1.314 +*
   1.315 +*  The canonical form of a partitioning equality is the following:
   1.316 +*
   1.317 +*     sum x[j] = 1,                                                  (1)
   1.318 +*   j in J
   1.319 +*
   1.320 +*  where all x[j] are binary variables.
   1.321 +*
   1.322 +*  In general case a partitioning equality may have one of the following
   1.323 +*  two forms:
   1.324 +*
   1.325 +*     sum  x[j] -  sum  x[j] = 1 - |J-|,                             (2)
   1.326 +*   j in J+      j in J-
   1.327 +*
   1.328 +*
   1.329 +*     sum  x[j] -  sum  x[j] = |J+| - 1.                             (3)
   1.330 +*   j in J+      j in J-
   1.331 +*
   1.332 +*  Obviously, the equality (2) can be transformed to the form (1) by
   1.333 +*  substitution x[j] = 1 - x'[j] for all j in J-, where x'[j] is the
   1.334 +*  negation of variable x[j]. And the equality (3) can be transformed
   1.335 +*  to (2) by multiplying both left- and right-hand sides by -1.
   1.336 +*
   1.337 +*  This routine returns one of the following codes:
   1.338 +*
   1.339 +*  0, if the specified row is not a partitioning equality;
   1.340 +*
   1.341 +*  1, if the specified row has the form (2);
   1.342 +*
   1.343 +*  2, if the specified row has the form (3). */
   1.344 +
   1.345 +int npp_sat_is_partn_eq(NPP *npp, NPPROW *row)
   1.346 +{     xassert(npp == npp);
   1.347 +      if (row->lb == row->ub)
   1.348 +      {  /* row is equality constraint */
   1.349 +         if (npp_sat_is_bin_comb(npp, row))
   1.350 +         {  /* row is a binary combination */
   1.351 +            if (row->lb == 1.0 - npp_sat_num_neg_coef(npp, row))
   1.352 +            {  /* row has the form (2) */
   1.353 +               return 1;
   1.354 +            }
   1.355 +            if (row->ub == npp_sat_num_pos_coef(npp, row) - 1.0)
   1.356 +            {  /* row has the form (3) */
   1.357 +               return 2;
   1.358 +            }
   1.359 +         }
   1.360 +      }
   1.361 +      /* row is not a partitioning equality */
   1.362 +      return 0;
   1.363 +}
   1.364 +
   1.365 +/***********************************************************************
   1.366 +*  npp_sat_reverse_row - multiply both sides of row by -1
   1.367 +*
   1.368 +*  This routines multiplies by -1 both left- and right-hand sides of
   1.369 +*  the specified row:
   1.370 +*
   1.371 +*     L <= sum x[j] <= U,
   1.372 +*
   1.373 +*  that results in the following row:
   1.374 +*
   1.375 +*     -U <= sum (-x[j]) <= -L.
   1.376 +*
   1.377 +*  If no integer overflow occured, the routine returns zero, otherwise
   1.378 +*  non-zero. */
   1.379 +
   1.380 +int npp_sat_reverse_row(NPP *npp, NPPROW *row)
   1.381 +{     NPPAIJ *aij;
   1.382 +      int temp, ret = 0;
   1.383 +      double old_lb, old_ub;
   1.384 +      xassert(npp == npp);
   1.385 +      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
   1.386 +      {  aij->val = -aij->val;
   1.387 +         temp = (int)aij->val;
   1.388 +         if ((double)temp != aij->val)
   1.389 +            ret = 1;
   1.390 +      }
   1.391 +      old_lb = row->lb, old_ub = row->ub;
   1.392 +      if (old_ub == +DBL_MAX)
   1.393 +         row->lb = -DBL_MAX;
   1.394 +      else
   1.395 +      {  row->lb = -old_ub;
   1.396 +         temp = (int)row->lb;
   1.397 +         if ((double)temp != row->lb)
   1.398 +            ret = 2;
   1.399 +      }
   1.400 +      if (old_lb == -DBL_MAX)
   1.401 +         row->ub = +DBL_MAX;
   1.402 +      else
   1.403 +      {  row->ub = -old_lb;
   1.404 +         temp = (int)row->ub;
   1.405 +         if ((double)temp != row->ub)
   1.406 +            ret = 3;
   1.407 +      }
   1.408 +      return ret;
   1.409 +}
   1.410 +
   1.411 +/***********************************************************************
   1.412 +*  npp_sat_split_pack - split packing inequality
   1.413 +*
   1.414 +*  Let there be given a packing inequality in canonical form:
   1.415 +*
   1.416 +*     sum  t[j] <= 1,                                                (1)
   1.417 +*   j in J
   1.418 +*
   1.419 +*  where t[j] = x[j] or t[j] = 1 - x[j], x[j] is a binary variable.
   1.420 +*  And let J = J1 U J2 is a partition of the set of literals. Then the
   1.421 +*  inequality (1) is obviously equivalent to the following two packing
   1.422 +*  inequalities:
   1.423 +*
   1.424 +*     sum  t[j] <= y       <-->   sum  t[j] + (1 - y) <= 1,          (2)
   1.425 +*   j in J1                     j in J1
   1.426 +*
   1.427 +*     sum  t[j] <= 1 - y   <-->   sum  t[j] + y       <= 1,          (3)
   1.428 +*   j in J2                     j in J2
   1.429 +*
   1.430 +*  where y is a new binary variable added to the transformed problem.
   1.431 +*
   1.432 +*  Assuming that the specified row is a packing inequality (1), this
   1.433 +*  routine constructs the set J1 by including there first nlit literals
   1.434 +*  (terms) from the specified row, and the set J2 = J \ J1. Then the
   1.435 +*  routine creates a new row, which corresponds to inequality (2), and
   1.436 +*  replaces the specified row with inequality (3). */
   1.437 +
   1.438 +NPPROW *npp_sat_split_pack(NPP *npp, NPPROW *row, int nlit)
   1.439 +{     NPPROW *rrr;
   1.440 +      NPPCOL *col;
   1.441 +      NPPAIJ *aij;
   1.442 +      int k;
   1.443 +      /* original row should be packing inequality (1) */
   1.444 +      xassert(npp_sat_is_pack_ineq(npp, row) == 1);
   1.445 +      /* and nlit should be less than the number of literals (terms)
   1.446 +         in the original row */
   1.447 +      xassert(0 < nlit && nlit < npp_row_nnz(npp, row));
   1.448 +      /* create new row corresponding to inequality (2) */
   1.449 +      rrr = npp_add_row(npp);
   1.450 +      rrr->lb = -DBL_MAX, rrr->ub = 1.0;
   1.451 +      /* move first nlit literals (terms) from the original row to the
   1.452 +         new row; the original row becomes inequality (3) */
   1.453 +      for (k = 1; k <= nlit; k++)
   1.454 +      {  aij = row->ptr;
   1.455 +         xassert(aij != NULL);
   1.456 +         /* add literal to the new row */
   1.457 +         npp_add_aij(npp, rrr, aij->col, aij->val);
   1.458 +         /* correct rhs */
   1.459 +         if (aij->val < 0.0)
   1.460 +            rrr->ub -= 1.0, row->ub += 1.0;
   1.461 +         /* remove literal from the original row */
   1.462 +         npp_del_aij(npp, aij);
   1.463 +      }
   1.464 +      /* create new binary variable y */
   1.465 +      col = npp_add_col(npp);
   1.466 +      col->is_int = 1, col->lb = 0.0, col->ub = 1.0;
   1.467 +      /* include literal (1 - y) in the new row */
   1.468 +      npp_add_aij(npp, rrr, col, -1.0);
   1.469 +      rrr->ub -= 1.0;
   1.470 +      /* include literal y in the original row */
   1.471 +      npp_add_aij(npp, row, col, +1.0);
   1.472 +      return rrr;
   1.473 +}
   1.474 +
   1.475 +/***********************************************************************
   1.476 +*  npp_sat_encode_pack - encode packing inequality
   1.477 +*
   1.478 +*  Given a packing inequality in canonical form:
   1.479 +*
   1.480 +*     sum  t[j] <= 1,                                                (1)
   1.481 +*   j in J
   1.482 +*
   1.483 +*  where t[j] = x[j] or t[j] = 1 - x[j], x[j] is a binary variable,
   1.484 +*  this routine translates it to CNF by replacing it with the following
   1.485 +*  equivalent set of edge packing inequalities:
   1.486 +*
   1.487 +*     t[j] + t[k] <= 1   for all j, k in J, j != k.                  (2)
   1.488 +*
   1.489 +*  Then the routine transforms each edge packing inequality (2) to
   1.490 +*  corresponding covering inequality (that encodes two-literal clause)
   1.491 +*  by multiplying both its part by -1:
   1.492 +*
   1.493 +*     - t[j] - t[k] >= -1   <-->   (1 - t[j]) + (1 - t[k]) >= 1.     (3)
   1.494 +*
   1.495 +*  On exit the routine removes the original row from the problem. */
   1.496 +
   1.497 +void npp_sat_encode_pack(NPP *npp, NPPROW *row)
   1.498 +{     NPPROW *rrr;
   1.499 +      NPPAIJ *aij, *aik;
   1.500 +      /* original row should be packing inequality (1) */
   1.501 +      xassert(npp_sat_is_pack_ineq(npp, row) == 1);
   1.502 +      /* create equivalent system of covering inequalities (3) */
   1.503 +      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
   1.504 +      {  /* due to symmetry only one of inequalities t[j] + t[k] <= 1
   1.505 +            and t[k] <= t[j] <= 1 can be considered */
   1.506 +         for (aik = aij->r_next; aik != NULL; aik = aik->r_next)
   1.507 +         {  /* create edge packing inequality (2) */
   1.508 +            rrr = npp_add_row(npp);
   1.509 +            rrr->lb = -DBL_MAX, rrr->ub = 1.0;
   1.510 +            npp_add_aij(npp, rrr, aij->col, aij->val);
   1.511 +            if (aij->val < 0.0)
   1.512 +               rrr->ub -= 1.0;
   1.513 +            npp_add_aij(npp, rrr, aik->col, aik->val);
   1.514 +            if (aik->val < 0.0)
   1.515 +               rrr->ub -= 1.0;
   1.516 +            /* and transform it to covering inequality (3) */
   1.517 +            npp_sat_reverse_row(npp, rrr);
   1.518 +            xassert(npp_sat_is_cover_ineq(npp, rrr) == 1);
   1.519 +         }
   1.520 +      }
   1.521 +      /* remove the original row from the problem */
   1.522 +      npp_del_row(npp, row);
   1.523 +      return;
   1.524 +}
   1.525 +
   1.526 +/***********************************************************************
   1.527 +*  npp_sat_encode_sum2 - encode 2-bit summation
   1.528 +*
   1.529 +*  Given a set containing two literals x and y this routine encodes
   1.530 +*  the equality
   1.531 +*
   1.532 +*     x + y = s + 2 * c,                                             (1)
   1.533 +*
   1.534 +*  where
   1.535 +*
   1.536 +*     s = (x + y) % 2                                                (2)
   1.537 +*
   1.538 +*  is a binary variable modeling the low sum bit, and
   1.539 +*
   1.540 +*     c = (x + y) / 2                                                (3)
   1.541 +*
   1.542 +*  is a binary variable modeling the high (carry) sum bit. */
   1.543 +
   1.544 +void npp_sat_encode_sum2(NPP *npp, NPPLSE *set, NPPSED *sed)
   1.545 +{     NPPROW *row;
   1.546 +      int x, y, s, c;
   1.547 +      /* the set should contain exactly two literals */
   1.548 +      xassert(set != NULL);
   1.549 +      xassert(set->next != NULL);
   1.550 +      xassert(set->next->next == NULL);
   1.551 +      sed->x = set->lit;
   1.552 +      xassert(sed->x.neg == 0 || sed->x.neg == 1);
   1.553 +      sed->y = set->next->lit;
   1.554 +      xassert(sed->y.neg == 0 || sed->y.neg == 1);
   1.555 +      sed->z.col = NULL, sed->z.neg = 0;
   1.556 +      /* perform encoding s = (x + y) % 2 */
   1.557 +      sed->s = npp_add_col(npp);
   1.558 +      sed->s->is_int = 1, sed->s->lb = 0.0, sed->s->ub = 1.0;
   1.559 +      for (x = 0; x <= 1; x++)
   1.560 +      {  for (y = 0; y <= 1; y++)
   1.561 +         {  for (s = 0; s <= 1; s++)
   1.562 +            {  if ((x + y) % 2 != s)
   1.563 +               {  /* generate CNF clause to disable infeasible
   1.564 +                     combination */
   1.565 +                  row = npp_add_row(npp);
   1.566 +                  row->lb = 1.0, row->ub = +DBL_MAX;
   1.567 +                  if (x == sed->x.neg)
   1.568 +                     npp_add_aij(npp, row, sed->x.col, +1.0);
   1.569 +                  else
   1.570 +                  {  npp_add_aij(npp, row, sed->x.col, -1.0);
   1.571 +                     row->lb -= 1.0;
   1.572 +                  }
   1.573 +                  if (y == sed->y.neg)
   1.574 +                     npp_add_aij(npp, row, sed->y.col, +1.0);
   1.575 +                  else
   1.576 +                  {  npp_add_aij(npp, row, sed->y.col, -1.0);
   1.577 +                     row->lb -= 1.0;
   1.578 +                  }
   1.579 +                  if (s == 0)
   1.580 +                     npp_add_aij(npp, row, sed->s, +1.0);
   1.581 +                  else
   1.582 +                  {  npp_add_aij(npp, row, sed->s, -1.0);
   1.583 +                     row->lb -= 1.0;
   1.584 +                  }
   1.585 +               }
   1.586 +            }
   1.587 +         }
   1.588 +      }
   1.589 +      /* perform encoding c = (x + y) / 2 */
   1.590 +      sed->c = npp_add_col(npp);
   1.591 +      sed->c->is_int = 1, sed->c->lb = 0.0, sed->c->ub = 1.0;
   1.592 +      for (x = 0; x <= 1; x++)
   1.593 +      {  for (y = 0; y <= 1; y++)
   1.594 +         {  for (c = 0; c <= 1; c++)
   1.595 +            {  if ((x + y) / 2 != c)
   1.596 +               {  /* generate CNF clause to disable infeasible
   1.597 +                     combination */
   1.598 +                  row = npp_add_row(npp);
   1.599 +                  row->lb = 1.0, row->ub = +DBL_MAX;
   1.600 +                  if (x == sed->x.neg)
   1.601 +                     npp_add_aij(npp, row, sed->x.col, +1.0);
   1.602 +                  else
   1.603 +                  {  npp_add_aij(npp, row, sed->x.col, -1.0);
   1.604 +                     row->lb -= 1.0;
   1.605 +                  }
   1.606 +                  if (y == sed->y.neg)
   1.607 +                     npp_add_aij(npp, row, sed->y.col, +1.0);
   1.608 +                  else
   1.609 +                  {  npp_add_aij(npp, row, sed->y.col, -1.0);
   1.610 +                     row->lb -= 1.0;
   1.611 +                  }
   1.612 +                  if (c == 0)
   1.613 +                     npp_add_aij(npp, row, sed->c, +1.0);
   1.614 +                  else
   1.615 +                  {  npp_add_aij(npp, row, sed->c, -1.0);
   1.616 +                     row->lb -= 1.0;
   1.617 +                  }
   1.618 +               }
   1.619 +            }
   1.620 +         }
   1.621 +      }
   1.622 +      return;
   1.623 +}
   1.624 +
   1.625 +/***********************************************************************
   1.626 +*  npp_sat_encode_sum3 - encode 3-bit summation
   1.627 +*
   1.628 +*  Given a set containing at least three literals this routine chooses
   1.629 +*  some literals x, y, z from that set and encodes the equality
   1.630 +*
   1.631 +*     x + y + z = s + 2 * c,                                         (1)
   1.632 +*
   1.633 +*  where
   1.634 +*
   1.635 +*     s = (x + y + z) % 2                                            (2)
   1.636 +*
   1.637 +*  is a binary variable modeling the low sum bit, and
   1.638 +*
   1.639 +*     c = (x + y + z) / 2                                            (3)
   1.640 +*
   1.641 +*  is a binary variable modeling the high (carry) sum bit. */
   1.642 +
   1.643 +void npp_sat_encode_sum3(NPP *npp, NPPLSE *set, NPPSED *sed)
   1.644 +{     NPPROW *row;
   1.645 +      int x, y, z, s, c;
   1.646 +      /* the set should contain at least three literals */
   1.647 +      xassert(set != NULL);
   1.648 +      xassert(set->next != NULL);
   1.649 +      xassert(set->next->next != NULL);
   1.650 +      sed->x = set->lit;
   1.651 +      xassert(sed->x.neg == 0 || sed->x.neg == 1);
   1.652 +      sed->y = set->next->lit;
   1.653 +      xassert(sed->y.neg == 0 || sed->y.neg == 1);
   1.654 +      sed->z = set->next->next->lit;
   1.655 +      xassert(sed->z.neg == 0 || sed->z.neg == 1);
   1.656 +      /* perform encoding s = (x + y + z) % 2 */
   1.657 +      sed->s = npp_add_col(npp);
   1.658 +      sed->s->is_int = 1, sed->s->lb = 0.0, sed->s->ub = 1.0;
   1.659 +      for (x = 0; x <= 1; x++)
   1.660 +      {  for (y = 0; y <= 1; y++)
   1.661 +         {  for (z = 0; z <= 1; z++)
   1.662 +            {  for (s = 0; s <= 1; s++)
   1.663 +               {  if ((x + y + z) % 2 != s)
   1.664 +                  {  /* generate CNF clause to disable infeasible
   1.665 +                        combination */
   1.666 +                     row = npp_add_row(npp);
   1.667 +                     row->lb = 1.0, row->ub = +DBL_MAX;
   1.668 +                     if (x == sed->x.neg)
   1.669 +                        npp_add_aij(npp, row, sed->x.col, +1.0);
   1.670 +                     else
   1.671 +                     {  npp_add_aij(npp, row, sed->x.col, -1.0);
   1.672 +                        row->lb -= 1.0;
   1.673 +                     }
   1.674 +                     if (y == sed->y.neg)
   1.675 +                        npp_add_aij(npp, row, sed->y.col, +1.0);
   1.676 +                     else
   1.677 +                     {  npp_add_aij(npp, row, sed->y.col, -1.0);
   1.678 +                        row->lb -= 1.0;
   1.679 +                     }
   1.680 +                     if (z == sed->z.neg)
   1.681 +                        npp_add_aij(npp, row, sed->z.col, +1.0);
   1.682 +                     else
   1.683 +                     {  npp_add_aij(npp, row, sed->z.col, -1.0);
   1.684 +                        row->lb -= 1.0;
   1.685 +                     }
   1.686 +                     if (s == 0)
   1.687 +                        npp_add_aij(npp, row, sed->s, +1.0);
   1.688 +                     else
   1.689 +                     {  npp_add_aij(npp, row, sed->s, -1.0);
   1.690 +                        row->lb -= 1.0;
   1.691 +                     }
   1.692 +                  }
   1.693 +               }
   1.694 +            }
   1.695 +         }
   1.696 +      }
   1.697 +      /* perform encoding c = (x + y + z) / 2 */
   1.698 +      sed->c = npp_add_col(npp);
   1.699 +      sed->c->is_int = 1, sed->c->lb = 0.0, sed->c->ub = 1.0;
   1.700 +      for (x = 0; x <= 1; x++)
   1.701 +      {  for (y = 0; y <= 1; y++)
   1.702 +         {  for (z = 0; z <= 1; z++)
   1.703 +            {  for (c = 0; c <= 1; c++)
   1.704 +               {  if ((x + y + z) / 2 != c)
   1.705 +                  {  /* generate CNF clause to disable infeasible
   1.706 +                        combination */
   1.707 +                     row = npp_add_row(npp);
   1.708 +                     row->lb = 1.0, row->ub = +DBL_MAX;
   1.709 +                     if (x == sed->x.neg)
   1.710 +                        npp_add_aij(npp, row, sed->x.col, +1.0);
   1.711 +                     else
   1.712 +                     {  npp_add_aij(npp, row, sed->x.col, -1.0);
   1.713 +                        row->lb -= 1.0;
   1.714 +                     }
   1.715 +                     if (y == sed->y.neg)
   1.716 +                        npp_add_aij(npp, row, sed->y.col, +1.0);
   1.717 +                     else
   1.718 +                     {  npp_add_aij(npp, row, sed->y.col, -1.0);
   1.719 +                        row->lb -= 1.0;
   1.720 +                     }
   1.721 +                     if (z == sed->z.neg)
   1.722 +                        npp_add_aij(npp, row, sed->z.col, +1.0);
   1.723 +                     else
   1.724 +                     {  npp_add_aij(npp, row, sed->z.col, -1.0);
   1.725 +                        row->lb -= 1.0;
   1.726 +                     }
   1.727 +                     if (c == 0)
   1.728 +                        npp_add_aij(npp, row, sed->c, +1.0);
   1.729 +                     else
   1.730 +                     {  npp_add_aij(npp, row, sed->c, -1.0);
   1.731 +                        row->lb -= 1.0;
   1.732 +                     }
   1.733 +                  }
   1.734 +               }
   1.735 +            }
   1.736 +         }
   1.737 +      }
   1.738 +      return;
   1.739 +}
   1.740 +
   1.741 +/***********************************************************************
   1.742 +*  npp_sat_encode_sum_ax - encode linear combination of 0-1 variables
   1.743 +*
   1.744 +*  PURPOSE
   1.745 +*
   1.746 +*  Given a linear combination of binary variables:
   1.747 +*
   1.748 +*     sum a[j] x[j],                                                 (1)
   1.749 +*      j
   1.750 +*
   1.751 +*  which is the linear form of the specified row, this routine encodes
   1.752 +*  (i.e. translates to CNF) the following equality:
   1.753 +*
   1.754 +*                        n
   1.755 +*     sum |a[j]| t[j] = sum 2**(k-1) * y[k],                         (2)
   1.756 +*      j                k=1
   1.757 +*
   1.758 +*  where t[j] = x[j] (if a[j] > 0) or t[j] = 1 - x[j] (if a[j] < 0),
   1.759 +*  and y[k] is either t[j] or a new literal created by the routine or
   1.760 +*  a constant zero. Note that the sum in the right-hand side of (2) can
   1.761 +*  be thought as a n-bit representation of the sum in the left-hand
   1.762 +*  side, which is a non-negative integer number.
   1.763 +*
   1.764 +*  ALGORITHM
   1.765 +*
   1.766 +*  First, the number of bits, n, sufficient to represent any value in
   1.767 +*  the left-hand side of (2) is determined. Obviously, n is the number
   1.768 +*  of bits sufficient to represent the sum (sum |a[j]|).
   1.769 +*
   1.770 +*  Let
   1.771 +*
   1.772 +*               n
   1.773 +*     |a[j]| = sum 2**(k-1) b[j,k],                                  (3)
   1.774 +*              k=1
   1.775 +*
   1.776 +*  where b[j,k] is k-th bit in a n-bit representation of |a[j]|. Then
   1.777 +*
   1.778 +*                          m            n
   1.779 +*     sum |a[j]| * t[j] = sum 2**(k-1) sum b[j,k] * t[j].            (4)
   1.780 +*      j                  k=1          j=1
   1.781 +*
   1.782 +*  Introducing the set
   1.783 +*
   1.784 +*     J[k] = { j : b[j,k] = 1 }                                      (5)
   1.785 +*
   1.786 +*  allows rewriting (4) as follows:
   1.787 +*
   1.788 +*                          n
   1.789 +*     sum |a[j]| * t[j] = sum 2**(k-1)  sum    t[j].                 (6)
   1.790 +*      j                  k=1         j in J[k]
   1.791 +*
   1.792 +*  Thus, our goal is to provide |J[k]| <= 1 for all k, in which case
   1.793 +*  we will have the representation (1).
   1.794 +*
   1.795 +*  Let |J[k]| = 2, i.e. J[k] has exactly two literals u and v. In this
   1.796 +*  case we can apply the following transformation:
   1.797 +*
   1.798 +*     u + v = s + 2 * c,                                             (7)
   1.799 +*
   1.800 +*  where s and c are, respectively, low (sum) and high (carry) bits of
   1.801 +*  the sum of two bits. This allows to replace two literals u and v in
   1.802 +*  J[k] by one literal s, and carry out literal c to J[k+1].
   1.803 +*
   1.804 +*  If |J[k]| >= 3, i.e. J[k] has at least three literals u, v, and w,
   1.805 +*  we can apply the following transformation:
   1.806 +*
   1.807 +*     u + v + w = s + 2 * c.                                         (8)
   1.808 +*
   1.809 +*  Again, literal s replaces literals u, v, and w in J[k], and literal
   1.810 +*  c goes into J[k+1].
   1.811 +*
   1.812 +*  On exit the routine stores each literal from J[k] in element y[k],
   1.813 +*  1 <= k <= n. If J[k] is empty, y[k] is set to constant false.
   1.814 +*
   1.815 +*  RETURNS
   1.816 +*
   1.817 +*  The routine returns n, the number of literals in the right-hand side
   1.818 +*  of (2), 0 <= n <= NBIT_MAX. If the sum (sum |a[j]|) is too large, so
   1.819 +*  more than NBIT_MAX (= 31) literals are needed to encode the original
   1.820 +*  linear combination, the routine returns a negative value. */
   1.821 +
   1.822 +#define NBIT_MAX 31
   1.823 +/* maximal number of literals in the right hand-side of (2) */
   1.824 +
   1.825 +static NPPLSE *remove_lse(NPP *npp, NPPLSE *set, NPPCOL *col)
   1.826 +{     /* remove specified literal from specified literal set */
   1.827 +      NPPLSE *lse, *prev = NULL;
   1.828 +      for (lse = set; lse != NULL; prev = lse, lse = lse->next)
   1.829 +         if (lse->lit.col == col) break;
   1.830 +      xassert(lse != NULL);
   1.831 +      if (prev == NULL)
   1.832 +         set = lse->next;
   1.833 +      else
   1.834 +         prev->next = lse->next;
   1.835 +      dmp_free_atom(npp->pool, lse, sizeof(NPPLSE));
   1.836 +      return set;
   1.837 +}
   1.838 +
   1.839 +int npp_sat_encode_sum_ax(NPP *npp, NPPROW *row, NPPLIT y[])
   1.840 +{     NPPAIJ *aij;
   1.841 +      NPPLSE *set[1+NBIT_MAX], *lse;
   1.842 +      NPPSED sed;
   1.843 +      int k, n, temp;
   1.844 +      double sum;
   1.845 +      /* compute the sum (sum |a[j]|) */
   1.846 +      sum = 0.0;
   1.847 +      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
   1.848 +         sum += fabs(aij->val);
   1.849 +      /* determine n, the number of bits in the sum */
   1.850 +      temp = (int)sum;
   1.851 +      if ((double)temp != sum)
   1.852 +         return -1; /* integer arithmetic error */
   1.853 +      for (n = 0; temp > 0; n++, temp >>= 1);
   1.854 +      xassert(0 <= n && n <= NBIT_MAX);
   1.855 +      /* build initial sets J[k], 1 <= k <= n; see (5) */
   1.856 +      /* set[k] is a pointer to the list of literals in J[k] */
   1.857 +      for (k = 1; k <= n; k++)
   1.858 +         set[k] = NULL;
   1.859 +      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
   1.860 +      {  temp = (int)fabs(aij->val);
   1.861 +         xassert((int)temp == fabs(aij->val));
   1.862 +         for (k = 1; temp > 0; k++, temp >>= 1)
   1.863 +         {  if (temp & 1)
   1.864 +            {  xassert(k <= n);
   1.865 +               lse = dmp_get_atom(npp->pool, sizeof(NPPLSE));
   1.866 +               lse->lit.col = aij->col;
   1.867 +               lse->lit.neg = (aij->val > 0.0 ? 0 : 1);
   1.868 +               lse->next = set[k];
   1.869 +               set[k] = lse;
   1.870 +            }
   1.871 +         }
   1.872 +      }
   1.873 +      /* main transformation loop */
   1.874 +      for (k = 1; k <= n; k++)
   1.875 +      {  /* reduce J[k] and set y[k] */
   1.876 +         for (;;)
   1.877 +         {  if (set[k] == NULL)
   1.878 +            {  /* J[k] is empty */
   1.879 +               /* set y[k] to constant false */
   1.880 +               y[k].col = NULL, y[k].neg = 0;
   1.881 +               break;
   1.882 +            }
   1.883 +            if (set[k]->next == NULL)
   1.884 +            {  /* J[k] contains one literal */
   1.885 +               /* set y[k] to that literal */
   1.886 +               y[k] = set[k]->lit;
   1.887 +               dmp_free_atom(npp->pool, set[k], sizeof(NPPLSE));
   1.888 +               break;
   1.889 +            }
   1.890 +            if (set[k]->next->next == NULL)
   1.891 +            {  /* J[k] contains two literals */
   1.892 +               /* apply transformation (7) */
   1.893 +               npp_sat_encode_sum2(npp, set[k], &sed);
   1.894 +            }
   1.895 +            else
   1.896 +            {  /* J[k] contains at least three literals */
   1.897 +               /* apply transformation (8) */
   1.898 +               npp_sat_encode_sum3(npp, set[k], &sed);
   1.899 +               /* remove third literal from set[k] */
   1.900 +               set[k] = remove_lse(npp, set[k], sed.z.col);
   1.901 +            }
   1.902 +            /* remove second literal from set[k] */
   1.903 +            set[k] = remove_lse(npp, set[k], sed.y.col);
   1.904 +            /* remove first literal from set[k] */
   1.905 +            set[k] = remove_lse(npp, set[k], sed.x.col);
   1.906 +            /* include new literal s to set[k] */
   1.907 +            lse = dmp_get_atom(npp->pool, sizeof(NPPLSE));
   1.908 +            lse->lit.col = sed.s, lse->lit.neg = 0;
   1.909 +            lse->next = set[k];
   1.910 +            set[k] = lse;
   1.911 +            /* include new literal c to set[k+1] */
   1.912 +            xassert(k < n); /* FIXME: can "overflow" happen? */
   1.913 +            lse = dmp_get_atom(npp->pool, sizeof(NPPLSE));
   1.914 +            lse->lit.col = sed.c, lse->lit.neg = 0;
   1.915 +            lse->next = set[k+1];
   1.916 +            set[k+1] = lse;
   1.917 +         }
   1.918 +      }
   1.919 +      return n;
   1.920 +}
   1.921 +
   1.922 +/***********************************************************************
   1.923 +*  npp_sat_normalize_clause - normalize clause
   1.924 +*
   1.925 +*  This routine normalizes the specified clause, which is a disjunction
   1.926 +*  of literals, by replacing multiple literals, which refer to the same
   1.927 +*  binary variable, with a single literal.
   1.928 +*
   1.929 +*  On exit the routine returns the number of literals in the resulting
   1.930 +*  clause. However, if the specified clause includes both a literal and
   1.931 +*  its negation, the routine returns a negative value meaning that the
   1.932 +*  clause is equivalent to the value true. */
   1.933 +
   1.934 +int npp_sat_normalize_clause(NPP *npp, int size, NPPLIT lit[])
   1.935 +{     int j, k, new_size;
   1.936 +      xassert(npp == npp);
   1.937 +      xassert(size >= 0);
   1.938 +      new_size = 0;
   1.939 +      for (k = 1; k <= size; k++)
   1.940 +      {  for (j = 1; j <= new_size; j++)
   1.941 +         {  if (lit[k].col == lit[j].col)
   1.942 +            {  /* lit[k] refers to the same variable as lit[j], which
   1.943 +                  is already included in the resulting clause */
   1.944 +               if (lit[k].neg == lit[j].neg)
   1.945 +               {  /* ignore lit[k] due to the idempotent law */
   1.946 +                  goto skip;
   1.947 +               }
   1.948 +               else
   1.949 +               {  /* lit[k] is NOT lit[j]; the clause is equivalent to
   1.950 +                     the value true */
   1.951 +                  return -1;
   1.952 +               }
   1.953 +            }
   1.954 +         }
   1.955 +         /* include lit[k] in the resulting clause */
   1.956 +         lit[++new_size] = lit[k];
   1.957 +skip:    ;
   1.958 +      }
   1.959 +      return new_size;
   1.960 +}
   1.961 +
   1.962 +/***********************************************************************
   1.963 +*  npp_sat_encode_clause - translate clause to cover inequality
   1.964 +*
   1.965 +*  Given a clause
   1.966 +*
   1.967 +*     OR  t[j],                                                      (1)
   1.968 +*   j in J
   1.969 +*
   1.970 +*  where t[j] is a literal, i.e. t[j] = x[j] or t[j] = NOT x[j], this
   1.971 +*  routine translates it to the following equivalent cover inequality,
   1.972 +*  which is added to the transformed problem:
   1.973 +*
   1.974 +*     sum t[j] >= 1,                                                 (2)
   1.975 +*   j in J
   1.976 +*
   1.977 +*  where t[j] = x[j] or t[j] = 1 - x[j].
   1.978 +*
   1.979 +*  If necessary, the clause should be normalized before a call to this
   1.980 +*  routine. */
   1.981 +
   1.982 +NPPROW *npp_sat_encode_clause(NPP *npp, int size, NPPLIT lit[])
   1.983 +{     NPPROW *row;
   1.984 +      int k;
   1.985 +      xassert(size >= 1);
   1.986 +      row = npp_add_row(npp);
   1.987 +      row->lb = 1.0, row->ub = +DBL_MAX;
   1.988 +      for (k = 1; k <= size; k++)
   1.989 +      {  xassert(lit[k].col != NULL);
   1.990 +         if (lit[k].neg == 0)
   1.991 +            npp_add_aij(npp, row, lit[k].col, +1.0);
   1.992 +         else if (lit[k].neg == 1)
   1.993 +         {  npp_add_aij(npp, row, lit[k].col, -1.0);
   1.994 +            row->lb -= 1.0;
   1.995 +         }
   1.996 +         else
   1.997 +            xassert(lit != lit);
   1.998 +      }
   1.999 +      return row;
  1.1000 +}
  1.1001 +
  1.1002 +/***********************************************************************
  1.1003 +*  npp_sat_encode_geq - encode "not less than" constraint
  1.1004 +*
  1.1005 +*  PURPOSE
  1.1006 +*
  1.1007 +*  This routine translates to CNF the following constraint:
  1.1008 +*
  1.1009 +*      n
  1.1010 +*     sum 2**(k-1) * y[k] >= b,                                      (1)
  1.1011 +*     k=1
  1.1012 +*
  1.1013 +*  where y[k] is either a literal (i.e. y[k] = x[k] or y[k] = 1 - x[k])
  1.1014 +*  or constant false (zero), b is a given lower bound.
  1.1015 +*
  1.1016 +*  ALGORITHM
  1.1017 +*
  1.1018 +*  If b < 0, the constraint is redundant, so assume that b >= 0. Let
  1.1019 +*
  1.1020 +*          n
  1.1021 +*     b = sum 2**(k-1) b[k],                                         (2)
  1.1022 +*         k=1
  1.1023 +*
  1.1024 +*  where b[k] is k-th binary digit of b. (Note that if b >= 2**n and
  1.1025 +*  therefore cannot be represented in the form (2), the constraint (1)
  1.1026 +*  is infeasible.) In this case the condition (1) is equivalent to the
  1.1027 +*  following condition:
  1.1028 +*
  1.1029 +*     y[n] y[n-1] ... y[2] y[1] >= b[n] b[n-1] ... b[2] b[1],        (3)
  1.1030 +*
  1.1031 +*  where ">=" is understood lexicographically.
  1.1032 +*
  1.1033 +*  Algorithmically the condition (3) can be tested as follows:
  1.1034 +*
  1.1035 +*     for (k = n; k >= 1; k--)
  1.1036 +*     {  if (y[k] < b[k])
  1.1037 +*           y is less than b;
  1.1038 +*        if (y[k] > b[k])
  1.1039 +*           y is greater than b;
  1.1040 +*     }
  1.1041 +*     y is equal to b;
  1.1042 +*
  1.1043 +*  Thus, y is less than b iff there exists k, 1 <= k <= n, for which
  1.1044 +*  the following condition is satisfied:
  1.1045 +*
  1.1046 +*     y[n] = b[n] AND ... AND y[k+1] = b[k+1] AND y[k] < b[k].       (4)
  1.1047 +*
  1.1048 +*  Negating the condition (4) we have that y is not less than b iff for
  1.1049 +*  all k, 1 <= k <= n, the following condition is satisfied:
  1.1050 +*
  1.1051 +*     y[n] != b[n] OR ... OR y[k+1] != b[k+1] OR y[k] >= b[k].       (5)
  1.1052 +*
  1.1053 +*  Note that if b[k] = 0, the literal y[k] >= b[k] is always true, in
  1.1054 +*  which case the entire clause (5) is true and can be omitted.
  1.1055 +*
  1.1056 +*  RETURNS
  1.1057 +*
  1.1058 +*  Normally the routine returns zero. However, if the constraint (1) is
  1.1059 +*  infeasible, the routine returns non-zero. */
  1.1060 +
  1.1061 +int npp_sat_encode_geq(NPP *npp, int n, NPPLIT y[], int rhs)
  1.1062 +{     NPPLIT lit[1+NBIT_MAX];
  1.1063 +      int j, k, size, temp, b[1+NBIT_MAX];
  1.1064 +      xassert(0 <= n && n <= NBIT_MAX);
  1.1065 +      /* if the constraint (1) is redundant, do nothing */
  1.1066 +      if (rhs < 0)
  1.1067 +         return 0;
  1.1068 +      /* determine binary digits of b according to (2) */
  1.1069 +      for (k = 1, temp = rhs; k <= n; k++, temp >>= 1)
  1.1070 +         b[k] = temp & 1;
  1.1071 +      if (temp != 0)
  1.1072 +      {  /* b >= 2**n; the constraint (1) is infeasible */
  1.1073 +         return 1;
  1.1074 +      }
  1.1075 +      /* main transformation loop */
  1.1076 +      for (k = 1; k <= n; k++)
  1.1077 +      {  /* build the clause (5) for current k */
  1.1078 +         size = 0; /* clause size = number of literals */
  1.1079 +         /* add literal y[k] >= b[k] */
  1.1080 +         if (b[k] == 0)
  1.1081 +         {  /* b[k] = 0 -> the literal is true */
  1.1082 +            goto skip;
  1.1083 +         }
  1.1084 +         else if (y[k].col == NULL)
  1.1085 +         {  /* y[k] = 0, b[k] = 1 -> the literal is false */
  1.1086 +            xassert(y[k].neg == 0);
  1.1087 +         }
  1.1088 +         else
  1.1089 +         {  /* add literal y[k] = 1 */
  1.1090 +            lit[++size] = y[k];
  1.1091 +         }
  1.1092 +         for (j = k+1; j <= n; j++)
  1.1093 +         {  /* add literal y[j] != b[j] */
  1.1094 +            if (y[j].col == NULL)
  1.1095 +            {  xassert(y[j].neg == 0);
  1.1096 +               if (b[j] == 0)
  1.1097 +               {  /* y[j] = 0, b[j] = 0 -> the literal is false */
  1.1098 +                  continue;
  1.1099 +               }
  1.1100 +               else
  1.1101 +               {  /* y[j] = 0, b[j] = 1 -> the literal is true */
  1.1102 +                  goto skip;
  1.1103 +               }
  1.1104 +            }
  1.1105 +            else
  1.1106 +            {  lit[++size] = y[j];
  1.1107 +               if (b[j] != 0)
  1.1108 +                  lit[size].neg = 1 - lit[size].neg;
  1.1109 +            }
  1.1110 +         }
  1.1111 +         /* normalize the clause */
  1.1112 +         size = npp_sat_normalize_clause(npp, size, lit);
  1.1113 +         if (size < 0)
  1.1114 +         {  /* the clause is equivalent to the value true */
  1.1115 +            goto skip;
  1.1116 +         }
  1.1117 +         if (size == 0)
  1.1118 +         {  /* the clause is equivalent to the value false; this means
  1.1119 +               that the constraint (1) is infeasible */
  1.1120 +            return 2;
  1.1121 +         }
  1.1122 +         /* translate the clause to corresponding cover inequality */
  1.1123 +         npp_sat_encode_clause(npp, size, lit);
  1.1124 +skip:    ;
  1.1125 +      }
  1.1126 +      return 0;
  1.1127 +}
  1.1128 +
  1.1129 +/***********************************************************************
  1.1130 +*  npp_sat_encode_leq - encode "not greater than" constraint
  1.1131 +*
  1.1132 +*  PURPOSE
  1.1133 +*
  1.1134 +*  This routine translates to CNF the following constraint:
  1.1135 +*
  1.1136 +*      n
  1.1137 +*     sum 2**(k-1) * y[k] <= b,                                      (1)
  1.1138 +*     k=1
  1.1139 +*
  1.1140 +*  where y[k] is either a literal (i.e. y[k] = x[k] or y[k] = 1 - x[k])
  1.1141 +*  or constant false (zero), b is a given upper bound.
  1.1142 +*
  1.1143 +*  ALGORITHM
  1.1144 +*
  1.1145 +*  If b < 0, the constraint is infeasible, so assume that b >= 0. Let
  1.1146 +*
  1.1147 +*          n
  1.1148 +*     b = sum 2**(k-1) b[k],                                         (2)
  1.1149 +*         k=1
  1.1150 +*
  1.1151 +*  where b[k] is k-th binary digit of b. (Note that if b >= 2**n and
  1.1152 +*  therefore cannot be represented in the form (2), the constraint (1)
  1.1153 +*  is redundant.) In this case the condition (1) is equivalent to the
  1.1154 +*  following condition:
  1.1155 +*
  1.1156 +*     y[n] y[n-1] ... y[2] y[1] <= b[n] b[n-1] ... b[2] b[1],        (3)
  1.1157 +*
  1.1158 +*  where "<=" is understood lexicographically.
  1.1159 +*
  1.1160 +*  Algorithmically the condition (3) can be tested as follows:
  1.1161 +*
  1.1162 +*     for (k = n; k >= 1; k--)
  1.1163 +*     {  if (y[k] < b[k])
  1.1164 +*           y is less than b;
  1.1165 +*        if (y[k] > b[k])
  1.1166 +*           y is greater than b;
  1.1167 +*     }
  1.1168 +*     y is equal to b;
  1.1169 +*
  1.1170 +*  Thus, y is greater than b iff there exists k, 1 <= k <= n, for which
  1.1171 +*  the following condition is satisfied:
  1.1172 +*
  1.1173 +*     y[n] = b[n] AND ... AND y[k+1] = b[k+1] AND y[k] > b[k].       (4)
  1.1174 +*
  1.1175 +*  Negating the condition (4) we have that y is not greater than b iff
  1.1176 +*  for all k, 1 <= k <= n, the following condition is satisfied:
  1.1177 +*
  1.1178 +*     y[n] != b[n] OR ... OR y[k+1] != b[k+1] OR y[k] <= b[k].       (5)
  1.1179 +*
  1.1180 +*  Note that if b[k] = 1, the literal y[k] <= b[k] is always true, in
  1.1181 +*  which case the entire clause (5) is true and can be omitted.
  1.1182 +*
  1.1183 +*  RETURNS
  1.1184 +*
  1.1185 +*  Normally the routine returns zero. However, if the constraint (1) is
  1.1186 +*  infeasible, the routine returns non-zero. */
  1.1187 +
  1.1188 +int npp_sat_encode_leq(NPP *npp, int n, NPPLIT y[], int rhs)
  1.1189 +{     NPPLIT lit[1+NBIT_MAX];
  1.1190 +      int j, k, size, temp, b[1+NBIT_MAX];
  1.1191 +      xassert(0 <= n && n <= NBIT_MAX);
  1.1192 +      /* check if the constraint (1) is infeasible */
  1.1193 +      if (rhs < 0)
  1.1194 +         return 1;
  1.1195 +      /* determine binary digits of b according to (2) */
  1.1196 +      for (k = 1, temp = rhs; k <= n; k++, temp >>= 1)
  1.1197 +         b[k] = temp & 1;
  1.1198 +      if (temp != 0)
  1.1199 +      {  /* b >= 2**n; the constraint (1) is redundant */
  1.1200 +         return 0;
  1.1201 +      }
  1.1202 +      /* main transformation loop */
  1.1203 +      for (k = 1; k <= n; k++)
  1.1204 +      {  /* build the clause (5) for current k */
  1.1205 +         size = 0; /* clause size = number of literals */
  1.1206 +         /* add literal y[k] <= b[k] */
  1.1207 +         if (b[k] == 1)
  1.1208 +         {  /* b[k] = 1 -> the literal is true */
  1.1209 +            goto skip;
  1.1210 +         }
  1.1211 +         else if (y[k].col == NULL)
  1.1212 +         {  /* y[k] = 0, b[k] = 0 -> the literal is true */
  1.1213 +            xassert(y[k].neg == 0);
  1.1214 +            goto skip;
  1.1215 +         }
  1.1216 +         else
  1.1217 +         {  /* add literal y[k] = 0 */
  1.1218 +            lit[++size] = y[k];
  1.1219 +            lit[size].neg = 1 - lit[size].neg;
  1.1220 +         }
  1.1221 +         for (j = k+1; j <= n; j++)
  1.1222 +         {  /* add literal y[j] != b[j] */
  1.1223 +            if (y[j].col == NULL)
  1.1224 +            {  xassert(y[j].neg == 0);
  1.1225 +               if (b[j] == 0)
  1.1226 +               {  /* y[j] = 0, b[j] = 0 -> the literal is false */
  1.1227 +                  continue;
  1.1228 +               }
  1.1229 +               else
  1.1230 +               {  /* y[j] = 0, b[j] = 1 -> the literal is true */
  1.1231 +                  goto skip;
  1.1232 +               }
  1.1233 +            }
  1.1234 +            else
  1.1235 +            {  lit[++size] = y[j];
  1.1236 +               if (b[j] != 0)
  1.1237 +                  lit[size].neg = 1 - lit[size].neg;
  1.1238 +            }
  1.1239 +         }
  1.1240 +         /* normalize the clause */
  1.1241 +         size = npp_sat_normalize_clause(npp, size, lit);
  1.1242 +         if (size < 0)
  1.1243 +         {  /* the clause is equivalent to the value true */
  1.1244 +            goto skip;
  1.1245 +         }
  1.1246 +         if (size == 0)
  1.1247 +         {  /* the clause is equivalent to the value false; this means
  1.1248 +               that the constraint (1) is infeasible */
  1.1249 +            return 2;
  1.1250 +         }
  1.1251 +         /* translate the clause to corresponding cover inequality */
  1.1252 +         npp_sat_encode_clause(npp, size, lit);
  1.1253 +skip:    ;
  1.1254 +      }
  1.1255 +      return 0;
  1.1256 +}
  1.1257 +
  1.1258 +/***********************************************************************
  1.1259 +*  npp_sat_encode_row - encode constraint (row) of general type
  1.1260 +*
  1.1261 +*  PURPOSE
  1.1262 +*
  1.1263 +*  This routine translates to CNF the following constraint (row):
  1.1264 +*
  1.1265 +*     L <= sum a[j] x[j] <= U,                                       (1)
  1.1266 +*           j
  1.1267 +*
  1.1268 +*  where all x[j] are binary variables.
  1.1269 +*
  1.1270 +*  ALGORITHM
  1.1271 +*
  1.1272 +*  First, the routine performs substitution x[j] = t[j] for j in J+
  1.1273 +*  and x[j] = 1 - t[j] for j in J-, where J+ = { j : a[j] > 0 } and
  1.1274 +*  J- = { j : a[j] < 0 }. This gives:
  1.1275 +*
  1.1276 +*     L <=  sum  a[j] t[j] +   sum  a[j] (1 - t[j]) <= U  ==>
  1.1277 +*         j in J+            j in J-
  1.1278 +*
  1.1279 +*     L' <= sum |a[j]| t[j] <= U',                                   (2)
  1.1280 +*            j
  1.1281 +*
  1.1282 +*  where
  1.1283 +*
  1.1284 +*     L' = L -   sum  a[j],   U' = U -   sum  a[j].                  (3)
  1.1285 +*              j in J-                 j in J-
  1.1286 +*
  1.1287 +*  (Actually only new bounds L' and U' are computed.)
  1.1288 +*
  1.1289 +*  Then the routine translates to CNF the following equality:
  1.1290 +*
  1.1291 +*                        n
  1.1292 +*     sum |a[j]| t[j] = sum 2**(k-1) * y[k],                         (4)
  1.1293 +*      j                k=1
  1.1294 +*
  1.1295 +*  where y[k] is either some t[j] or a new literal or a constant zero
  1.1296 +*  (see the routine npp_sat_encode_sum_ax).
  1.1297 +*
  1.1298 +*  Finally, the routine translates to CNF the following conditions:
  1.1299 +*
  1.1300 +*      n
  1.1301 +*     sum 2**(k-1) * y[k] >= L'                                      (5)
  1.1302 +*     k=1
  1.1303 +*
  1.1304 +*  and
  1.1305 +*
  1.1306 +*      n
  1.1307 +*     sum 2**(k-1) * y[k] <= U'                                      (6)
  1.1308 +*     k=1
  1.1309 +*
  1.1310 +*  (see the routines npp_sat_encode_geq and npp_sat_encode_leq).
  1.1311 +*
  1.1312 +*  All resulting clauses are encoded as cover inequalities and included
  1.1313 +*  into the transformed problem.
  1.1314 +*
  1.1315 +*  Note that on exit the routine removes the specified constraint (row)
  1.1316 +*  from the original problem.
  1.1317 +*
  1.1318 +*  RETURNS
  1.1319 +*
  1.1320 +*  The routine returns one of the following codes:
  1.1321 +*
  1.1322 +*  0 - translation was successful;
  1.1323 +*  1 - constraint (1) was found infeasible;
  1.1324 +*  2 - integer arithmetic error occured. */
  1.1325 +
  1.1326 +int npp_sat_encode_row(NPP *npp, NPPROW *row)
  1.1327 +{     NPPAIJ *aij;
  1.1328 +      NPPLIT y[1+NBIT_MAX];
  1.1329 +      int n, rhs;
  1.1330 +      double lb, ub;
  1.1331 +      /* the row should not be free */
  1.1332 +      xassert(!(row->lb == -DBL_MAX && row->ub == +DBL_MAX));
  1.1333 +      /* compute new bounds L' and U' (3) */
  1.1334 +      lb = row->lb;
  1.1335 +      ub = row->ub;
  1.1336 +      for (aij = row->ptr; aij != NULL; aij = aij->r_next)
  1.1337 +      {  if (aij->val < 0.0)
  1.1338 +         {  if (lb != -DBL_MAX)
  1.1339 +               lb -= aij->val;
  1.1340 +            if (ub != -DBL_MAX)
  1.1341 +               ub -= aij->val;
  1.1342 +         }
  1.1343 +      }
  1.1344 +      /* encode the equality (4) */
  1.1345 +      n = npp_sat_encode_sum_ax(npp, row, y);
  1.1346 +      if (n < 0)
  1.1347 +         return 2; /* integer arithmetic error */
  1.1348 +      /* encode the condition (5) */
  1.1349 +      if (lb != -DBL_MAX)
  1.1350 +      {  rhs = (int)lb;
  1.1351 +         if ((double)rhs != lb)
  1.1352 +            return 2; /* integer arithmetic error */
  1.1353 +         if (npp_sat_encode_geq(npp, n, y, rhs) != 0)
  1.1354 +            return 1; /* original constraint is infeasible */
  1.1355 +      }
  1.1356 +      /* encode the condition (6) */
  1.1357 +      if (ub != +DBL_MAX)
  1.1358 +      {  rhs = (int)ub;
  1.1359 +         if ((double)rhs != ub)
  1.1360 +            return 2; /* integer arithmetic error */
  1.1361 +         if (npp_sat_encode_leq(npp, n, y, rhs) != 0)
  1.1362 +            return 1; /* original constraint is infeasible */
  1.1363 +      }
  1.1364 +      /* remove the specified row from the problem */
  1.1365 +      npp_del_row(npp, row);
  1.1366 +      return 0;
  1.1367 +}
  1.1368 +
  1.1369 +/***********************************************************************
  1.1370 +*  npp_sat_encode_prob - encode 0-1 feasibility problem
  1.1371 +*
  1.1372 +*  This routine translates the specified 0-1 feasibility problem to an
  1.1373 +*  equivalent SAT-CNF problem.
  1.1374 +*
  1.1375 +*  N.B. Currently this is a very crude implementation.
  1.1376 +*
  1.1377 +*  RETURNS
  1.1378 +*
  1.1379 +*  0           success;
  1.1380 +*
  1.1381 +*  GLP_ENOPFS  primal/integer infeasibility detected;
  1.1382 +*
  1.1383 +*  GLP_ERANGE  integer overflow occured. */
  1.1384 +
  1.1385 +int npp_sat_encode_prob(NPP *npp)
  1.1386 +{     NPPROW *row, *next_row, *prev_row;
  1.1387 +      NPPCOL *col, *next_col;
  1.1388 +      int cover = 0, pack = 0, partn = 0, ret;
  1.1389 +      /* process and remove free rows */
  1.1390 +      for (row = npp->r_head; row != NULL; row = next_row)
  1.1391 +      {  next_row = row->next;
  1.1392 +         if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
  1.1393 +            npp_sat_free_row(npp, row);
  1.1394 +      }
  1.1395 +      /* process and remove fixed columns */
  1.1396 +      for (col = npp->c_head; col != NULL; col = next_col)
  1.1397 +      {  next_col = col->next;
  1.1398 +         if (col->lb == col->ub)
  1.1399 +            xassert(npp_sat_fixed_col(npp, col) == 0);
  1.1400 +      }
  1.1401 +      /* only binary variables should remain */
  1.1402 +      for (col = npp->c_head; col != NULL; col = col->next)
  1.1403 +         xassert(col->is_int && col->lb == 0.0 && col->ub == 1.0);
  1.1404 +      /* new rows may be added to the end of the row list, so we walk
  1.1405 +         from the end to beginning of the list */
  1.1406 +      for (row = npp->r_tail; row != NULL; row = prev_row)
  1.1407 +      {  prev_row = row->prev;
  1.1408 +         /* process special cases */
  1.1409 +         ret = npp_sat_is_cover_ineq(npp, row);
  1.1410 +         if (ret != 0)
  1.1411 +         {  /* row is covering inequality */
  1.1412 +            cover++;
  1.1413 +            /* since it already encodes a clause, just transform it to
  1.1414 +               canonical form */
  1.1415 +            if (ret == 2)
  1.1416 +            {  xassert(npp_sat_reverse_row(npp, row) == 0);
  1.1417 +               ret = npp_sat_is_cover_ineq(npp, row);
  1.1418 +            }
  1.1419 +            xassert(ret == 1);
  1.1420 +            continue;
  1.1421 +         }
  1.1422 +         ret = npp_sat_is_partn_eq(npp, row);
  1.1423 +         if (ret != 0)
  1.1424 +         {  /* row is partitioning equality */
  1.1425 +            NPPROW *cov;
  1.1426 +            NPPAIJ *aij;
  1.1427 +            partn++;
  1.1428 +            /* transform it to canonical form */
  1.1429 +            if (ret == 2)
  1.1430 +            {  xassert(npp_sat_reverse_row(npp, row) == 0);
  1.1431 +               ret = npp_sat_is_partn_eq(npp, row);
  1.1432 +            }
  1.1433 +            xassert(ret == 1);
  1.1434 +            /* and split it into covering and packing inequalities,
  1.1435 +               both in canonical forms */
  1.1436 +            cov = npp_add_row(npp);
  1.1437 +            cov->lb = row->lb, cov->ub = +DBL_MAX;
  1.1438 +            for (aij = row->ptr; aij != NULL; aij = aij->r_next)
  1.1439 +               npp_add_aij(npp, cov, aij->col, aij->val);
  1.1440 +            xassert(npp_sat_is_cover_ineq(npp, cov) == 1);
  1.1441 +            /* the cover inequality already encodes a clause and do
  1.1442 +               not need any further processing */
  1.1443 +            row->lb = -DBL_MAX;
  1.1444 +            xassert(npp_sat_is_pack_ineq(npp, row) == 1);
  1.1445 +            /* the packing inequality will be processed below */
  1.1446 +            pack--;
  1.1447 +         }
  1.1448 +         ret = npp_sat_is_pack_ineq(npp, row);
  1.1449 +         if (ret != 0)
  1.1450 +         {  /* row is packing inequality */
  1.1451 +            NPPROW *rrr;
  1.1452 +            int nlit, desired_nlit = 4;
  1.1453 +            pack++;
  1.1454 +            /* transform it to canonical form */
  1.1455 +            if (ret == 2)
  1.1456 +            {  xassert(npp_sat_reverse_row(npp, row) == 0);
  1.1457 +               ret = npp_sat_is_pack_ineq(npp, row);
  1.1458 +            }
  1.1459 +            xassert(ret == 1);
  1.1460 +            /* process the packing inequality */
  1.1461 +            for (;;)
  1.1462 +            {  /* determine the number of literals in the remaining
  1.1463 +                  inequality */
  1.1464 +               nlit = npp_row_nnz(npp, row);
  1.1465 +               if (nlit <= desired_nlit)
  1.1466 +                  break;
  1.1467 +               /* split the current inequality into one having not more
  1.1468 +                  than desired_nlit literals and remaining one */
  1.1469 +               rrr = npp_sat_split_pack(npp, row, desired_nlit-1);
  1.1470 +               /* translate the former inequality to CNF and remove it
  1.1471 +                  from the original problem */
  1.1472 +               npp_sat_encode_pack(npp, rrr);
  1.1473 +            }
  1.1474 +            /* translate the remaining inequality to CNF and remove it
  1.1475 +               from the original problem */
  1.1476 +            npp_sat_encode_pack(npp, row);
  1.1477 +            continue;
  1.1478 +         }
  1.1479 +         /* translate row of general type to CNF and remove it from the
  1.1480 +            original problem */
  1.1481 +         ret = npp_sat_encode_row(npp, row);
  1.1482 +         if (ret == 0)
  1.1483 +            ;
  1.1484 +         else if (ret == 1)
  1.1485 +            ret = GLP_ENOPFS;
  1.1486 +         else if (ret == 2)
  1.1487 +            ret = GLP_ERANGE;
  1.1488 +         else
  1.1489 +            xassert(ret != ret);
  1.1490 +         if (ret != 0)
  1.1491 +            goto done;
  1.1492 +      }
  1.1493 +      ret = 0;
  1.1494 +      if (cover != 0)
  1.1495 +         xprintf("%d covering inequalities\n", cover);
  1.1496 +      if (pack != 0)
  1.1497 +         xprintf("%d packing inequalities\n", pack);
  1.1498 +      if (partn != 0)
  1.1499 +         xprintf("%d partitioning equalities\n", partn);
  1.1500 +done: return ret;
  1.1501 +}
  1.1502 +
  1.1503 +/* eof */