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