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