Pip/traiter.c

   1 /*--------------------------------------------------------------------*/
   2 /*                             T R A I T E R . C                      */
   3 /* Copyright Paul Feautrier, 1988, 1993, 1994, 1996                   */
   4 /* This file is part of the PIP software                              */
   5 /* PIP is NOT public domain. It can be                                */
   6 /* used freely for research and teaching                              */
   7 /* but NOT for commercial advantage.                                  */
   8 /*--------------------------------------------------------------------*/
   9
  10 #include <stdio.h>
  11 #include <stdlib.h>
  12
  13 #include <polylib/arithmetique.h>
  14 #include <polylib/polylib.h>
  15
  16 #include "type.h"
  17 #include "tab.h"
  18 #include "funcall.h"
  19
  20 #define MAX_BUFFER      1024
  21
  22 extern long int cross_product, limit;
  23 extern int verbose;
  24 extern FILE *dump;
  25 extern int profondeur;
  26 extern float clock;
  27 extern int compa_count;
  28
  29 int llog(Entier x)
  30 {int n = 0;
  31  while(x) x >>= 1, n++;
  32  return(n);
  33 }
  34
  35 int chercher(Tableau *p, int masque, int n)
  36 {int i;
  37  for(i = 0; i<n; i++)
  38      if(p->row[i].flags & masque) break;
  39  return(i);
  40 }
  41
  42 /* il est convenu que traiter ne doit modifier ni le tableau, ni le contexte;
  43    le tableau peut grandir en cas de coupure (+1 en hauteur et +1 en largeur
  44    si nparm != 0) et en cas de partage (+1 en hauteur)(seulement si nparm != 0).
  45    le contexte peut grandir en cas de coupure (+2 en hauteur et +1 en largeur)
  46    (seulement si nparm !=0) et en cas de partage (+1 en hauteur)(nparm !=0).
  47    On estime le nombre de coupures a llog(D) et le nombre de partages a
  48    ni.
  49 */
  50
  51 Tableau *expanser(Tableau *tp, int virt, int reel, int ncol,
  52                                int off, int dh, int dw)
  53 {
  54  int i, j, ff;
  55  char *q; Entier *pq;
  56  Entier *pp, *qq;
  57  Tableau *rp;
  58  rp = tab_alloc(reel+dh, ncol+dw, virt);
  59  q = (char *)rp + 2 * sizeof(int) + (virt + reel + dh) * sizeof(struct L);
  60  pq = (Entier *) q;
  61  for(i = off; i<virt + reel; i++)
  62      {ff = Flag(rp, i) = Flag(tp, i-off);
  63       Denom(rp, i) = Denom(tp, i-off);
  64       if(ff & Unit) rp->row[i].objet.unit = tp->row[i-off].objet.unit;
  65       else {
  66           rp->row[i].objet.val = pq;
  67           pq +=(ncol + dw);
  68           pp = tp->row[i-off].objet.val;
  69           qq = rp->row[i].objet.val;
  70           for(j = 0; j<ncol; j++) *qq++ = *pp++;
  71           }
  72       }
  73  return(rp);
  74 }
  75
  76 int exam_coef(Tableau *tp, int nvar, int ncol, int bigparm)
  77 {int i, j;
  78  int ff, fff;
  79  Entier x, *p;
  80  for(i = 0; i<tp->height; i++)
  81      {ff = Flag(tp,i);
  82       if(ff == 0) break;
  83       if(ff == Unknown)
  84           {if(bigparm >= 0 && (x = Index(tp,i, bigparm)))
  85               {if(x<0) {Flag(tp, i) = Minus;
  86                         return(i);
  87                        }
  88                else    Flag(tp, i) = Plus;
  89                continue;
  90               }
  91            ff = Zero;
  92            p = tp->row[i].objet.val+nvar+1;
  93            for(j = nvar+1; j<ncol; j++)
  94                {x = *p++;
  95                 if(x<0) fff = Minus;
  96                 else if (x>0) fff = Plus;
  97                 else fff = Zero;
  98                 if(fff != Zero && fff != ff) {
  99                     if(ff == Zero) ff = fff;
 100                     else {ff = Unknown;
 101                           break;
 102                          } }
 103                }
 104 /* bug de'tecte' par [paf], 16/2/93 !
 105    Si tous les coefficients des parame`tres sont ne'gatifs
 106    et si le terme constant est nul, le signe est inconnu!!
 107    On traite donc spe'cialement le terme constant. */
 108            x = Index(tp, i, nvar);
 109            if(x<0) fff = Minus;
 110            else if(x>0) fff = Plus;
 111            else fff = Zero;
 112 /* ici on a le signe du terme constant */
 113            switch(ff){
 114 /* le signe est inconnu si les coefficients sont positifs et
 115    le terme constant ne'gatif */
 116            case Plus: if(fff == Minus) ff = Unknown; break;
 117 /* si les coefficients sont tous nuls, le signe est celui
 118    du terme constant */
 119            case Zero: ff = fff; break;
 120 /* le signe est inconnu si les coefficients sont ne'gatifs,
 121    sauf si le terme constant est ‚galement n‚gatif. */
 122            case Minus: if(fff != Minus) ff = Unknown; break;
 123 /* enfin, il n'y a rien a` dire si le signe des coefficients est inconnu */
 124            }
 125            Flag(tp, i) = ff;
 126            if(ff == Minus) return(i);
 127           }
 128       }
 129  return(i);
 130 }
 131
 132 void compa_test(Tableau *tp, Tableau *context,
 133                 int ni, int nvar, int nparm, int nc)
 134 {
 135  Entier discr[MAXPARM];
 136  int i, j;
 137  int ff;
 138  int cPlus, cMinus, isCritic;
 139  int verbold;
 140  Tableau *tPlus, *tMinus;
 141  int p;
 142  struct high_water_mark q;
 143
 144  if(nparm == 0) return;
 145  if(nparm >= MAXPARM) {
 146      fprintf(stderr, "Too much parameters : %d\n", nparm);
 147      exit(1);
 148      }
 149  q = tab_hwm();
 150  for(i = 0; i<ni + nvar; i++)
 151      {ff = Flag(tp,i);
 152       if(ff & (Critic | Unknown))
 153           {isCritic = True;
 154            for(j = 0; j<nvar; j++) if(Index(tp, i, j) > 0)
 155                  {isCritic = False;
 156                   break;
 157                  }
 158            compa_count++;
 159            for(j = 0; j < nparm; j++) discr[j] = Index(tp, i, j+nvar+1);
 160            discr[nparm] = Index(tp, i, nvar)- (isCritic ? 0 : 1);
 161            tPlus = expanser(context, nparm, nc, nparm+1, nparm, 1, 0);
 162            Flag(tPlus, nparm+nc) = Unknown;
 163            for(j = 0; j<=nparm; j++)Index(tPlus, nparm+nc, j) = discr[j];
 164            Denom(tPlus, nparm+nc) = UN;
 165            verbold = verbose; verbose = 0;
 166            if(verbold && limit <= compa_count && compa_count <= limit+10) {
 167                tab_display(tPlus, dump);
 168              }
 169            p = sol_hwm();
 170            traiter(tPlus, NULL, True, UN, nparm, 0, nc+1, 0, -1);
 171            cPlus = is_not_Nil(p);
 172            sol_reset(p);
 173            for(j = 0; j<nparm+1; j++) discr[j] = -discr[j];
 174            discr[nparm] = discr[nparm] - (isCritic ? 1 : 2);
 175            tMinus = expanser(context, nparm, nc, nparm+1, nparm, 1, 0);
 176            Flag(tMinus, nparm+nc) = Unknown;
 177            for(j = 0; j<= nparm; j++)Index(tMinus, nparm+nc, j) = discr[j];
 178            Denom(tMinus, nparm+nc) = UN;
 179            if(verbold && limit <= compa_count && compa_count <= limit+10) {
 180                tab_display(tMinus, dump);
 181              }
 182            traiter(tMinus, NULL, True, UN, nparm, 0, nc+1, 0, -1);
 183            cMinus = is_not_Nil(p);
 184            verbose = verbold;
 185            sol_reset(p);
 186            if (cPlus && cMinus) {
 187                Flag(tp,i) = isCritic ? Critic : Unknown;
 188              }
 189            else if (cMinus)
 190               {Flag(tp,i) = Minus;
 191                break;
 192               }
 193            else {
 194              Flag(tp,i) = cPlus ? Plus : Zero;
 195            }
 196           }
 197      }
 198  tab_reset(q);
 199  return;
 200 }
 201
 202 Entier valeur(Tableau *tp, int i, int j, Entier D)
 203 {
 204  if(Flag(tp, i) & Unit)
 205      return(tp->row[i].objet.unit == j ? Denom(tp,i) : 0);
 206  else return(Index(tp, i, j));
 207 }
 208
 209 void solution(Tableau *tp, int nvar, int nparm, Entier D)
 210 {int i, j;
 211  int ncol = nvar + nparm + 1;
 212  sol_list(nvar);
 213  for(i = 0; i<nvar; i++)
 214      {sol_forme(nparm+1);
 215       for(j = nvar+1; j<ncol; j++)
 216          sol_val(valeur(tp, i, j, D), Denom(tp,i));
 217       sol_val(valeur(tp, i, nvar, D), Denom(tp,i));
 218      }
 219 }
 220
 221 int choisir_piv(Tableau *tp, int pivi, int nvar, int nligne, Entier D)
 222 {
 223  int j, k;
 224  long pivot, foo, x;
 225  int pivj = -1;
 226  for(j = 0; j<nvar; j++)
 227      {if((foo = Index(tp, pivi, j)) <= 0) continue;
 228       if(pivj < 0)
 229           {pivj = j;
 230            pivot = foo;
 231            continue;
 232           }
 233       for(k = 0; k<nligne; k++)
 234           {x = pivot * valeur(tp, k, j, D) - valeur(tp, k, pivj, D) * foo;
 235            cross_product++;
 236            if(x) break;
 237           }
 238       if(x < 0)
 239           {pivj = j;
 240            pivot = foo;
 241           }
 242      }
 243  return(pivj);
 244 }
 245
 246 Entier pivoter(Tableau *tp, int pivi, Entier D, int nvar,
 247                int nparm, int ni, int iq)
 248 {int pivj;
 249  int ncol = nvar + nparm + 1;
 250  int nligne = nvar + ni;
 251  int j, k;
 252  Entier x, d, gcd, dpiv;
 253  int ff, fff;
 254  Entier pivot, foo, z;
 255  Entier new[MAXCOL], *p, *q;
 256
 257  if(ncol >= MAXCOL) {
 258      fprintf(stderr, "Too much variables\n");
 259      exit(1);
 260      }
 261  if(0 > pivi || pivi >= nligne || Flag(tp, pivi) == Unit) {
 262      fprintf(stderr, "Syserr : pivoter : wrong pivot row\n");
 263      exit(1);
 264      }
 265  pivj = choisir_piv(tp, pivi, nvar, nligne, D);
 266  if(pivj < 0) return(-1);
 267  if(pivj >= nvar) {
 268      fprintf(stderr, "Syserr : pivoter : wrong pivot\n");
 269      exit(1);
 270      }
 271  pivot = Index(tp, pivi, pivj);
 272  dpiv = Denom(tp, pivi);
 273  if(D%dpiv) {
 274      char format[MAX_BUFFER];
 275      sprintf(format,"Integer overflow : %s\n",FORMAT);
 276      fprintf(stderr, format, D);
 277      tab_display(tp, stdout);
 278      exit(1);
 279      }
 280  D = (D/dpiv)*pivot;
 281  for(j = 0; j<ncol; j++)
 282       new[j] = (j == pivj ? dpiv : -Index(tp, pivi, j));
 283  for(k = 0; k<nligne; k++)
 284      {Entier lpiv;
 285       if(Flag(tp,k) & Unit)continue;
 286       if(k == pivi)continue;
 287       foo = Index(tp, k, pivj);
 288       d = pgcd_entier(pivot, foo);
 289       lpiv = pivot/d;
 290       foo /= d;
 291       d = Denom(tp,k);
 292       gcd = lpiv * d;
 293       Denom(tp, k) = gcd;
 294       p = tp->row[k].objet.val;
 295       q = tp->row[pivi].objet.val;
 296       for(j = 0; j<ncol; j++)
 297            {if(j == pivj)
 298                 z = dpiv * foo;
 299             else
 300                 z = (*p) * lpiv - (*q) * foo;
 301             q++;
 302             cross_product++;
 303             *p++ = z;
 304             if(gcd != 1) gcd = pgcd_entier(gcd, z);
 305            }
 306        if(gcd != 1) {
 307            p = tp->row[k].objet.val;
 308            for(j = 0; j<ncol; j++)
 309                *p++ /= gcd;
 310            }
 311        Denom(tp,k) = Denom(tp,k)/gcd;
 312       }
 313  p = tp->row[pivi].objet.val;
 314  for(k = 0; k<nligne; k++)
 315      if((Flag(tp, k) & Unit) && tp->row[k].objet.unit == pivj) break;
 316  Flag(tp, k) = Plus;
 317  tp->row[k].objet.val = p;
 318  for(j = 0; j<ncol; j++)
 319       *p++ = new[j];
 320  Denom(tp, k) = pivot;
 321  Flag(tp, pivi) = Unit | Zero;
 322  Denom(tp, pivi) = UN;
 323  tp->row[pivi].objet.unit = pivj;
 324
 325  for(k = 0; k<nligne; k++)
 326     {ff = Flag(tp, k);
 327      if(ff & Unit) continue;
 328      x = Index(tp, k, pivj);
 329      if(x < 0) fff = Minus;
 330      else if(x == 0) fff = Zero;
 331      else fff = Plus;
 332      if(fff != Zero && fff != ff) {
 333          if(ff == Zero) ff = (fff == Minus ? Unknown : fff);
 334          else ff = Unknown; }
 335      Flag(tp, k) = ff;
 336     }
 337  return(D);
 338 }
 339
 340 /* dans cette version, "traiter" modifie ineq; par contre
 341    le contexte est immediatement recopie' */
 342
 343 static Entier discr[MAXPARM];
 344
 345 Entier traiter(tp, ctxt, iq, D, nvar, nparm, ni, nc, bigparm)
 346 Tableau *tp, *ctxt;
 347 int iq, nvar, nparm, ni, nc, bigparm;
 348 Entier D;
 349 {
 350  int j;
 351  int pivi, nligne, ncol;
 352  struct high_water_mark x;
 353  Tableau *context;
 354  int dch, dcw;
 355  dcw = llog(D);
 356  dch = 2 * dcw + 1;
 357  x = tab_hwm();
 358
 359  context = expanser(ctxt, 0, nc, nparm+1, 0, dch, dcw);
 360  for(;;) {
 361      nligne = nvar+ni; ncol = nvar+nparm+1;
 362      if(nligne > tp->height || ncol > tp->width) {
 363          fprintf(stderr, "Syserr : traiter : tableau too small\n");
 364          exit(1);
 365          }
 366      pivi = chercher(tp, Minus, nligne);
 367      if(pivi < nligne) goto pirouette;         /* There is a negative row   */
 368
 369      pivi = exam_coef(tp, nvar, ncol, bigparm);
 370      if(pivi < nligne) goto pirouette;
 371                         /* There is a row whose coefficients are negative */
 372      compa_test(tp, context, ni, nvar, nparm, nc);
 373      pivi = chercher(tp, Minus, nligne);
 374      if(pivi < nligne) goto pirouette;
 375                       /* The compatibility test has found a negative row */
 376      pivi = chercher(tp, Critic, nligne);
 377      if(pivi >= nligne)pivi = chercher(tp, Unknown, nligne);
 378                                   /* Here, the problem tree splits        */
 379      if(pivi < nligne) {
 380          Tableau * ntp;
 381          Entier com_dem;
 382          struct high_water_mark q;
 383          if(nc >= context->height) {
 384              dcw = llog(D);
 385              dch = 2 * dcw + 1;
 386              context = expanser(context, 0, nc, nparm+1, 0, dch, dcw);
 387              }
 388          if(nparm >= MAXPARM) {
 389              char format[MAX_BUFFER];
 390              sprintf(format,"Too much parameters : %s\n",FORMAT);
 391              fprintf(stderr, format, nparm);
 392              exit(2);
 393              }
 394          q = tab_hwm();
 395          if(verbose){
 396              char format[MAX_BUFFER];
 397              sprintf(format,"profondeur %s %%lx\n",FORMAT);
 398              fprintf(dump,format, profondeur, q.top);
 399              }
 400          ntp = expanser(tp, nvar, ni, ncol, 0, 0, 0);
 401          sol_if();
 402          sol_forme(nparm+1);
 403          com_dem = 0;
 404          for(j = 0; j<nparm; j++) {
 405              discr[j] = Index(tp, pivi, j + nvar +1);
 406              com_dem = pgcd_entier(com_dem, discr[j]);
 407              }
 408          discr[nparm] = Index(tp, pivi, nvar);
 409          com_dem = pgcd_entier(com_dem, discr[nparm]);
 410          for(j = 0; j<=nparm; j++) {
 411              discr[j] /= com_dem;
 412              Index(context, nc, j) = discr[j];
 413              sol_val(discr[j], UN);
 414              }
 415          Flag(context, nc) = Unknown;
 416          Denom(context, nc) = UN;
 417          Flag(ntp, pivi) = Plus;
 418          profondeur++;
 419          traiter(ntp, context, iq, D, nvar, nparm, ni, nc+1, bigparm);
 420          profondeur--;
 421          tab_reset(q);
 422          if(verbose){
 423              char format[MAX_BUFFER];
 424              sprintf(format,"descente %s %%lx\n",FORMAT);
 425              fprintf(dump, format, profondeur, tab_hwm().top);
 426              }
 427          for(j = 0; j<nparm; j++)
 428              Index(context, nc, j) = - Index(context, nc, j);
 429          Index(context, nc, nparm) = - Index(context, nc, nparm) -1;
 430          Flag(tp, pivi) = Minus;
 431          Denom(context, nc) = UN;
 432          nc++;
 433          goto pirouette;
 434          }
 435 /* Here, all rows are positive. Do we need an integral solution?      */
 436     if(!iq) {
 437         solution(tp, nvar, nparm, D);
 438         break;
 439         }
 440 /* Yes we do! */
 441     pivi = integrer(&tp, &context, D, &nvar, &nparm, &ni, &nc);
 442     if(pivi > 0) goto pirouette;
 443                     /* A cut has been inserted and is always negative */
 444 /* Here, either there is an integral solution, */
 445     if(pivi == 0) solution(tp, nvar, nparm, D);
 446 /* or no solution exists */
 447     else sol_nil();
 448     break;
 449
 450 /* Here, a negative row has been found. The call to <<pivoter>> executes
 451       a pivoting step                                                 */
 452
 453 pirouette :
 454     if((D = pivoter(tp, pivi, D, nvar, nparm, ni, iq)) < 0) {
 455         sol_nil();
 456         break;
 457         }
 458     }
 459 /* Danger : a premature return would induce memory leaks   */
 460  tab_reset(x);
 461  return D;
 462 }
 463
 464
 465
 466
 467
 468
 469