source/integrerMP.c

   1 /********************************************************/
   2 /* Part of MultiPrecision PIP port by Zbigniew Chamski  */
   3 /* and Paul Feautrier.                                  */
   4 /* Based on straight PIP E.1 version by Paul Feautrier  */
   5 /* <Paul.Feautrier@inria.fr>                            */
   6 /*                                                      */
   7 /* and a previous port (C) Zbigniew CHAMSKI, 1993.      */
   8 /* <Zbigniew.Chamski@philips.com>                       */
   9 /*                                                      */
  10 /* Copying subject to the terms and conditions of the   */
  11 /* GNU General Public License.                          */
  12 /*                                                      */
  13 /* Send questions, bug reports and fixes to:            */
  14 /* <Paul.Feautrier@inria.fr>                            */
  15 /********************************************************/
  16
  17 #include <stdlib.h>
  18 #include <stdio.h>
  19
  20 #include <piplib/piplib.h>
  21
  22 /*  The routines in this file are used to build a Gomory cut from
  23     a non-integral row of the problem tableau                             */
  24
  25 extern long int cross_product, limit;
  26 extern int verbose;
  27 extern FILE * dump;
  28
  29 /* mod(x,y) computes the remainder of x when divided by y. The difference
  30    with x%y is that the result is guaranteed to be positive, which is not
  31    always true for x%y
  32
  33 This function is replaced by mpz_fdiv_r.
  34 */
  35
  36
  37 Tableau *expanser();
  38
  39 /* integrer(.....) add a cut to the problem tableau, or return 0 when an
  40    integral solution has been found, or -1 when no integral solution
  41    exists.
  42
  43    Since integrer may add rows and columns to the problem tableau, its
  44    arguments are pointers rather than values. If a cut is constructed,
  45    ni increases by 1. If the cut is parametric, nparm increases by 1 and
  46    nc increases by 2.
  47                                                                          */
  48
  49 int integrer(ptp, pcontext, pnvar, pnparm, pni, pnc)
  50 Tableau **ptp, **pcontext;
  51 int *pnvar, *pnparm, *pni, *pnc;
  52
  53 {int ncol = *pnvar+*pnparm+1;
  54  int nligne = *pnvar + *pni;
  55  int nparm = *pnparm;
  56  int nvar = *pnvar;
  57  int ni = *pni;
  58  int nc = *pnc;
  59  Entier coupure[MAXCOL];
  60  int i, j, k, ff;
  61  Entier x, d, D;
  62  int ok_var, ok_const, ok_parm;
  63  Entier discrp[MAXPARM], discrm[MAXPARM];
  64
  65
  66  for(i=0; i<=ncol; i++)
  67    mpz_init(coupure[i]);
  68
  69  for(i=0; i<=nparm+1; i++){
  70    mpz_init(discrp[i]);
  71    mpz_init(discrm[i]);
  72  }
  73
  74  mpz_init(x); mpz_init(d); mpz_init(D);
  75
  76  if(ncol+1 >= MAXCOL) {
  77       fprintf(stdout, "Too much variables : %d\n", ncol);
  78       exit(3);
  79       }
  80
  81 /* search for a non-integral row */
  82
  83  for(i = 0; i<nvar; i++) {
  84    ff = Flag(*ptp, i);
  85    if(ff & Unit)continue;
  86 /* If the row is a Unit, it is integral        */
  87
  88 /* D = Denom(*ptp,i); */
  89    mpz_set(D, Denom(*ptp, i));
  90
  91 /* if(D == 1) continue; */
  92    if(mpz_cmp_ui(D, 1) == 0) continue;
  93
  94 /*                          If the common denominator of the row is 1
  95                             the row is integral                         */
  96
  97 /*                          Here a potential candidate has been found.
  98                             Build the cut by reducing each coefficient
  99                             modulo D, the common denominator            */
 100    ok_var = False;
 101    for(j = 0; j<nvar; j++) {
 102      /* x = coupure[j] = mod(Index(*ptp, i, j), D); */
 103      mpz_fdiv_r(x, Index(*ptp, i, j), D);
 104      mpz_set(coupure[j], x);
 105
 106      if(x > 0) ok_var = True;
 107    }
 108 /*                          Done for the coefficient of the variables.  */
 109
 110    /* x = coupure[nvar] = - mod(-Index(*ptp, i, nvar), D); */
 111    mpz_neg(x, Index(*ptp, i, nvar));
 112    mpz_fdiv_r(x, x, D);
 113    mpz_neg(x, x);
 114    mpz_set(coupure[nvar], x);
 115    /* ok_const = (x != 0); */
 116    ok_const = mpz_cmp_ui(x, 0);
 117 /*                          This is the constant term                   */
 118    ok_parm = False;
 119    for(j = nvar+1; j<ncol; j++) {
 120      /* x = coupure[j] = - mod(- Index(*ptp, i, j), D); */       /* (1) */
 121      mpz_neg(x, Index(*ptp, i, j));
 122      mpz_fdiv_r(x, x, D);
 123      mpz_neg(x, x);
 124      mpz_set(coupure[j], x);
 125      /* if(x != 0) ok_parm = True; */
 126      if(mpz_cmp_ui(x, 0) != 0) ok_parm = True;
 127    }
 128 /*                          These are the parametric terms              */
 129
 130    /* coupure[ncol] = D;     Just in case                             */
 131    mpz_set(coupure[ncol], D);
 132
 133 /* The question now is whether the cut is valid. The answer is given
 134 by the following decision table:
 135
 136 ok_var   ok_parm   ok_const
 137
 138   F        F         F       (a) continue, integral row
 139   F        F         T       (b) return -1, no solution
 140   F        T         F
 141                              (c) if the <<constant>> part is not divisible
 142                              by D then bottom else ....
 143   F        T         T
 144   T        F         F       (a) continue, integral row
 145   T        F         T       (d) constant cut
 146   T        T         F
 147                              (e) parametric cut
 148   T        T         T
 149
 150                                                                 case (a)  */
 151
 152    if(!ok_parm && !ok_const) continue;
 153    if(!ok_parm)
 154      if(ok_var) {                                   /*     case (d)  */
 155        if(nligne >= (*ptp)->height) {
 156          int d, dth, dtw;
 157          d = mpz_sizeinbase(D, 2);
 158          dth = d;
 159          *ptp = expanser(*ptp, nvar, ni, ncol, 0, dth, 0);
 160        }
 161        /* The cut has a negative <<constant>> part      */
 162        Flag(*ptp, nligne) = Minus;
 163        /* Denom(*ptp, nligne) = D; */
 164        mpz_set(Denom(*ptp, nligne), D);
 165        /* Insert the cut */
 166        for(j = 0; j<ncol; j++)
 167          /* Index(*ptp, nligne, j) = coupure[j]; */
 168          mpz_set(Index(*ptp, nligne, j), coupure[j]);
 169                       /* A new row has been added to the problem tableau. */
 170        (*pni)++;
 171        /*       return(nligne);   */
 172        goto clear;
 173      }
 174    /* else return -1;  */
 175       else {
 176             nligne = -1;                              /*   case (b)    */
 177             goto clear;
 178       }
 179 /* In cases (c) and (e), one has to introduce a new parameter and
 180    introduce its defining inequalities into the context.
 181
 182    Let the cut be    sum_{j=0}^{nvar-1} c_j x_j + c_{nvar} +             (2)
 183                      sum_{j=0}^{nparm-1} c_{nvar + 1 + j} p_j >= 0.       */
 184
 185
 186    if(nparm >= MAXPARM) {
 187      fprintf(stdout, "Too much parameters : %d\n", *pnparm);
 188      exit(4);
 189    }
 190 /*        Build the definition of the new parameter into the solution :
 191           p_{nparm} = -(sum_{j=0}^{nparm-1} c_{nvar + 1 + j} p_j
 192                      + c_{nvar})/D                             (3)
 193          The minus sign is there to compensate the one in (1)     */
 194
 195
 196    sol_new(nparm);
 197    sol_div();
 198    sol_forme(nparm+1);
 199    for(j = 0; j<nparm; j++) {
 200      /* sol_val(-coupure[j+nvar+1], UN); */
 201      mpz_neg(x, coupure[j+nvar+1]);
 202      sol_val(x, UN);
 203    }
 204    /* sol_val(-coupure[*pnvar], UN); */
 205    mpz_neg(x, coupure[*pnvar]);
 206    sol_val(x, UN);
 207    sol_val(D, UN);                        /* The divisor                */
 208    if(verbose>0) fflush(dump);
 209
 210 /* The value of the new parameter is specified by applying the definition of
 211    Euclidean division to (3) :
 212
 213  0<= - sum_{j=0}^{nparm-1} c_{nvar+1+j} p_j - c_{nvar} - D * p_{nparm} < D (4)
 214
 215    This formula gives two inequalities which are stored in the context    */
 216
 217    for(j = 0; j<nparm; j++) {
 218      /* x = coupure[j+nvar+1]; */
 219      mpz_set(x, coupure[j+nvar+1]);
 220      /* discrp[j] = -x; */
 221      mpz_neg(discrp[j], x);
 222      /* discrm[j] = x; */
 223      mpz_set(discrm[j], x);
 224    }
 225    /* discrp[nparm] = -D; */
 226    mpz_neg(discrp[nparm], D);
 227    /* discrm[nparm] = D; */
 228    mpz_set(discrm[nparm], D);
 229    /* x = coupure[nvar]; */
 230    mpz_set(x, coupure[nvar]);
 231    if(verbose > 0){
 232      fprintf(dump, "Take care %d %d\n", nparm, nvar);
 233      fflush(dump);
 234    }
 235    /* discrp[(nparm)+1] = -x; */
 236    mpz_neg(discrp[nparm+1], x);
 237    /* discrm[(nparm)+1] = x + D -1; */
 238    mpz_sub_ui(x, x, 1);
 239    mpz_add(discrm[nparm+1], x, D);
 240
 241    if(nc+2 > (*pcontext)->height || nparm+1 > (*pcontext)->width) {
 242      int dcw, dch;
 243      dcw = mpz_sizeinbase(D, 2);
 244      dch = 2 * dcw + *pni;
 245      *pcontext = expanser(*pcontext, 0, nc, nparm+1, 0, dch, dcw);
 246    }
 247    /*   Flag(*pcontext, *pnc) = 0; Probably useless see line A */
 248 /* Since a new parameter is to be added, the constant term has to be moved
 249    right and a zero has to be inserted in all rows of the old context    */
 250
 251    for(k = 0; k < nc; k++) {
 252      /* Index(*pcontext, k, nparm+1) = Index(*pcontext, k, nparm); */
 253      mpz_set(Index(*pcontext, k, nparm+1), Index(*pcontext, k, nparm));
 254      /* Index(*pcontext, k, nparm) = 0; */
 255      mpz_set_ui(Index(*pcontext, k, nparm), 0);
 256    }
 257 /* Now, insert the new rows                                              */
 258
 259    for(j = 0; j <= nparm+1; j++) {
 260      /* Index(*pcontext, nc, j) = discrp[j]; */
 261      mpz_set(Index(*pcontext, nc, j), discrp[j]);
 262      /* Index(*pcontext, nc+1, j) = discrm[j]; */
 263      mpz_set(Index(*pcontext, nc+1, j), discrm[j]);
 264    }
 265    Flag(*pcontext, nc) = Unknown;                            /* A */
 266    /* Denom(*pcontext, nc) = UN; */
 267    mpz_set(Denom(*pcontext, nc), UN);
 268    Flag(*pcontext, nc+1) = Unknown;
 269    /* Denom(*pcontext, nc+1) = UN; */
 270    mpz_set(Denom(*pcontext, nc+1), UN);
 271    (*pnparm)++;
 272    (*pnc) += 2;
 273    if(verbose > 0){
 274      fprintf(dump, "enlarged context %d x %d\n", *pnparm, *pnc);
 275      fflush(dump);
 276    }
 277                          /* end of the construction of the new parameter */
 278
 279    if(ok_var) {                                 /*   case (e)         */
 280      if(nligne >= (*ptp)->height || ncol >= (*ptp)->width) {
 281        int d, dth, dtw;
 282        d = mpz_sizeinbase(D, 2);
 283        dth = d + ni;
 284        dtw = d;
 285        *ptp = expanser(*ptp, nvar, ni, ncol, 0, dth, dtw);
 286      }
 287                          /* Zeroing out the new column seems to be useless
 288                             since <<expanser>> does it anyway            */
 289
 290                          /* The cut has a negative <<constant>> part    */
 291      Flag(*ptp, nligne) = Minus;
 292      /* Denom(*ptp, nligne) = D; */
 293      mpz_set(Denom(*ptp, nligne), D);
 294                          /* Insert the cut */
 295      for(j = 0; j<ncol+1; j++)
 296        /* Index(*ptp, nligne, j) = coupure[j]; */
 297        mpz_set(Index(*ptp, nligne, j), coupure[j]);
 298                  /* A new row has been added to the problem tableau.    */
 299      (*pni)++;
 300      /*     return(nligne);    */
 301      goto clear;
 302    }
 303                                                   /*  case (c)          */
 304                         /* The new parameter has already been defined as a
 305                            quotient. It remains to express that the
 306                            remainder of that division is zero           */
 307    sol_if();
 308    sol_forme(nparm + 2);
 309    for (j = 0; j < nparm+1 ; j++)
 310      sol_val(discrm[j], UN);
 311    mpz_neg(x, UN);
 312    sol_val(x, UN);
 313    sol_nil();    /* No solution if the division is not even      */
 314                         /* Add a new column */
 315    if(ncol+1 >= (*ptp)-> width) {
 316      int dtw;
 317      dtw = mpz_sizeinbase(D, 2);
 318      *ptp = expanser(*ptp, *pnvar, *pni, ncol, 0, 0, dtw);
 319    }
 320           /* The new column is zeroed out by <<expanser>>          */
 321 /* Let c be the coefficient of parameter p in the i row. In <<coupure>>,
 322    this parameter has coefficient  - mod(-c, D). In <<discrp>>, this same
 323    parameter has coefficient mod(-c, D). The sum c + mod(-c, D) is obviously
 324    divisible by D.                                                      */
 325
 326    for (j = 0; j <= nparm; j++)
 327      /* Index(*ptp, i, j + nvar + 1) += discrp[j]; */
 328      mpz_add(Index(*ptp, i, j + nvar + 1),
 329              Index(*ptp, i, j + nvar + 1), discrp[j]);
 330    continue;
 331  }
 332  /* The solution is integral.                              */
 333  /* return 0;   */
 334  nligne = 0;
 335 clear :
 336  for(i=0; i <= ncol; i++)
 337    mpz_clear(coupure[i]);
 338  for(i=0; i <= nparm+1; i++){
 339    mpz_clear(discrp[i]);
 340    mpz_clear(discrm[i]);
 341  }
 342  mpz_clear(x); mpz_clear(d); mpz_clear(D);
 343  return(nligne);
 344 }
 345