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(Tableau **ptp, Tableau **pcontext,
 147              int *pnvar, int *pnparm, int *pni, int *pnc)
 148 {int ncol = *pnvar+*pnparm+1;
 149  int nligne = *pnvar + *pni;
 150  int nparm = *pnparm;
 151  int nvar = *pnvar;
 152  int ni = *pni;
 153  int nc = *pnc;
 154  Entier coupure[MAXCOL];
 155  int i, j, k, ff;
 156  Entier x, d;
 157  int ok_var, ok_const, ok_parm;
 158  Entier discrp[MAXPARM], discrm[MAXPARM];
 159  int llog();
 160  Entier D;
 161
 162  Entier t, delta, tau, lambda;
 163
 164  #if defined(LINEAR_VALUE_IS_MP)
 165  for(i=0; i<=ncol; i++)
 166    mpz_init(coupure[i]);
 167
 168  for(i=0; i<=nparm+1; i++){
 169    mpz_init(discrp[i]);
 170    mpz_init(discrm[i]);
 171  }
 172
 173  mpz_init(x); mpz_init(d); mpz_init(D);
 174  mpz_init(t); mpz_init(delta); mpz_init(tau); mpz_init(lambda);
 175  #endif
 176
 177
 178  if(ncol+1 >= MAXCOL) {
 179       fprintf(stderr, "Too much variables : %d\n", ncol);
 180       exit(3);
 181       }
 182
 183 /* search for a non-integral row */
 184  for(i = 0; i<nvar; i++) {
 185       #if defined(LINEAR_VALUE_IS_MP)
 186       mpz_set(D, Denom(*ptp, i));
 187       if(mpz_cmp_ui(D, 1) == 0) continue;
 188       #else
 189       D = Denom(*ptp,i);
 190       if(D == 1) continue;
 191       #endif
 192 /*                          If the common denominator of the row is 1
 193                             the row is integral                         */
 194       ff = Flag(*ptp, i);
 195       if(ff & Unit)continue;
 196 /*                          If the row is a Unit, it is integral        */
 197
 198 /*                          Here a portential candidate has been found.
 199                             Build the cut by reducing each coefficient
 200                             modulo D, the common denominator            */
 201       ok_var = Pip_False;
 202       for(j = 0; j<nvar; j++) {
 203          #if defined(LINEAR_VALUE_IS_MP)
 204          mpz_fdiv_r(x, Index(*ptp, i, j), D);
 205          mpz_set(coupure[j], x);
 206          #else
 207          x = coupure[j] = mod(Index(*ptp, i, j), D);
 208          #endif
 209           if(x > 0) ok_var = Pip_True;
 210           }
 211 /*                          Done for the coefficient of the variables.  */
 212
 213       #if defined(LINEAR_VALUE_IS_MP)
 214       mpz_neg(x, Index(*ptp, i, nvar));
 215       mpz_fdiv_r(x, x, D);
 216       mpz_neg(x, x);
 217       mpz_set(coupure[nvar], x);
 218       ok_const = mpz_cmp_ui(x, 0);
 219       #else
 220       x = coupure[nvar] = - mod(-Index(*ptp, i, nvar), D);
 221       ok_const = (x != 0);
 222       #endif
 223 /*                          This is the constant term                   */
 224       ok_parm = Pip_False;
 225       for(j = nvar+1; j<ncol; j++) {
 226          value_oppose(x, Index(*ptp, i, j));
 227          value_pmodulus(x, x, D);
 228          value_oppose(coupure[j], x);
 229          if (value_notzero_p(coupure[j]))
 230             ok_parm = Pip_True;
 231       }
 232 /*                          These are the parametric terms              */
 233
 234       #if defined(LINEAR_VALUE_IS_MP)
 235       mpz_set(coupure[ncol], D);
 236       #else
 237       coupure[ncol] = D;    /* Just in case                             */
 238       #endif
 239
 240 /* The question now is whether the cut is valid. The answer is given
 241 by the following decision table:
 242
 243 ok_var   ok_parm   ok_const
 244
 245   F        F         F       (a) continue, integral row
 246   F        F         T       (b) return -1, no solution
 247   F        T         F
 248                              (c) if the <<constant>> part is not divisible
 249                              by D then bottom else ....
 250   F        T         T
 251   T        F         F       (a) continue, integral row
 252   T        F         T       (d) constant cut
 253   T        T         F
 254                              (e) parametric cut
 255   T        T         T
 256
 257                                                                 case (a)  */
 258
 259       if(!ok_parm && !ok_const) continue;
 260       if(!ok_parm)
 261           if(ok_var) {                                   /*     case (d)  */
 262               if(nligne >= (*ptp)->height) {
 263                   int d, dth, dtw;
 264                   #if defined(LINEAR_VALUE_IS_MP)
 265                   d = mpz_sizeinbase(D, 2);
 266                   #else
 267                   d = llog(D);
 268                   #endif
 269                   dth = d;
 270                   *ptp = expanser(*ptp, nvar, ni, ncol, 0, dth, 0);
 271                   }
 272               /* Find the deepest cut*/
 273               if(deepest_cut){
 274 #if defined(LINEAR_VALUE_IS_MP)
 275               mpz_neg(t, coupure[nvar]);
 276               mpz_gcd(delta, t, D);
 277               mpz_divexact(tau, t, delta);
 278               mpz_divexact(d, D, delta);
 279               mpz_sub_ui(t, d, 1);
 280               bezout(t, tau, d, &lambda);
 281               mpz_gcd(t, lambda, D);
 282               while(mpz_cmp_ui(t, 1) != 0){
 283                 mpz_add(lambda, lambda, d);
 284                 mpz_gcd(t, lambda, D);
 285               }
 286               for(j=0; j<nvar; j++){
 287                 mpz_mul(t, lambda, coupure[j]);
 288                 mpz_fdiv_r(coupure[j], t, D);
 289               }
 290               mpz_mul(t, coupure[nvar], lambda);
 291               mpz_mod(t, t, D);
 292               mpz_sub(t, D, t);
 293               mpz_neg(coupure[nvar], t);
 294 #else
 295               t = -coupure[nvar];
 296               delta = pgcd(t,D);
 297               tau = t/delta;
 298               d = D/delta;
 299               bezout(d-1, tau, d, &lambda);
 300               while(pgcd(lambda, D) != 1)
 301                 lambda += d;
 302               for(j=0; j<nvar; j++)
 303                 coupure[j] = mod(lambda*coupure[j], D);
 304               coupure[nvar] = -mod(-lambda*coupure[nvar], D);
 305 #endif
 306               }
 307                          /* The cut has a negative <<constant>> part      */
 308               Flag(*ptp, nligne) = Minus;
 309               #if defined(LINEAR_VALUE_IS_MP)
 310               mpz_set(Denom(*ptp, nligne), D);
 311               #else
 312               Denom(*ptp, nligne) = D;
 313               #endif
 314                          /* Insert the cut */
 315               for(j = 0; j<ncol; j++)
 316                   #if defined(LINEAR_VALUE_IS_MP)
 317                   mpz_set(Index(*ptp, nligne, j), coupure[j]);
 318                   #else
 319                   Index(*ptp, nligne, j) = coupure[j];
 320                   #endif
 321                       /* A new row has been added to the problem tableau. */
 322               (*pni)++;
 323               if(verbose > 0) {
 324                 fprintf(dump, "just cut ");
 325                 if(deepest_cut){
 326                   fprintf(dump, "Bezout multiplier ");
 327 #if defined(LINEAR_VALUE_IS_MP)
 328                   mpz_out_str(dump, 10, lambda);
 329 #else
 330                   fprintf(dump, FORMAT, lambda);
 331 #endif
 332                 }
 333                 fprintf(dump, "\n");
 334                 k=0;
 335                 for(i=0; i<nvar; i++){
 336                   if(Flag(*ptp, i) & Unit){
 337 #if defined(LINEAR_VALUE_IS_MP)
 338                     fprintf(dump, "0 ");
 339 #else
 340                     sprintf(compose+k, "0 ");
 341 #endif
 342                     k += 2;
 343                   }
 344                   else {
 345 #if defined(LINEAR_VALUE_IS_MP)
 346                     k += mpz_out_str(dump, 10, Index(*ptp, i, nvar));
 347                     fprintf(dump, "/");
 348                     k++;
 349                     k += mpz_out_str(dump, 10, Denom(*ptp, i));
 350                     fprintf(dump, " ");
 351                     k++;
 352                     if(k > 60){
 353                       putc('\n', dump);
 354                       k = 0;
 355                     }
 356 #else
 357                     sprintf(compose+k, FORMAT, Index(*ptp, i, nvar));
 358                     k = strlen(compose);
 359                     sprintf(compose+k, "/");
 360                     k++;
 361                     sprintf(compose+k, FORMAT, Denom(*ptp, i));
 362                     k = strlen(compose);
 363                     sprintf(compose+k, " ");
 364                     k++;
 365                     if(k>60)  {
 366                       fputs(compose, dump);
 367                       putc('\n', dump);
 368                       k=0;
 369                     }
 370 #endif
 371                   }
 372                 }
 373                 fputs(compose, dump);
 374                 putc('\n', dump);
 375               }
 376               if(verbose > 2) tab_display(*ptp, dump);
 377 #if defined(LINEAR_VALUE_IS_MP)
 378               goto clear;
 379 #else
 380               return(nligne);
 381 #endif
 382               }
 383           else                                                               /*   case (b)    */
 384 #if defined(LINEAR_VALUE_IS_MP)
 385           { nligne = -1;
 386             goto clear;
 387           }
 388 #else
 389           return -1;
 390 #endif
 391 /* In cases (c) and (e), one has to introduce a new parameter and
 392    introduce its defining inequalities into the context.
 393
 394    Let the cut be    sum_{j=0}^{nvar-1} c_j x_j + c_{nvar} +             (2)
 395                      sum_{j=0}^{nparm-1} c_{nvar + 1 + j} p_j >= 0.       */
 396
 397
 398       if(nparm >= MAXPARM) {
 399           fprintf(stderr, "Too much parameters : %d\n", *pnparm);
 400           exit(4);
 401           }
 402 /*        Build the definition of the new parameter into the solution :
 403       p_{nparm} = -(sum_{j=0}^{nparm-1} c_{nvar + 1 + j} p_j
 404                      + c_{nvar})/D                             (3)
 405          The minus sign is there to compensate the one in (1)     */
 406
 407       sol_new(nparm);
 408       sol_div();
 409       sol_forme(nparm+1);
 410       for(j = 0; j<nparm; j++)
 411       #if defined(LINEAR_VALUE_IS_MP)
 412       { mpz_neg(x, coupure[j+nvar+1]);
 413         sol_val(x, UN);
 414       }
 415       mpz_neg(x, coupure[*pnvar]);
 416       sol_val(x, UN);
 417       #else
 418       sol_val(-coupure[j+nvar+1], UN); /* loop body. */
 419       sol_val(-coupure[*pnvar], UN);
 420       #endif
 421       sol_val(D, UN);                     /* The divisor                */
 422
 423 /* The value of the new parameter is specified by applying the definition of
 424    Euclidean division to (3) :
 425
 426  0<= - sum_{j=0}^{nparm-1} c_{nvar+1+j} p_j - c_{nvar} - D * p_{nparm} < D (4)
 427
 428    This formula gives two inequalities which are stored in the context    */
 429
 430       for(j = 0; j<nparm; j++) {
 431           #if defined(LINEAR_VALUE_IS_MP)
 432           mpz_set(x, coupure[j+nvar+1]);
 433           mpz_neg(discrp[j], x);
 434           mpz_set(discrm[j], x);
 435           #else
 436           x = coupure[j+nvar+1];
 437           discrp[j] = -x;
 438           discrm[j] = x;
 439           #endif
 440           }
 441       #if defined(LINEAR_VALUE_IS_MP)
 442       mpz_neg(discrp[nparm], D);
 443       mpz_set(discrm[nparm], D);
 444       mpz_set(x, coupure[nvar]);
 445       mpz_neg(discrp[nparm+1], x);
 446       mpz_sub_ui(x, x, 1);
 447       mpz_add(discrm[nparm+1], x, D);
 448       #else
 449       discrp[nparm] = -D;
 450       discrm[nparm] = D;
 451       x = coupure[nvar];
 452       discrp[(nparm)+1] = -x;
 453       discrm[(nparm)+1] = x + D -1;
 454       #endif
 455       if(nc+2 > (*pcontext)->height || nparm+1 > (*pcontext)->width) {
 456           int dcw, dch;
 457           #if defined(LINEAR_VALUE_IS_MP)
 458           dcw = mpz_sizeinbase(D, 2);
 459           #else
 460           dcw = llog(D);
 461           #endif
 462           dch = 2 * dcw + *pni;
 463           *pcontext = expanser(*pcontext, 0, nc, nparm+1, 0, dch, dcw);
 464           }
 465       /* Flag(*pcontext, *pnc) = 0; Probably useless see line A */
 466
 467 /* Since a new parameter is to be added, the constant term has to be moved
 468    right and a zero has to be inserted in all rows of the old context    */
 469
 470       for(k = 0; k < nc; k++) {
 471           #if defined(LINEAR_VALUE_IS_MP)
 472           mpz_set(Index(*pcontext, k, nparm+1), Index(*pcontext, k, nparm));
 473           mpz_set_ui(Index(*pcontext, k, nparm), 0);
 474           #else
 475           Index(*pcontext, k, nparm+1) = Index(*pcontext, k, nparm);
 476           Index(*pcontext, k, nparm) = 0;
 477           #endif
 478           }
 479 /* Now, insert the new rows                                              */
 480
 481       for(j = 0; j <= nparm+1; j++) {
 482           #if defined(LINEAR_VALUE_IS_MP)
 483           mpz_set(Index(*pcontext, nc, j), discrp[j]);
 484           mpz_set(Index(*pcontext, nc+1, j), discrm[j]);
 485           #else
 486           Index(*pcontext, nc, j) = discrp[j];
 487           Index(*pcontext, nc+1, j) = discrm[j];
 488           #endif
 489           }
 490       Flag(*pcontext, nc) = Unknown;                                /* A */
 491       Flag(*pcontext, nc+1) = Unknown;
 492       #if defined(LINEAR_VALUE_IS_MP)
 493       mpz_set(Denom(*pcontext, nc), UN);
 494       mpz_set(Denom(*pcontext, nc+1), UN);
 495       #else
 496       Denom(*pcontext, nc) = UN;
 497       Denom(*pcontext, nc+1) = UN;
 498       #endif
 499       (*pnparm)++;
 500       (*pnc) += 2;
 501       if(verbose > 0){
 502         fprintf(dump, "enlarged context %d x %d\n", *pnparm, *pnc);
 503         fflush(dump);
 504       }
 505                          /* end of the construction of the new parameter */
 506
 507       if(ok_var) {                                 /*   case (e)         */
 508           if(nligne >= (*ptp)->height || ncol >= (*ptp)->width) {
 509               int d, dth, dtw;
 510              #if defined(LINEAR_VALUE_IS_MP)
 511              d = mpz_sizeinbase(D, 2);
 512              #else
 513              d = llog(D);
 514              #endif
 515               dth = d + ni;
 516               dtw = d;
 517               *ptp = expanser(*ptp, nvar, ni, ncol, 0, dth, dtw);
 518               }
 519                          /* Zeroing out the new column seems to be useless
 520                             since <<expanser>> does it anyway            */
 521
 522                          /* The cut has a negative <<constant>> part    */
 523           Flag(*ptp, nligne) = Minus;
 524           #if defined(LINEAR_VALUE_IS_MP)
 525           mpz_set(Denom(*ptp, nligne), D);
 526           #else
 527           Denom(*ptp, nligne) = D;
 528           #endif
 529                  /* Insert the cut */
 530           for(j = 0; j<ncol+1; j++)
 531               #if defined(LINEAR_VALUE_IS_MP)
 532               mpz_set(Index(*ptp, nligne, j), coupure[j]);
 533               #else
 534               Index(*ptp, nligne, j) = coupure[j];
 535               #endif
 536                  /* A new row has been added to the problem tableau.    */
 537           (*pni)++;
 538           #if defined(LINEAR_VALUE_IS_MP)
 539           goto clear;
 540           #else
 541           return(nligne);
 542           #endif
 543           }
 544                                                   /*  case (c)          */
 545                         /* The new parameter has already been defined as a
 546                            quotient. It remains to express that the
 547                            remainder of that division is zero           */
 548       sol_if();
 549       sol_forme(nparm + 2);
 550       for (j = 0; j < nparm+1 ; j++)
 551           sol_val(discrm[j], UN);
 552           #if defined(LINEAR_VALUE_IS_MP)
 553           mpz_neg(x, UN);
 554           sol_val(x, UN);
 555           #else
 556           sol_val(-UN, UN);
 557           #endif
 558       sol_nil();    /* No solution if the division is not even      */
 559                         /* Add a new column */
 560       if(ncol+1 >= (*ptp)-> width) {
 561           int dtw;
 562           #if defined(LINEAR_VALUE_IS_MP)
 563           dtw = mpz_sizeinbase(D, 2);
 564           #else
 565           dtw = llog(D);
 566           #endif
 567           *ptp = expanser(*ptp, *pnvar, *pni, ncol, 0, 0, dtw);
 568           }
 569           /* The new column is zeroed out by <<expanser>>          */
 570 /* Let c be the coefficient of parameter p in the i row. In <<coupure>>,
 571    this parameter has coefficient  - mod(-c, D). In <<discrp>>, this same
 572    parameter has coefficient mod(-c, D). The sum c + mod(-c, D) is obviously
 573    divisible by D.                                                      */
 574
 575       for (j = 0; j <= nparm; j++)
 576           #if defined(LINEAR_VALUE_IS_MP)
 577           mpz_add(Index(*ptp, i, j + nvar + 1),
 578                   Index(*ptp, i, j + nvar + 1), discrp[j]);
 579           #else
 580           Index(*ptp, i, j + nvar + 1) += discrp[j];
 581           #endif
 582       tab_display(*ptp, stderr);
 583       exit(0);
 584       continue;
 585       }
 586  /* The solution is integral.                              */
 587  #if defined(LINEAR_VALUE_IS_MP)
 588  nligne = 0;
 589  clear :
 590    for(i=0; i <= ncol; i++)
 591      mpz_clear(coupure[i]);
 592    for(i=0; i <= nparm+1; i++){
 593      mpz_clear(discrp[i]);
 594      mpz_clear(discrm[i]);
 595    }
 596    mpz_clear(x); mpz_clear(d); mpz_clear(D);
 597    mpz_clear(t); mpz_clear(tau); mpz_clear(lambda); mpz_clear(delta);
 598  return(nligne);
 599  #else
 600    return 0;
 601  #endif
 602 }
 603