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