source/traiter.c

   1 /******************************************************************************
   2  *                     PIP : Parametric Integer Programming                   *
   3  ******************************************************************************
   4  *                                 traiter.c                                  *
   5  ******************************************************************************
   6  *                                                                            *
   7  * Copyright Paul Feautrier, 1988-2005                                        *
   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
  30 #include "pip.h"
  31
  32 #define max(x,y) ((x) > (y)? (x) : (y))
  33
  34 extern long int cross_product, limit;
  35 extern int verbose;
  36 extern FILE *dump;
  37 extern int profondeur;
  38 extern int compa_count;
  39
  40 #if !defined(LINEAR_VALUE_IS_MP)
  41 int llog(Entier x)
  42 {int n = 0;
  43 /* x must be positive, you dummy */
  44  if(x<0) x=-x;
  45  while(x) x >>= 1, n++;
  46  return(n);
  47 }
  48 #endif
  49
  50 int chercher(Tableau *p, int masque, int n)
  51 {int i;
  52  for(i = 0; i<n; i++)
  53      if(p->row[i].flags & masque) break;
  54  return(i);
  55 }
  56
  57 /* il est convenu que traiter ne doit modifier ni le tableau, ni le contexte;
  58    le tableau peut grandir en cas de coupure (+1 en hauteur et +1 en largeur
  59    si nparm != 0) et en cas de partage (+1 en hauteur)(seulement si nparm != 0).
  60    le contexte peut grandir en cas de coupure (+2 en hauteur et +1 en largeur)
  61    (seulement si nparm !=0) et en cas de partage (+1 en hauteur)(nparm !=0).
  62    On estime le nombre de coupures a llog(D) et le nombre de partages a
  63    ni.
  64 */
  65
  66 Tableau *expanser(Tableau *tp, int virt, int reel, int ncol,
  67                                int off, int dh, int dw)
  68 {
  69  int i, j, ff;
  70  char *q; Entier *pq;
  71  Entier *pp, *qq;
  72  Tableau *rp;
  73  if(tp == NULL) return(NULL);
  74  rp = tab_alloc(reel+dh, ncol+dw, virt);
  75
  76  #if defined(LINEAR_VALUE_IS_MP)
  77  mpz_set(rp->determinant, tp->determinant);
  78  #else
  79  rp->l_determinant = tp->l_determinant;
  80  for(i=0; i<tp->l_determinant; i++)
  81      rp->determinant[i] = tp->determinant[i];
  82  #endif
  83  pq = (Entier *) & (rp->row[virt+reel+dh]);
  84  for(i = off; i<virt + reel; i++)
  85      {ff = Flag(rp, i) = Flag(tp, i-off);
  86       #if defined(LINEAR_VALUE_IS_MP)
  87       mpz_set(Denom(rp, i), Denom(tp, i-off));
  88       #else
  89       Denom(rp, i) = Denom(tp, i-off);
  90       #endif
  91       if(ff & Unit) rp->row[i].objet.unit = tp->row[i-off].objet.unit;
  92       else {
  93           rp->row[i].objet.val = pq;
  94           pq +=(ncol + dw);
  95           pp = tp->row[i-off].objet.val;
  96           qq = rp->row[i].objet.val;
  97           for(j = 0; j<ncol; j++)
  98              #if defined(LINEAR_VALUE_IS_MP)
  99              mpz_set(*qq++, *pp++);
 100              #else
 101              *qq++ = *pp++;
 102              #endif
 103           }
 104       }
 105  return(rp);
 106 }
 107
 108 /* Check for "obvious" signs of the parametric constant terms
 109  * of the inequalities.  As soon as a negative sign is found
 110  * we return from this function and handle this constraint
 111  * in the calling function.  The signs of the other constraints
 112  * are then mostly irrelevant.
 113  * If any of the negative signs is due to the "big parameter",
 114  * then we want to use this constraint first.
 115  * We therefore check for signs determined by the coefficient
 116  * of the big parameter first.
 117  */
 118 int exam_coef(Tableau *tp, int nvar, int ncol, int bigparm)
 119 {int i, j ;
 120  int ff, fff;
 121  #if defined(LINEAR_VALUE_IS_MP)
 122  int x;
 123  #else
 124  Entier x;
 125  #endif
 126  Entier *p;
 127
 128  if (bigparm >= 0)
 129     for (i = 0; i<tp->height; i++) {
 130         if (Flag(tp, i) != Unknown)
 131             continue;
 132         x = entier_sgn(Index(tp,i, bigparm));
 133         if (x < 0) {
 134             Flag(tp, i) = Minus;
 135             return i;
 136         } else if (x > 0)
 137             Flag(tp, i) = Plus;
 138     }
 139
 140  for(i = 0; i<tp->height; i++)
 141      {ff = Flag(tp,i);
 142       if(ff == 0) break;
 143       if(ff == Unknown) {
 144            ff = Zero;
 145            p = &(tp->row[i].objet.val[nvar+1]);
 146            for(j = nvar+1; j<ncol; j++) {
 147                 #if defined(LINEAR_VALUE_IS_MP)
 148                 x = mpz_sgn(*p); p++ ;
 149                 #else
 150                 x = *p++;
 151                 #endif
 152                 if(x<0) fff = Minus;
 153                 else if (x>0) fff = Plus;
 154                 else fff = Zero;
 155                 if(fff != Zero && fff != ff)
 156                     if(ff == Zero) ff = fff;
 157                     else {ff = Unknown;
 158                           break;
 159                          }
 160                }
 161 /* bug de'tecte' par [paf], 16/2/93 !
 162    Si tous les coefficients des parame`tres sont ne'gatifs
 163    et si le terme constant est nul, le signe est inconnu!!
 164    On traite donc spe'cialement le terme constant. */
 165            #if defined(LINEAR_VALUE_IS_MP)
 166            x = mpz_sgn(Index(tp, i, nvar));
 167            #else
 168            x = Index(tp, i, nvar);
 169            #endif
 170            if(x<0) fff = Minus;
 171            else if(x>0) fff = Plus;
 172            else fff = Zero;
 173 /* ici on a le signe du terme constant */
 174            switch(ff){
 175 /* le signe est inconnu si les coefficients sont positifs et
 176    le terme constant ne'gatif */
 177            case Plus: if(fff == Minus) ff = Unknown; break;
 178 /* si les coefficients sont tous nuls, le signe est celui
 179    du terme constant */
 180            case Zero: ff = fff; break;
 181 /* le signe est inconnu si les coefficients sont ne'gatifs,
 182    sauf si le terme constant est egalement negatif. */
 183            case Minus: if(fff != Minus) ff = Unknown; break;
 184 /* enfin, il n'y a rien a` dire si le signe des coefficients est inconnu */
 185            }
 186            Flag(tp, i) = ff;
 187            if(ff == Minus) return(i);
 188           }
 189       }
 190  return(i);
 191 }
 192
 193 void compa_test(Tableau *tp, Tableau *context,
 194                 int ni, int nvar, int nparm, int nc)
 195 {
 196  Entier discr[MAXPARM];
 197  int i, j;
 198  int ff;
 199  int cPlus, cMinus, isCritic;
 200  int verbold;
 201  Tableau *tPlus, *tMinus;
 202  int p;
 203  struct high_water_mark q;
 204
 205  if(nparm == 0) return;
 206  if(nparm >= MAXPARM) {
 207      fprintf(stderr, "Too much parameters : %d\n", nparm);
 208      exit(1);
 209      }
 210  q = tab_hwm();
 211  #if defined(LINEAR_VALUE_IS_MP)
 212  for(i=0; i<=nparm; i++)
 213    mpz_init(discr[i]);
 214  #endif
 215
 216  for(i = 0; i<ni + nvar; i++)
 217      {ff = Flag(tp,i);
 218       if(ff & (Critic | Unknown))
 219           {isCritic = Pip_True;
 220            for(j = 0; j<nvar; j++)
 221                  #if defined(LINEAR_VALUE_IS_MP)
 222                  if(mpz_sgn(Index(tp, i, j)) > 0)
 223                  #else
 224                  if(Index(tp, i, j) > 0)
 225                  #endif
 226                  {isCritic = Pip_False;
 227                   break;
 228                  }
 229            compa_count++;
 230            for(j = 0; j < nparm; j++)
 231            #if defined(LINEAR_VALUE_IS_MP)
 232            mpz_set(discr[j], Index(tp, i, j+nvar+1));   /* loop body. */
 233            mpz_set(discr[nparm], Index(tp, i, nvar));
 234            mpz_sub_ui(discr[nparm], discr[nparm], (isCritic ? 0 : 1));
 235            #else
 236            discr[j] = Index(tp, i, j+nvar+1);           /* loop body. */
 237            discr[nparm] = Index(tp, i, nvar)- (isCritic ? 0 : 1);
 238            #endif
 239            tPlus = expanser(context, nparm, nc, nparm+1, nparm, 1, 0);
 240            Flag(tPlus, nparm+nc) = Unknown;
 241            for(j = 0; j<=nparm; j++)
 242            #if defined(LINEAR_VALUE_IS_MP)
 243            mpz_set(Index(tPlus, nparm+nc, j),discr[j]); /* loop body. */
 244            mpz_set(Denom(tPlus, nparm+nc), UN);
 245            #else
 246            Index(tPlus, nparm+nc, j) = discr[j];        /* loop body. */
 247            Denom(tPlus, nparm+nc) = UN;
 248            #endif
 249
 250            p = sol_hwm();
 251            traiter(tPlus, NULL, nparm, 0, nc+1, 0, -1, TRAITER_INT);
 252            cPlus = is_not_Nil(p);
 253            if(verbose>0){
 254              fprintf(dump, "\nThe positive case has been found ");
 255              fprintf(dump, cPlus? "possible\n": "impossible\n");
 256              fflush(dump);
 257            }
 258
 259            sol_reset(p);
 260            for(j = 0; j<nparm+1; j++)
 261            #if defined(LINEAR_VALUE_IS_MP)
 262            mpz_neg(discr[j], discr[j]);                 /* loop body. */
 263            mpz_sub_ui(discr[nparm], discr[nparm], (isCritic ? 1 : 2));
 264            #else
 265            discr[j] = -discr[j];                        /* loop body. */
 266            discr[nparm] = discr[nparm] - (isCritic ? 1 : 2);
 267            #endif
 268            tMinus = expanser(context, nparm, nc, nparm+1, nparm, 1, 0);
 269            Flag(tMinus, nparm+nc) = Unknown;
 270            for(j = 0; j<= nparm; j++)
 271            #if defined(LINEAR_VALUE_IS_MP)
 272            mpz_set(Index(tMinus, nparm+nc, j),discr[j]);/* loop body. */
 273            mpz_set(Denom(tMinus, nparm+nc), UN);
 274            #else
 275            Index(tMinus, nparm+nc, j) = discr[j];       /* loop body. */
 276            Denom(tMinus, nparm+nc) = UN;
 277            #endif
 278            traiter(tMinus, NULL, nparm, 0, nc+1, 0, -1, TRAITER_INT);
 279            cMinus = is_not_Nil(p);
 280            if(verbose>0){
 281              fprintf(dump, "\nThe negative case has been found ");
 282              fprintf(dump, cMinus? "possible\n": "impossible\n");
 283              fflush(dump);
 284            }
 285
 286            sol_reset(p);
 287            if (cPlus && cMinus) {
 288                Flag(tp,i) = isCritic ? Critic : Unknown;
 289              }
 290            else if (cMinus)
 291               {Flag(tp,i) = Minus;
 292                break;
 293               }
 294            else {
 295              Flag(tp,i) = cPlus ? Plus : Zero;
 296            }
 297           }
 298      }
 299  tab_reset(q);
 300
 301  #if defined(LINEAR_VALUE_IS_MP)
 302  for(i=0; i<=nparm; i++)
 303    mpz_clear(discr[i]);
 304  #endif
 305
 306  return;
 307 }
 308
 309 Entier *valeur(Tableau *tp, int i, int j)
 310 {
 311  if(Flag(tp, i) & Unit)
 312      return(tp->row[i].objet.unit == j ? &Denom(tp,i) : &ZERO);
 313  else return(&Index(tp, i, j));
 314 }
 315
 316 void solution(Tableau *tp, int nvar, int nparm)
 317 {int i, j;
 318  int ncol = nvar + nparm + 1;
 319
 320  sol_list(nvar);
 321  for(i = 0; i<nvar; i++)
 322      {sol_forme(nparm+1);
 323       for(j = nvar+1; j<ncol; j++)
 324          sol_val(*valeur(tp, i, j), Denom(tp,i));
 325       sol_val(*valeur(tp, i, nvar), Denom(tp,i));
 326      }
 327 }
 328
 329 static void solution_dual(Tableau *tp, int nvar, int nparm, int *pos)
 330 {
 331     int i;
 332
 333     sol_list(tp->height - nvar);
 334     for (i = 0; i < tp->height - nvar; ++i) {
 335         sol_forme(1);
 336         if (Flag(tp, pos[i]) & Unit)
 337             sol_val(*valeur(tp, 0, tp->row[pos[i]].objet.unit), Denom(tp, 0));
 338         else
 339             sol_val(ZERO, UN);
 340     }
 341 }
 342
 343 int choisir_piv(Tableau *tp, int pivi, int nvar, int nligne)
 344 {
 345  int j, k;
 346  Entier pivot, foo, x, y;
 347  int sgn_x, pivj = -1;
 348
 349  #if defined(LINEAR_VALUE_IS_MP)
 350  mpz_init(pivot); mpz_init(foo); mpz_init(x); mpz_init(y);
 351  #endif
 352
 353  for(j = 0; j<nvar; j++) {
 354     #if defined(LINEAR_VALUE_IS_MP)
 355     mpz_set(foo, Index(tp, pivi, j));
 356     if(mpz_sgn(foo) <= 0) continue;
 357     if(pivj < 0)
 358         {pivj = j;
 359          mpz_set(pivot, foo);
 360          continue;
 361         }
 362     for(k = 0; k<nligne; k++)
 363         {mpz_mul(x, pivot, *valeur(tp, k, j));
 364          mpz_mul(y, *valeur(tp, k, pivj), foo);
 365          mpz_sub(x, x, y);
 366          cross_product++;
 367          sgn_x = mpz_sgn(x);
 368          if(sgn_x) break;
 369         }
 370     if(sgn_x < 0)
 371         {pivj = j;
 372          mpz_set(pivot, foo);
 373         }
 374     #else
 375     if((foo = Index(tp, pivi, j)) <= 0) continue;
 376     if(pivj < 0)
 377         {pivj = j;
 378          pivot = foo;
 379          continue;
 380         }
 381     for(k = 0; k<nligne; k++)
 382         {x = pivot * (*valeur(tp, k, j)) - (*valeur(tp, k, pivj)) * foo;
 383          cross_product++;
 384          if(x) break;
 385         }
 386     if(x < 0)
 387         {pivj = j;
 388          pivot = foo;
 389         }
 390     #endif
 391  }
 392
 393  #if defined(LINEAR_VALUE_IS_MP)
 394  mpz_clear(pivot); mpz_clear(foo); mpz_clear(x); mpz_clear(y);
 395  #endif
 396
 397  return(pivj);
 398 }
 399
 400
 401 int pivoter(Tableau *tp, int pivi, int nvar, int nparm, int ni)
 402
 403 {int pivj;
 404  int ncol = nvar + nparm + 1;
 405  int nligne = nvar + ni;
 406  int i, j, k;
 407  Entier x, y, d, gcd, u, dpiv;
 408  int ff, fff;
 409  Entier pivot, foo, z;
 410  Entier ppivot, dppiv;
 411  Entier new[MAXCOL], *p, *q;
 412  Entier lpiv;
 413  int sgn_x;
 414  #if !defined(LINEAR_VALUE_IS_MP)
 415  char format_format[32];
 416
 417  sprintf(format_format, "\nPivot %s/%s\n", FORMAT, FORMAT);
 418  #endif
 419
 420  if(ncol >= MAXCOL) {
 421    fprintf(stdout, "Too much variables\n");
 422    exit(1);
 423  }
 424  if(0 > pivi || pivi >= nligne || Flag(tp, pivi) == Unit) {
 425    fprintf(stdout, "Syserr : pivoter : wrong pivot row\n");
 426    exit(1);
 427  }
 428
 429  pivj = choisir_piv(tp, pivi, nvar, nligne);
 430  if(pivj < 0) return(-1);
 431  if(pivj >= nvar) {
 432    fprintf(stdout, "Syserr : pivoter : wrong pivot\n");
 433    exit(1);
 434  }
 435
 436  #if defined(LINEAR_VALUE_IS_MP)
 437  mpz_init(x); mpz_init(y); mpz_init(d);
 438  mpz_init(gcd); mpz_init(u); mpz_init(dpiv);
 439  mpz_init(lpiv); mpz_init(pivot); mpz_init(foo);
 440  mpz_init(z); mpz_init(ppivot); mpz_init(dppiv);
 441
 442  for(i=0; i<ncol; i++)
 443    mpz_init(new[i]);
 444
 445  mpz_set(pivot, Index(tp, pivi, pivj));
 446  mpz_set(dpiv, Denom(tp, pivi));
 447  mpz_gcd(d, pivot, dpiv);
 448  mpz_divexact(ppivot, pivot, d);
 449  mpz_divexact(dppiv, dpiv, d);
 450  #else
 451  pivot = Index(tp, pivi, pivj);
 452  dpiv = Denom(tp, pivi);
 453  d = pgcd(pivot, dpiv);
 454  ppivot = pivot/d;
 455  dppiv = dpiv/d;
 456  #endif
 457
 458  if(verbose>1){
 459    #if defined(LINEAR_VALUE_IS_MP)
 460    fprintf(dump, "Pivot ");
 461    mpz_out_str(dump, 10, ppivot);
 462    putc('/', dump);
 463    mpz_out_str(dump, 10, dppiv);
 464    putc('\n', dump);
 465    #else
 466    fprintf(dump, format_format, ppivot, dppiv);
 467    #endif
 468    fprintf(dump, "%d x %d\n", pivi, pivj);
 469  }
 470
 471  #if defined(LINEAR_VALUE_IS_MP)
 472  mpz_fdiv_qr(x, y, tp->determinant, dppiv);
 473  #else
 474  for(i=0; i< tp->l_determinant; i++){
 475      d=pgcd(tp->determinant[i], dppiv);
 476      tp->determinant[i] /= d;
 477      dppiv /= d;
 478      }
 479  #endif
 480
 481  #if defined(LINEAR_VALUE_IS_MP)
 482  if(mpz_sgn(y) != 0){
 483  #else
 484  if(dppiv != 1) {
 485  #endif
 486    fprintf(stderr, "Integer overflow\n");
 487    if(verbose>0) fflush(dump);
 488    exit(1);
 489  }
 490
 491  #if defined(LINEAR_VALUE_IS_MP)
 492  mpz_mul(tp->determinant, x, ppivot);
 493  #else
 494  for(i=0; i<tp->l_determinant; i++)
 495      if(llog(tp->determinant[i]) + llog(ppivot) < 8*sizeof(Entier)){
 496          tp->determinant[i] *= ppivot;
 497          break;
 498          }
 499  if(i >= tp->l_determinant){
 500      tp->l_determinant++;
 501      if(tp->l_determinant >= MAX_DETERMINANT){
 502          fprintf(stderr, "Integer overflow : %d\n", tp->l_determinant);
 503          exit(1);
 504          }
 505      tp->determinant[i] = ppivot;
 506      }
 507  #endif
 508
 509  if(verbose>1){
 510    fprintf(dump, "determinant ");
 511    #if defined(LINEAR_VALUE_IS_MP)
 512    mpz_out_str(dump, 10, tp->determinant);
 513    #else
 514    for(i=0; i<tp->l_determinant; i++)
 515         fprintf(dump, FORMAT, tp->determinant[i]);
 516    #endif
 517    fprintf(dump, "\n");
 518  }
 519
 520
 521  for(j = 0; j<ncol; j++)
 522    #if defined(LINEAR_VALUE_IS_MP)
 523    if(j==pivj)
 524      mpz_set(new[j], dpiv);
 525    else
 526      mpz_neg(new[j], Index(tp, pivi, j));
 527    #else
 528    new[j] = (j == pivj ? dpiv : -Index(tp, pivi, j));
 529    #endif
 530
 531  for(k = 0; k<nligne; k++){
 532    if(Flag(tp,k) & Unit)continue;
 533    if(k == pivi)continue;
 534    #if defined(LINEAR_VALUE_IS_MP)
 535    mpz_set(foo, Index(tp, k, pivj));
 536    mpz_gcd(d, pivot, foo);
 537    mpz_divexact(lpiv, pivot, d);
 538    mpz_divexact(foo, foo, d);
 539    mpz_set(d, Denom(tp,k));
 540    mpz_mul(gcd, lpiv, d);
 541    mpz_set(Denom(tp, k), gcd);
 542    #else
 543    foo = Index(tp, k, pivj);
 544    d = pgcd(pivot, foo);
 545    lpiv = pivot/d;
 546    foo /= d;
 547    d = Denom(tp,k);
 548    gcd = lpiv * d;
 549    Denom(tp, k) = gcd;
 550    #endif
 551    p = tp->row[k].objet.val;
 552    q = tp->row[pivi].objet.val;
 553    for(j = 0; j<ncol; j++){
 554      if(j == pivj)
 555      #if defined(LINEAR_VALUE_IS_MP)
 556        mpz_mul(z, dpiv, foo);
 557      #else
 558        z = dpiv * foo;
 559      #endif
 560      else {
 561      #if defined(LINEAR_VALUE_IS_MP)
 562        mpz_mul(z, *p, lpiv);
 563        mpz_mul(y, *q, foo);
 564        mpz_sub(z, z, y);
 565      #else
 566        z = (*p) * lpiv - (*q) * foo;
 567      #endif
 568      }
 569      q++;
 570      cross_product++;
 571      #if defined(LINEAR_VALUE_IS_MP)
 572      mpz_set(*p, z);
 573      p++;
 574      if(mpz_cmp_ui(gcd, 1) != 0)
 575        mpz_gcd(gcd, gcd, z);
 576      #else
 577      *p++ = z;
 578      if(gcd != 1)
 579        gcd = pgcd(gcd, z);
 580      #endif
 581    }
 582    #if defined(LINEAR_VALUE_IS_MP)
 583    if(mpz_cmp_ui(gcd, 1) != 0){
 584      p = tp->row[k].objet.val;
 585      for(j = 0; j<ncol; j++){
 586        mpz_divexact(*p, *p, gcd);
 587        p++;
 588      }
 589    }
 590    mpz_divexact(Denom(tp,k), Denom(tp,k), gcd);
 591    #else
 592    if(gcd != 1) {
 593     p = tp->row[k].objet.val;
 594     for(j = 0; j<ncol; j++)
 595       *p++ /= gcd;
 596       Denom(tp,k) = Denom(tp,k)/gcd;
 597    }
 598    #endif
 599  }
 600  p = tp->row[pivi].objet.val;
 601  for(k = 0; k<nligne; k++)
 602    if((Flag(tp, k) & Unit) && tp->row[k].objet.unit == pivj) break;
 603  Flag(tp, k) = Plus;
 604  tp->row[k].objet.val = p;
 605  for(j = 0; j<ncol; j++)
 606    #if defined(LINEAR_VALUE_IS_MP)
 607    mpz_set(*p++, new[j]);
 608    #else
 609    *p++ = new[j];
 610    #endif
 611
 612  #if defined(LINEAR_VALUE_IS_MP)
 613  mpz_set(Denom(tp, k), pivot);
 614  Flag(tp, pivi) = Unit | Zero;
 615  mpz_set(Denom(tp, pivi), UN);
 616  #else
 617  Denom(tp, k) = pivot;
 618  Flag(tp, pivi) = Unit | Zero;
 619  Denom(tp, pivi) = UN;
 620  #endif
 621  tp->row[pivi].objet.unit = pivj;
 622
 623  for(k = 0; k<nligne; k++){
 624    ff = Flag(tp, k);
 625    if(ff & Unit) continue;
 626    #if defined(LINEAR_VALUE_IS_MP)
 627    sgn_x = mpz_sgn(Index(tp, k, pivj));
 628    #else
 629    sgn_x = Index(tp, k, pivj);
 630    #endif
 631    if(sgn_x < 0) fff = Minus;
 632    else if(sgn_x == 0) fff = Zero;
 633    else fff = Plus;
 634    if(fff != Zero && fff != ff)
 635      if(ff == Zero) ff = (fff == Minus ? Unknown : fff);
 636      else ff = Unknown;
 637    Flag(tp, k) = ff;
 638  }
 639
 640  if(verbose>2){
 641    fprintf(dump, "just pivoted\n");
 642    tab_display(tp, dump);
 643  }
 644
 645  #if defined(LINEAR_VALUE_IS_MP)
 646  mpz_clear(x); mpz_clear(y); mpz_clear(d); mpz_clear(gcd);
 647  mpz_clear(u); mpz_clear(dpiv); mpz_clear(lpiv);
 648  mpz_clear(pivot); mpz_clear(foo); mpz_clear(z);
 649  mpz_clear(ppivot); mpz_clear(dppiv);
 650
 651  for(i=0; i<ncol; i++)
 652    mpz_clear(new[i]);
 653  #endif
 654
 655  return(0);
 656 }
 657
 658 /*
 659  * Sort the rows in increasing order of the largest coefficient
 660  * and (if TRAITER_DUAL is set) return the new position of the
 661  * original constraints.
 662  */
 663 static int *tab_sort_rows(Tableau *tp, int nvar, int nligne, int flags)
 664 {
 665     int i, j;
 666     int pivi;
 667     double s, t, d, smax = 0;
 668     struct L temp;
 669     int *pos = NULL, *ineq = NULL;
 670
 671     if (flags & TRAITER_DUAL) {
 672         ineq = malloc(tp->height * sizeof(int));
 673         pos = malloc((tp->height-nvar) * sizeof(int));
 674         if (!ineq || !pos) {
 675             fprintf(stderr, "Memory Overflow.\n") ;
 676             exit(1) ;
 677         }
 678     }
 679
 680     for (i = nvar; i < nligne; i++) {
 681         if (Flag(tp,i) & Unit)
 682             continue;
 683         s = 0;
 684         d = ENTIER_TO_DOUBLE(Denom(tp, i));
 685         for (j = 0; j < nvar; j++) {
 686             t = ENTIER_TO_DOUBLE(Index(tp,i,j))/d;
 687             s = max(s, abs(t));
 688         }
 689         tp->row[i].size = s;
 690         smax = max(s, smax);
 691         if (flags & TRAITER_DUAL)
 692             ineq[i] = i-nvar;
 693     }
 694
 695     for (i = nvar; i < nligne; i++) {
 696         if (Flag(tp,i) & Unit)
 697             continue;
 698         s = smax;
 699         pivi = i;
 700         for (j = i; j < nligne; j++) {
 701             if (Flag(tp,j) & Unit)
 702                 continue;
 703             if (tp->row[j].size < s) {
 704                 s = tp->row[j].size;
 705                 pivi = j;
 706             }
 707         }
 708         if (pivi != i) {
 709             temp = tp->row[pivi];
 710             tp->row[pivi] = tp->row[i];
 711             tp->row[i] = temp;
 712             if (flags & TRAITER_DUAL) {
 713                 j = ineq[i];
 714                 ineq[i] = ineq[pivi];
 715                 ineq[pivi] = j;
 716              }
 717         }
 718     }
 719
 720     if (flags & TRAITER_DUAL) {
 721         for (i = nvar; i < nligne; i++)
 722             pos[ineq[i]] = i;
 723         free(ineq);
 724     }
 725
 726     return pos;
 727 }
 728
 729 /* dans cette version, "traiter" modifie ineq; par contre
 730    le contexte est immediatement recopie' */
 731
 732 void traiter(Tableau *tp, Tableau *ctxt, int nvar, int nparm, int ni, int nc,
 733              int bigparm, int flags)
 734 {
 735  int j;
 736  int pivi, nligne, ncol;
 737  struct high_water_mark x;
 738  Tableau *context;
 739  int dch, dcw;
 740  int i;
 741  Entier discr[MAXPARM];
 742  int *pos;
 743
 744  #if !defined(LINEAR_VALUE_IS_MP)
 745  Entier D = UN;
 746  #endif
 747
 748  #if defined(LINEAR_VALUE_IS_MP)
 749  for(i=0; i<MAXPARM; i++)
 750    mpz_init(discr[i]);
 751  dcw = mpz_sizeinbase(tp->determinant, 2);
 752  #else
 753  dcw = 0;
 754  for(i=0; i<tp->l_determinant; i++)
 755    dcw += llog(tp->determinant[i]);
 756  #endif
 757  dch = 2 * dcw + 1;
 758  x = tab_hwm();
 759  nligne = nvar+ni;
 760
 761  context = expanser(ctxt, 0, nc, nparm+1, 0, dch, dcw);
 762
 763  pos = tab_sort_rows(tp, nvar, nligne, flags);
 764
 765  for(;;) {
 766    if(verbose>2){
 767      fprintf(dump, "debut for\n");
 768      tab_display(tp, dump);
 769      fflush(dump);
 770    }
 771    nligne = nvar+ni; ncol = nvar+nparm+1;
 772    if(nligne > tp->height || ncol > tp->width) {
 773      fprintf(stdout, "Syserr : traiter : tableau too small\n");
 774      exit(1);
 775    }
 776    pivi = chercher(tp, Minus, nligne);
 777    if(pivi < nligne) goto pirouette;           /* There is a negative row   */
 778
 779    pivi = exam_coef(tp, nvar, ncol, bigparm);
 780
 781    if(verbose>2){
 782      fprintf(dump, "coefs examined\n");
 783      tab_display(tp, dump);
 784      fflush(dump);
 785    }
 786
 787    if(pivi < nligne) goto pirouette;
 788    /* There is a row whose coefficients are negative */
 789    compa_test(tp, context, ni, nvar, nparm, nc);
 790    if(verbose>2){
 791      fprintf(dump, "compatibility tested\n");
 792      tab_display(tp, dump);
 793      fflush(dump);
 794    }
 795
 796    pivi = chercher(tp, Minus, nligne);
 797    if(pivi < nligne) goto pirouette;
 798    /* The compatibility test has found a negative row */
 799    pivi = chercher(tp, Critic, nligne);
 800    if(pivi >= nligne)pivi = chercher(tp, Unknown, nligne);
 801    /* Here, the problem tree splits        */
 802    if(pivi < nligne) {
 803      Tableau * ntp;
 804      Entier com_dem;
 805      struct high_water_mark q;
 806      if(nc >= context->height) {
 807        #if defined(LINEAR_VALUE_IS_MP)
 808        dcw = mpz_sizeinbase(context->determinant,2);
 809        #else
 810        dcw = 0;
 811        for(i=0; i<tp->l_determinant; i++)
 812        dcw += llog(tp->determinant[i]);
 813        #endif
 814        dch = 2 * dcw + 1;
 815        context = expanser(context, 0, nc, nparm+1, 0, dch, dcw);
 816      }
 817      if(nparm >= MAXPARM) {
 818        fprintf(stdout, "Too much parameters : %d\n", nparm);
 819        exit(2);
 820      }
 821      q = tab_hwm();
 822      if(verbose>1)
 823        fprintf(stdout,"profondeur %d %lx\n", profondeur, q.top);
 824      ntp = expanser(tp, nvar, ni, ncol, 0, 0, 0);
 825      fflush(stdout);
 826      sol_if();
 827      sol_forme(nparm+1);
 828      entier_init_zero(com_dem);
 829      for (j = 0; j < nparm; j++)
 830        entier_gcd(com_dem, com_dem, Index(tp, pivi, j + nvar +1));
 831      if (!(flags & TRAITER_INT))
 832          entier_gcd(com_dem, com_dem, Index(tp, pivi, nvar));
 833      for (j = 0; j < nparm; j++) {
 834        entier_divexact(Index(context, nc, j), Index(tp, pivi, j + nvar + 1), com_dem);
 835        sol_val(Index(context, nc, j), UN);
 836      }
 837      if (!(flags & TRAITER_INT))
 838          entier_divexact(Index(context, nc, nparm), Index(tp, pivi, nvar), com_dem);
 839      else
 840          entier_pdivision(Index(context, nc, nparm), Index(tp, pivi, nvar), com_dem);
 841      sol_val(Index(context, nc, nparm), UN);
 842      entier_clear(com_dem);
 843      Flag(context, nc) = Unknown;
 844      entier_set_si(Denom(context, nc), 1);
 845      Flag(ntp, pivi) = Plus;
 846      profondeur++;
 847      fflush(stdout);
 848      if(verbose > 0) fflush(dump);
 849      #if defined(LINEAR_VALUE_IS_MP)
 850      traiter(ntp, context, nvar, nparm, ni, nc+1, bigparm, flags);
 851      profondeur--;
 852      tab_reset(q);
 853      if(verbose>1)
 854        fprintf(stdout, "descente %d %lx\n", profondeur, tab_hwm().top);
 855      for(j = 0; j<nparm; j++)
 856        mpz_neg(Index(context, nc, j), Index(context, nc, j));
 857      mpz_add_ui(Index(context, nc, nparm), Index(context, nc, nparm), 1);
 858      mpz_neg(Index(context, nc, nparm), Index(context, nc, nparm));
 859      Flag(tp, pivi) = Minus;
 860      mpz_set(Denom(context, nc), UN);
 861      #else
 862      traiter(ntp, context, nvar, nparm, ni, nc+1, bigparm, flags);
 863      profondeur--;
 864      tab_reset(q);
 865      if(verbose>1)
 866        fprintf(stderr, "descente %d %lx\n", profondeur, tab_hwm().top);
 867      for(j = 0; j<nparm; j++)
 868        Index(context, nc, j) = - Index(context, nc, j);
 869      Index(context, nc, nparm) = - Index(context, nc, nparm) -1;
 870      Flag(tp, pivi) = Minus;
 871      Denom(context, nc) = UN;
 872      #endif
 873      nc++;
 874      goto pirouette;
 875    }
 876 /* Here, all rows are positive. Do we need an integral solution?      */
 877    if (!(flags & TRAITER_INT)) {
 878      solution(tp, nvar, nparm);
 879      if (flags & TRAITER_DUAL)
 880         solution_dual(tp, nvar, nparm, pos);
 881      break;
 882    }
 883 /* Yes we do! */
 884    pivi = integrer(&tp, &context, &nvar, &nparm, &ni, &nc, bigparm);
 885    if(pivi > 0) goto pirouette;
 886                     /* A cut has been inserted and is always negative */
 887 /* Here, either there is an integral solution, */
 888    if(pivi == 0) solution(tp, nvar, nparm);
 889 /* or no solution exists */
 890    else sol_nil();
 891    break;
 892
 893 /* Here, a negative row has been found. The call to <<pivoter>> executes
 894       a pivoting step                                                 */
 895
 896 pirouette :
 897      if (pivoter(tp, pivi, nvar, nparm, ni) < 0) {
 898        sol_nil();
 899        break;
 900      }
 901  }
 902 /* Danger : a premature return would induce memory leaks   */
 903  tab_reset(x);
 904  #if defined(LINEAR_VALUE_IS_MP)
 905  for(i=0; i<MAXPARM; i++)
 906    mpz_clear(discr[i]);
 907  #endif
 908  free(pos);
 909  return;
 910 }