source/tabMP.c

   1 /********************************************************/
   2 /* Part of MultiPrecision PIP port (C) Zbigniew Chamski */
   3 /* and Paul Feautrier, 2001.                            */
   4 /* Based on straight PIP E.1 version by Paul Feautrier  */
   5 /* <Paul.Feautrier@inria.fr>                            */
   6 /*                                                      */
   7 /* and a previous port (C) Zbigniew CHAMSKI, 1993.      */
   8 /* <Zbigniew.Chamski@philips.com>                       */
   9 /*                                                      */
  10 /* Copying subject to the terms and conditions of the   */
  11 /* GNU General Public License.                          */
  12 /*                                                      */
  13 /* Send questions, bug reports and fixes to:            */
  14 /* <Paul.Feautrier@inria.fr>                            */
  15 /********************************************************/
  16
  17
  18 #include <stdio.h>
  19 #include <stdlib.h>
  20
  21 #include <piplib/piplib.h>
  22
  23 #define TAB_CHUNK 4096*sizeof(Entier)
  24
  25 static char *tab_free, *tab_top;
  26 static struct A *tab_base;
  27
  28 extern int allocation;
  29 extern long int cross_product, limit;
  30 static int chunk_count;
  31
  32 extern FILE * dump;
  33
  34 /*   Structure of the heap after execution of tab_init.
  35
  36      tab_base         tab_free         tab_top
  37         |                |                |
  38         |                |                |
  39         V                |                |
  40      struct A {          V                |
  41 NULL <---- precedent;    .  <-------------.
  42            bout;
  43            }
  44 */
  45
  46 void tab_init(void)
  47 {
  48  tab_free = malloc(sizeof (struct A));
  49  if(tab_free == NULL)
  50      {fprintf(stdout, "Your computer doesn't have memory\n");
  51       exit(1);
  52      }
  53  allocation = 1;
  54  tab_top = tab_free + sizeof (struct A);
  55  tab_base = (struct A *)tab_free;
  56  tab_free += sizeof(struct A);
  57  tab_base->precedent = NULL;
  58  tab_base->bout = tab_top;
  59  chunk_count = 1;
  60 }
  61
  62 struct high_water_mark tab_hwm(void)
  63 {struct high_water_mark p;
  64  p.chunk = chunk_count;
  65  p.top = tab_free;
  66  return p;
  67 }
  68
  69 /* the clear_tab routine clears the GMP objects which may be referenced
  70    in the given Tableau.
  71 */
  72
  73 void tab_clear(Tableau *tp)
  74 {
  75   int i, j;
  76   /* clear the determinant */
  77   mpz_clear(tp->determinant);
  78
  79   for(i=0; i<tp->height; i++){
  80     /* clear the denominator */
  81     mpz_clear(Denom(tp, i));
  82     if((Flag(tp, i) & Unit) == 0)
  83       for(j=0; j<tp->width; j++)
  84         mpz_clear(Index(tp,i,j));
  85   }
  86 }
  87
  88 void tab_reset(struct high_water_mark by_the_mark)
  89
  90 {struct A *g;
  91  char *p;
  92  while(chunk_count > by_the_mark.chunk)
  93      {
  94       g = tab_base->precedent;
  95 /* Before actually freeing the memory, one has to clear the included Tableaux.
  96    If this is not done, the GMP objects referenced in the Tableaux will be
  97    orphaned.
  98 */
  99
 100 /* Enumerate the included tableaux.
 101  */
 102       p = (char *)tab_base + sizeof(struct A);
 103       while(p < tab_free){
 104         Tableau *pt;
 105         pt = (Tableau *) p;
 106         tab_clear(pt);
 107         p += pt->taille;
 108       }
 109       free(tab_base);
 110       tab_base = g;
 111       tab_free = tab_base->bout;
 112       chunk_count--;
 113      }
 114  if(chunk_count > 0) tab_free = by_the_mark.top;
 115  else {
 116      fprintf(stdout, "Syserr: tab_reset : error in memory allocation\n");
 117      exit(1);
 118      }
 119 }
 120
 121 Tableau * tab_alloc(int h, int w, int n)
 122
 123 /* h : le nombre de ligne reelles;
 124    n : le nombre de lignes virtuelles
 125 */
 126 {
 127  char *p; Tableau *tp;
 128  Entier *q;
 129  unsigned long the_taille;
 130  int i, j;
 131  the_taille = sizeof(struct T) + (h+n-1) * sizeof (struct L)
 132           + h * w * sizeof (Entier);
 133  if(tab_free + the_taille >= tab_top)
 134      {struct A * g;
 135       unsigned long d;
 136       d = the_taille + sizeof(struct A);
 137       if(d < TAB_CHUNK) d = TAB_CHUNK;
 138       tab_free = malloc(d);
 139       if(tab_free == NULL)
 140           {printf("Memory overflow\n");
 141            exit(23);
 142           }
 143       chunk_count++;
 144       g = (struct A *)tab_free;
 145       g->precedent = tab_base;
 146       tab_top = tab_free + d;
 147       tab_free += sizeof(struct A);
 148       tab_base = g;
 149       g->bout = tab_top;
 150      }
 151  p = tab_free;
 152  tab_free += the_taille;
 153  tp = (Tableau *)p;
 154  q = (Entier *)(p +  sizeof(struct T) + (h+n-1) * sizeof (struct L));
 155  mpz_init_set_ui(tp->determinant, 1);
 156  for(i = 0; i<n ; i++){
 157    tp->row[i].flags = Unit;
 158    tp->row[i].objet.unit = i;
 159    mpz_init_set_ui(Denom(tp, i), 1);
 160  }
 161  for(i = n; i < (h+n); i++){
 162    tp->row[i].flags = 0;
 163    tp->row[i].objet.val = q;
 164    for(j = 0; j < w; j++)
 165      mpz_init_set_ui(*q++, 0);
 166    mpz_init_set_ui(Denom(tp, i), 0);
 167  }
 168  tp->height = h + n; tp->width = w; tp->taille = the_taille;
 169  return(tp);
 170 }
 171
 172 Tableau * tab_get(foo, h, w, n)
 173 FILE * foo;
 174 int h, w, n;
 175 {
 176  Tableau *p;
 177  int i, j, c;
 178  int x;
 179  p = tab_alloc(h, w, n);
 180  while((c = dgetc(foo)) != EOF)
 181       if(c == '(')break;
 182  for(i = n; i<h+n; i++){
 183    p->row[i].flags = Unknown;
 184    mpz_set_ui(Denom(p, i), 1);
 185    while((c = dgetc(foo)) != EOF)if(c == '[')break;
 186    for(j = 0; j<w; j++){
 187      if(dscanf(foo, FORMAT, &x) < 0)
 188        return NULL;
 189      else mpz_set_si(p->row[i].objet.val[j], x);
 190    }
 191  }
 192  while((c = dgetc(foo)) != EOF)if(c == ']')break;
 193
 194  return(p);
 195 }
 196
 197
 198 char *Attr[] = {"Unit", "+", "-", "0", "*", "?"};
 199
 200 void tab_display(p, foo)
 201 FILE *foo;
 202 Tableau *p;
 203 {
 204
 205  int i, j, ff, fff, n;
 206  Entier d;
 207  mpz_init(d);
 208  fprintf(foo, "%ld/[%d * %d]\n", cross_product, p->height, p->width);
 209  for(i = 0; i<p->height; i++){
 210    fff = ff = p->row[i].flags;
 211    if(fff == 0) continue;
 212    mpz_set(d, Denom(p, i));
 213    n = 0;
 214    while(fff){
 215      if(fff & 1) fprintf(foo, "%s ",Attr[n]);
 216      n++; fff >>= 1;
 217    }
 218    fprintf(foo, "%f #[", p->row[i].size);
 219    if(ff & Unit)
 220      for(j = 0; j<p->width; j++)
 221        fprintf(foo, " /%d/",(j == p->row[i].objet.unit)? 1: 0);
 222    else
 223      for(j = 0; j<p->width; j++){
 224        mpz_out_str(foo, 10, Index(p, i, j));
 225        putc(' ', foo);
 226      }
 227    fprintf(foo, "]/");
 228    mpz_out_str(foo, 10, d);
 229    putc('\n', foo);
 230  }
 231  mpz_clear(d);
 232 }
 233
 234
 235 /* Fonction tab_Matrix2Tableau :
 236  * Cette fonction effectue la conversion du format de matrice de la polylib
 237  * vers le format de traitement de Pip. matrix est la matrice a convertir.
 238  * Nineq est le nombre d'inequations necessaires (dans le format de la
 239  * polylib, le premier element d'une ligne indique si l'equation decrite
 240  * est une inequation ou une egalite. Pip ne gere que les inequations. On
 241  * compte donc le nombre d'inequations total pour reserver la place
 242  * necessaire, et on scinde toute egalite p(x)=0 en p(x)>=0 et -p(x)>=0).
 243  * Nv est le nombre de variables dans la premiere serie de variables (c'est
 244  * a dire que si les premiers coefficients dans les lignes de la matrice
 245  * sont ceux des inconnues, Nv est le nombre d'inconnues, resp. parametres).
 246  * n est le nombre de lignes 'virtuelles' contenues dans la matrice (c'est
 247  * a dire en fait le nombre d'inconnues).
 248  * 27 juillet 2001 : Premiere version, Ced.
 249  * 30 juillet 2001 : Nombreuses modifications. Le calcul du nombre total
 250  *                   d'inequations (Nineq) se fait a present a l'exterieur.
 251  *  3 octobre 2001 : Pas mal d'ameliorations.
 252  * 24 octobre 2002 : Premiere version MP.
 253  */
 254 Tableau * tab_Matrix2Tableau(PipMatrix * matrix, int Nineq, int Nv, int n)
 255 { Tableau * p ;
 256   unsigned i, j, end, current, new, nb_columns, decal=0 ;
 257   Entier * entier, inequality ;
 258
 259   mpz_init(inequality) ;
 260   nb_columns = matrix->NbColumns - 1 ;
 261   p = tab_alloc(Nineq,nb_columns,n) ;
 262
 263   /* La variable decal sert a prendre en compte les lignes supplementaires
 264    * issues des egalites.
 265    */
 266   end = matrix->NbRows + n ;
 267   for (i=n;i<end;i++)
 268   { current = i + decal ;
 269     Flag(p,current) = Unknown ;
 270     mpz_set_ui(Denom(p,current),1) ;
 271     entier = *(matrix->p + i - n) ;
 272     /* Pour passer l'indicateur d'egalite/inegalite. */
 273     mpz_set(inequality,*entier) ;
 274     entier ++ ;
 275
 276     /* Dans le format de la polylib, l'element constant est place en
 277      * dernier. Dans le format de Pip, il se trouve apres la premiere
 278      * serie de variables (inconnues ou parametres). On remet donc les
 279      * choses dans l'ordre de Pip. Ici pour p(x) >= 0.
 280      */
 281     for (j=0;j<Nv;j++)
 282     mpz_set(*(p->row[current].objet.val + j),*entier++) ;
 283     for (j=Nv+1;j<nb_columns;j++)
 284     mpz_set(*(p->row[current].objet.val + j),*entier++) ;
 285     mpz_set(*(p->row[current].objet.val + Nv),*entier) ;
 286
 287     /* Et ici lors de l'ajout de -p(x) >= 0 quand on traite une egalite. */
 288     if (!inequality)
 289     { decal ++ ;
 290       new = current + 1 ;
 291       Flag(p,new)= Unknown ;
 292       mpz_set(Denom(p,new),UN) ;
 293
 294       for (j=0;j<nb_columns;j++)
 295       mpz_neg(*(p->row[new].objet.val + j),*(p->row[current].objet.val + j)) ;
 296     }
 297   }
 298   mpz_clear(inequality);
 299   return(p);
 300 }
 301