source/tab.c

   1 /******************************************************************************
   2  *                     PIP : Parametric Integer Programming                   *
   3  ******************************************************************************
   4  *                                   tab.h                                    *
   5  ******************************************************************************
   6  *                                                                            *
   7  * Copyright Paul Feautrier, 1988, 1993, 1994, 1996, 2002                     *
   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 and Cedric Bastoul                               *
  24  *                                                                            *
  25  *****************************************************************************/
  26
  27
  28 #include <stdio.h>
  29 #include <stdlib.h>
  30
  31 #include "pip.h"
  32
  33 #define TAB_CHUNK 4096*sizeof(Entier)
  34
  35 static char *tab_free, *tab_top;
  36 static struct A *tab_base;
  37
  38 extern int allocation;
  39 extern long int cross_product, limit;
  40 static int chunk_count;
  41
  42 int dgetc(FILE *);
  43 #if defined(LINEAR_VALUE_IS_MP)
  44 int dscanf(FILE *, Entier);
  45 #else
  46 int dscanf(FILE *, Entier *);
  47 #endif
  48
  49 extern FILE * dump;
  50
  51 #define sizeof_struct_A ((sizeof(struct A) % sizeof(Entier)) ?              \
  52                          (sizeof(struct A) + sizeof(Entier)                 \
  53                                 - (sizeof(struct A) % sizeof(Entier))) :    \
  54                           sizeof(struct A))
  55
  56 void tab_init(void)
  57 {
  58  tab_free = malloc(sizeof_struct_A);
  59  if(tab_free == NULL)
  60      {fprintf(stderr, "Your computer doesn't have enough memory\n");
  61       exit(1);
  62      }
  63  allocation = 1;
  64  tab_top = tab_free + sizeof_struct_A;
  65  tab_base = (struct A *)tab_free;
  66  tab_free += sizeof_struct_A;
  67  tab_base->precedent = NULL;
  68  tab_base->bout = tab_top;
  69  tab_base->free = tab_free;
  70  chunk_count = 1;
  71 }
  72
  73
  74 void tab_close(void)
  75 {
  76   if (tab_base) free(tab_base);
  77 }
  78
  79
  80 struct high_water_mark tab_hwm(void)
  81 {struct high_water_mark p;
  82  p.chunk = chunk_count;
  83  p.top = tab_free;
  84  return p;
  85 }
  86
  87
  88 #if defined(LINEAR_VALUE_IS_MP)
  89 /* the clear_tab routine clears the GMP objects which may be referenced
  90    in the given Tableau.
  91 */
  92 void tab_clear(Tableau *tp)
  93 {
  94   int i, j;
  95   /* clear the determinant */
  96   mpz_clear(tp->determinant);
  97
  98   for(i=0; i<tp->height; i++){
  99     /* clear the denominator */
 100     mpz_clear(Denom(tp, i));
 101     if((Flag(tp, i) & Unit) == 0)
 102       for(j=0; j<tp->width; j++)
 103         mpz_clear(Index(tp,i,j));
 104   }
 105 }
 106 #endif
 107
 108 void tab_reset(struct high_water_mark by_the_mark)
 109
 110 {struct A *g;
 111  char *p;
 112  while(chunk_count > by_the_mark.chunk)
 113      {
 114       g = tab_base->precedent;
 115
 116       #if defined(LINEAR_VALUE_IS_MP)
 117       /* Before actually freeing the memory, one has to clear the
 118        * included Tableaux. If this is not done, the GMP objects
 119        * referenced in the Tableaux will be orphaned.
 120        */
 121
 122       /* Enumerate the included tableaux. */
 123       p = (char *)tab_base + sizeof_struct_A;
 124       while(p < tab_base->free){
 125         Tableau *pt;
 126         pt = (Tableau *) p;
 127         tab_clear(pt);
 128         p += pt->taille;
 129       }
 130       #endif
 131
 132       free(tab_base);
 133       tab_base = g;
 134       tab_top = tab_base->bout;
 135       chunk_count--;
 136      }
 137  if(chunk_count > 0) {
 138      #if defined(LINEAR_VALUE_IS_MP)
 139      /* Do not forget to clear the tables in the current chunk above the
 140         high water mark */
 141      p = (char *)by_the_mark.top;
 142      while(p < tab_base->free) {
 143         Tableau *pt;
 144         pt = (Tableau *) p;
 145         tab_clear(pt);
 146         p += pt->taille;
 147         }
 148      #endif
 149      tab_free = by_the_mark.top;
 150      tab_base->free = tab_free;
 151      }
 152  else {
 153      fprintf(stderr, "Syserr: tab_reset : error in memory allocation\n");
 154      exit(1);
 155      }
 156 }
 157
 158 Tableau * tab_alloc(int h, int w, int n)
 159
 160 /* h : le nombre de ligne reelles;
 161    n : le nombre de lignes virtuelles
 162 */
 163 {
 164  char *p; Tableau *tp;
 165  Entier *q;
 166  unsigned long taille;
 167  int i, j;
 168  taille = sizeof(struct T) + (h+n-1) * sizeof (struct L)
 169           + h * w * sizeof (Entier);
 170  if(tab_free + taille >= tab_top)
 171      {struct A * g;
 172       unsigned long d;
 173       d = taille + sizeof_struct_A;
 174       if(d < TAB_CHUNK) d = TAB_CHUNK;
 175       tab_free = malloc(d);
 176       if(tab_free == NULL)
 177           {printf("Memory overflow\n");
 178            exit(23);
 179           }
 180       chunk_count++;
 181       g = (struct A *)tab_free;
 182       g->precedent = tab_base;
 183       tab_top = tab_free + d;
 184       tab_free += sizeof_struct_A;
 185       tab_base = g;
 186       g->bout = tab_top;
 187      }
 188  p = tab_free;
 189  tab_free += taille;
 190  tab_base->free = tab_free;
 191  tp = (Tableau *)p;
 192  q = (Entier *)(p +  sizeof(struct T) + (h+n-1) * sizeof (struct L));
 193  #if defined(LINEAR_VALUE_IS_MP)
 194  mpz_init_set_ui(tp->determinant,1);
 195  #else
 196  tp->determinant[0] = (Entier) 1;
 197  tp->l_determinant = 1;
 198  #endif
 199  for(i = 0; i<n ; i++){
 200    tp->row[i].flags = Unit;
 201    tp->row[i].objet.unit = i;
 202    #if defined(LINEAR_VALUE_IS_MP)
 203    mpz_init_set_ui(Denom(tp, i), 1);
 204    #else
 205    Denom(tp, i) = UN ;
 206    #endif
 207  }
 208  for(i = n; i < (h+n); i++){
 209    tp->row[i].flags = 0;
 210    tp->row[i].objet.val = q;
 211    for(j = 0; j < w; j++)
 212    #if defined(LINEAR_VALUE_IS_MP)
 213    mpz_init_set_ui(*q++, 0); /* loop body. */
 214    mpz_init_set_ui(Denom(tp, i), 0);
 215    #else
 216    *q++ = 0;                 /* loop body. */
 217    Denom(tp, i) = ZERO ;
 218    #endif
 219  }
 220  tp->height = h + n; tp->width = w;
 221  #if defined(LINEAR_VALUE_IS_MP)
 222  tp->taille = taille ;
 223  #endif
 224
 225  return(tp);
 226 }
 227
 228 Tableau * tab_get(foo, h, w, n)
 229 FILE * foo;
 230 int h, w, n;
 231 {
 232  Tableau *p;
 233  int i, j, c;
 234  Entier x;
 235  #if defined(LINEAR_VALUE_IS_MP)
 236  mpz_init(x);
 237  #endif
 238
 239  p = tab_alloc(h, w, n);
 240  while((c = dgetc(foo)) != EOF)
 241       if(c == '(')break;
 242  for(i = n; i<h+n; i++)
 243      {p->row[i].flags = Unknown;
 244       #if defined(LINEAR_VALUE_IS_MP)
 245       mpz_set_ui(Denom(p, i), 1);
 246       #else
 247       Denom(p, i) = UN;
 248       #endif
 249       while((c = dgetc(foo)) != EOF)if(c == '[')break;
 250       for(j = 0; j<w; j++){
 251         #if defined(LINEAR_VALUE_IS_MP)
 252         if(dscanf(foo, x) < 0)
 253           return NULL;
 254         else
 255           mpz_set(p->row[i].objet.val[j], x);
 256         #else
 257         if(dscanf(foo, &x) < 0)
 258           return NULL;
 259         else
 260           p->row[i].objet.val[j] = x;
 261         #endif
 262         }
 263       }
 264       while((c = dgetc(foo)) != EOF)if(c == ']')break;
 265
 266  #if defined(LINEAR_VALUE_IS_MP)
 267  mpz_clear(x);
 268  #endif
 269
 270  return(p);
 271 }
 272
 273
 274 /* Fonction tab_Matrix2Tableau :
 275  * Cette fonction effectue la conversion du format de matrice de la polylib
 276  * vers le format de traitement de Pip. matrix est la matrice a convertir.
 277  * Nineq est le nombre d'inequations necessaires (dans le format de la
 278  * polylib, le premier element d'une ligne indique si l'equation decrite
 279  * est une inequation ou une egalite. Pip ne gere que les inequations. On
 280  * compte donc le nombre d'inequations total pour reserver la place
 281  * necessaire, et on scinde toute egalite p(x)=0 en p(x)>=0 et -p(x)>=0).
 282  * Nv est le nombre de variables dans la premiere serie de variables (c'est
 283  * a dire que si les premiers coefficients dans les lignes de la matrice
 284  * sont ceux des inconnues, Nv est le nombre d'inconnues, resp. parametres).
 285  * n est le nombre de lignes 'virtuelles' contenues dans la matrice (c'est
 286  * a dire en fait le nombre d'inconnues). Si Shift vaut 0, on va rechercher
 287  * le minimum lexicographique non-negatif, sinon on recherche le maximum
 288  * (Shift = 1) ou bien le minimum tout court (Shift = -1). La fonction
 289  * met alors en place le bignum s'il n'y est pas deja et prepare les
 290  * contraintes au calcul du maximum lexicographique.
 291  *
 292  * This function is called both for both the context (only parameters)
 293  * and the actual domain (variables + parameters).
 294  * Let Np be the number of parameters and Nn the number of variables.
 295  *
 296  * For the context, the columns in matrix are
 297  *              1 Np 1
 298  * while the result has
 299  *              Np Bg Urs_parms 1
 300  * Nv = Np + Bg; n = -1
 301  *
 302  * For the domain, matrix has
 303  *              1 Nn Np 1
 304  * while the result has
 305  *              Nn 1 Np Bg Urs_parms
 306  * Nv = Nn; n >= 0
 307  *
 308  * 27 juillet 2001 : Premiere version, Ced.
 309  * 30 juillet 2001 : Nombreuses modifications. Le calcul du nombre total
 310  *                   d'inequations (Nineq) se fait a present a l'exterieur.
 311  *  3 octobre 2001 : Pas mal d'ameliorations.
 312  * 18 octobre 2003 : Mise en place de la possibilite de calculer le
 313  *                   maximum lexicographique (parties 'if (Max)').
 314  */
 315 Tableau * tab_Matrix2Tableau(matrix, Nineq, Nv, n, Shift, Bg, Urs_parms)
 316 PipMatrix * matrix ;
 317 int Nineq, Nv, n, Shift, Bg, Urs_parms;
 318 { Tableau * p ;
 319   unsigned i, j, k, current, new, nb_columns, decal=0, bignum_is_new ;
 320   unsigned cst;
 321   int inequality, ctx;
 322   Entier bignum;
 323
 324   /* Are we dealing with the context? */
 325   ctx = n == -1;
 326   if (ctx)
 327     n = 0;
 328   entier_init(bignum);
 329   nb_columns = matrix->NbColumns - 1 ;
 330   /* S'il faut un BigNum et qu'il n'existe pas, on lui reserve sa place. */
 331   bignum_is_new = Shift && (Bg+ctx > (matrix->NbColumns - 2));
 332   if (bignum_is_new)
 333     nb_columns++;
 334   if (ctx) {
 335     Shift = 0;
 336     cst = Nv + Urs_parms;
 337   } else
 338     cst = Nv;
 339
 340   p = tab_alloc(Nineq,nb_columns+Urs_parms,n) ;
 341
 342   /* La variable decal sert a prendre en compte les lignes supplementaires
 343    * issues des egalites.
 344    */
 345   for (i = 0; i < matrix->NbRows; i++) {
 346     current = i + n + decal;
 347     Flag(p,current) = Unknown ;
 348     entier_set_si(Denom(p,current), 1);
 349     if (Shift)
 350       entier_set_si(bignum, 0);
 351     /* Pour passer l'indicateur d'egalite/inegalite. */
 352     inequality = entier_notzero_p(matrix->p[i][0]);
 353
 354     /* Dans le format de la polylib, l'element constant est place en
 355      * dernier. Dans le format de Pip, il se trouve apres la premiere
 356      * serie de variables (inconnues ou parametres). On remet donc les
 357      * choses dans l'ordre de Pip. Ici pour p(x) >= 0.
 358      */
 359     for (j=0;j<Nv;j++) {
 360       if (bignum_is_new && j == Bg)
 361         continue;
 362       if (Shift)
 363         entier_addto(bignum, bignum, matrix->p[i][1+j]);
 364       if (Shift > 0)
 365         entier_oppose(p->row[current].objet.val[j], matrix->p[i][1+j]);
 366       else
 367         entier_assign(p->row[current].objet.val[j], matrix->p[i][1+j]);
 368     }
 369     for (k=j=Nv+1;j<nb_columns;j++) {
 370         if (bignum_is_new && j == Bg)
 371           continue;
 372         entier_assign(p->row[current].objet.val[j], matrix->p[i][k++]);
 373     }
 374     for (j=0; j < Urs_parms; ++j) {
 375         int pos_n = nb_columns - ctx + j;
 376         int pos = pos_n - Urs_parms;
 377         if (pos <= Bg)
 378             --pos;
 379         entier_oppose(p->row[current].objet.val[pos_n],
 380                      p->row[current].objet.val[pos]);
 381     }
 382     entier_assign(p->row[current].objet.val[cst],
 383                  matrix->p[i][matrix->NbColumns-1]);
 384     if (Shift) {
 385       if (Shift < 0)
 386         entier_oppose(bignum, bignum);
 387
 388       if (bignum_is_new)
 389         entier_assign(p->row[current].objet.val[Bg], bignum);
 390       else
 391         entier_addto(p->row[current].objet.val[Bg],
 392                     p->row[current].objet.val[Bg], bignum);
 393     }
 394
 395     /* Et ici lors de l'ajout de -p(x) >= 0 quand on traite une egalite. */
 396     if (!inequality) {
 397       decal ++ ;
 398       new = current + 1 ;
 399       Flag(p,new)= Unknown ;
 400       entier_set_si(Denom(p,new), 1);
 401
 402       for (j=0;j<nb_columns+Urs_parms;j++)
 403         entier_oppose(p->row[new].objet.val[j], p->row[current].objet.val[j]);
 404     }
 405   }
 406   entier_clear(bignum);
 407
 408   return(p);
 409 }
 410
 411
 412 int tab_simplify(Tableau *tp, int cst)
 413 {
 414     int i, j;
 415     Entier gcd;
 416
 417     entier_init(gcd);
 418     for (i = 0; i < tp->height; ++i) {
 419         if (Flag(tp, i) & Unit)
 420             continue;
 421         entier_set_si(gcd, 0);
 422         for (j = 0; j < tp->width; ++j) {
 423             if (j == cst)
 424                 continue;
 425             entier_gcd(gcd, gcd, Index(tp, i, j));
 426             if (entier_one_p(gcd))
 427                 break;
 428         }
 429         if (entier_zero_p(gcd))
 430             continue;
 431         if (entier_one_p(gcd))
 432             continue;
 433         for (j = 0; j < tp->width; ++j) {
 434             if (j == cst)
 435                 entier_pdivision(Index(tp, i, j), Index(tp, i, j), gcd);
 436             else
 437                 entier_divexact(Index(tp, i, j), Index(tp, i, j), gcd);
 438         }
 439     }
 440     entier_clear(gcd);
 441
 442     return 0;
 443 }
 444
 445
 446 char *Attr[] = {"Unit", "+", "-", "0", "*", "?"};
 447
 448 void tab_display(p, foo)
 449 FILE *foo;
 450 Tableau *p;
 451 {
 452
 453  int i, j, ff, fff, n;
 454  Entier x, d;
 455  #if defined(LINEAR_VALUE_IS_MP)
 456  mpz_init(d);
 457  #endif
 458
 459  fprintf(foo, "%ld/[%d * %d]\n", cross_product, p->height, p->width);
 460  for(i = 0; i<p->height; i++){
 461    fff = ff = p->row[i].flags;
 462    /* if(fff ==0) continue; */
 463    #if defined(LINEAR_VALUE_IS_MP)
 464    mpz_set(d, Denom(p, i));
 465    #else
 466    d = Denom(p, i);
 467    #endif
 468    n = 0;
 469    while(fff){
 470      if(fff & 1) fprintf(foo, "%s ",Attr[n]);
 471      n++; fff >>= 1;
 472    }
 473    fprintf(foo, "%f #[", p->row[i].size);
 474    if(ff & Unit)
 475      for(j = 0; j<p->width; j++)
 476        fprintf(foo, " /%d/",(j == p->row[i].objet.unit)? 1: 0);
 477    else
 478      for(j = 0; j<p->width; j++){
 479        #if defined(LINEAR_VALUE_IS_MP)
 480        mpz_out_str(foo, 10, Index(p, i, j));
 481        putc(' ', foo);
 482        #else
 483        x = Index(p,i,j);
 484        fprintf(foo, FORMAT, x);
 485        fprintf(foo, " ");
 486        #endif
 487      }
 488    fprintf(foo, "]/");
 489    #if defined(LINEAR_VALUE_IS_MP)
 490    mpz_out_str(foo, 10, d);
 491    #else
 492    fprintf(foo, "%d", (int)d);
 493    #endif
 494    putc('\n', foo);
 495  }
 496  #if defined(LINEAR_VALUE_IS_MP)
 497  mpz_clear(d);
 498  #endif
 499 }