source/integrer.c

   1 /******************************************************************************
   2  *                     PIP : Parametric Integer Programming                   *
   3  ******************************************************************************
   4  *                                integrer.c                                  *
   5  ******************************************************************************
   6  *                                                                            *
   7  * Copyright Paul Feautrier, 1988, 1993, 1994, 1996, 2002                     *
   8  *                                                                            *
   9  * This library is free software; you can redistribute it and/or modify it    *
  10  * under the terms of the GNU Lesser General Public License as published by   *
  11  * the Free Software Foundation; either version 2.1 of the License, or (at    *
  12  * your option) any later version.                                            *
  13  *                                                                            *
  14  * This software is distributed in the hope that it will be useful, but       *
  15  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
  16  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
  17  * for more details.                                                          *
  18  *                                                                            *
  19  * You should have received a copy of the GNU Lesser General Public License   *
  20  * along with this library; if not, write to the Free Software Foundation,    *
  21  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA         *
  22  *                                                                            *
  23  * Written by Paul Feautrier                                                  *
  24  *                                                                            *
  25  *****************************************************************************/
  26
  27 #include <stdio.h>
  28 #include <stdlib.h>
  29 #include <assert.h>
  30 #include <string.h>
  31
  32 #include "pip.h"
  33
  34 /*  The routines in this file are used to build a Gomory cut from
  35     a non-integral row of the problem tableau                             */
  36
  37 extern long int cross_product, limit;
  38 extern int verbose;
  39 extern int deepest_cut;
  40 extern FILE * dump;
  41 char compose[256];
  42
  43 /* mod(x,y) computes the remainder of x when divided by y. The difference
  44    with x%y is that the result is guaranteed to be positive, which is not
  45    always true for x%y.  This function is replaced by mpz_fdiv_r for MP. */
  46
  47 #if !defined(LINEAR_VALUE_IS_MP)
  48 Entier mod(x, y)
  49 Entier x, y;
  50 {Entier r;
  51  r = x % y;
  52  if(r<0) r += y;
  53  return(r);
  54 }
  55 #endif
  56
  57 /* this routine is useless at present                                     */
  58
  59 int non_borne(tp, nvar, D, bigparm)
  60 Tableau *tp;
  61 int nvar, bigparm;
  62 Entier D;
  63 {int i, ff;
  64  for(i = 0; i<nvar; i++)
  65      {ff = Flag(tp, i);
  66       if(bigparm > 0)
  67          {if(ff & Unit)return(Pip_True);
  68           if(Index(tp, i, bigparm) != D) return(Pip_True);
  69          }
  70       }
  71  return(Pip_False);
  72 }
  73
  74 /* This routine solve for z in the equation z.y = x (mod delta), provided
  75    y and delta are mutually prime. Remember that for multiple precision
  76    operation, the responsibility of creating and destroying <<z>> is the
  77    caller's.                                                                */
  78
  79 void bezout(Entier x, Entier y, Entier delta, Entier *z){
  80   Entier a, b, c, d, e, f, u, v, q, r;
  81 #if defined(LINEAR_VALUE_IS_MP)
  82   mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
  83   mpz_init(e); mpz_init(f); mpz_init(u); mpz_init(v);
  84   mpz_init(q); mpz_init(r);
  85   mpz_set_ui(a, 1); mpz_set_ui(b, 0); mpz_set_ui(c, 0);
  86   mpz_set_ui(d, 1); mpz_set(u, y); mpz_set(v, delta);
  87 #else
  88   a = 1; b = 0; c = 0; d = 1;
  89   u = y; v = delta;
  90 #endif
  91   for(;;){
  92 #if defined(LINEAR_VALUE_IS_MP)
  93     mpz_fdiv_qr(q, r, u, v);
  94     if(mpz_cmp_ui(r, 0) == 0) break;
  95     mpz_set(u, v);
  96     mpz_set(v, r);
  97     mpz_mul(e, q, c);
  98     mpz_sub(e, a, e);
  99     mpz_mul(f, q, d);
 100     mpz_sub(f, b, f);
 101     mpz_set(a, c);
 102     mpz_set(b, d);
 103     mpz_set(c, e);
 104     mpz_set(d, f);
 105 #else
 106     q = u / v;
 107     r = mod(u, v);
 108     if(r == 0) break;
 109     u = v;
 110     v = r;
 111     e = a - q*c;
 112     f = b - q*d;
 113     a = c;
 114     b = d;
 115     c = e;
 116     d = f;
 117 #endif
 118   }
 119 #if defined(LINEAR_VALUE_IS_MP)
 120   if(mpz_cmp_ui(v, 1) != 0)
 121     mpz_set_ui(*z, 0);
 122   else {
 123     mpz_mul(a, c, x);
 124     mpz_mod(*z, a, delta);
 125   }
 126   mpz_clear(a); mpz_clear(b); mpz_clear(c); mpz_clear(d);
 127   mpz_clear(e); mpz_clear(f); mpz_clear(u); mpz_clear(v);
 128   mpz_clear(q); mpz_clear(r);
 129
 130 #else
 131   if(v != 1) *z = 0; /* y and delta are not mutually prime */
 132   else *z = mod(c*x, delta);
 133 #endif
 134 }
 135
 136 Tableau *expanser();
 137
 138 /* cut: constant parameters denominator */
 139 static int add_parm(Tableau **pcontext, int nr, int *pnparm, int *pni, int *pnc,
 140                     Entier *cut)
 141 {
 142     int nparm = *pnparm;
 143     int i, j, k;
 144     Entier x;
 145
 146     entier_init(x);
 147
 148 /*        Build the definition of the new parameter into the solution :
 149       p_{nparm} = -(sum_{j=0}^{nparm-1} c_{nvar + 1 + j} p_j
 150                      + c_{nvar})/D                             (3)
 151          The minus sign is there to compensate the one in (1)     */
 152
 153     sol_new(nparm);
 154     sol_div();
 155     sol_forme(nparm+1);
 156     for (j = 0; j < nparm; j++) {
 157         entier_oppose(x, cut[1+j]);
 158         sol_val(x, UN);
 159     }
 160     entier_oppose(x, cut[0]);
 161     sol_val(x, UN);
 162     sol_val(cut[1+nparm], UN);              /* The divisor                */
 163
 164     if (nr+2 > (*pcontext)->height || nparm+1+1 > (*pcontext)->width) {
 165         int dcw, dch;
 166         dcw = entier_llog(cut[1+nparm]);
 167         dch = 2 * dcw + *pni;
 168         *pcontext = expanser(*pcontext, 0, nr, nparm+1, 0, dch, dcw);
 169     }
 170
 171 /* Since a new parameter is to be added, the constant term has to be moved
 172    right and a zero has to be inserted in all rows of the old context    */
 173
 174     for (k = 0; k < nr; k++) {
 175         entier_assign(Index(*pcontext, k, nparm+1), Index(*pcontext, k, nparm));
 176         entier_set_si(Index(*pcontext, k, nparm), 0);
 177     }
 178
 179 /* The value of the new parameter is specified by applying the definition of
 180    Euclidean division to (3) :
 181
 182  0<= - sum_{j=0}^{nparm-1} c_{nvar+1+j} p_j - c_{nvar} - D * p_{nparm} < D (4)
 183
 184    This formula gives two inequalities which are stored in the context    */
 185
 186     for (j = 0; j < nparm; j++) {
 187         entier_oppose(Index(*pcontext, nr, j), cut[1+j]);
 188         entier_assign(Index(*pcontext, nr+1, j), cut[1+j]);
 189     }
 190     entier_oppose(Index(*pcontext, nr, nparm), cut[1+nparm]);
 191     entier_assign(Index(*pcontext, nr+1, nparm), cut[1+nparm]);
 192     entier_assign(x, cut[0]);
 193     entier_oppose(Index(*pcontext, nr, nparm+1), x);
 194     entier_decrement(x, x);
 195     entier_addto(Index(*pcontext, nr+1, nparm+1), x, cut[1+nparm]);
 196
 197     Flag(*pcontext, nr) = Unknown;
 198     Flag(*pcontext, nr+1) = Unknown;
 199     entier_set_si(Denom(*pcontext, nr), 1);
 200     entier_set_si(Denom(*pcontext, nr+1), 1);
 201     (*pnparm)++;
 202     (*pnc) += 2;
 203     if (verbose > 0) {
 204         fprintf(dump, "enlarged context %d x %d\n", *pnparm, *pnc);
 205         fflush(dump);
 206     }
 207
 208     entier_clear(x);
 209 }
 210
 211 static int has_cut(Tableau *context, int nr, int nparm, int p, Entier *cut)
 212 {
 213     int row, col;
 214
 215     for (row = 0; row < nr; ++row) {
 216         if (entier_ne(Index(context, row, p), cut[1+nparm]))
 217             continue;
 218         if (entier_ne(Index(context, row, nparm), cut[0]))
 219             continue;
 220         for (col = p+1; col < nparm; ++col)
 221             if (entier_notzero_p(Index(context, row, col)))
 222                 break;
 223         if (col < nparm)
 224             continue;
 225         for (col = 0; col < p; ++col)
 226             if (entier_ne(Index(context, row, col), cut[1+col]))
 227                 break;
 228         if (col < p)
 229             continue;
 230         return 1;
 231     }
 232     return 0;
 233 }
 234
 235 /* cut: constant parameters denominator */
 236 static int find_parm(Tableau *context, int nr, int nparm, Entier *cut)
 237 {
 238     int p;
 239     int col;
 240     int found;
 241
 242     if (entier_notzero_p(cut[1+nparm-1]))
 243         return -1;
 244
 245     entier_addto(cut[0], cut[0], cut[1+nparm]);
 246     entier_decrement(cut[0], cut[0]);
 247     for (p = nparm-1; p >= 0; --p) {
 248         if (entier_notzero_p(cut[1+p]))
 249             break;
 250         if (!has_cut(context, nr, nparm, p, cut))
 251             continue;
 252         entier_increment(cut[0], cut[0]);
 253         entier_subtract(cut[0], cut[0], cut[1+nparm]);
 254         for (col = 0; col < 1+nparm+1; ++col)
 255             entier_oppose(cut[col], cut[col]);
 256         found = has_cut(context, nr, nparm, p, cut);
 257         for (col = 0; col < 1+nparm+1; ++col)
 258             entier_oppose(cut[col], cut[col]);
 259         if (found)
 260             return p;
 261         entier_addto(cut[0], cut[0], cut[1+nparm]);
 262         entier_decrement(cut[0], cut[0]);
 263     }
 264     entier_increment(cut[0], cut[0]);
 265     entier_subtract(cut[0], cut[0], cut[1+nparm]);
 266     return -1;
 267 }
 268
 269 /* integrer(.....) add a cut to the problem tableau, or return 0 when an
 270    integral solution has been found, or -1 when no integral solution
 271    exists.
 272
 273    Since integrer may add rows and columns to the problem tableau, its
 274    arguments are pointers rather than values. If a cut is constructed,
 275    ni increases by 1. If the cut is parametric, nparm increases by 1 and
 276    nc increases by 2.
 277                                                                          */
 278
 279 int integrer(Tableau **ptp, Tableau **pcontext,
 280              int *pnvar, int *pnparm, int *pni, int *pnc, int bigparm)
 281 {int ncol = *pnvar+*pnparm+1;
 282  int nligne = *pnvar + *pni;
 283  int nparm = *pnparm;
 284  int nvar = *pnvar;
 285  int ni = *pni;
 286  int nc = *pnc;
 287  Entier coupure[MAXCOL];
 288  int i, j, k, ff;
 289  Entier x, d;
 290  int ok_var, ok_const, ok_parm;
 291  Entier D;
 292     int parm;
 293
 294  Entier t, delta, tau, lambda;
 295
 296     if (ncol >= MAXCOL) {
 297         fprintf(stderr, "Too many variables: %d\n", ncol);
 298         exit(3);
 299     }
 300
 301  #if defined(LINEAR_VALUE_IS_MP)
 302  for(i=0; i<=ncol; i++)
 303    mpz_init(coupure[i]);
 304
 305  mpz_init(x); mpz_init(d); mpz_init(D);
 306  mpz_init(t); mpz_init(delta); mpz_init(tau); mpz_init(lambda);
 307  #endif
 308
 309
 310 /* search for a non-integral row */
 311  for(i = 0; i<nvar; i++) {
 312       #if defined(LINEAR_VALUE_IS_MP)
 313       mpz_set(D, Denom(*ptp, i));
 314       if(mpz_cmp_ui(D, 1) == 0) continue;
 315       #else
 316       D = Denom(*ptp,i);
 317       if(D == 1) continue;
 318       #endif
 319 /*                          If the common denominator of the row is 1
 320                             the row is integral                         */
 321       ff = Flag(*ptp, i);
 322       if(ff & Unit)continue;
 323 /*                          If the row is a Unit, it is integral        */
 324
 325 /*                          Here a portential candidate has been found.
 326                             Build the cut by reducing each coefficient
 327                             modulo D, the common denominator            */
 328       ok_var = Pip_False;
 329       for(j = 0; j<nvar; j++) {
 330          #if defined(LINEAR_VALUE_IS_MP)
 331          mpz_fdiv_r(x, Index(*ptp, i, j), D);
 332          mpz_set(coupure[j], x);
 333          #else
 334          x = coupure[j] = mod(Index(*ptp, i, j), D);
 335          #endif
 336             if (entier_pos_p(x))
 337                 ok_var = Pip_True;
 338           }
 339 /*                          Done for the coefficient of the variables.  */
 340
 341       #if defined(LINEAR_VALUE_IS_MP)
 342       mpz_neg(x, Index(*ptp, i, nvar));
 343       mpz_fdiv_r(x, x, D);
 344       mpz_neg(x, x);
 345       mpz_set(coupure[nvar], x);
 346       ok_const = mpz_cmp_ui(x, 0);
 347       #else
 348       x = coupure[nvar] = - mod(-Index(*ptp, i, nvar), D);
 349       ok_const = (x != 0);
 350       #endif
 351 /*                          This is the constant term                   */
 352       ok_parm = Pip_False;
 353       for(j = nvar+1; j<ncol; j++) {
 354          /* We assume that the big parameter is divisible by any number. */
 355          if (j == bigparm) {
 356             entier_set_si(coupure[j], 0);
 357             continue;
 358          }
 359          entier_oppose(x, Index(*ptp, i, j));
 360          entier_pmodulus(x, x, D);
 361          entier_oppose(coupure[j], x);
 362          if (entier_notzero_p(coupure[j]))
 363             ok_parm = Pip_True;
 364       }
 365 /*                          These are the parametric terms              */
 366
 367       #if defined(LINEAR_VALUE_IS_MP)
 368       mpz_set(coupure[ncol], D);
 369       #else
 370       coupure[ncol] = D;    /* Just in case                             */
 371       #endif
 372
 373 /* The question now is whether the cut is valid. The answer is given
 374 by the following decision table:
 375
 376 ok_var   ok_parm   ok_const
 377
 378   F        F         F       (a) continue, integral row
 379   F        F         T       (b) return -1, no solution
 380   F        T         F
 381                              (c) if the <<constant>> part is not divisible
 382                              by D then bottom else ....
 383   F        T         T
 384   T        F         F       (a) continue, integral row
 385   T        F         T       (d) constant cut
 386   T        T         F
 387                              (e) parametric cut
 388   T        T         T
 389
 390                                                                 case (a)  */
 391
 392       if(!ok_parm && !ok_const) continue;
 393       if(!ok_parm)
 394           if(ok_var) {                                   /*     case (d)  */
 395               if(nligne >= (*ptp)->height) {
 396                   int d, dth, dtw;
 397                   #if defined(LINEAR_VALUE_IS_MP)
 398                   d = mpz_sizeinbase(D, 2);
 399                   #else
 400                   d = llog(D);
 401                   #endif
 402                   dth = d;
 403                   *ptp = expanser(*ptp, nvar, ni, ncol, 0, dth, 0);
 404                   }
 405               /* Find the deepest cut*/
 406               if(deepest_cut){
 407 #if defined(LINEAR_VALUE_IS_MP)
 408               mpz_neg(t, coupure[nvar]);
 409               mpz_gcd(delta, t, D);
 410               mpz_divexact(tau, t, delta);
 411               mpz_divexact(d, D, delta);
 412               mpz_sub_ui(t, d, 1);
 413               bezout(t, tau, d, &lambda);
 414               mpz_gcd(t, lambda, D);
 415               while(mpz_cmp_ui(t, 1) != 0){
 416                 mpz_add(lambda, lambda, d);
 417                 mpz_gcd(t, lambda, D);
 418               }
 419               for(j=0; j<nvar; j++){
 420                 mpz_mul(t, lambda, coupure[j]);
 421                 mpz_fdiv_r(coupure[j], t, D);
 422               }
 423               mpz_mul(t, coupure[nvar], lambda);
 424               mpz_mod(t, t, D);
 425               mpz_sub(t, D, t);
 426               mpz_neg(coupure[nvar], t);
 427 #else
 428               t = -coupure[nvar];
 429               delta = pgcd(t,D);
 430               tau = t/delta;
 431               d = D/delta;
 432               bezout(d-1, tau, d, &lambda);
 433               while(pgcd(lambda, D) != 1)
 434                 lambda += d;
 435               for(j=0; j<nvar; j++)
 436                 coupure[j] = mod(lambda*coupure[j], D);
 437               coupure[nvar] = -mod(-lambda*coupure[nvar], D);
 438 #endif
 439               }
 440                          /* The cut has a negative <<constant>> part      */
 441               Flag(*ptp, nligne) = Minus;
 442               #if defined(LINEAR_VALUE_IS_MP)
 443               mpz_set(Denom(*ptp, nligne), D);
 444               #else
 445               Denom(*ptp, nligne) = D;
 446               #endif
 447                          /* Insert the cut */
 448               for(j = 0; j<ncol; j++)
 449                   #if defined(LINEAR_VALUE_IS_MP)
 450                   mpz_set(Index(*ptp, nligne, j), coupure[j]);
 451                   #else
 452                   Index(*ptp, nligne, j) = coupure[j];
 453                   #endif
 454                       /* A new row has been added to the problem tableau. */
 455               (*pni)++;
 456               if(verbose > 0) {
 457                 fprintf(dump, "just cut ");
 458                 if(deepest_cut){
 459                   fprintf(dump, "Bezout multiplier ");
 460 #if defined(LINEAR_VALUE_IS_MP)
 461                   mpz_out_str(dump, 10, lambda);
 462 #else
 463                   fprintf(dump, FORMAT, lambda);
 464 #endif
 465                 }
 466                 fprintf(dump, "\n");
 467                 k=0;
 468                 for(i=0; i<nvar; i++){
 469                   if(Flag(*ptp, i) & Unit){
 470 #if defined(LINEAR_VALUE_IS_MP)
 471                     fprintf(dump, "0 ");
 472 #else
 473                     sprintf(compose+k, "0 ");
 474 #endif
 475                     k += 2;
 476                   }
 477                   else {
 478 #if defined(LINEAR_VALUE_IS_MP)
 479                     k += mpz_out_str(dump, 10, Index(*ptp, i, nvar));
 480                     fprintf(dump, "/");
 481                     k++;
 482                     k += mpz_out_str(dump, 10, Denom(*ptp, i));
 483                     fprintf(dump, " ");
 484                     k++;
 485                     if(k > 60){
 486                       putc('\n', dump);
 487                       k = 0;
 488                     }
 489 #else
 490                     sprintf(compose+k, FORMAT, Index(*ptp, i, nvar));
 491                     k = strlen(compose);
 492                     sprintf(compose+k, "/");
 493                     k++;
 494                     sprintf(compose+k, FORMAT, Denom(*ptp, i));
 495                     k = strlen(compose);
 496                     sprintf(compose+k, " ");
 497                     k++;
 498                     if(k>60)  {
 499                       fputs(compose, dump);
 500                       putc('\n', dump);
 501                       k=0;
 502                     }
 503 #endif
 504                   }
 505                 }
 506                 fputs(compose, dump);
 507                 putc('\n', dump);
 508               }
 509               if(verbose > 2) tab_display(*ptp, dump);
 510               goto clear;
 511               }
 512           else {                                         /*   case (b)    */
 513             nligne = -1;
 514             goto clear;
 515           }
 516 /* In case (e), one has to introduce a new parameter and
 517    introduce its defining inequalities into the context.
 518
 519    Let the cut be    sum_{j=0}^{nvar-1} c_j x_j + c_{nvar} +             (2)
 520                      sum_{j=0}^{nparm-1} c_{nvar + 1 + j} p_j >= 0.       */
 521
 522         parm = find_parm(*pcontext, nc, nparm, coupure+nvar);
 523         if (parm == -1) {
 524             add_parm(pcontext, nc, pnparm, pni, pnc, coupure+nvar);
 525             parm = nparm;
 526         }
 527
 528         assert(ok_var);
 529           if(nligne >= (*ptp)->height || ncol >= (*ptp)->width) {
 530               int d, dth, dtw;
 531              #if defined(LINEAR_VALUE_IS_MP)
 532              d = mpz_sizeinbase(D, 2);
 533              #else
 534              d = llog(D);
 535              #endif
 536               dth = d + ni;
 537               dtw = d;
 538               *ptp = expanser(*ptp, nvar, ni, ncol, 0, dth, dtw);
 539               }
 540                          /* Zeroing out the new column seems to be useless
 541                             since <<expanser>> does it anyway            */
 542
 543                          /* The cut has a negative <<constant>> part    */
 544           Flag(*ptp, nligne) = Minus;
 545           #if defined(LINEAR_VALUE_IS_MP)
 546           mpz_set(Denom(*ptp, nligne), D);
 547           #else
 548           Denom(*ptp, nligne) = D;
 549           #endif
 550                  /* Insert the cut */
 551         for (j = 0; j < ncol; j++)
 552               #if defined(LINEAR_VALUE_IS_MP)
 553               mpz_set(Index(*ptp, nligne, j), coupure[j]);
 554               #else
 555               Index(*ptp, nligne, j) = coupure[j];
 556               #endif
 557         entier_addto(Index(*ptp, nligne, nvar+1+parm),
 558                     Index(*ptp, nligne, nvar+1+parm), coupure[ncol]);
 559                  /* A new row has been added to the problem tableau.    */
 560           (*pni)++;
 561           goto clear;
 562       }
 563  /* The solution is integral.                              */
 564     nligne = 0;
 565 clear:
 566    for(i=0; i <= ncol; i++)
 567         entier_clear(coupure[i]);
 568     entier_clear(x); entier_clear(d); entier_clear(D);
 569     entier_clear(t); entier_clear(tau); entier_clear(lambda); entier_clear(delta);
 570     return nligne;
 571 }
 572