source/integrer.c

   1 /******************************************************************************
   2  *                     PIP : Parametric Integer Programming                   *
   3  ******************************************************************************
   4  *                                                                            *
   5  * Copyright Paul Feautrier, 1988, 1993, 1994, 1996                           *
   6  *                                                                            *
   7  * This is free software; you can redistribute it and/or modify it under the  *
   8  * terms of the GNU General Public License as published by the Free Software  *
   9  * Foundation; either version 2 of the License, or (at your option) any later *
  10  * version.                                                                   *
  11  *                                                                            *
  12  * This software is distributed in the hope that it will be useful, but       *
  13  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
  14  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
  15  * for more details.                                                          *
  16  *                                                                            *
  17  * You should have received a copy of the GNU General Public License along    *
  18  * with software; if not, write to the Free Software Foundation, Inc.,        *
  19  * 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA                     *
  20  *                                                                            *
  21  * Written by Paul Feautrier                                                  *
  22  *                                                                            *
  23  ******************************************************************************/
  24
  25 #include <stdio.h>
  26 #include <stdlib.h>
  27
  28 #include <piplib/piplib.h>
  29
  30 /*  The routines in this file are used to build a Gomory cut from
  31     a non-integral row of the problem tableau                             */
  32
  33 extern long int cross_product, limit;
  34 extern int verbose;
  35 extern FILE * dump;
  36
  37 /* mod(x,y) computes the remainder of x when divided by y. The difference
  38    with x%y is that the result is guaranteed to be positive, which is not
  39    always true for x%y                                                    */
  40
  41 Entier mod(x, y)
  42 Entier x, y;
  43 {Entier r;
  44  r = x % y;
  45  if(r<0) r += y;
  46  return(r);
  47 }
  48
  49 /* this routine is useless at present                                     */
  50
  51 int non_borne(tp, nvar, D, bigparm)
  52 Tableau *tp;
  53 int nvar, bigparm;
  54 Entier D;
  55 {int i, ff;
  56  for(i = 0; i<nvar; i++)
  57      {ff = Flag(tp, i);
  58       if(bigparm > 0)
  59          {if(ff & Unit)return(True);
  60           if(Index(tp, i, bigparm) != D) return(True);
  61          }
  62       }
  63  return(False);
  64 }
  65
  66 Tableau *expanser();
  67
  68 /* integrer(.....) add a cut to the problem tableau, or return 0 when an
  69    integral solution has been found, or -1 when no integral solution
  70    exists.
  71
  72    Since integrer may add rows and columns to the problem tableau, its
  73    arguments are pointers rather than values. If a cut is constructed,
  74    ni increases by 1. If the cut is parametric, nparm increases by 1 and
  75    nc increases by 2.
  76                                                                          */
  77
  78 int integrer(ptp, pcontext, D, pnvar, pnparm, pni, pnc)
  79 Tableau **ptp, **pcontext;
  80 int *pnvar, *pnparm, *pni, *pnc;
  81 Entier D;
  82 {int ncol = *pnvar+*pnparm+1;
  83  int nligne = *pnvar + *pni;
  84  int nparm = *pnparm;
  85  int nvar = *pnvar;
  86  int ni = *pni;
  87  int nc = *pnc;
  88  Entier coupure[MAXCOL];
  89  int i, j, k, ff;
  90  Entier x, d;
  91  int ok_var, ok_const, ok_parm;
  92  Entier discrp[MAXPARM], discrm[MAXPARM];
  93  int llog();
  94
  95  if(ncol+1 >= MAXCOL) {
  96       fprintf(stderr, "Too much variables : %d\n", ncol);
  97       exit(3);
  98       }
  99
 100 /* search for a non-integral row */
 101  for(i = 0; i<nvar; i++) {
 102       D = Denom(*ptp,i);
 103       if(D == 1) continue;
 104 /*                          If the common denominator of the row is 1
 105                             the row is integral                         */
 106       ff = Flag(*ptp, i);
 107       if(ff & Unit)continue;
 108 /*                          If the row is a Unit, it is integral        */
 109
 110 /*                          Here a portential candidate has been found.
 111                             Build the cut by reducing each coefficient
 112                             modulo D, the common denominator            */
 113       ok_var = False;
 114       for(j = 0; j<nvar; j++) {
 115           x = coupure[j] = mod(Index(*ptp, i, j), D);
 116           if(x > 0) ok_var = True;
 117           }
 118 /*                          Done for the coefficient of the variables.  */
 119
 120       x = coupure[nvar] = - mod(-Index(*ptp, i, nvar), D);
 121       ok_const = (x != 0);
 122 /*                          This is the constant term                   */
 123       ok_parm = False;
 124       for(j = nvar+1; j<ncol; j++) {
 125           x = coupure[j] = - mod(- Index(*ptp, i, j), D);          /* (1) */
 126           if(x != 0) ok_parm = True;
 127           }
 128 /*                          These are the parametric terms              */
 129
 130       coupure[ncol] = D;    /* Just in case                             */
 131
 132 /* The question now is whether the cut is valid. The answer is given
 133 by the following decision table:
 134
 135 ok_var   ok_parm   ok_const
 136
 137   F        F         F       (a) continue, integral row
 138   F        F         T       (b) return -1, no solution
 139   F        T         F
 140                              (c) if the <<constant>> part is not divisible
 141                              by D then bottom else ....
 142   F        T         T
 143   T        F         F       (a) continue, integral row
 144   T        F         T       (d) constant cut
 145   T        T         F
 146                              (e) parametric cut
 147   T        T         T
 148
 149                                                                 case (a)  */
 150
 151       if(!ok_parm && !ok_const) continue;
 152       if(!ok_parm)
 153           if(ok_var) {                                   /*     case (d)  */
 154               if(nligne >= (*ptp)->height) {
 155                   int d, dth, dtw;
 156                   d = llog(D);
 157                   dth = d;
 158                   *ptp = expanser(*ptp, nvar, ni, ncol, 0, dth, 0);
 159                   }
 160                          /* The cut has a negative <<constant>> part      */
 161               Flag(*ptp, nligne) = Minus;
 162               Denom(*ptp, nligne) = D;
 163                          /* Insert the cut */
 164               for(j = 0; j<ncol; j++)
 165                   Index(*ptp, nligne, j) = coupure[j];
 166                       /* A new row has been added to the problem tableau. */
 167               (*pni)++;
 168               return(nligne);
 169               }
 170           else return -1;                                /*   case (b)    */
 171
 172 /* In cases (c) and (e), one has to introduce a new parameter and
 173    introduce its defining inequalities into the context.
 174
 175    Let the cut be    sum_{j=0}^{nvar-1} c_j x_j + c_{nvar} +             (2)
 176                      sum_{j=0}^{nparm-1} c_{nvar + 1 + j} p_j >= 0.       */
 177
 178
 179       if(nparm >= MAXPARM) {
 180           fprintf(stderr, "Too much parameters : %d\n", *pnparm);
 181           exit(4);
 182           }
 183 /*        Build the definition of the new parameter into the solution :
 184       p_{nparm} = -(sum_{j=0}^{nparm-1} c_{nvar + 1 + j} p_j
 185                      + c_{nvar})/D                             (3)
 186          The minus sign is there to compensate the one in (1)     */
 187
 188       sol_new(nparm);
 189       sol_div();
 190       sol_forme(nparm+1);
 191       for(j = 0; j<nparm; j++) sol_val(-coupure[j+nvar+1], UN);
 192       sol_val(-coupure[*pnvar], UN);
 193       sol_val(D, UN);                     /* The divisor                */
 194
 195 /* The value of the new parameter is specified by applying the definition of
 196    Euclidean division to (3) :
 197
 198  0<= - sum_{j=0}^{nparm-1} c_{nvar+1+j} p_j - c_{nvar} - D * p_{nparm} < D (4)
 199
 200    This formula gives two inequalities which are stored in the context    */
 201
 202       for(j = 0; j<nparm; j++) {
 203           x = coupure[j+nvar+1];
 204           discrp[j] = -x;
 205           discrm[j] = x;
 206           }
 207       discrp[nparm] = -D;
 208       discrm[nparm] = D;
 209       x = coupure[nvar];
 210       discrp[(nparm)+1] = -x;
 211       discrm[(nparm)+1] = x + D -1;
 212       if(nc+2 > (*pcontext)->height || nparm+1 > (*pcontext)->width) {
 213           int dcw, dch;
 214           dcw = llog(D);
 215           dch = 2 * dcw + *pni;
 216           *pcontext = expanser(*pcontext, 0, nc, nparm+1, 0, dch, dcw);
 217           }
 218       /* Flag(*pcontext, *pnc) = 0; Probably useless see line A */
 219
 220 /* Since a new parameter is to be added, the constant term has to be moved
 221    right and a zero has to be inserted in all rows of the old context    */
 222
 223       for(k = 0; k < nc; k++) {
 224           Index(*pcontext, k, nparm+1) = Index(*pcontext, k, nparm);
 225           Index(*pcontext, k, nparm) = 0;
 226           }
 227 /* Now, insert the new rows                                              */
 228
 229       for(j = 0; j <= nparm+1; j++) {
 230           Index(*pcontext, nc, j) = discrp[j];
 231           Index(*pcontext, nc+1, j) = discrm[j];
 232           }
 233       Flag(*pcontext, nc) = Unknown;                                /* A */
 234       Denom(*pcontext, nc) = UN;
 235       Flag(*pcontext, nc+1) = Unknown;
 236       Denom(*pcontext, nc+1) = UN;
 237       (*pnparm)++;
 238       (*pnc) += 2;
 239                          /* end of the construction of the new parameter */
 240
 241       if(ok_var) {                                 /*   case (e)         */
 242           if(nligne >= (*ptp)->height || ncol >= (*ptp)->width) {
 243               int d, dth, dtw;
 244               d = llog(D);
 245               dth = d + ni;
 246               dtw = d;
 247               *ptp = expanser(*ptp, nvar, ni, ncol, 0, dth, dtw);
 248               }
 249                          /* Zeroing out the new column seems to be useless
 250                             since <<expanser>> does it anyway            */
 251
 252                          /* The cut has a negative <<constant>> part    */
 253           Flag(*ptp, nligne) = Minus;
 254           Denom(*ptp, nligne) = D;
 255                  /* Insert the cut */
 256           for(j = 0; j<ncol+1; j++)
 257               Index(*ptp, nligne, j) = coupure[j];
 258                  /* A new row has been added to the problem tableau.    */
 259           (*pni)++;
 260           return(nligne);
 261           }
 262                                                   /*  case (c)          */
 263                         /* The new parameter has already been defined as a
 264                            quotient. It remains to express that the
 265                            remainder of that division is zero           */
 266       sol_if();
 267       sol_forme(nparm + 2);
 268       for (j = 0; j < nparm+1 ; j++)
 269           sol_val(discrm[j], UN);
 270       sol_val(-UN, UN);
 271       sol_nil();    /* No solution if the division is not even      */
 272                         /* Add a new column */
 273       if(ncol+1 >= (*ptp)-> width) {
 274           int dtw;
 275           dtw = llog(D);
 276           *ptp = expanser(*ptp, *pnvar, *pni, ncol, 0, 0, dtw);
 277           }
 278           /* The new column is zeroed out by <<expanser>>          */
 279 /* Let c be the coefficient of parameter p in the i row. In <<coupure>>,
 280    this parameter has coefficient  - mod(-c, D). In <<discrp>>, this same
 281    parameter has coefficient mod(-c, D). The sum c + mod(-c, D) is obviously
 282    divisible by D.                                                      */
 283
 284       for (j = 0; j <= nparm; j++)
 285           Index(*ptp, i, j + nvar + 1) += discrp[j];
 286       tab_display(*ptp, stderr);
 287       exit(0);
 288       continue;
 289       }
 290  /* The solution is integral.                              */
 291  return 0;
 292 }
 293