source/traiter.c

   1 /******************************************************************************
   2  *                     PIP : Parametric Integer Programming                   *
   3  ******************************************************************************
   4  *                                 traiter.c                                  *
   5  ******************************************************************************
   6  *                                                                            *
   7  * Copyright Paul Feautrier, 1988-2005                                        *
   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 #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            /* NdCed : Attention au contexte == NULL ! */
 235            tPlus = expanser(context, nparm, nc, nparm+1, nparm, 1, 0);
 236            Flag(tPlus, nparm+nc) = Unknown;
 237            for(j = 0; j<=nparm; j++)
 238            #if defined(LINEAR_VALUE_IS_MP)
 239            mpz_set(Index(tPlus, nparm+nc, j),discr[j]); /* loop body. */
 240            mpz_set(Denom(tPlus, nparm+nc), UN);
 241            #else
 242            Index(tPlus, nparm+nc, j) = discr[j];        /* loop body. */
 243            Denom(tPlus, nparm+nc) = UN;
 244            #endif
 245
 246            p = sol_hwm();
 247            traiter(tPlus, NULL, Pip_True, nparm, 0, nc+1, 0, -1);
 248            cPlus = is_not_Nil(p);
 249            if(verbose>0){
 250              fprintf(dump, "\nThe positive case has been found ");
 251              fprintf(dump, cPlus? "possible\n": "impossible\n");
 252              fflush(dump);
 253            }
 254
 255            sol_reset(p);
 256            for(j = 0; j<nparm+1; j++)
 257            #if defined(LINEAR_VALUE_IS_MP)
 258            mpz_neg(discr[j], discr[j]);                 /* loop body. */
 259            mpz_sub_ui(discr[nparm], discr[nparm], (isCritic ? 1 : 2));
 260            #else
 261            discr[j] = -discr[j];                        /* loop body. */
 262            discr[nparm] = discr[nparm] - (isCritic ? 1 : 2);
 263            #endif
 264            tMinus = expanser(context, nparm, nc, nparm+1, nparm, 1, 0);
 265            Flag(tMinus, nparm+nc) = Unknown;
 266            for(j = 0; j<= nparm; j++)
 267            #if defined(LINEAR_VALUE_IS_MP)
 268            mpz_set(Index(tMinus, nparm+nc, j),discr[j]);/* loop body. */
 269            mpz_set(Denom(tMinus, nparm+nc), UN);
 270            #else
 271            Index(tMinus, nparm+nc, j) = discr[j];       /* loop body. */
 272            Denom(tMinus, nparm+nc) = UN;
 273            #endif
 274            traiter(tMinus, NULL, Pip_True, nparm, 0, nc+1, 0, -1);
 275            cMinus = is_not_Nil(p);
 276            if(verbose>0){
 277              fprintf(dump, "\nThe negative case has been found ");
 278              fprintf(dump, cMinus? "possible\n": "impossible\n");
 279              fflush(dump);
 280            }
 281
 282            sol_reset(p);
 283            if (cPlus && cMinus) {
 284                Flag(tp,i) = isCritic ? Critic : Unknown;
 285              }
 286            else if (cMinus)
 287               {Flag(tp,i) = Minus;
 288                break;
 289               }
 290            else {
 291              Flag(tp,i) = cPlus ? Plus : Zero;
 292            }
 293           }
 294      }
 295  tab_reset(q);
 296
 297  #if defined(LINEAR_VALUE_IS_MP)
 298  for(i=0; i<=nparm; i++)
 299    mpz_clear(discr[i]);
 300  #endif
 301
 302  return;
 303 }
 304
 305 Entier *valeur(Tableau *tp, int i, int j)
 306 {
 307  if(Flag(tp, i) & Unit)
 308      return(tp->row[i].objet.unit == j ? &Denom(tp,i) : &ZERO);
 309  else return(&Index(tp, i, j));
 310 }
 311
 312 void solution(Tableau *tp, int nvar, int nparm)
 313 {int i, j;
 314  int ncol = nvar + nparm + 1;
 315
 316  sol_list(nvar);
 317  for(i = 0; i<nvar; i++)
 318      {sol_forme(nparm+1);
 319       for(j = nvar+1; j<ncol; j++)
 320          sol_val(*valeur(tp, i, j), Denom(tp,i));
 321       sol_val(*valeur(tp, i, nvar), Denom(tp,i));
 322      }
 323 }
 324
 325 int choisir_piv(Tableau *tp, int pivi, int nvar, int nligne)
 326 {
 327  int j, k;
 328  Entier pivot, foo, x, y;
 329  int sgn_x, pivj = -1;
 330
 331  #if defined(LINEAR_VALUE_IS_MP)
 332  mpz_init(pivot); mpz_init(foo); mpz_init(x); mpz_init(y);
 333  #endif
 334
 335  for(j = 0; j<nvar; j++) {
 336     #if defined(LINEAR_VALUE_IS_MP)
 337     mpz_set(foo, Index(tp, pivi, j));
 338     if(mpz_sgn(foo) <= 0) continue;
 339     if(pivj < 0)
 340         {pivj = j;
 341          mpz_set(pivot, foo);
 342          continue;
 343         }
 344     for(k = 0; k<nligne; k++)
 345         {mpz_mul(x, pivot, *valeur(tp, k, j));
 346          mpz_mul(y, *valeur(tp, k, pivj), foo);
 347          mpz_sub(x, x, y);
 348          cross_product++;
 349          sgn_x = mpz_sgn(x);
 350          if(sgn_x) break;
 351         }
 352     if(sgn_x < 0)
 353         {pivj = j;
 354          mpz_set(pivot, foo);
 355         }
 356     #else
 357     if((foo = Index(tp, pivi, j)) <= 0) continue;
 358     if(pivj < 0)
 359         {pivj = j;
 360          pivot = foo;
 361          continue;
 362         }
 363     for(k = 0; k<nligne; k++)
 364         {x = pivot * (*valeur(tp, k, j)) - (*valeur(tp, k, pivj)) * foo;
 365          cross_product++;
 366          if(x) break;
 367         }
 368     if(x < 0)
 369         {pivj = j;
 370          pivot = foo;
 371         }
 372     #endif
 373  }
 374
 375  #if defined(LINEAR_VALUE_IS_MP)
 376  mpz_clear(pivot); mpz_clear(foo); mpz_clear(x); mpz_clear(y);
 377  #endif
 378
 379  return(pivj);
 380 }
 381
 382
 383 int pivoter(Tableau *tp, int pivi, int nvar, int nparm, int ni, int iq)
 384
 385 {int pivj;
 386  int ncol = nvar + nparm + 1;
 387  int nligne = nvar + ni;
 388  int i, j, k;
 389  Entier x, y, d, gcd, u, dpiv;
 390  int ff, fff;
 391  Entier pivot, foo, z;
 392  Entier ppivot, dppiv;
 393  Entier new[MAXCOL], *p, *q;
 394  Entier lpiv;
 395  int sgn_x;
 396  #if !defined(LINEAR_VALUE_IS_MP)
 397  char format_format[32];
 398
 399  sprintf(format_format, "\nPivot %s/%s\n", FORMAT, FORMAT);
 400  #endif
 401
 402  if(ncol >= MAXCOL) {
 403    fprintf(stdout, "Too much variables\n");
 404    exit(1);
 405  }
 406  if(0 > pivi || pivi >= nligne || Flag(tp, pivi) == Unit) {
 407    fprintf(stdout, "Syserr : pivoter : wrong pivot row\n");
 408    exit(1);
 409  }
 410
 411  pivj = choisir_piv(tp, pivi, nvar, nligne);
 412  if(pivj < 0) return(-1);
 413  if(pivj >= nvar) {
 414    fprintf(stdout, "Syserr : pivoter : wrong pivot\n");
 415    exit(1);
 416  }
 417
 418  #if defined(LINEAR_VALUE_IS_MP)
 419  mpz_init(x); mpz_init(y); mpz_init(d);
 420  mpz_init(gcd); mpz_init(u); mpz_init(dpiv);
 421  mpz_init(lpiv); mpz_init(pivot); mpz_init(foo);
 422  mpz_init(z); mpz_init(ppivot); mpz_init(dppiv);
 423
 424  for(i=0; i<ncol; i++)
 425    mpz_init(new[i]);
 426
 427  mpz_set(pivot, Index(tp, pivi, pivj));
 428  mpz_set(dpiv, Denom(tp, pivi));
 429  mpz_gcd(d, pivot, dpiv);
 430  mpz_divexact(ppivot, pivot, d);
 431  mpz_divexact(dppiv, dpiv, d);
 432  #else
 433  pivot = Index(tp, pivi, pivj);
 434  dpiv = Denom(tp, pivi);
 435  d = pgcd(pivot, dpiv);
 436  ppivot = pivot/d;
 437  dppiv = dpiv/d;
 438  #endif
 439
 440  if(verbose>1){
 441    #if defined(LINEAR_VALUE_IS_MP)
 442    fprintf(dump, "Pivot ");
 443    mpz_out_str(dump, 10, ppivot);
 444    putc('/', dump);
 445    mpz_out_str(dump, 10, dppiv);
 446    putc('\n', dump);
 447    #else
 448    fprintf(dump, format_format, ppivot, dppiv);
 449    #endif
 450    fprintf(dump, "%d x %d\n", pivi, pivj);
 451  }
 452
 453  #if defined(LINEAR_VALUE_IS_MP)
 454  mpz_fdiv_qr(x, y, tp->determinant, dppiv);
 455  #else
 456  for(i=0; i< tp->l_determinant; i++){
 457      d=pgcd(tp->determinant[i], dppiv);
 458      tp->determinant[i] /= d;
 459      dppiv /= d;
 460      }
 461  #endif
 462
 463  #if defined(LINEAR_VALUE_IS_MP)
 464  if(mpz_sgn(y) != 0){
 465  #else
 466  if(dppiv != 1) {
 467  #endif
 468    fprintf(stderr, "Integer overflow\n");
 469    if(verbose>0) fflush(dump);
 470    exit(1);
 471  }
 472
 473  #if defined(LINEAR_VALUE_IS_MP)
 474  mpz_mul(tp->determinant, x, ppivot);
 475  #else
 476  for(i=0; i<tp->l_determinant; i++)
 477      if(llog(tp->determinant[i]) + llog(ppivot) < 8*sizeof(Entier)){
 478          tp->determinant[i] *= ppivot;
 479          break;
 480          }
 481  if(i >= tp->l_determinant){
 482      tp->l_determinant++;
 483      if(tp->l_determinant >= MAX_DETERMINANT){
 484          fprintf(stderr, "Integer overflow : %d\n", tp->l_determinant);
 485          exit(1);
 486          }
 487      tp->determinant[i] = ppivot;
 488      }
 489  #endif
 490
 491  if(verbose>1){
 492    fprintf(dump, "determinant ");
 493    #if defined(LINEAR_VALUE_IS_MP)
 494    mpz_out_str(dump, 10, tp->determinant);
 495    #else
 496    for(i=0; i<tp->l_determinant; i++)
 497         fprintf(dump, FORMAT, tp->determinant[i]);
 498    #endif
 499    fprintf(dump, "\n");
 500  }
 501
 502
 503  for(j = 0; j<ncol; j++)
 504    #if defined(LINEAR_VALUE_IS_MP)
 505    if(j==pivj)
 506      mpz_set(new[j], dpiv);
 507    else
 508      mpz_neg(new[j], Index(tp, pivi, j));
 509    #else
 510    new[j] = (j == pivj ? dpiv : -Index(tp, pivi, j));
 511    #endif
 512
 513  for(k = 0; k<nligne; k++){
 514    if(Flag(tp,k) & Unit)continue;
 515    if(k == pivi)continue;
 516    #if defined(LINEAR_VALUE_IS_MP)
 517    mpz_set(foo, Index(tp, k, pivj));
 518    mpz_gcd(d, pivot, foo);
 519    mpz_divexact(lpiv, pivot, d);
 520    mpz_divexact(foo, foo, d);
 521    mpz_set(d, Denom(tp,k));
 522    mpz_mul(gcd, lpiv, d);
 523    mpz_set(Denom(tp, k), gcd);
 524    #else
 525    foo = Index(tp, k, pivj);
 526    d = pgcd(pivot, foo);
 527    lpiv = pivot/d;
 528    foo /= d;
 529    d = Denom(tp,k);
 530    gcd = lpiv * d;
 531    Denom(tp, k) = gcd;
 532    #endif
 533    p = tp->row[k].objet.val;
 534    q = tp->row[pivi].objet.val;
 535    for(j = 0; j<ncol; j++){
 536      if(j == pivj)
 537      #if defined(LINEAR_VALUE_IS_MP)
 538        mpz_mul(z, dpiv, foo);
 539      #else
 540        z = dpiv * foo;
 541      #endif
 542      else {
 543      #if defined(LINEAR_VALUE_IS_MP)
 544        mpz_mul(z, *p, lpiv);
 545        mpz_mul(y, *q, foo);
 546        mpz_sub(z, z, y);
 547      #else
 548        z = (*p) * lpiv - (*q) * foo;
 549      #endif
 550      }
 551      q++;
 552      cross_product++;
 553      #if defined(LINEAR_VALUE_IS_MP)
 554      mpz_set(*p, z);
 555      p++;
 556      if(mpz_cmp_ui(gcd, 1) != 0)
 557        mpz_gcd(gcd, gcd, z);
 558      #else
 559      *p++ = z;
 560      if(gcd != 1)
 561        gcd = pgcd(gcd, z);
 562      #endif
 563    }
 564    #if defined(LINEAR_VALUE_IS_MP)
 565    if(mpz_cmp_ui(gcd, 1) != 0){
 566      p = tp->row[k].objet.val;
 567      for(j = 0; j<ncol; j++){
 568        mpz_divexact(*p, *p, gcd);
 569        p++;
 570      }
 571    }
 572    mpz_divexact(Denom(tp,k), Denom(tp,k), gcd);
 573    #else
 574    if(gcd != 1) {
 575     p = tp->row[k].objet.val;
 576     for(j = 0; j<ncol; j++)
 577       *p++ /= gcd;
 578       Denom(tp,k) = Denom(tp,k)/gcd;
 579    }
 580    #endif
 581  }
 582  p = tp->row[pivi].objet.val;
 583  for(k = 0; k<nligne; k++)
 584    if((Flag(tp, k) & Unit) && tp->row[k].objet.unit == pivj) break;
 585  Flag(tp, k) = Plus;
 586  tp->row[k].objet.val = p;
 587  for(j = 0; j<ncol; j++)
 588    #if defined(LINEAR_VALUE_IS_MP)
 589    mpz_set(*p++, new[j]);
 590    #else
 591    *p++ = new[j];
 592    #endif
 593
 594  #if defined(LINEAR_VALUE_IS_MP)
 595  mpz_set(Denom(tp, k), pivot);
 596  Flag(tp, pivi) = Unit | Zero;
 597  mpz_set(Denom(tp, pivi), UN);
 598  #else
 599  Denom(tp, k) = pivot;
 600  Flag(tp, pivi) = Unit | Zero;
 601  Denom(tp, pivi) = UN;
 602  #endif
 603  tp->row[pivi].objet.unit = pivj;
 604
 605  for(k = 0; k<nligne; k++){
 606    ff = Flag(tp, k);
 607    if(ff & Unit) continue;
 608    #if defined(LINEAR_VALUE_IS_MP)
 609    sgn_x = mpz_sgn(Index(tp, k, pivj));
 610    #else
 611    sgn_x = Index(tp, k, pivj);
 612    #endif
 613    if(sgn_x < 0) fff = Minus;
 614    else if(sgn_x == 0) fff = Zero;
 615    else fff = Plus;
 616    if(fff != Zero && fff != ff)
 617      if(ff == Zero) ff = (fff == Minus ? Unknown : fff);
 618      else ff = Unknown;
 619    Flag(tp, k) = ff;
 620  }
 621
 622  if(verbose>2){
 623    fprintf(dump, "just pivoted\n");
 624    tab_display(tp, dump);
 625  }
 626
 627  #if defined(LINEAR_VALUE_IS_MP)
 628  mpz_clear(x); mpz_clear(y); mpz_clear(d); mpz_clear(gcd);
 629  mpz_clear(u); mpz_clear(dpiv); mpz_clear(lpiv);
 630  mpz_clear(pivot); mpz_clear(foo); mpz_clear(z);
 631  mpz_clear(ppivot); mpz_clear(dppiv);
 632
 633  for(i=0; i<ncol; i++)
 634    mpz_clear(new[i]);
 635  #endif
 636
 637  return(0);
 638 }
 639
 640 /* dans cette version, "traiter" modifie ineq; par contre
 641    le contexte est immediatement recopie' */
 642
 643 void traiter(tp, ctxt, iq, nvar, nparm, ni, nc, bigparm)
 644 Tableau *tp, *ctxt;
 645 int iq, nvar, nparm, ni, nc, bigparm;
 646 {
 647  int j;
 648  int pivi, nligne, ncol;
 649  struct high_water_mark x;
 650  Tableau *context;
 651  int dch, dcw;
 652  double s, t, d, smax;
 653  int i;
 654  struct L temp;
 655  Entier discr[MAXPARM];
 656
 657  #if !defined(LINEAR_VALUE_IS_MP)
 658  Entier D = UN;
 659  #endif
 660
 661  #if defined(LINEAR_VALUE_IS_MP)
 662  for(i=0; i<MAXPARM; i++)
 663    mpz_init(discr[i]);
 664  dcw = mpz_sizeinbase(tp->determinant, 2);
 665  #else
 666  dcw = 0;
 667  for(i=0; i<tp->l_determinant; i++)
 668    dcw += llog(tp->determinant[i]);
 669  #endif
 670  dch = 2 * dcw + 1;
 671  x = tab_hwm();
 672  nligne = nvar+ni;
 673
 674  context = expanser(ctxt, 0, nc, nparm+1, 0, dch, dcw);
 675 /*
 676  sort the rows in increasing order of the largest coefficient
 677 */
 678
 679  smax = 0.;
 680
 681  for(i=nvar; i<nligne; i++){
 682    if(Flag(tp,i) & Unit) continue;
 683    s = 0.;
 684    #if defined(LINEAR_VALUE_IS_MP)
 685    d = mpz_get_d(Denom(tp,i));
 686    for(j=0; j<nvar; j++){
 687      t = mpz_get_d(Index(tp,i,j))/d;
 688      s = max(s, abs(t));
 689    }
 690    #else
 691    d = (float) Denom(tp,i);
 692    for(j=0; j<nvar; j++){
 693        t = Index(tp,i,j)/d;
 694        s = max(s, abs(t));
 695        }
 696    #endif
 697    tp->row[i].size = s;
 698    smax = max(s, smax);
 699  }
 700
 701  for(i=nvar; i<nligne; i++){
 702    if(Flag(tp,i) & Unit) continue;
 703    s = smax;
 704    pivi = i;
 705    for(j=i; j<nligne; j++){
 706      if(Flag(tp,j) & Unit) continue;
 707      if(tp->row[j].size < s){
 708        s = tp->row[i].size;
 709        pivi = j;
 710      }
 711    }
 712    if(pivi != i) {
 713      temp = tp->row[pivi];
 714      tp->row[pivi] = tp->row[i];
 715      tp->row[i]=temp;
 716    }
 717  }
 718
 719
 720  for(;;) {
 721    if(verbose>2){
 722      fprintf(dump, "debut for\n");
 723      tab_display(tp, dump);
 724      fflush(dump);
 725    }
 726    nligne = nvar+ni; ncol = nvar+nparm+1;
 727    if(nligne > tp->height || ncol > tp->width) {
 728      fprintf(stdout, "Syserr : traiter : tableau too small\n");
 729      exit(1);
 730    }
 731    pivi = chercher(tp, Minus, nligne);
 732    if(pivi < nligne) goto pirouette;           /* There is a negative row   */
 733
 734    pivi = exam_coef(tp, nvar, ncol, bigparm);
 735
 736    if(verbose>2){
 737      fprintf(dump, "coefs examined\n");
 738      tab_display(tp, dump);
 739      fflush(dump);
 740    }
 741
 742    if(pivi < nligne) goto pirouette;
 743    /* There is a row whose coefficients are negative */
 744    compa_test(tp, context, ni, nvar, nparm, nc);
 745    if(verbose>2){
 746      fprintf(dump, "compatibility tested\n");
 747      tab_display(tp, dump);
 748      fflush(dump);
 749    }
 750
 751    pivi = chercher(tp, Minus, nligne);
 752    if(pivi < nligne) goto pirouette;
 753    /* The compatibility test has found a negative row */
 754    pivi = chercher(tp, Critic, nligne);
 755    if(pivi >= nligne)pivi = chercher(tp, Unknown, nligne);
 756    /* Here, the problem tree splits        */
 757    if(pivi < nligne) {
 758      Tableau * ntp;
 759      Entier com_dem;
 760      struct high_water_mark q;
 761      if(nc >= context->height) {
 762        #if defined(LINEAR_VALUE_IS_MP)
 763        dcw = mpz_sizeinbase(context->determinant,2);
 764        #else
 765        dcw = 0;
 766        for(i=0; i<tp->l_determinant; i++)
 767        dcw += llog(tp->determinant[i]);
 768        #endif
 769        dch = 2 * dcw + 1;
 770        context = expanser(context, 0, nc, nparm+1, 0, dch, dcw);
 771      }
 772      if(nparm >= MAXPARM) {
 773        fprintf(stdout, "Too much parameters : %d\n", nparm);
 774        exit(2);
 775      }
 776      q = tab_hwm();
 777      if(verbose>1)
 778        fprintf(stdout,"profondeur %d %lx\n", profondeur, q.top);
 779      ntp = expanser(tp, nvar, ni, ncol, 0, 0, 0);
 780      fflush(stdout);
 781      sol_if();
 782      sol_forme(nparm+1);
 783      #if defined(LINEAR_VALUE_IS_MP)
 784      mpz_init_set_ui(com_dem, 0);
 785      for(j = 0; j<nparm; j++) {
 786        mpz_set(discr[j], Index(tp, pivi, j + nvar +1));
 787        mpz_gcd(com_dem, com_dem, discr[j]);
 788      }
 789      mpz_set(discr[nparm], Index(tp, pivi, nvar));
 790      mpz_gcd(com_dem, com_dem, discr[nparm]);
 791      for(j = 0; j<=nparm; j++) {
 792        mpz_divexact(discr[j], discr[j], com_dem);
 793        mpz_set(Index(context, nc, j), discr[j]);
 794        sol_val(discr[j], UN);
 795      }
 796      mpz_clear(com_dem);
 797      Flag(context, nc) = Unknown;
 798      mpz_set(Denom(context, nc), UN);
 799      #else
 800      com_dem = 0;
 801      for(j = 0; j<nparm; j++) {
 802        discr[j] = Index(tp, pivi, j + nvar +1);
 803        com_dem = pgcd(com_dem, discr[j]);
 804      }
 805      discr[nparm] = Index(tp, pivi, nvar);
 806      com_dem = pgcd(com_dem, discr[nparm]);
 807      for(j = 0; j<=nparm; j++) {
 808        discr[j] /= com_dem;
 809        Index(context, nc, j) = discr[j];
 810        sol_val(discr[j], UN);
 811      }
 812      Flag(context, nc) = Unknown;
 813      Denom(context, nc) = UN;
 814      #endif
 815      Flag(ntp, pivi) = Plus;
 816      profondeur++;
 817      fflush(stdout);
 818      if(verbose > 0) fflush(dump);
 819      #if defined(LINEAR_VALUE_IS_MP)
 820      traiter(ntp, context, iq, nvar, nparm, ni, nc+1, bigparm);
 821      profondeur--;
 822      tab_reset(q);
 823      if(verbose>1)
 824        fprintf(stdout, "descente %d %lx\n", profondeur, tab_hwm().top);
 825      for(j = 0; j<nparm; j++)
 826        mpz_neg(Index(context, nc, j), Index(context, nc, j));
 827      mpz_add_ui(Index(context, nc, nparm), Index(context, nc, nparm), 1);
 828      mpz_neg(Index(context, nc, nparm), Index(context, nc, nparm));
 829      Flag(tp, pivi) = Minus;
 830      mpz_set(Denom(context, nc), UN);
 831      #else
 832      traiter(ntp, context, iq, nvar, nparm, ni, nc+1, bigparm);
 833      profondeur--;
 834      tab_reset(q);
 835      if(verbose>1)
 836        fprintf(stderr, "descente %d %lx\n", profondeur, tab_hwm().top);
 837      for(j = 0; j<nparm; j++)
 838        Index(context, nc, j) = - Index(context, nc, j);
 839      Index(context, nc, nparm) = - Index(context, nc, nparm) -1;
 840      Flag(tp, pivi) = Minus;
 841      Denom(context, nc) = UN;
 842      #endif
 843      nc++;
 844      goto pirouette;
 845    }
 846 /* Here, all rows are positive. Do we need an integral solution?      */
 847    if(!iq) {
 848      solution(tp, nvar, nparm);
 849      break;
 850    }
 851 /* Yes we do! */
 852    pivi = integrer(&tp, &context, &nvar, &nparm, &ni, &nc);
 853    if(pivi > 0) goto pirouette;
 854                     /* A cut has been inserted and is always negative */
 855 /* Here, either there is an integral solution, */
 856    if(pivi == 0) solution(tp, nvar, nparm);
 857 /* or no solution exists */
 858    else sol_nil();
 859    break;
 860
 861 /* Here, a negative row has been found. The call to <<pivoter>> executes
 862       a pivoting step                                                 */
 863
 864 pirouette :
 865      if(pivoter(tp, pivi, nvar, nparm, ni, iq) < 0) {
 866        sol_nil();
 867        break;
 868      }
 869  }
 870 /* Danger : a premature return would induce memory leaks   */
 871  tab_reset(x);
 872  #if defined(LINEAR_VALUE_IS_MP)
 873  for(i=0; i<MAXPARM; i++)
 874    mpz_clear(discr[i]);
 875  #endif
 876  return;
 877 }