ev_operations.c

   1 #include <assert.h>
   2 #include <math.h>
   3 #include <stdlib.h>
   4
   5 #include "ev_operations.h"
   6 #include "util.h"
   7
   8 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
   9
  10 void evalue_set_si(evalue *ev, int n, int d) {
  11     value_set_si(ev->d, d);
  12     value_init(ev->x.n);
  13     value_set_si(ev->x.n, n);
  14 }
  15
  16 void evalue_set(evalue *ev, Value n, Value d) {
  17     value_assign(ev->d, d);
  18     value_init(ev->x.n);
  19     value_assign(ev->x.n, n);
  20 }
  21
  22 void aep_evalue(evalue *e, int *ref) {
  23
  24     enode *p;
  25     int i;
  26
  27     if (value_notzero_p(e->d))
  28         return;         /* a rational number, its already reduced */
  29     if(!(p = e->x.p))
  30         return;         /* hum... an overflow probably occured */
  31
  32     /* First check the components of p */
  33     for (i=0;i<p->size;i++)
  34         aep_evalue(&p->arr[i],ref);
  35
  36     /* Then p itself */
  37     switch (p->type) {
  38     case polynomial:
  39     case periodic:
  40     case evector:
  41         p->pos = ref[p->pos-1]+1;
  42     }
  43     return;
  44 } /* aep_evalue */
  45
  46 /** Comments */
  47 void addeliminatedparams_evalue(evalue *e,Matrix *CT) {
  48
  49     enode *p;
  50     int i, j;
  51     int *ref;
  52
  53     if (value_notzero_p(e->d))
  54         return;          /* a rational number, its already reduced */
  55     if(!(p = e->x.p))
  56         return;          /* hum... an overflow probably occured */
  57
  58     /* Compute ref */
  59     ref = (int *)malloc(sizeof(int)*(CT->NbRows-1));
  60     for(i=0;i<CT->NbRows-1;i++)
  61         for(j=0;j<CT->NbColumns;j++)
  62             if(value_notzero_p(CT->p[i][j])) {
  63                 ref[i] = j;
  64                 break;
  65             }
  66
  67     /* Transform the references in e, using ref */
  68     aep_evalue(e,ref);
  69     free( ref );
  70     return;
  71 } /* addeliminatedparams_evalue */
  72
  73 static int mod_rational_smaller(evalue *e1, evalue *e2)
  74 {
  75     int r;
  76     Value m;
  77     value_init(m);
  78
  79     assert(value_notzero_p(e1->d));
  80     assert(value_notzero_p(e2->d));
  81     value_multiply(m, e1->x.n, e2->d);
  82     value_division(m, m, e1->d);
  83     if (value_lt(m, e2->x.n))
  84         r = 1;
  85     else if (value_gt(m, e2->x.n))
  86         r = 0;
  87     else
  88         r = -1;
  89     value_clear(m);
  90
  91     return r;
  92 }
  93
  94 static int mod_term_smaller_r(evalue *e1, evalue *e2)
  95 {
  96     if (value_notzero_p(e1->d)) {
  97         if (value_zero_p(e2->d))
  98             return 1;
  99         int r = mod_rational_smaller(e1, e2);
 100         return r == -1 ? 0 : r;
 101     }
 102     if (value_notzero_p(e2->d))
 103         return 0;
 104     if (e1->x.p->pos < e2->x.p->pos)
 105         return 1;
 106     else if (e1->x.p->pos > e2->x.p->pos)
 107         return 0;
 108     else {
 109         int r = mod_rational_smaller(&e1->x.p->arr[1], &e2->x.p->arr[1]);
 110         return r == -1
 111                  ? mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0])
 112                  : r;
 113     }
 114 }
 115
 116 static int mod_term_smaller(evalue *e1, evalue *e2)
 117 {
 118     assert(value_zero_p(e1->d));
 119     assert(value_zero_p(e2->d));
 120     assert(e1->x.p->type == fractional);
 121     assert(e2->x.p->type == fractional);
 122     return mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0]);
 123 }
 124
 125 /* Negative pos means inequality */
 126 /* s is negative of substitution if m is not zero */
 127 struct fixed_param {
 128     int     pos;
 129     evalue  s;
 130     Value   d;
 131     Value   m;
 132 };
 133
 134 struct subst {
 135     struct fixed_param *fixed;
 136     int                 n;
 137     int                 max;
 138 };
 139
 140 static int relations_depth(evalue *e)
 141 {
 142     int d;
 143
 144     for (d = 0;
 145          value_zero_p(e->d) && e->x.p->type == relation;
 146          e = &e->x.p->arr[1], ++d);
 147     return d;
 148 }
 149
 150 static void Lcm3(Value i, Value j, Value *res)
 151 {
 152     Value aux;
 153
 154     value_init(aux);
 155     Gcd(i,j,&aux);
 156     value_multiply(*res,i,j);
 157     value_division(*res, *res, aux);
 158     value_clear(aux);
 159 }
 160
 161 static void poly_denom(evalue *p, Value *d)
 162 {
 163     value_set_si(*d, 1);
 164
 165     while (value_zero_p(p->d)) {
 166         assert(p->x.p->type == polynomial);
 167         assert(p->x.p->size == 2);
 168         assert(value_notzero_p(p->x.p->arr[1].d));
 169         Lcm3(*d, p->x.p->arr[1].d, d);
 170         p = &p->x.p->arr[0];
 171     }
 172     Lcm3(*d, p->d, d);
 173 }
 174
 175 #define EVALUE_IS_ONE(ev)       (value_pos_p((ev).d) && value_one_p((ev).x.n))
 176
 177 static void realloc_substitution(struct subst *s, int d)
 178 {
 179     struct fixed_param *n;
 180     int i;
 181     NALLOC(n, s->max+d);
 182     for (i = 0; i < s->n; ++i)
 183         n[i] = s->fixed[i];
 184     free(s->fixed);
 185     s->fixed = n;
 186     s->max += d;
 187 }
 188
 189 static int add_modulo_substitution(struct subst *s, evalue *r)
 190 {
 191     evalue *p;
 192     evalue *f;
 193     evalue *m;
 194
 195     assert(value_zero_p(r->d) && r->x.p->type == relation);
 196     m = &r->x.p->arr[0];
 197
 198     /* May have been reduced already */
 199     if (value_notzero_p(m->d))
 200         return 0;
 201
 202     assert(value_zero_p(m->d) && m->x.p->type == fractional);
 203     assert(m->x.p->size == 3);
 204
 205     /* fractional was inverted during reduction
 206      * invert it back and move constant in
 207      */
 208     if (!EVALUE_IS_ONE(m->x.p->arr[2])) {
 209         assert(value_pos_p(m->x.p->arr[2].d));
 210         assert(value_mone_p(m->x.p->arr[2].x.n));
 211         value_set_si(m->x.p->arr[2].x.n, 1);
 212         value_increment(m->x.p->arr[1].x.n, m->x.p->arr[1].x.n);
 213         assert(value_eq(m->x.p->arr[1].x.n, m->x.p->arr[1].d));
 214         value_set_si(m->x.p->arr[1].x.n, 1);
 215         eadd(&m->x.p->arr[1], &m->x.p->arr[0]);
 216         value_set_si(m->x.p->arr[1].x.n, 0);
 217         value_set_si(m->x.p->arr[1].d, 1);
 218     }
 219
 220     /* Oops.  Nested identical relations. */
 221     if (!EVALUE_IS_ZERO(m->x.p->arr[1]))
 222         return 0;
 223
 224     if (s->n >= s->max) {
 225         int d = relations_depth(r);
 226         realloc_substitution(s, d);
 227     }
 228
 229     p = &m->x.p->arr[0];
 230     assert(value_zero_p(p->d) && p->x.p->type == polynomial);
 231     assert(p->x.p->size == 2);
 232     f = &p->x.p->arr[1];
 233
 234     assert(value_pos_p(f->x.n));
 235
 236     value_init(s->fixed[s->n].m);
 237     value_assign(s->fixed[s->n].m, f->d);
 238     s->fixed[s->n].pos = p->x.p->pos;
 239     value_init(s->fixed[s->n].d);
 240     value_assign(s->fixed[s->n].d, f->x.n);
 241     value_init(s->fixed[s->n].s.d);
 242     evalue_copy(&s->fixed[s->n].s, &p->x.p->arr[0]);
 243     ++s->n;
 244
 245     return 1;
 246 }
 247
 248 void _reduce_evalue (evalue *e, struct subst *s, int fract) {
 249
 250     enode *p;
 251     int i, j, k;
 252     int add;
 253
 254     if (value_notzero_p(e->d)) {
 255         if (fract)
 256             mpz_fdiv_r(e->x.n, e->x.n, e->d);
 257         return; /* a rational number, its already reduced */
 258     }
 259
 260     if(!(p = e->x.p))
 261         return; /* hum... an overflow probably occured */
 262
 263     /* First reduce the components of p */
 264     add = p->type == relation;
 265     for (i=0; i<p->size; i++) {
 266         if (add && i == 1)
 267             add = add_modulo_substitution(s, e);
 268
 269         if (i == 0 && p->type==fractional)
 270             _reduce_evalue(&p->arr[i], s, 1);
 271         else
 272             _reduce_evalue(&p->arr[i], s, fract);
 273
 274         if (add && i == p->size-1) {
 275             --s->n;
 276             value_clear(s->fixed[s->n].m);
 277             value_clear(s->fixed[s->n].d);
 278             free_evalue_refs(&s->fixed[s->n].s);
 279         } else if (add && i == 1)
 280             s->fixed[s->n-1].pos *= -1;
 281     }
 282
 283     if (p->type==periodic) {
 284
 285         /* Try to reduce the period */
 286         for (i=1; i<=(p->size)/2; i++) {
 287             if ((p->size % i)==0) {
 288
 289                 /* Can we reduce the size to i ? */
 290                 for (j=0; j<i; j++)
 291                     for (k=j+i; k<e->x.p->size; k+=i)
 292                         if (!eequal(&p->arr[j], &p->arr[k])) goto you_lose;
 293
 294                 /* OK, lets do it */
 295                 for (j=i; j<p->size; j++) free_evalue_refs(&p->arr[j]);
 296                 p->size = i;
 297                 break;
 298
 299 you_lose:       /* OK, lets not do it */
 300                 continue;
 301             }
 302         }
 303
 304         /* Try to reduce its strength */
 305         if (p->size == 1) {
 306             value_clear(e->d);
 307             memcpy(e,&p->arr[0],sizeof(evalue));
 308             free(p);
 309         }
 310     }
 311     else if (p->type==polynomial) {
 312         for (k = 0; s && k < s->n; ++k) {
 313             if (s->fixed[k].pos == p->pos) {
 314                 int divide = value_notone_p(s->fixed[k].d);
 315                 evalue d;
 316
 317                 if (value_notzero_p(s->fixed[k].m)) {
 318                     if (!fract)
 319                         continue;
 320                     assert(p->size == 2);
 321                     if (divide && value_ne(s->fixed[k].d, p->arr[1].x.n))
 322                         continue;
 323                     if (!mpz_divisible_p(s->fixed[k].m, p->arr[1].d))
 324                         continue;
 325                     divide = 1;
 326                 }
 327
 328                 if (divide) {
 329                     value_init(d.d);
 330                     value_assign(d.d, s->fixed[k].d);
 331                     value_init(d.x.n);
 332                     if (value_notzero_p(s->fixed[k].m))
 333                         value_oppose(d.x.n, s->fixed[k].m);
 334                     else
 335                         value_set_si(d.x.n, 1);
 336                 }
 337
 338                 for (i=p->size-1;i>=1;i--) {
 339                     emul(&s->fixed[k].s, &p->arr[i]);
 340                     if (divide)
 341                         emul(&d, &p->arr[i]);
 342                     eadd(&p->arr[i], &p->arr[i-1]);
 343                     free_evalue_refs(&(p->arr[i]));
 344                 }
 345                 p->size = 1;
 346                 _reduce_evalue(&p->arr[0], s, fract);
 347
 348                 if (divide)
 349                     free_evalue_refs(&d);
 350
 351                 break;
 352             }
 353         }
 354
 355         /* Try to reduce the degree */
 356         for (i=p->size-1;i>=1;i--) {
 357             if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
 358                 break;
 359             /* Zero coefficient */
 360             free_evalue_refs(&(p->arr[i]));
 361         }
 362         if (i+1<p->size)
 363             p->size = i+1;
 364
 365         /* Try to reduce its strength */
 366         if (p->size == 1) {
 367             value_clear(e->d);
 368             memcpy(e,&p->arr[0],sizeof(evalue));
 369             free(p);
 370         }
 371     }
 372     else if (p->type==fractional) {
 373         int reorder = 0;
 374         evalue v;
 375
 376         if (value_notzero_p(p->arr[0].d)) {
 377             value_init(v.d);
 378             value_assign(v.d, p->arr[0].d);
 379             value_init(v.x.n);
 380             mpz_fdiv_r(v.x.n, p->arr[0].x.n,  p->arr[0].d);
 381
 382             reorder = 1;
 383         } else {
 384             evalue *f, *base;
 385             evalue *pp = &p->arr[0];
 386             assert(value_zero_p(pp->d) && pp->x.p->type == polynomial);
 387             assert(pp->x.p->size == 2);
 388
 389             /* search for exact duplicate among the modulo inequalities */
 390             do {
 391                 f = &pp->x.p->arr[1];
 392                 for (k = 0; s && k < s->n; ++k) {
 393                     if (-s->fixed[k].pos == pp->x.p->pos &&
 394                             value_eq(s->fixed[k].d, f->x.n) &&
 395                             value_eq(s->fixed[k].m, f->d) &&
 396                             eequal(&s->fixed[k].s, &pp->x.p->arr[0]))
 397                         break;
 398                 }
 399                 if (k < s->n) {
 400                     Value g;
 401                     value_init(g);
 402
 403                     /* replace { E/m } by { (E-1)/m } + 1/m */
 404                     poly_denom(pp, &g);
 405                     if (reorder) {
 406                         evalue extra;
 407                         value_init(extra.d);
 408                         evalue_set_si(&extra, 1, 1);
 409                         value_assign(extra.d, g);
 410                         eadd(&extra, &v.x.p->arr[1]);
 411                         free_evalue_refs(&extra);
 412
 413                         /* We've been going in circles; stop now */
 414                         if (value_ge(v.x.p->arr[1].x.n, v.x.p->arr[1].d)) {
 415                             free_evalue_refs(&v);
 416                             value_init(v.d);
 417                             evalue_set_si(&v, 0, 1);
 418                             break;
 419                         }
 420                     } else {
 421                         value_init(v.d);
 422                         value_set_si(v.d, 0);
 423                         v.x.p = new_enode(fractional, 3, -1);
 424                         evalue_set_si(&v.x.p->arr[1], 1, 1);
 425                         value_assign(v.x.p->arr[1].d, g);
 426                         evalue_set_si(&v.x.p->arr[2], 1, 1);
 427                         evalue_copy(&v.x.p->arr[0], &p->arr[0]);
 428                     }
 429
 430                     for (f = &v.x.p->arr[0]; value_zero_p(f->d);
 431                                              f = &f->x.p->arr[0])
 432                         ;
 433                     value_division(f->d, g, f->d);
 434                     value_multiply(f->x.n, f->x.n, f->d);
 435                     value_assign(f->d, g);
 436                     value_decrement(f->x.n, f->x.n);
 437                     mpz_fdiv_r(f->x.n, f->x.n, f->d);
 438
 439                     Gcd(f->d, f->x.n, &g);
 440                     value_division(f->d, f->d, g);
 441                     value_division(f->x.n, f->x.n, g);
 442
 443                     value_clear(g);
 444                     pp = &v.x.p->arr[0];
 445
 446                     reorder = 1;
 447                 }
 448             } while (k < s->n);
 449
 450             /* reduction may have made this fractional arg smaller */
 451             i = reorder ? p->size : 1;
 452             for ( ; i < p->size; ++i)
 453                 if (value_zero_p(p->arr[i].d) &&
 454                         p->arr[i].x.p->type == fractional &&
 455                         !mod_term_smaller(e, &p->arr[i]))
 456                     break;
 457             if (i < p->size) {
 458                 value_init(v.d);
 459                 value_set_si(v.d, 0);
 460                 v.x.p = new_enode(fractional, 3, -1);
 461                 evalue_set_si(&v.x.p->arr[1], 0, 1);
 462                 evalue_set_si(&v.x.p->arr[2], 1, 1);
 463                 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
 464
 465                 reorder = 1;
 466             }
 467
 468             if (!reorder) {
 469                 int invert = 0;
 470                 Value twice;
 471                 value_init(twice);
 472
 473                 for (pp = &p->arr[0]; value_zero_p(pp->d);
 474                                       pp = &pp->x.p->arr[0]) {
 475                     f = &pp->x.p->arr[1];
 476                     assert(value_pos_p(f->d));
 477                     mpz_mul_ui(twice, f->x.n, 2);
 478                     if (value_lt(twice, f->d))
 479                         break;
 480                     if (value_eq(twice, f->d))
 481                         continue;
 482                     invert = 1;
 483                     break;
 484                 }
 485
 486                 if (invert) {
 487                     value_init(v.d);
 488                     value_set_si(v.d, 0);
 489                     v.x.p = new_enode(fractional, 3, -1);
 490                     evalue_set_si(&v.x.p->arr[1], 0, 1);
 491                     poly_denom(&p->arr[0], &twice);
 492                     value_assign(v.x.p->arr[1].d, twice);
 493                     value_decrement(v.x.p->arr[1].x.n, twice);
 494                     evalue_set_si(&v.x.p->arr[2], -1, 1);
 495                     evalue_copy(&v.x.p->arr[0], &p->arr[0]);
 496
 497                     for (pp = &v.x.p->arr[0]; value_zero_p(pp->d);
 498                                               pp = &pp->x.p->arr[0]) {
 499                         f = &pp->x.p->arr[1];
 500                         value_oppose(f->x.n, f->x.n);
 501                         mpz_fdiv_r(f->x.n, f->x.n,  f->d);
 502                     }
 503                     value_division(pp->d, twice, pp->d);
 504                     value_multiply(pp->x.n, pp->x.n, pp->d);
 505                     value_assign(pp->d, twice);
 506                     value_oppose(pp->x.n, pp->x.n);
 507                     value_decrement(pp->x.n, pp->x.n);
 508                     mpz_fdiv_r(pp->x.n, pp->x.n, pp->d);
 509
 510                     reorder = 1;
 511                 }
 512
 513                 value_clear(twice);
 514             }
 515         }
 516
 517         if (reorder) {
 518             for (i=p->size-1;i>=2;i--) {
 519                 emul(&v, &p->arr[i]);
 520                 eadd(&p->arr[i], &p->arr[i-1]);
 521                 free_evalue_refs(&(p->arr[i]));
 522             }
 523             p->size = 2;
 524             free_evalue_refs(&v);
 525             _reduce_evalue(&p->arr[1], s, fract);
 526         }
 527
 528         /* Try to reduce the degree */
 529         for (i=p->size-1;i>=2;i--) {
 530             if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
 531                 break;
 532             /* Zero coefficient */
 533             free_evalue_refs(&(p->arr[i]));
 534         }
 535         if (i+1<p->size)
 536             p->size = i+1;
 537
 538         /* Try to reduce its strength */
 539         if (p->size == 2) {
 540             value_clear(e->d);
 541             memcpy(e,&p->arr[1],sizeof(evalue));
 542             free_evalue_refs(&(p->arr[0]));
 543             free(p);
 544         }
 545     }
 546     else if (p->type == relation) {
 547         if (p->size == 3 && eequal(&p->arr[1], &p->arr[2])) {
 548             free_evalue_refs(&(p->arr[2]));
 549             free_evalue_refs(&(p->arr[0]));
 550             p->size = 2;
 551             value_clear(e->d);
 552             *e = p->arr[1];
 553             free(p);
 554             return;
 555         }
 556         if (p->size == 3 && EVALUE_IS_ZERO(p->arr[2])) {
 557             free_evalue_refs(&(p->arr[2]));
 558             p->size = 2;
 559         }
 560         if (p->size == 2 && EVALUE_IS_ZERO(p->arr[1])) {
 561             free_evalue_refs(&(p->arr[1]));
 562             free_evalue_refs(&(p->arr[0]));
 563             evalue_set_si(e, 0, 1);
 564             free(p);
 565         } else {
 566             int reduced = 0;
 567             evalue *m;
 568             m = &p->arr[0];
 569
 570             /* Relation was reduced by means of an identical
 571              * inequality => remove
 572              */
 573             if (value_zero_p(m->d) && !EVALUE_IS_ZERO(m->x.p->arr[1]))
 574                 reduced = 1;
 575
 576             if (reduced || value_notzero_p(p->arr[0].d)) {
 577                 if (!reduced && value_zero_p(p->arr[0].x.n)) {
 578                     value_clear(e->d);
 579                     memcpy(e,&p->arr[1],sizeof(evalue));
 580                     if (p->size == 3)
 581                         free_evalue_refs(&(p->arr[2]));
 582                 } else {
 583                     if (p->size == 3) {
 584                         value_clear(e->d);
 585                         memcpy(e,&p->arr[2],sizeof(evalue));
 586                     } else
 587                         evalue_set_si(e, 0, 1);
 588                     free_evalue_refs(&(p->arr[1]));
 589                 }
 590                 free_evalue_refs(&(p->arr[0]));
 591                 free(p);
 592             }
 593         }
 594     }
 595     return;
 596 } /* reduce_evalue */
 597
 598 static void add_substitution(struct subst *s, Value *row, unsigned dim)
 599 {
 600     int k, l;
 601     evalue *r;
 602
 603     for (k = 0; k < dim; ++k)
 604         if (value_notzero_p(row[k+1]))
 605             break;
 606
 607     Vector_Normalize_Positive(row+1, dim+1, k);
 608     assert(s->n < s->max);
 609     value_init(s->fixed[s->n].d);
 610     value_init(s->fixed[s->n].m);
 611     value_assign(s->fixed[s->n].d, row[k+1]);
 612     s->fixed[s->n].pos = k+1;
 613     value_set_si(s->fixed[s->n].m, 0);
 614     r = &s->fixed[s->n].s;
 615     value_init(r->d);
 616     for (l = k+1; l < dim; ++l)
 617         if (value_notzero_p(row[l+1])) {
 618             value_set_si(r->d, 0);
 619             r->x.p = new_enode(polynomial, 2, l + 1);
 620             value_init(r->x.p->arr[1].x.n);
 621             value_oppose(r->x.p->arr[1].x.n, row[l+1]);
 622             value_set_si(r->x.p->arr[1].d, 1);
 623             r = &r->x.p->arr[0];
 624         }
 625     value_init(r->x.n);
 626     value_oppose(r->x.n, row[dim+1]);
 627     value_set_si(r->d, 1);
 628     ++s->n;
 629 }
 630
 631 void reduce_evalue (evalue *e) {
 632     struct subst s = { NULL, 0, 0 };
 633
 634     if (value_notzero_p(e->d))
 635         return; /* a rational number, its already reduced */
 636
 637     if (e->x.p->type == partition) {
 638         int i;
 639         unsigned dim = -1;
 640         for (i = 0; i < e->x.p->size/2; ++i) {
 641             s.n = 0;
 642             Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
 643             /* This shouldn't really happen;
 644              * Empty domains should not be added.
 645              */
 646             if (emptyQ(D))
 647                 goto discard;
 648
 649             dim = D->Dimension;
 650             if (D->next)
 651                 D = DomainConvex(D, 0);
 652             if (!D->next && D->NbEq) {
 653                 int j, k;
 654                 if (s.max < dim) {
 655                     if (s.max != 0)
 656                         realloc_substitution(&s, dim);
 657                     else {
 658                         int d = relations_depth(&e->x.p->arr[2*i+1]);
 659                         s.max = dim+d;
 660                         NALLOC(s.fixed, s.max);
 661                     }
 662                 }
 663                 for (j = 0; j < D->NbEq; ++j)
 664                     add_substitution(&s, D->Constraint[j], dim);
 665             }
 666             if (D != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 667                 Domain_Free(D);
 668             _reduce_evalue(&e->x.p->arr[2*i+1], &s, 0);
 669             if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
 670 discard:
 671                 free_evalue_refs(&e->x.p->arr[2*i+1]);
 672                 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
 673                 value_clear(e->x.p->arr[2*i].d);
 674                 e->x.p->size -= 2;
 675                 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
 676                 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
 677                 --i;
 678             }
 679             if (s.n != 0) {
 680                 int j;
 681                 for (j = 0; j < s.n; ++j) {
 682                     value_clear(s.fixed[j].d);
 683                     free_evalue_refs(&s.fixed[j].s);
 684                 }
 685             }
 686         }
 687         if (e->x.p->size == 0) {
 688             free(e->x.p);
 689             evalue_set_si(e, 0, 1);
 690         }
 691         if (s.max != 0)
 692             free(s.fixed);
 693     } else
 694         _reduce_evalue(e, &s, 0);
 695 }
 696
 697 void print_evalue(FILE *DST,evalue *e,char **pname) {
 698
 699   if(value_notzero_p(e->d)) {
 700     if(value_notone_p(e->d)) {
 701       value_print(DST,VALUE_FMT,e->x.n);
 702       fprintf(DST,"/");
 703       value_print(DST,VALUE_FMT,e->d);
 704     }
 705     else {
 706       value_print(DST,VALUE_FMT,e->x.n);
 707     }
 708   }
 709   else
 710     print_enode(DST,e->x.p,pname);
 711   return;
 712 } /* print_evalue */
 713
 714 void print_enode(FILE *DST,enode *p,char **pname) {
 715
 716   int i;
 717
 718   if (!p) {
 719     fprintf(DST, "NULL");
 720     return;
 721   }
 722   switch (p->type) {
 723   case evector:
 724     fprintf(DST, "{ ");
 725     for (i=0; i<p->size; i++) {
 726       print_evalue(DST, &p->arr[i], pname);
 727       if (i!=(p->size-1))
 728         fprintf(DST, ", ");
 729     }
 730     fprintf(DST, " }\n");
 731     break;
 732   case polynomial:
 733     fprintf(DST, "( ");
 734     for (i=p->size-1; i>=0; i--) {
 735       print_evalue(DST, &p->arr[i], pname);
 736       if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
 737       else if (i>1)
 738         fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
 739     }
 740     fprintf(DST, " )\n");
 741     break;
 742   case periodic:
 743     fprintf(DST, "[ ");
 744     for (i=0; i<p->size; i++) {
 745       print_evalue(DST, &p->arr[i], pname);
 746       if (i!=(p->size-1)) fprintf(DST, ", ");
 747     }
 748     fprintf(DST," ]_%s", pname[p->pos-1]);
 749     break;
 750   case fractional:
 751     fprintf(DST, "( ");
 752     for (i=p->size-1; i>=1; i--) {
 753       print_evalue(DST, &p->arr[i], pname);
 754       if (i >= 2) {
 755         fprintf(DST, " * ");
 756         fprintf(DST, "{");
 757         print_evalue(DST, &p->arr[0], pname);
 758         fprintf(DST, "}");
 759         if (i>2)
 760           fprintf(DST, "^%d + ", i-1);
 761         else
 762           fprintf(DST, " + ");
 763       }
 764     }
 765     fprintf(DST, " )\n");
 766     break;
 767   case relation:
 768     fprintf(DST, "[ ");
 769     print_evalue(DST, &p->arr[0], pname);
 770     fprintf(DST, "= 0 ] * \n");
 771     print_evalue(DST, &p->arr[1], pname);
 772     if (p->size > 2) {
 773         fprintf(DST, " +\n [ ");
 774         print_evalue(DST, &p->arr[0], pname);
 775         fprintf(DST, "!= 0 ] * \n");
 776         print_evalue(DST, &p->arr[2], pname);
 777     }
 778     break;
 779   case partition:
 780     for (i=0; i<p->size/2; i++) {
 781         Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), pname);
 782         print_evalue(DST, &p->arr[2*i+1], pname);
 783     }
 784     break;
 785   default:
 786     assert(0);
 787   }
 788   return;
 789 } /* print_enode */
 790
 791 static void eadd_rev(evalue *e1, evalue *res)
 792 {
 793     evalue ev;
 794     value_init(ev.d);
 795     evalue_copy(&ev, e1);
 796     eadd(res, &ev);
 797     free_evalue_refs(res);
 798     *res = ev;
 799 }
 800
 801 static void eadd_rev_cst (evalue *e1, evalue *res)
 802 {
 803     evalue ev;
 804     value_init(ev.d);
 805     evalue_copy(&ev, e1);
 806     eadd(res, &ev.x.p->arr[ev.x.p->type==fractional]);
 807     free_evalue_refs(res);
 808     *res = ev;
 809 }
 810
 811 struct section { Polyhedron * D; evalue E; };
 812
 813 void eadd_partitions (evalue *e1,evalue *res)
 814 {
 815     int n, i, j;
 816     Polyhedron *d, *fd;
 817     struct section *s;
 818     s = (struct section *)
 819             malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
 820                    sizeof(struct section));
 821     assert(s);
 822
 823     n = 0;
 824     for (j = 0; j < e1->x.p->size/2; ++j) {
 825         assert(res->x.p->size >= 2);
 826         fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
 827                               EVALUE_DOMAIN(res->x.p->arr[0]), 0);
 828         if (!emptyQ(fd))
 829             for (i = 1; i < res->x.p->size/2; ++i) {
 830                 Polyhedron *t = fd;
 831                 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
 832                 Domain_Free(t);
 833                 if (emptyQ(fd))
 834                     break;
 835             }
 836         if (emptyQ(fd)) {
 837             Domain_Free(fd);
 838             continue;
 839         }
 840         value_init(s[n].E.d);
 841         evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
 842         s[n].D = fd;
 843         ++n;
 844     }
 845     for (i = 0; i < res->x.p->size/2; ++i) {
 846         fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
 847         for (j = 0; j < e1->x.p->size/2; ++j) {
 848             Polyhedron *t;
 849             d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
 850                                    EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
 851             if (emptyQ(d)) {
 852                 Domain_Free(d);
 853                 continue;
 854             }
 855             t = fd;
 856             fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
 857             if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
 858                 Domain_Free(t);
 859             value_init(s[n].E.d);
 860             evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
 861             eadd(&e1->x.p->arr[2*j+1], &s[n].E);
 862             s[n].D = d;
 863             ++n;
 864         }
 865         if (!emptyQ(fd)) {
 866             s[n].E = res->x.p->arr[2*i+1];
 867             s[n].D = fd;
 868             ++n;
 869         } else {
 870             free_evalue_refs(&res->x.p->arr[2*i+1]);
 871             Domain_Free(fd);
 872         }
 873         if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
 874             Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
 875         value_clear(res->x.p->arr[2*i].d);
 876     }
 877
 878     free(res->x.p);
 879     res->x.p = new_enode(partition, 2*n, -1);
 880     for (j = 0; j < n; ++j) {
 881         s[j].D = DomainConstraintSimplify(s[j].D, 0);
 882         EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
 883         value_clear(res->x.p->arr[2*j+1].d);
 884         res->x.p->arr[2*j+1] = s[j].E;
 885     }
 886
 887     free(s);
 888 }
 889
 890 static void explicit_complement(evalue *res)
 891 {
 892     enode *rel = new_enode(relation, 3, 0);
 893     assert(rel);
 894     value_clear(rel->arr[0].d);
 895     rel->arr[0] = res->x.p->arr[0];
 896     value_clear(rel->arr[1].d);
 897     rel->arr[1] = res->x.p->arr[1];
 898     value_set_si(rel->arr[2].d, 1);
 899     value_init(rel->arr[2].x.n);
 900     value_set_si(rel->arr[2].x.n, 0);
 901     free(res->x.p);
 902     res->x.p = rel;
 903 }
 904
 905 void eadd(evalue *e1,evalue *res) {
 906
 907  int i;
 908     if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
 909          /* Add two rational numbers */
 910          Value g,m1,m2;
 911          value_init(g);
 912          value_init(m1);
 913          value_init(m2);
 914
 915          value_multiply(m1,e1->x.n,res->d);
 916          value_multiply(m2,res->x.n,e1->d);
 917          value_addto(res->x.n,m1,m2);
 918          value_multiply(res->d,e1->d,res->d);
 919          Gcd(res->x.n,res->d,&g);
 920          if (value_notone_p(g)) {
 921               value_division(res->d,res->d,g);
 922               value_division(res->x.n,res->x.n,g);
 923          }
 924          value_clear(g); value_clear(m1); value_clear(m2);
 925          return ;
 926      }
 927      else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
 928           switch (res->x.p->type) {
 929           case polynomial:
 930               /* Add the constant to the constant term of a polynomial*/
 931                eadd(e1, &res->x.p->arr[0]);
 932                return ;
 933           case periodic:
 934               /* Add the constant to all elements of a periodic number */
 935               for (i=0; i<res->x.p->size; i++) {
 936                   eadd(e1, &res->x.p->arr[i]);
 937               }
 938               return ;
 939           case evector:
 940               fprintf(stderr, "eadd: cannot add const with vector\n");
 941               return;
 942           case fractional:
 943                eadd(e1, &res->x.p->arr[1]);
 944                return ;
 945           case partition:
 946                 assert(EVALUE_IS_ZERO(*e1));
 947                 break;                          /* Do nothing */
 948           case relation:
 949                 /* Create (zero) complement if needed */
 950                 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
 951                     explicit_complement(res);
 952                 for (i = 1; i < res->x.p->size; ++i)
 953                     eadd(e1, &res->x.p->arr[i]);
 954                 break;
 955           default:
 956                 assert(0);
 957           }
 958      }
 959      /* add polynomial or periodic to constant
 960       * you have to exchange e1 and res, before doing addition */
 961
 962      else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
 963           eadd_rev(e1, res);
 964           return;
 965      }
 966      else {   // ((e1->d==0) && (res->d==0))
 967         assert(!((e1->x.p->type == partition) ^
 968                  (res->x.p->type == partition)));
 969         if (e1->x.p->type == partition) {
 970             eadd_partitions(e1, res);
 971             return;
 972         }
 973         if (e1->x.p->type == relation &&
 974             (res->x.p->type != relation ||
 975              mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
 976                 eadd_rev(e1, res);
 977                 return;
 978         }
 979         if (res->x.p->type == relation) {
 980             if (e1->x.p->type == relation &&
 981                 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
 982                     if (res->x.p->size < 3 && e1->x.p->size == 3)
 983                         explicit_complement(res);
 984                     if (e1->x.p->size < 3 && res->x.p->size == 3)
 985                         explicit_complement(e1);
 986                     for (i = 1; i < res->x.p->size; ++i)
 987                         eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
 988                     return;
 989             }
 990             if (res->x.p->size < 3)
 991                 explicit_complement(res);
 992             for (i = 1; i < res->x.p->size; ++i)
 993                 eadd(e1, &res->x.p->arr[i]);
 994             return;
 995         }
 996                  if ((e1->x.p->type != res->x.p->type) ) {
 997                       /* adding to evalues of different type. two cases are possible
 998                        * res is periodic and e1 is polynomial, you have to exchange
 999                        * e1 and res then to add e1 to the constant term of res */
1000                      if (e1->x.p->type == polynomial) {
1001                           eadd_rev_cst(e1, res);
1002                      }
1003                      else if (res->x.p->type == polynomial) {
1004                           /* res is polynomial and e1 is periodic,
1005                             add e1 to the constant term of res */
1006
1007                           eadd(e1,&res->x.p->arr[0]);
1008                      } else
1009                         assert(0);
1010
1011                      return;
1012                  }
1013                  else if (e1->x.p->pos != res->x.p->pos ||
1014                             (res->x.p->type == fractional &&
1015                              !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1016                  /* adding evalues of different position (i.e function of different unknowns
1017                   * to case are possible  */
1018
1019                         switch (res->x.p->type) {
1020                         case fractional:
1021                             if(mod_term_smaller(res, e1))
1022                                 eadd(e1,&res->x.p->arr[1]);
1023                             else
1024                                 eadd_rev_cst(e1, res);
1025                             return;
1026                         case polynomial: //  res and e1 are polynomials
1027                                //  add e1 to the constant term of res
1028
1029                             if(res->x.p->pos < e1->x.p->pos)
1030                                eadd(e1,&res->x.p->arr[0]);
1031                             else
1032                                 eadd_rev_cst(e1, res);
1033                               // value_clear(g); value_clear(m1); value_clear(m2);
1034                                return;
1035                         case periodic:  // res and e1 are pointers to periodic numbers
1036                                   //add e1 to all elements of res
1037
1038                             if(res->x.p->pos < e1->x.p->pos)
1039                                   for (i=0;i<res->x.p->size;i++) {
1040                                        eadd(e1,&res->x.p->arr[i]);
1041                                   }
1042                             else
1043                                 eadd_rev(e1, res);
1044                             return;
1045                         }
1046                  }
1047
1048
1049                  //same type , same pos  and same size
1050                  if (e1->x.p->size == res->x.p->size) {
1051                       // add any element in e1 to the corresponding element in res
1052                       if (res->x.p->type == fractional)
1053                         assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1054                       i = res->x.p->type == fractional ? 1 : 0;
1055                       for (; i<res->x.p->size; i++) {
1056                             eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1057                       }
1058                       return ;
1059                 }
1060
1061                 /* Sizes are different */
1062                 switch(res->x.p->type) {
1063                 case polynomial:
1064                 case fractional:
1065                     /* VIN100: if e1-size > res-size you have to copy e1 in a   */
1066                     /* new enode and add res to that new node. If you do not do */
1067                     /* that, you lose the the upper weight part of e1 !         */
1068
1069                      if(e1->x.p->size > res->x.p->size)
1070                           eadd_rev(e1, res);
1071                      else {
1072
1073                         if (res->x.p->type == fractional)
1074                             assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1075                         i = res->x.p->type == fractional ? 1 : 0;
1076                         for (; i<e1->x.p->size ; i++) {
1077                              eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1078                         }
1079
1080                         return ;
1081                      }
1082                      break;
1083
1084     /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1085                 case periodic:
1086                 {
1087                       /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1088                        of the sizes of e1 and res, then to copy res periodicaly in ne, after
1089                        to add periodicaly elements of e1 to elements of ne, and finaly to
1090                        return ne. */
1091                        int x,y,p;
1092                        Value ex, ey ,ep;
1093                        evalue *ne;
1094                        value_init(ex); value_init(ey);value_init(ep);
1095                        x=e1->x.p->size;
1096                        y= res->x.p->size;
1097                        value_set_si(ex,e1->x.p->size);
1098                        value_set_si(ey,res->x.p->size);
1099                        value_assign (ep,*Lcm(ex,ey));
1100                        p=(int)mpz_get_si(ep);
1101                        ne= (evalue *) malloc (sizeof(evalue));
1102                        value_init(ne->d);
1103                        value_set_si( ne->d,0);
1104
1105                        ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
1106                        for(i=0;i<p;i++)  {
1107                               value_assign(ne->x.p->arr[i].d, res->x.p->arr[i%y].d);
1108                               if (value_notzero_p(ne->x.p->arr[i].d))   {
1109                                 value_init(ne->x.p->arr[i].x.n);
1110                                 value_assign(ne->x.p->arr[i].x.n, res->x.p->arr[i%y].x.n);
1111                               }
1112                               else {
1113                                   ne->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%y].x.p);
1114                               }
1115                        }
1116                        for(i=0;i<p;i++)  {
1117                             eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
1118                        }
1119
1120                       value_assign(res->d, ne->d);
1121                       res->x.p=ne->x.p;
1122
1123                         return ;
1124                 }
1125                 case evector:
1126                      fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
1127                      return ;
1128                 }
1129      }
1130      return ;
1131  }/* eadd  */
1132
1133 static void emul_rev (evalue *e1, evalue *res)
1134 {
1135     evalue ev;
1136     value_init(ev.d);
1137     evalue_copy(&ev, e1);
1138     emul(res, &ev);
1139     free_evalue_refs(res);
1140     *res = ev;
1141 }
1142
1143 static void emul_poly (evalue *e1, evalue *res)
1144 {
1145     int i, j, o = res->x.p->type == fractional;
1146     evalue tmp;
1147     int size=(e1->x.p->size + res->x.p->size - o - 1);
1148     value_init(tmp.d);
1149     value_set_si(tmp.d,0);
1150     tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
1151     if (o)
1152         evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
1153     for (i=o; i < e1->x.p->size; i++) {
1154         evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
1155         emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
1156     }
1157     for (; i<size; i++)
1158         evalue_set_si(&tmp.x.p->arr[i], 0, 1);
1159     for (i=o+1; i<res->x.p->size; i++)
1160         for (j=o; j<e1->x.p->size; j++) {
1161             evalue ev;
1162             value_init(ev.d);
1163             evalue_copy(&ev, &e1->x.p->arr[j]);
1164             emul(&res->x.p->arr[i], &ev);
1165             eadd(&ev, &tmp.x.p->arr[i+j-o]);
1166             free_evalue_refs(&ev);
1167         }
1168     free_evalue_refs(res);
1169     *res = tmp;
1170 }
1171
1172 void emul_partitions (evalue *e1,evalue *res)
1173 {
1174     int n, i, j, k;
1175     Polyhedron *d;
1176     struct section *s;
1177     s = (struct section *)
1178             malloc((e1->x.p->size/2) * (res->x.p->size/2) *
1179                    sizeof(struct section));
1180     assert(s);
1181
1182     n = 0;
1183     for (i = 0; i < res->x.p->size/2; ++i) {
1184         for (j = 0; j < e1->x.p->size/2; ++j) {
1185             d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1186                                    EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1187             if (emptyQ(d)) {
1188                 Domain_Free(d);
1189                 continue;
1190             }
1191
1192             /* This code is only needed because the partitions
1193                are not true partitions.
1194              */
1195             for (k = 0; k < n; ++k) {
1196                 if (DomainIncludes(s[k].D, d))
1197                     break;
1198                 if (DomainIncludes(d, s[k].D)) {
1199                     Domain_Free(s[k].D);
1200                     free_evalue_refs(&s[k].E);
1201                     if (n > k)
1202                         s[k] = s[--n];
1203                     --k;
1204                 }
1205             }
1206             if (k < n) {
1207                 Domain_Free(d);
1208                 continue;
1209             }
1210
1211             value_init(s[n].E.d);
1212             evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1213             emul(&e1->x.p->arr[2*j+1], &s[n].E);
1214             s[n].D = d;
1215             ++n;
1216         }
1217         Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1218         value_clear(res->x.p->arr[2*i].d);
1219         free_evalue_refs(&res->x.p->arr[2*i+1]);
1220     }
1221
1222     free(res->x.p);
1223     res->x.p = new_enode(partition, 2*n, -1);
1224     for (j = 0; j < n; ++j) {
1225         s[j].D = DomainConstraintSimplify(s[j].D, 0);
1226         EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1227         value_clear(res->x.p->arr[2*j+1].d);
1228         res->x.p->arr[2*j+1] = s[j].E;
1229     }
1230
1231     free(s);
1232 }
1233
1234 #define value_two_p(val)        (mpz_cmp_si(val,2) == 0)
1235
1236 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1237  * do a copy of "res" befor calling this function if you nead it after. The vector type of
1238  * evalues is not treated here */
1239
1240 void emul (evalue *e1, evalue *res ){
1241     int i,j;
1242
1243 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
1244     fprintf(stderr, "emul: do not proced on evector type !\n");
1245     return;
1246 }
1247
1248     if (value_zero_p(e1->d) && e1->x.p->type == partition) {
1249         if (value_zero_p(res->d) && res->x.p->type == partition)
1250             emul_partitions(e1, res);
1251         else
1252             emul_rev(e1, res);
1253     } else if (value_zero_p(res->d) && res->x.p->type == partition) {
1254         for (i = 0; i < res->x.p->size/2; ++i)
1255             emul(e1, &res->x.p->arr[2*i+1]);
1256     } else
1257    if (value_zero_p(res->d) && res->x.p->type == relation) {
1258         if (value_zero_p(e1->d) && e1->x.p->type == relation &&
1259             eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1260                 if (res->x.p->size < 3 && e1->x.p->size == 3)
1261                     explicit_complement(res);
1262                 if (e1->x.p->size < 3 && res->x.p->size == 3)
1263                     explicit_complement(e1);
1264                 for (i = 1; i < res->x.p->size; ++i)
1265                     emul(&e1->x.p->arr[i], &res->x.p->arr[i]);
1266                 return;
1267         }
1268         for (i = 1; i < res->x.p->size; ++i)
1269             emul(e1, &res->x.p->arr[i]);
1270    } else
1271    if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
1272        switch(e1->x.p->type) {
1273        case polynomial:
1274            switch(res->x.p->type) {
1275            case polynomial:
1276                if(e1->x.p->pos == res->x.p->pos) {
1277                /* Product of two polynomials of the same variable */
1278                     emul_poly(e1, res);
1279                     return;
1280                }
1281                else {
1282                   /* Product of two polynomials of different variables */
1283
1284                 if(res->x.p->pos < e1->x.p->pos)
1285                      for( i=0; i<res->x.p->size ; i++)
1286                             emul(e1, &res->x.p->arr[i]);
1287                 else
1288                     emul_rev(e1, res);
1289
1290                  return ;
1291                }
1292            case periodic:
1293            case fractional:
1294                 /* Product of a polynomial and a periodic or fractional */
1295                 emul_rev(e1, res);
1296                 return;
1297            }
1298        case periodic:
1299            switch(res->x.p->type) {
1300            case periodic:
1301                 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
1302                  /* Product of two periodics of the same parameter and period */
1303
1304                      for(i=0; i<res->x.p->size;i++)
1305                            emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
1306
1307                      return;
1308                 }
1309                 else{
1310                   if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
1311                    /* Product of two periodics of the same parameter and different periods */
1312                     evalue *newp;
1313                     Value x,y,z;
1314                     int ix,iy,lcm;
1315                     value_init(x); value_init(y);value_init(z);
1316                     ix=e1->x.p->size;
1317                     iy=res->x.p->size;
1318                     value_set_si(x,e1->x.p->size);
1319                     value_set_si(y,res->x.p->size);
1320                     value_assign (z,*Lcm(x,y));
1321                     lcm=(int)mpz_get_si(z);
1322                     newp= (evalue *) malloc (sizeof(evalue));
1323                     value_init(newp->d);
1324                     value_set_si( newp->d,0);
1325                     newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
1326                         for(i=0;i<lcm;i++)  {
1327                            value_assign(newp->x.p->arr[i].d, res->x.p->arr[i%iy].d);
1328                            if (value_notzero_p(newp->x.p->arr[i].d))   {
1329                            value_assign(newp->x.p->arr[i].x.n, res->x.p->arr[i%iy].x.n);
1330                            }
1331                            else {
1332                            newp->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%iy].x.p);
1333                            }
1334
1335                         }
1336                         for(i=0;i<lcm;i++)
1337                             emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
1338
1339                         value_assign(res->d,newp->d);
1340                         res->x.p=newp->x.p;
1341
1342                           value_clear(x); value_clear(y);value_clear(z);
1343                           return ;
1344                   }
1345                   else {
1346                      /* Product of two periodics of different parameters */
1347
1348                         for(i=0; i<res->x.p->size; i++)
1349                             emul(e1, &(res->x.p->arr[i]));
1350
1351                         return;
1352                   }
1353                 }
1354            case polynomial:
1355                   /* Product of a periodic and a polynomial */
1356
1357                        for(i=0; i<res->x.p->size ; i++)
1358                             emul(e1, &(res->x.p->arr[i]));
1359
1360                        return;
1361
1362            }
1363        case fractional:
1364             switch(res->x.p->type) {
1365             case polynomial:
1366                 for(i=0; i<res->x.p->size ; i++)
1367                     emul(e1, &(res->x.p->arr[i]));
1368                 return;
1369             case periodic:
1370                 assert(0);
1371             case fractional:
1372                 if (e1->x.p->pos == res->x.p->pos &&
1373                             eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1374                     evalue d;
1375                     value_init(d.d);
1376                     poly_denom(&e1->x.p->arr[0], &d.d);
1377                     if (!value_two_p(d.d))
1378                         emul_poly(e1, res);
1379                     else {
1380                         value_init(d.x.n);
1381                         value_set_si(d.x.n, 1);
1382                         evalue tmp;
1383                         /* { x }^2 == { x }/2 */
1384                         /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1385                         assert(e1->x.p->size == 3);
1386                         assert(res->x.p->size == 3);
1387                         value_init(tmp.d);
1388                         evalue_copy(&tmp, &res->x.p->arr[2]);
1389                         emul(&d, &tmp);
1390                         eadd(&res->x.p->arr[1], &tmp);
1391                         emul(&e1->x.p->arr[2], &tmp);
1392                         emul(&e1->x.p->arr[1], res);
1393                         eadd(&tmp, &res->x.p->arr[2]);
1394                         free_evalue_refs(&tmp);
1395                         value_clear(d.x.n);
1396                     }
1397                     value_clear(d.d);
1398                 } else {
1399                     if(mod_term_smaller(res, e1))
1400                         for(i=1; i<res->x.p->size ; i++)
1401                             emul(e1, &(res->x.p->arr[i]));
1402                     else
1403                         emul_rev(e1, res);
1404                     return;
1405                 }
1406             }
1407             break;
1408        case relation:
1409             emul_rev(e1, res);
1410             break;
1411        default:
1412             assert(0);
1413        }
1414    }
1415    else {
1416        if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
1417            /* Product of two rational numbers */
1418
1419             Value g;
1420             value_init(g);
1421             value_multiply(res->d,e1->d,res->d);
1422             value_multiply(res->x.n,e1->x.n,res->x.n );
1423             Gcd(res->x.n, res->d,&g);
1424            if (value_notone_p(g)) {
1425                value_division(res->d,res->d,g);
1426                value_division(res->x.n,res->x.n,g);
1427            }
1428            value_clear(g);
1429            return ;
1430        }
1431        else {
1432              if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
1433               /* Product of an expression (polynomial or peririodic) and a rational number */
1434
1435                 emul_rev(e1, res);
1436                  return ;
1437              }
1438              else {
1439                /* Product of a rationel number and an expression (polynomial or peririodic) */
1440
1441                    i = res->x.p->type == fractional ? 1 : 0;
1442                    for (; i<res->x.p->size; i++)
1443                       emul(e1, &res->x.p->arr[i]);
1444
1445                  return ;
1446              }
1447        }
1448    }
1449
1450    return ;
1451 }
1452
1453 /* Frees mask ! */
1454 void emask(evalue *mask, evalue *res) {
1455     int n, i, j;
1456     Polyhedron *d, *fd;
1457     struct section *s;
1458     evalue mone;
1459
1460     if (EVALUE_IS_ZERO(*res))
1461         return;
1462
1463     assert(value_zero_p(mask->d));
1464     assert(mask->x.p->type == partition);
1465     assert(value_zero_p(res->d));
1466     assert(res->x.p->type == partition);
1467
1468     s = (struct section *)
1469             malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1470                    sizeof(struct section));
1471     assert(s);
1472
1473     value_init(mone.d);
1474     evalue_set_si(&mone, -1, 1);
1475
1476     n = 0;
1477     for (j = 0; j < res->x.p->size/2; ++j) {
1478         assert(mask->x.p->size >= 2);
1479         fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1480                               EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1481         if (!emptyQ(fd))
1482             for (i = 1; i < mask->x.p->size/2; ++i) {
1483                 Polyhedron *t = fd;
1484                 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1485                 Domain_Free(t);
1486                 if (emptyQ(fd))
1487                     break;
1488             }
1489         if (emptyQ(fd)) {
1490             Domain_Free(fd);
1491             continue;
1492         }
1493         value_init(s[n].E.d);
1494         evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1495         s[n].D = fd;
1496         ++n;
1497     }
1498     for (i = 0; i < mask->x.p->size/2; ++i) {
1499         if (EVALUE_IS_ONE(mask->x.p->arr[2*i+1]))
1500             continue;
1501
1502         fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1503         eadd(&mone, &mask->x.p->arr[2*i+1]);
1504         emul(&mone, &mask->x.p->arr[2*i+1]);
1505         for (j = 0; j < res->x.p->size/2; ++j) {
1506             Polyhedron *t;
1507             d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1508                                    EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1509             if (emptyQ(d)) {
1510                 Domain_Free(d);
1511                 continue;
1512             }
1513             t = fd;
1514             fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1515             if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1516                 Domain_Free(t);
1517             value_init(s[n].E.d);
1518             evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1519             emul(&mask->x.p->arr[2*i+1], &s[n].E);
1520             s[n].D = d;
1521             ++n;
1522         }
1523
1524         if (!emptyQ(fd)) {
1525             /* Just ignore; this may have been previously masked off */
1526         }
1527     }
1528
1529     free_evalue_refs(&mone);
1530     free_evalue_refs(mask);
1531     free_evalue_refs(res);
1532     value_init(res->d);
1533     if (n == 0)
1534         evalue_set_si(res, 0, 1);
1535     else {
1536         res->x.p = new_enode(partition, 2*n, -1);
1537         for (j = 0; j < n; ++j) {
1538             EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1539             value_clear(res->x.p->arr[2*j+1].d);
1540             res->x.p->arr[2*j+1] = s[j].E;
1541         }
1542     }
1543
1544     free(s);
1545 }
1546
1547 void evalue_copy(evalue *dst, evalue *src)
1548 {
1549     value_assign(dst->d, src->d);
1550     if(value_notzero_p(src->d)) {
1551          value_init(dst->x.n);
1552          value_assign(dst->x.n, src->x.n);
1553     } else
1554          dst->x.p = ecopy(src->x.p);
1555 }
1556
1557 enode *new_enode(enode_type type,int size,int pos) {
1558
1559   enode *res;
1560   int i;
1561
1562   if(size == 0) {
1563     fprintf(stderr, "Allocating enode of size 0 !\n" );
1564     return NULL;
1565   }
1566   res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1567   res->type = type;
1568   res->size = size;
1569   res->pos = pos;
1570   for(i=0; i<size; i++) {
1571     value_init(res->arr[i].d);
1572     value_set_si(res->arr[i].d,0);
1573     res->arr[i].x.p = 0;
1574   }
1575   return res;
1576 } /* new_enode */
1577
1578 enode *ecopy(enode *e) {
1579
1580   enode *res;
1581   int i;
1582
1583   res = new_enode(e->type,e->size,e->pos);
1584   for(i=0;i<e->size;++i) {
1585     value_assign(res->arr[i].d,e->arr[i].d);
1586     if(value_zero_p(res->arr[i].d))
1587       res->arr[i].x.p = ecopy(e->arr[i].x.p);
1588     else if (EVALUE_IS_DOMAIN(res->arr[i]))
1589       EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1590     else {
1591       value_init(res->arr[i].x.n);
1592       value_assign(res->arr[i].x.n,e->arr[i].x.n);
1593     }
1594   }
1595   return(res);
1596 } /* ecopy */
1597
1598 int ecmp(const evalue *e1, const evalue *e2)
1599 {
1600     enode *p1, *p2;
1601     int i;
1602     int r;
1603
1604     if (value_notzero_p(e1->d) && value_notzero_p(e2->d)) {
1605         Value m, m2;
1606         value_init(m);
1607         value_init(m2);
1608         value_multiply(m, e1->x.n, e2->d);
1609         value_multiply(m2, e2->x.n, e1->d);
1610
1611         if (value_lt(m, m2))
1612             r = -1;
1613         else if (value_gt(m, m2))
1614             r = 1;
1615         else
1616             r = 0;
1617
1618         value_clear(m);
1619         value_clear(m2);
1620
1621         return r;
1622     }
1623     if (value_notzero_p(e1->d))
1624         return -1;
1625     if (value_notzero_p(e2->d))
1626         return 1;
1627
1628     p1 = e1->x.p;
1629     p2 = e2->x.p;
1630
1631     if (p1->type != p2->type)
1632         return p1->type - p2->type;
1633     if (p1->pos != p2->pos)
1634         return p1->pos - p2->pos;
1635     if (p1->size != p2->size)
1636         return p1->size - p2->size;
1637
1638     for (i = p1->size-1; i >= 0; --i)
1639         if ((r = ecmp(&p1->arr[i], &p2->arr[i])) != 0)
1640             return r;
1641
1642     return 0;
1643 }
1644
1645 int eequal(evalue *e1,evalue *e2) {
1646
1647     int i;
1648     enode *p1, *p2;
1649
1650     if (value_ne(e1->d,e2->d))
1651         return 0;
1652
1653     /* e1->d == e2->d */
1654     if (value_notzero_p(e1->d)) {
1655         if (value_ne(e1->x.n,e2->x.n))
1656             return 0;
1657
1658         /* e1->d == e2->d != 0  AND e1->n == e2->n */
1659         return 1;
1660     }
1661
1662     /* e1->d == e2->d == 0 */
1663     p1 = e1->x.p;
1664     p2 = e2->x.p;
1665     if (p1->type != p2->type) return 0;
1666     if (p1->size != p2->size) return 0;
1667     if (p1->pos  != p2->pos) return 0;
1668     for (i=0; i<p1->size; i++)
1669         if (!eequal(&p1->arr[i], &p2->arr[i]) )
1670             return 0;
1671     return 1;
1672 } /* eequal */
1673
1674 void free_evalue_refs(evalue *e) {
1675
1676   enode *p;
1677   int i;
1678
1679   if (EVALUE_IS_DOMAIN(*e)) {
1680     Domain_Free(EVALUE_DOMAIN(*e));
1681     value_clear(e->d);
1682     return;
1683   } else if (value_pos_p(e->d)) {
1684
1685     /* 'e' stores a constant */
1686     value_clear(e->d);
1687     value_clear(e->x.n);
1688     return;
1689   }
1690   assert(value_zero_p(e->d));
1691   value_clear(e->d);
1692   p = e->x.p;
1693   if (!p) return;       /* null pointer */
1694   for (i=0; i<p->size; i++) {
1695     free_evalue_refs(&(p->arr[i]));
1696   }
1697   free(p);
1698   return;
1699 } /* free_evalue_refs */
1700
1701 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
1702                         Vector * val, evalue *res)
1703 {
1704     unsigned nparam = periods->Size;
1705
1706     if (p == nparam) {
1707         double d = compute_evalue(e, val->p);
1708         d *= VALUE_TO_DOUBLE(m);
1709         if (d > 0)
1710             d += .25;
1711         else
1712             d -= .25;
1713         value_assign(res->d, m);
1714         value_init(res->x.n);
1715         value_set_double(res->x.n, d);
1716         mpz_fdiv_r(res->x.n, res->x.n, m);
1717         return;
1718     }
1719     if (value_one_p(periods->p[p]))
1720         mod2table_r(e, periods, m, p+1, val, res);
1721     else {
1722         Value tmp;
1723         value_init(tmp);
1724
1725         value_assign(tmp, periods->p[p]);
1726         value_set_si(res->d, 0);
1727         res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
1728         do {
1729             value_decrement(tmp, tmp);
1730             value_assign(val->p[p], tmp);
1731             mod2table_r(e, periods, m, p+1, val,
1732                         &res->x.p->arr[VALUE_TO_INT(tmp)]);
1733         } while (value_pos_p(tmp));
1734
1735         value_clear(tmp);
1736     }
1737 }
1738
1739 static void rel2table(evalue *e, int zero)
1740 {
1741     if (value_pos_p(e->d)) {
1742         if (value_zero_p(e->x.n) == zero)
1743             value_set_si(e->x.n, 1);
1744         else
1745             value_set_si(e->x.n, 0);
1746         value_set_si(e->d, 1);
1747     } else {
1748         int i;
1749         for (i = 0; i < e->x.p->size; ++i)
1750             rel2table(&e->x.p->arr[i], zero);
1751     }
1752 }
1753
1754 void evalue_mod2table(evalue *e, int nparam)
1755 {
1756   enode *p;
1757   int i;
1758
1759   if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
1760     return;
1761   p = e->x.p;
1762   for (i=0; i<p->size; i++) {
1763     evalue_mod2table(&(p->arr[i]), nparam);
1764   }
1765   if (p->type == relation) {
1766     evalue copy;
1767     evalue *ev;
1768
1769     if (p->size > 2) {
1770       value_init(copy.d);
1771       evalue_copy(&copy, &p->arr[0]);
1772     }
1773     rel2table(&p->arr[0], 1);
1774     emul(&p->arr[0], &p->arr[1]);
1775     if (p->size > 2) {
1776       rel2table(&copy, 0);
1777       emul(&copy, &p->arr[2]);
1778       eadd(&p->arr[2], &p->arr[1]);
1779       free_evalue_refs(&p->arr[2]);
1780       free_evalue_refs(&copy);
1781     }
1782     free_evalue_refs(&p->arr[0]);
1783     ev = &p->arr[1];
1784     free(p);
1785     value_clear(e->d);
1786     *e = *ev;
1787   } else if (p->type == fractional) {
1788     Vector *periods = Vector_Alloc(nparam);
1789     Vector *val = Vector_Alloc(nparam);
1790     Value tmp;
1791     evalue *ev;
1792     evalue EP, res;
1793
1794     value_init(tmp);
1795     value_set_si(tmp, 1);
1796     Vector_Set(periods->p, 1, nparam);
1797     Vector_Set(val->p, 0, nparam);
1798     for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
1799       enode *p = ev->x.p;
1800
1801       assert(p->type == polynomial);
1802       assert(p->size == 2);
1803       value_assign(periods->p[p->pos-1], p->arr[1].d);
1804       Lcm3(tmp, p->arr[1].d, &tmp);
1805     }
1806     value_init(EP.d);
1807     mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
1808
1809     value_init(res.d);
1810     evalue_set_si(&res, 0, 1);
1811     /* Compute the polynomial using Horner's rule */
1812     for (i=p->size-1;i>1;i--) {
1813       eadd(&p->arr[i], &res);
1814       emul(&EP, &res);
1815     }
1816     eadd(&p->arr[1], &res);
1817
1818     free_evalue_refs(e);
1819     free_evalue_refs(&EP);
1820     *e = res;
1821
1822     value_clear(tmp);
1823     Vector_Free(val);
1824     Vector_Free(periods);
1825   }
1826 } /* evalue_mod2table */
1827
1828 /********************************************************/
1829 /* function in domain                                   */
1830 /*    check if the parameters in list_args              */
1831 /*    verifies the constraints of Domain P              */
1832 /********************************************************/
1833 int in_domain(Polyhedron *P, Value *list_args) {
1834
1835   int col,row;
1836   Value v; /* value of the constraint of a row when
1837                parameters are instanciated*/
1838   Value tmp;
1839
1840   value_init(v);
1841   value_init(tmp);
1842
1843   /*P->Constraint constraint matrice of polyhedron P*/
1844   for(row=0;row<P->NbConstraints;row++) {
1845     value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
1846     for(col=1;col<P->Dimension+1;col++) {
1847       value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
1848       value_addto(v,v,tmp);
1849     }
1850     if (value_notzero_p(P->Constraint[row][0])) {
1851
1852       /*if v is not >=0 then this constraint is not respected */
1853       if (value_neg_p(v)) {
1854 next:
1855         value_clear(v);
1856         value_clear(tmp);
1857         return P->next ? in_domain(P->next, list_args) : 0;
1858       }
1859     }
1860     else {
1861
1862       /*if v is not = 0 then this constraint is not respected */
1863       if (value_notzero_p(v))
1864         goto next;
1865     }
1866   }
1867
1868   /*if not return before this point => all
1869     the constraints are respected */
1870   value_clear(v);
1871   value_clear(tmp);
1872   return 1;
1873 } /* in_domain */
1874
1875 /****************************************************/
1876 /* function compute enode                           */
1877 /*     compute the value of enode p with parameters */
1878 /*     list "list_args                              */
1879 /*     compute the polynomial or the periodic       */
1880 /****************************************************/
1881
1882 static double compute_enode(enode *p, Value *list_args) {
1883
1884   int i;
1885   Value m, param;
1886   double res=0.0;
1887
1888   if (!p)
1889     return(0.);
1890
1891   value_init(m);
1892   value_init(param);
1893
1894   if (p->type == polynomial) {
1895     if (p->size > 1)
1896         value_assign(param,list_args[p->pos-1]);
1897
1898     /* Compute the polynomial using Horner's rule */
1899     for (i=p->size-1;i>0;i--) {
1900       res +=compute_evalue(&p->arr[i],list_args);
1901       res *=VALUE_TO_DOUBLE(param);
1902     }
1903     res +=compute_evalue(&p->arr[0],list_args);
1904   }
1905   else if (p->type == fractional) {
1906     double d = compute_evalue(&p->arr[0], list_args);
1907     d -= floor(d+1e-10);
1908
1909     /* Compute the polynomial using Horner's rule */
1910     for (i=p->size-1;i>1;i--) {
1911       res +=compute_evalue(&p->arr[i],list_args);
1912       res *=d;
1913     }
1914     res +=compute_evalue(&p->arr[1],list_args);
1915   }
1916   else if (p->type == periodic) {
1917     value_assign(param,list_args[p->pos-1]);
1918
1919     /* Choose the right element of the periodic */
1920     value_absolute(m,param);
1921     value_set_si(param,p->size);
1922     value_modulus(m,m,param);
1923     res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
1924   }
1925   else if (p->type == relation) {
1926     if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
1927       res = compute_evalue(&p->arr[1], list_args);
1928     else if (p->size > 2)
1929       res = compute_evalue(&p->arr[2], list_args);
1930   }
1931   else if (p->type == partition) {
1932     for (i = 0; i < p->size/2; ++i)
1933       if (in_domain(EVALUE_DOMAIN(p->arr[2*i]), list_args)) {
1934         res = compute_evalue(&p->arr[2*i+1], list_args);
1935         break;
1936       }
1937   }
1938   value_clear(m);
1939   value_clear(param);
1940   return res;
1941 } /* compute_enode */
1942
1943 /*************************************************/
1944 /* return the value of Ehrhart Polynomial        */
1945 /* It returns a double, because since it is      */
1946 /* a recursive function, some intermediate value */
1947 /* might not be integral                         */
1948 /*************************************************/
1949
1950 double compute_evalue(evalue *e,Value *list_args) {
1951
1952   double res;
1953
1954   if (value_notzero_p(e->d)) {
1955     if (value_notone_p(e->d))
1956       res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
1957     else
1958       res = VALUE_TO_DOUBLE(e->x.n);
1959   }
1960   else
1961     res = compute_enode(e->x.p,list_args);
1962   return res;
1963 } /* compute_evalue */
1964
1965
1966 /****************************************************/
1967 /* function compute_poly :                          */
1968 /* Check for the good validity domain               */
1969 /* return the number of point in the Polyhedron     */
1970 /* in allocated memory                              */
1971 /* Using the Ehrhart pseudo-polynomial              */
1972 /****************************************************/
1973 Value *compute_poly(Enumeration *en,Value *list_args) {
1974
1975   Value *tmp;
1976   /*    double d; int i; */
1977
1978   tmp = (Value *) malloc (sizeof(Value));
1979   assert(tmp != NULL);
1980   value_init(*tmp);
1981   value_set_si(*tmp,0);
1982
1983   if(!en)
1984     return(tmp);        /* no ehrhart polynomial */
1985   if(en->ValidityDomain) {
1986     if(!en->ValidityDomain->Dimension) { /* no parameters */
1987       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1988       return(tmp);
1989     }
1990   }
1991   else
1992     return(tmp);  /* no Validity Domain */
1993   while(en) {
1994     if(in_domain(en->ValidityDomain,list_args)) {
1995
1996 #ifdef EVAL_EHRHART_DEBUG
1997       Print_Domain(stdout,en->ValidityDomain);
1998       print_evalue(stdout,&en->EP);
1999 #endif
2000
2001       /*                        d = compute_evalue(&en->EP,list_args);
2002                                 i = d;
2003                                 printf("(double)%lf = %d\n", d, i ); */
2004       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2005       return(tmp);
2006     }
2007     else
2008       en=en->next;
2009   }
2010   value_set_si(*tmp,0);
2011   return(tmp); /* no compatible domain with the arguments */
2012 } /* compute_poly */
2013
2014 size_t value_size(Value v) {
2015     return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2016             * sizeof(v[0]._mp_d[0]);
2017 }
2018
2019 size_t domain_size(Polyhedron *D)
2020 {
2021     int i, j;
2022     size_t s = sizeof(*D);
2023
2024     for (i = 0; i < D->NbConstraints; ++i)
2025         for (j = 0; j < D->Dimension+2; ++j)
2026             s += value_size(D->Constraint[i][j]);
2027
2028     for (i = 0; i < D->NbRays; ++i)
2029         for (j = 0; j < D->Dimension+2; ++j)
2030             s += value_size(D->Ray[i][j]);
2031
2032     return D->next ? s+domain_size(D->next) : s;
2033 }
2034
2035 size_t enode_size(enode *p) {
2036     size_t s = sizeof(*p) - sizeof(p->arr[0]);
2037     int i;
2038
2039     if (p->type == partition)
2040         for (i = 0; i < p->size/2; ++i) {
2041             s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2042             s += evalue_size(&p->arr[2*i+1]);
2043         }
2044     else
2045         for (i = 0; i < p->size; ++i) {
2046             s += evalue_size(&p->arr[i]);
2047         }
2048     return s;
2049 }
2050
2051 size_t evalue_size(evalue *e)
2052 {
2053     size_t s = sizeof(*e);
2054     s += value_size(e->d);
2055     if (value_notzero_p(e->d))
2056         s += value_size(e->x.n);
2057     else
2058         s += enode_size(e->x.p);
2059     return s;
2060 }
2061
2062 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2063 {
2064     evalue *found = NULL;
2065     evalue offset;
2066     evalue copy;
2067     int i;
2068
2069     if (value_pos_p(e->d) || e->x.p->type != fractional)
2070         return NULL;
2071
2072     value_init(offset.d);
2073     value_init(offset.x.n);
2074     poly_denom(&e->x.p->arr[0], &offset.d);
2075     Lcm3(m, offset.d, &offset.d);
2076     value_set_si(offset.x.n, 1);
2077
2078     value_init(copy.d);
2079     evalue_copy(&copy, cst);
2080
2081     eadd(&offset, cst);
2082     mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2083
2084     if (eequal(base, &e->x.p->arr[0]))
2085         found = &e->x.p->arr[0];
2086     else {
2087         value_set_si(offset.x.n, -2);
2088
2089         eadd(&offset, cst);
2090         mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2091
2092         if (eequal(base, &e->x.p->arr[0]))
2093             found = base;
2094     }
2095     free_evalue_refs(cst);
2096     free_evalue_refs(&offset);
2097     *cst = copy;
2098
2099     for (i = 1; !found && i < e->x.p->size; ++i)
2100         found = find_second(base, cst, &e->x.p->arr[i], m);
2101
2102     return found;
2103 }
2104
2105 static evalue *find_relation_pair(evalue *e)
2106 {
2107     int i;
2108     evalue *found = NULL;
2109
2110     if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2111         return NULL;
2112
2113     if (e->x.p->type == fractional) {
2114         Value m;
2115         evalue *cst;
2116
2117         value_init(m);
2118         poly_denom(&e->x.p->arr[0], &m);
2119
2120         for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2121                                     cst = &cst->x.p->arr[0])
2122             ;
2123
2124         for (i = 1; !found && i < e->x.p->size; ++i)
2125             found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2126
2127         value_clear(m);
2128     }
2129
2130     i = e->x.p->type == relation;
2131     for (; !found && i < e->x.p->size; ++i)
2132         found = find_relation_pair(&e->x.p->arr[i]);
2133
2134     return found;
2135 }
2136
2137 void evalue_mod2relation(evalue *e) {
2138     evalue *d;
2139
2140     if (value_zero_p(e->d) && e->x.p->type == partition) {
2141         int i;
2142
2143         for (i = 0; i < e->x.p->size/2; ++i)
2144             evalue_mod2relation(&e->x.p->arr[2*i+1]);
2145
2146         return;
2147     }
2148
2149     while ((d = find_relation_pair(e)) != NULL) {
2150         evalue split;
2151         evalue *ev;
2152
2153         value_init(split.d);
2154         value_set_si(split.d, 0);
2155         split.x.p = new_enode(relation, 3, 0);
2156         evalue_set_si(&split.x.p->arr[1], 1, 1);
2157         evalue_set_si(&split.x.p->arr[2], 1, 1);
2158
2159         ev = &split.x.p->arr[0];
2160         value_set_si(ev->d, 0);
2161         ev->x.p = new_enode(fractional, 3, -1);
2162         evalue_set_si(&ev->x.p->arr[1], 0, 1);
2163         evalue_set_si(&ev->x.p->arr[2], 1, 1);
2164         evalue_copy(&ev->x.p->arr[0], d);
2165
2166         emul(&split, e);
2167
2168         reduce_evalue(e);
2169
2170         free_evalue_refs(&split);
2171     }
2172 }
2173
2174 static int evalue_comp(const void * a, const void * b)
2175 {
2176     const evalue *e1 = *(const evalue **)a;
2177     const evalue *e2 = *(const evalue **)b;
2178     return ecmp(e1, e2);
2179 }
2180
2181 void evalue_combine(evalue *e)
2182 {
2183     evalue **evs;
2184     int i, k;
2185     enode *p;
2186     evalue tmp;
2187
2188     if (value_notzero_p(e->d) || e->x.p->type != partition)
2189         return;
2190
2191     NALLOC(evs, e->x.p->size/2);
2192     for (i = 0; i < e->x.p->size/2; ++i)
2193         evs[i] = &e->x.p->arr[2*i+1];
2194     qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2195     p = new_enode(partition, e->x.p->size, -1);
2196     for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2197         if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2198             p->arr[2*k] = *(evs[i]-1);
2199             p->arr[2*k+1] = *(evs[i]);
2200             ++k;
2201         } else {
2202             Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2203             Polyhedron *L = D;
2204
2205             while (L->next)
2206                 L = L->next;
2207             L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2208             EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2209             free_evalue_refs(evs[i]);
2210         }
2211     }
2212     free(e->x.p);
2213     p->size = 2*k;
2214     e->x.p = p;
2215
2216     for (i = 0; i < e->x.p->size/2; ++i) {
2217         Polyhedron *H;
2218         if (value_notzero_p(e->x.p->arr[2*i+1].d))
2219             continue;
2220         H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2221         if (H == NULL)
2222             continue;
2223         for (k = 0; k < e->x.p->size/2; ++k) {
2224             Polyhedron *D, *N, **P;
2225             if (i == k)
2226                 continue;
2227             P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2228             D = *P;
2229             if (D == NULL)
2230                 continue;
2231             for (; D; D = N) {
2232                 *P = D;
2233                 N = D->next;
2234                 if (D->NbEq <= H->NbEq) {
2235                     P = &D->next;
2236                     continue;
2237                 }
2238
2239                 value_init(tmp.d);
2240                 tmp.x.p = new_enode(partition, 2, -1);
2241                 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2242                 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2243                 reduce_evalue(&tmp);
2244                 if (value_notzero_p(tmp.d) ||
2245                         ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2246                     P = &D->next;
2247                 else {
2248                     D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2249                     EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2250                     *P = NULL;
2251                 }
2252                 free_evalue_refs(&tmp);
2253             }
2254         }
2255     }
2256
2257     for (i = 0; i < e->x.p->size/2; ++i) {
2258         Polyhedron *H, *E;
2259         Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2260         if (!D) {
2261             free_evalue_refs(&e->x.p->arr[2*i+1]);
2262             e->x.p->size -= 2;
2263             if (2*i < e->x.p->size) {
2264                 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2265                 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2266             }
2267             --i;
2268             continue;
2269         }
2270         if (!D->next)
2271             continue;
2272         H = DomainConvex(D, 0);
2273         E = DomainDifference(H, D, 0);
2274         Domain_Free(D);
2275         D = DomainDifference(H, E, 0);
2276         Domain_Free(H);
2277         Domain_Free(E);
2278         EVALUE_SET_DOMAIN(p->arr[2*i], D);
2279     }
2280 }