ev_operations.c

   1 #include <assert.h>
   2 #include <math.h>
   3
   4 #include "ev_operations.h"
   5
   6 void evalue_set_si(evalue *ev, int n, int d) {
   7     value_set_si(ev->d, d);
   8     value_init(ev->x.n);
   9     value_set_si(ev->x.n, n);
  10 }
  11
  12 void evalue_set(evalue *ev, Value n, Value d) {
  13     value_assign(ev->d, d);
  14     value_init(ev->x.n);
  15     value_assign(ev->x.n, n);
  16 }
  17
  18 void aep_evalue(evalue *e, int *ref) {
  19
  20     enode *p;
  21     int i;
  22
  23     if (value_notzero_p(e->d))
  24         return;         /* a rational number, its already reduced */
  25     if(!(p = e->x.p))
  26         return;         /* hum... an overflow probably occured */
  27
  28     /* First check the components of p */
  29     for (i=0;i<p->size;i++)
  30         aep_evalue(&p->arr[i],ref);
  31
  32     /* Then p itself */
  33     if (p->type != modulo)
  34         p->pos = ref[p->pos-1]+1;
  35     return;
  36 } /* aep_evalue */
  37
  38 /** Comments */
  39 void addeliminatedparams_evalue(evalue *e,Matrix *CT) {
  40
  41     enode *p;
  42     int i, j;
  43     int *ref;
  44
  45     if (value_notzero_p(e->d))
  46         return;          /* a rational number, its already reduced */
  47     if(!(p = e->x.p))
  48         return;          /* hum... an overflow probably occured */
  49
  50     /* Compute ref */
  51     ref = (int *)malloc(sizeof(int)*(CT->NbRows-1));
  52     for(i=0;i<CT->NbRows-1;i++)
  53         for(j=0;j<CT->NbColumns;j++)
  54             if(value_notzero_p(CT->p[i][j])) {
  55                 ref[i] = j;
  56                 break;
  57             }
  58
  59     /* Transform the references in e, using ref */
  60     aep_evalue(e,ref);
  61     free( ref );
  62     return;
  63 } /* addeliminatedparams_evalue */
  64
  65 void reduce_evalue (evalue *e) {
  66
  67     enode *p;
  68     int i, j, k;
  69
  70     if (value_notzero_p(e->d))
  71         return; /* a rational number, its already reduced */
  72     if(!(p = e->x.p))
  73         return; /* hum... an overflow probably occured */
  74
  75     /* First reduce the components of p */
  76     for (i=0; i<p->size; i++)
  77         reduce_evalue(&p->arr[i]);
  78
  79     if (p->type==periodic) {
  80
  81         /* Try to reduce the period */
  82         for (i=1; i<=(p->size)/2; i++) {
  83             if ((p->size % i)==0) {
  84
  85                 /* Can we reduce the size to i ? */
  86                 for (j=0; j<i; j++)
  87                     for (k=j+i; k<e->x.p->size; k+=i)
  88                         if (!eequal(&p->arr[j], &p->arr[k])) goto you_lose;
  89
  90                 /* OK, lets do it */
  91                 for (j=i; j<p->size; j++) free_evalue_refs(&p->arr[j]);
  92                 p->size = i;
  93                 break;
  94
  95 you_lose:       /* OK, lets not do it */
  96                 continue;
  97             }
  98         }
  99
 100         /* Try to reduce its strength */
 101         if (p->size == 1) {
 102             value_clear(e->d);
 103             memcpy(e,&p->arr[0],sizeof(evalue));
 104             free(p);
 105         }
 106     }
 107     else if (p->type==polynomial) {
 108
 109         /* Try to reduce the degree */
 110         for (i=p->size-1;i>=1;i--) {
 111             if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
 112                 break;
 113             /* Zero coefficient */
 114             free_evalue_refs(&(p->arr[i]));
 115         }
 116         if (i+1<p->size)
 117             p->size = i+1;
 118
 119         /* Try to reduce its strength */
 120         if (p->size == 1) {
 121             value_clear(e->d);
 122             memcpy(e,&p->arr[0],sizeof(evalue));
 123             free(p);
 124         }
 125     }
 126     else if (p->type==modulo) {
 127
 128         /* Try to reduce the degree */
 129         for (i=p->size-1;i>=2;i--) {
 130             if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
 131                 break;
 132             /* Zero coefficient */
 133             free_evalue_refs(&(p->arr[i]));
 134         }
 135         if (i+1<p->size)
 136             p->size = i+1;
 137
 138         /* Try to reduce its strength */
 139         if (p->size == 2) {
 140             value_clear(e->d);
 141             memcpy(e,&p->arr[1],sizeof(evalue));
 142             free_evalue_refs(&(p->arr[0]));
 143             free(p);
 144         }
 145     }
 146 } /* reduce_evalue */
 147
 148 void print_evalue(FILE *DST,evalue *e,char **pname) {
 149
 150   if(value_notzero_p(e->d)) {
 151     if(value_notone_p(e->d)) {
 152       value_print(DST,VALUE_FMT,e->x.n);
 153       fprintf(DST,"/");
 154       value_print(DST,VALUE_FMT,e->d);
 155     }
 156     else {
 157       value_print(DST,VALUE_FMT,e->x.n);
 158     }
 159   }
 160   else
 161     print_enode(DST,e->x.p,pname);
 162   return;
 163 } /* print_evalue */
 164
 165 void print_enode(FILE *DST,enode *p,char **pname) {
 166
 167   int i;
 168
 169   if (!p) {
 170     fprintf(DST, "NULL");
 171     return;
 172   }
 173   switch (p->type) {
 174   case evector:
 175     fprintf(DST, "{ ");
 176     for (i=0; i<p->size; i++) {
 177       print_evalue(DST, &p->arr[i], pname);
 178       if (i!=(p->size-1))
 179         fprintf(DST, ", ");
 180     }
 181     fprintf(DST, " }\n");
 182     break;
 183   case polynomial:
 184     fprintf(DST, "( ");
 185     for (i=p->size-1; i>=0; i--) {
 186       print_evalue(DST, &p->arr[i], pname);
 187       if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
 188       else if (i>1)
 189         fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
 190     }
 191     fprintf(DST, " )\n");
 192     break;
 193   case periodic:
 194     fprintf(DST, "[ ");
 195     for (i=0; i<p->size; i++) {
 196       print_evalue(DST, &p->arr[i], pname);
 197       if (i!=(p->size-1)) fprintf(DST, ", ");
 198     }
 199     fprintf(DST," ]_%s", pname[p->pos-1]);
 200     break;
 201   case modulo:
 202     fprintf(DST, "( ");
 203     for (i=p->size-1; i>=1; i--) {
 204       print_evalue(DST, &p->arr[i], pname);
 205       if (i >= 2) {
 206         fprintf(DST, " * ");
 207         if (i > 2)
 208           fprintf(DST, "(");
 209         fprintf(DST, "(");
 210         print_evalue(DST, &p->arr[0], pname);
 211         fprintf(DST, ") mod %d", p->pos);
 212         if (i>2)
 213           fprintf(DST, ")^%d + ", i-1);
 214         else
 215           fprintf(DST, " + ", i-1);
 216       }
 217     }
 218     fprintf(DST, " )\n");
 219     break;
 220   case relation:
 221     fprintf(DST, "[ ");
 222     print_evalue(DST, &p->arr[0], pname);
 223     fprintf(DST, "= 0 ] * \n");
 224     print_evalue(DST, &p->arr[1], pname);
 225     if (p->size > 2) {
 226         fprintf(DST, " +\n [ ");
 227         print_evalue(DST, &p->arr[0], pname);
 228         fprintf(DST, "!= 0 ] * \n");
 229         print_evalue(DST, &p->arr[2], pname);
 230     }
 231     break;
 232   case partition:
 233     for (i=0; i<p->size/2; i++) {
 234         Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), pname);
 235         print_evalue(DST, &p->arr[2*i+1], pname);
 236     }
 237     break;
 238   default:
 239     assert(0);
 240   }
 241   return;
 242 } /* print_enode */
 243
 244 static int mod_term_smaller(evalue *e1, evalue *e2)
 245 {
 246     if (value_notzero_p(e1->d)) {
 247         if (value_zero_p(e2->d))
 248             return 1;
 249         return value_lt(e1->x.n, e2->x.n);
 250     }
 251     if (value_notzero_p(e2->d))
 252         return 0;
 253     if (e1->x.p->pos < e2->x.p->pos)
 254         return 1;
 255     else if (e1->x.p->pos > e2->x.p->pos)
 256         return 0;
 257     else
 258         return mod_term_smaller(&e1->x.p->arr[0], &e2->x.p->arr[0]);
 259 }
 260
 261 static void eadd_rev(evalue *e1, evalue *res)
 262 {
 263     evalue ev;
 264     value_init(ev.d);
 265     evalue_copy(&ev, e1);
 266     eadd(res, &ev);
 267     free_evalue_refs(res);
 268     *res = ev;
 269 }
 270
 271 static void eadd_rev_cst (evalue *e1, evalue *res)
 272 {
 273     evalue ev;
 274     value_init(ev.d);
 275     evalue_copy(&ev, e1);
 276     eadd(res, &ev.x.p->arr[ev.x.p->type==modulo]);
 277     free_evalue_refs(res);
 278     *res = ev;
 279 }
 280
 281 struct section { Polyhedron * D; evalue E; };
 282
 283 void eadd_partitions (evalue *e1,evalue *res)
 284 {
 285     int n, i, j;
 286     Polyhedron *d, *fd;
 287     struct section *s;
 288     s = (struct section *)
 289             malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
 290                    sizeof(struct section));
 291     assert(s);
 292
 293     n = 0;
 294     for (j = 0; j < e1->x.p->size/2; ++j) {
 295         assert(res->x.p->size >= 2);
 296         fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
 297                               EVALUE_DOMAIN(res->x.p->arr[0]), 0);
 298         if (!emptyQ(fd))
 299             for (i = 1; i < res->x.p->size/2; ++i) {
 300                 Polyhedron *t = fd;
 301                 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
 302                 Domain_Free(t);
 303                 if (emptyQ(fd))
 304                     break;
 305             }
 306         if (emptyQ(fd)) {
 307             Domain_Free(fd);
 308             continue;
 309         }
 310         value_init(s[n].E.d);
 311         evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
 312         s[n].D = fd;
 313         ++n;
 314     }
 315     for (i = 0; i < res->x.p->size/2; ++i) {
 316         fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
 317         for (j = 0; j < e1->x.p->size/2; ++j) {
 318             Polyhedron *t;
 319             d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
 320                                    EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
 321             if (emptyQ(d)) {
 322                 Domain_Free(d);
 323                 continue;
 324             }
 325             t = fd;
 326             fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
 327             if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
 328                 Domain_Free(t);
 329             value_init(s[n].E.d);
 330             evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
 331             eadd(&e1->x.p->arr[2*j+1], &s[n].E);
 332             s[n].D = d;
 333             ++n;
 334         }
 335         if (!emptyQ(fd)) {
 336             s[n].E = res->x.p->arr[2*i+1];
 337             s[n].D = fd;
 338             ++n;
 339         } else
 340             free_evalue_refs(&res->x.p->arr[2*i+1]);
 341         if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
 342             Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
 343         value_clear(res->x.p->arr[2*i].d);
 344     }
 345
 346     free(res->x.p);
 347     res->x.p = new_enode(partition, 2*n, -1);
 348     for (j = 0; j < n; ++j) {
 349         EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
 350         value_clear(res->x.p->arr[2*j+1].d);
 351         res->x.p->arr[2*j+1] = s[j].E;
 352     }
 353
 354     free(s);
 355 }
 356
 357 static void explicit_complement(evalue *res)
 358 {
 359     enode *rel = new_enode(relation, 3, 0);
 360     assert(rel);
 361     value_clear(rel->arr[0].d);
 362     rel->arr[0] = res->x.p->arr[0];
 363     value_clear(rel->arr[1].d);
 364     rel->arr[1] = res->x.p->arr[1];
 365     value_set_si(rel->arr[2].d, 1);
 366     value_init(rel->arr[2].x.n);
 367     value_set_si(rel->arr[2].x.n, 0);
 368     free(res->x.p);
 369     res->x.p = rel;
 370 }
 371
 372 void eadd(evalue *e1,evalue *res) {
 373
 374  int i;
 375     if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
 376          /* Add two rational numbers */
 377          Value g,m1,m2;
 378          value_init(g);
 379          value_init(m1);
 380          value_init(m2);
 381
 382          value_multiply(m1,e1->x.n,res->d);
 383          value_multiply(m2,res->x.n,e1->d);
 384          value_addto(res->x.n,m1,m2);
 385          value_multiply(res->d,e1->d,res->d);
 386          Gcd(res->x.n,res->d,&g);
 387          if (value_notone_p(g)) {
 388               value_division(res->d,res->d,g);
 389               value_division(res->x.n,res->x.n,g);
 390          }
 391          value_clear(g); value_clear(m1); value_clear(m2);
 392          return ;
 393      }
 394      else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
 395           switch (res->x.p->type) {
 396           case polynomial:
 397               /* Add the constant to the constant term of a polynomial*/
 398                eadd(e1, &res->x.p->arr[0]);
 399                return ;
 400           case periodic:
 401               /* Add the constant to all elements of a periodic number */
 402               for (i=0; i<res->x.p->size; i++) {
 403                   eadd(e1, &res->x.p->arr[i]);
 404               }
 405               return ;
 406           case evector:
 407               fprintf(stderr, "eadd: cannot add const with vector\n");
 408               return;
 409           case modulo:
 410                eadd(e1, &res->x.p->arr[1]);
 411                return ;
 412           case partition:
 413                 assert(EVALUE_IS_ZERO(*e1));
 414                 break;                          /* Do nothing */
 415           case relation:
 416                 /* Create (zero) complement if needed */
 417                 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
 418                     explicit_complement(res);
 419                 for (i = 1; i < res->x.p->size; ++i)
 420                     eadd(e1, &res->x.p->arr[i]);
 421                 break;
 422           default:
 423                 assert(0);
 424           }
 425      }
 426      /* add polynomial or periodic to constant
 427       * you have to exchange e1 and res, before doing addition */
 428
 429      else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
 430           eadd_rev(e1, res);
 431           return;
 432      }
 433      else {   // ((e1->d==0) && (res->d==0))
 434         assert(!((e1->x.p->type == partition) ^
 435                  (res->x.p->type == partition)));
 436         if (e1->x.p->type == partition) {
 437             eadd_partitions(e1, res);
 438             return;
 439         }
 440         if (e1->x.p->type == relation &&
 441             (res->x.p->type != relation ||
 442              mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
 443                 eadd_rev(e1, res);
 444                 return;
 445         }
 446         if (res->x.p->type == relation) {
 447             if (e1->x.p->type == relation &&
 448                 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
 449                     if (res->x.p->size < 3 && e1->x.p->size == 3)
 450                         explicit_complement(res);
 451                     if (e1->x.p->size < 3 && res->x.p->size == 3)
 452                         explicit_complement(e1);
 453                     for (i = 1; i < res->x.p->size; ++i)
 454                         eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
 455                     return;
 456             }
 457             if (res->x.p->size < 3)
 458                 explicit_complement(res);
 459             for (i = 1; i < res->x.p->size; ++i)
 460                 eadd(e1, &res->x.p->arr[i]);
 461             return;
 462         }
 463                  if ((e1->x.p->type != res->x.p->type) ) {
 464                       /* adding to evalues of different type. two cases are possible
 465                        * res is periodic and e1 is polynomial, you have to exchange
 466                        * e1 and res then to add e1 to the constant term of res */
 467                      if (e1->x.p->type == polynomial) {
 468                           eadd_rev_cst(e1, res);
 469                      }
 470                      else if (res->x.p->type == polynomial) {
 471                           /* res is polynomial and e1 is periodic,
 472                             add e1 to the constant term of res */
 473
 474                           eadd(e1,&res->x.p->arr[0]);
 475                      } else
 476                         assert(0);
 477
 478                      return;
 479                  }
 480                  else if (e1->x.p->pos != res->x.p->pos ||
 481                             (res->x.p->type == modulo &&
 482                              !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
 483                  /* adding evalues of different position (i.e function of different unknowns
 484                   * to case are possible  */
 485
 486                         switch (res->x.p->type) {
 487                         case modulo:
 488                             if(mod_term_smaller(res, e1))
 489                                 eadd(e1,&res->x.p->arr[1]);
 490                             else
 491                                 eadd_rev_cst(e1, res);
 492                             return;
 493                         case polynomial: //  res and e1 are polynomials
 494                                //  add e1 to the constant term of res
 495
 496                             if(res->x.p->pos < e1->x.p->pos)
 497                                eadd(e1,&res->x.p->arr[0]);
 498                             else
 499                                 eadd_rev_cst(e1, res);
 500                               // value_clear(g); value_clear(m1); value_clear(m2);
 501                                return;
 502                         case periodic:  // res and e1 are pointers to periodic numbers
 503                                   //add e1 to all elements of res
 504
 505                             if(res->x.p->pos < e1->x.p->pos)
 506                                   for (i=0;i<res->x.p->size;i++) {
 507                                        eadd(e1,&res->x.p->arr[i]);
 508                                   }
 509                             else
 510                                 eadd_rev(e1, res);
 511                             return;
 512                         }
 513                  }
 514
 515
 516                  //same type , same pos  and same size
 517                  if (e1->x.p->size == res->x.p->size) {
 518                       // add any element in e1 to the corresponding element in res
 519                       if (res->x.p->type == modulo)
 520                         assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
 521                       i = res->x.p->type == modulo ? 1 : 0;
 522                       for (; i<res->x.p->size; i++) {
 523                             eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
 524                       }
 525                       return ;
 526                 }
 527
 528                 /* Sizes are different */
 529                 switch(res->x.p->type) {
 530                 case polynomial:
 531                 case modulo:
 532                     /* VIN100: if e1-size > res-size you have to copy e1 in a   */
 533                     /* new enode and add res to that new node. If you do not do */
 534                     /* that, you lose the the upper weight part of e1 !         */
 535
 536                      if(e1->x.p->size > res->x.p->size)
 537                           eadd_rev(e1, res);
 538                      else {
 539
 540                         if (res->x.p->type == modulo)
 541                             assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
 542                         i = res->x.p->type == modulo ? 1 : 0;
 543                         for (; i<e1->x.p->size ; i++) {
 544                              eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
 545                         }
 546
 547                         return ;
 548                      }
 549                      break;
 550
 551     /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
 552                 case periodic:
 553                 {
 554                       /* you have to create a new evalue 'ne' in whitch size equals to the lcm
 555                        of the sizes of e1 and res, then to copy res periodicaly in ne, after
 556                        to add periodicaly elements of e1 to elements of ne, and finaly to
 557                        return ne. */
 558                        int x,y,p;
 559                        Value ex, ey ,ep;
 560                        evalue *ne;
 561                        value_init(ex); value_init(ey);value_init(ep);
 562                        x=e1->x.p->size;
 563                        y= res->x.p->size;
 564                        value_set_si(ex,e1->x.p->size);
 565                        value_set_si(ey,res->x.p->size);
 566                        value_assign (ep,*Lcm(ex,ey));
 567                        p=(int)mpz_get_si(ep);
 568                        ne= (evalue *) malloc (sizeof(evalue));
 569                        value_init(ne->d);
 570                        value_set_si( ne->d,0);
 571
 572                        ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
 573                        for(i=0;i<p;i++)  {
 574                               value_assign(ne->x.p->arr[i].d, res->x.p->arr[i%y].d);
 575                               if (value_notzero_p(ne->x.p->arr[i].d))   {
 576                                 value_init(ne->x.p->arr[i].x.n);
 577                                 value_assign(ne->x.p->arr[i].x.n, res->x.p->arr[i%y].x.n);
 578                               }
 579                               else {
 580                                   ne->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%y].x.p);
 581                               }
 582                        }
 583                        for(i=0;i<p;i++)  {
 584                             eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
 585                        }
 586
 587                       value_assign(res->d, ne->d);
 588                       res->x.p=ne->x.p;
 589
 590                         return ;
 591                 }
 592                 case evector:
 593                      fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
 594                      return ;
 595                 }
 596      }
 597      return ;
 598  }/* eadd  */
 599
 600 static void emul_rev (evalue *e1, evalue *res)
 601 {
 602     evalue ev;
 603     value_init(ev.d);
 604     evalue_copy(&ev, e1);
 605     emul(res, &ev);
 606     free_evalue_refs(res);
 607     *res = ev;
 608 }
 609
 610 static void emul_poly (evalue *e1, evalue *res)
 611 {
 612     int i, j, o = res->x.p->type == modulo;
 613     evalue tmp;
 614     int size=(e1->x.p->size + res->x.p->size - o - 1);
 615     value_init(tmp.d);
 616     value_set_si(tmp.d,0);
 617     tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
 618     if (o)
 619         evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
 620     for (i=o; i < e1->x.p->size; i++) {
 621         evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
 622         emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
 623     }
 624     for (; i<size; i++)
 625         evalue_set_si(&tmp.x.p->arr[i], 0, 1);
 626     for (i=o+1; i<res->x.p->size; i++)
 627         for (j=o; j<e1->x.p->size; j++) {
 628             evalue ev;
 629             value_init(ev.d);
 630             evalue_copy(&ev, &e1->x.p->arr[j]);
 631             emul(&res->x.p->arr[i], &ev);
 632             eadd(&ev, &tmp.x.p->arr[i+j-o]);
 633             free_evalue_refs(&ev);
 634         }
 635     free_evalue_refs(res);
 636     *res = tmp;
 637 }
 638
 639 void emul_partitions (evalue *e1,evalue *res)
 640 {
 641     int n, i, j, k;
 642     Polyhedron *d;
 643     struct section *s;
 644     s = (struct section *)
 645             malloc((e1->x.p->size/2) * (res->x.p->size/2) *
 646                    sizeof(struct section));
 647     assert(s);
 648
 649     n = 0;
 650     for (i = 0; i < res->x.p->size/2; ++i) {
 651         for (j = 0; j < e1->x.p->size/2; ++j) {
 652             d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
 653                                    EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
 654             if (emptyQ(d)) {
 655                 Domain_Free(d);
 656                 continue;
 657             }
 658
 659             /* This code is only needed because the partitions
 660                are not true partitions.
 661              */
 662             for (k = 0; k < n; ++k) {
 663                 if (DomainIncludes(s[k].D, d))
 664                     break;
 665                 if (DomainIncludes(d, s[k].D)) {
 666                     Domain_Free(s[k].D);
 667                     free_evalue_refs(&s[k].E);
 668                     if (n > k)
 669                         s[k] = s[--n];
 670                     --k;
 671                 }
 672             }
 673             if (k < n) {
 674                 Domain_Free(d);
 675                 continue;
 676             }
 677
 678             value_init(s[n].E.d);
 679             evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
 680             emul(&e1->x.p->arr[2*j+1], &s[n].E);
 681             s[n].D = d;
 682             ++n;
 683         }
 684         Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
 685         value_clear(res->x.p->arr[2*i].d);
 686         free_evalue_refs(&res->x.p->arr[2*i+1]);
 687     }
 688
 689     free(res->x.p);
 690     res->x.p = new_enode(partition, 2*n, -1);
 691     for (j = 0; j < n; ++j) {
 692         EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
 693         value_clear(res->x.p->arr[2*j+1].d);
 694         res->x.p->arr[2*j+1] = s[j].E;
 695     }
 696
 697     free(s);
 698 }
 699
 700 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
 701  * do a copy of "res" befor calling this function if you nead it after. The vector type of
 702  * evalues is not treated here */
 703
 704 void emul (evalue *e1, evalue *res ){
 705     int i,j;
 706
 707 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
 708     fprintf(stderr, "emul: do not proced on evector type !\n");
 709     return;
 710 }
 711
 712     if (value_zero_p(e1->d) && e1->x.p->type == partition) {
 713         if (value_zero_p(res->d) && res->x.p->type == partition)
 714             emul_partitions(e1, res);
 715         else
 716             emul_rev(e1, res);
 717     } else if (value_zero_p(res->d) && res->x.p->type == partition) {
 718         for (i = 0; i < res->x.p->size/2; ++i)
 719             emul(e1, &res->x.p->arr[2*i+1]);
 720     } else
 721    if (value_zero_p(res->d) && res->x.p->type == relation) {
 722         for (i = 1; i < res->x.p->size; ++i)
 723             emul(e1, &res->x.p->arr[i]);
 724    } else
 725    if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
 726        switch(e1->x.p->type) {
 727        case polynomial:
 728            switch(res->x.p->type) {
 729            case polynomial:
 730                if(e1->x.p->pos == res->x.p->pos) {
 731                /* Product of two polynomials of the same variable */
 732                     emul_poly(e1, res);
 733                     return;
 734                }
 735                else {
 736                   /* Product of two polynomials of different variables */
 737
 738                 if(res->x.p->pos < e1->x.p->pos)
 739                      for( i=0; i<res->x.p->size ; i++)
 740                             emul(e1, &res->x.p->arr[i]);
 741                 else
 742                     emul_rev(e1, res);
 743
 744                  return ;
 745                }
 746            case periodic:
 747            case modulo:
 748                 /* Product of a polynomial and a periodic or modulo */
 749                 emul_rev(e1, res);
 750                 return;
 751            }
 752        case periodic:
 753            switch(res->x.p->type) {
 754            case periodic:
 755                 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
 756                  /* Product of two periodics of the same parameter and period */
 757
 758                      for(i=0; i<res->x.p->size;i++)
 759                            emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
 760
 761                      return;
 762                 }
 763                 else{
 764                   if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
 765                    /* Product of two periodics of the same parameter and different periods */
 766                     evalue *newp;
 767                     Value x,y,z;
 768                     int ix,iy,lcm;
 769                     value_init(x); value_init(y);value_init(z);
 770                     ix=e1->x.p->size;
 771                     iy=res->x.p->size;
 772                     value_set_si(x,e1->x.p->size);
 773                     value_set_si(y,res->x.p->size);
 774                     value_assign (z,*Lcm(x,y));
 775                     lcm=(int)mpz_get_si(z);
 776                     newp= (evalue *) malloc (sizeof(evalue));
 777                     value_init(newp->d);
 778                     value_set_si( newp->d,0);
 779                     newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
 780                         for(i=0;i<lcm;i++)  {
 781                            value_assign(newp->x.p->arr[i].d, res->x.p->arr[i%iy].d);
 782                            if (value_notzero_p(newp->x.p->arr[i].d))   {
 783                            value_assign(newp->x.p->arr[i].x.n, res->x.p->arr[i%iy].x.n);
 784                            }
 785                            else {
 786                            newp->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%iy].x.p);
 787                            }
 788
 789                         }
 790                         for(i=0;i<lcm;i++)
 791                             emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
 792
 793                         value_assign(res->d,newp->d);
 794                         res->x.p=newp->x.p;
 795
 796                           value_clear(x); value_clear(y);value_clear(z);
 797                           return ;
 798                   }
 799                   else {
 800                      /* Product of two periodics of different parameters */
 801
 802                         for(i=0; i<res->x.p->size; i++)
 803                             emul(e1, &(res->x.p->arr[i]));
 804
 805                         return;
 806                   }
 807                 }
 808            case polynomial:
 809                   /* Product of a periodic and a polynomial */
 810
 811                        for(i=0; i<res->x.p->size ; i++)
 812                             emul(e1, &(res->x.p->arr[i]));
 813
 814                        return;
 815
 816            }
 817        case modulo:
 818             switch(res->x.p->type) {
 819             case polynomial:
 820                 for(i=0; i<res->x.p->size ; i++)
 821                     emul(e1, &(res->x.p->arr[i]));
 822                 return;
 823             case periodic:
 824                 assert(0);
 825             case modulo:
 826                 if (e1->x.p->pos == res->x.p->pos &&
 827                             eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
 828                     if (e1->x.p->pos != 2)
 829                         emul_poly(e1, res);
 830                     else {
 831                         evalue tmp;
 832                         /* x mod 2 == (x mod 2)^2 */
 833                         /* a0 b0 + (a0 b1 + a1 b0 + a1 b1) (x mod 2) */
 834                         assert(e1->x.p->size == 3);
 835                         assert(res->x.p->size == 3);
 836                         value_init(tmp.d);
 837                         evalue_copy(&tmp, &res->x.p->arr[1]);
 838                         eadd(&res->x.p->arr[2], &tmp);
 839                         emul(&e1->x.p->arr[2], &tmp);
 840                         emul(&e1->x.p->arr[1], res);
 841                         eadd(&tmp, &res->x.p->arr[2]);
 842                         free_evalue_refs(&tmp);
 843                     }
 844                 } else {
 845                     if(mod_term_smaller(res, e1))
 846                         for(i=1; i<res->x.p->size ; i++)
 847                             emul(e1, &(res->x.p->arr[i]));
 848                     else
 849                         emul_rev(e1, res);
 850                     return;
 851                 }
 852             }
 853             break;
 854        case relation:
 855             emul_rev(e1, res);
 856             break;
 857        default:
 858             assert(0);
 859        }
 860    }
 861    else {
 862        if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
 863            /* Product of two rational numbers */
 864
 865             Value g;
 866             value_init(g);
 867             value_multiply(res->d,e1->d,res->d);
 868             value_multiply(res->x.n,e1->x.n,res->x.n );
 869             Gcd(res->x.n, res->d,&g);
 870            if (value_notone_p(g)) {
 871                value_division(res->d,res->d,g);
 872                value_division(res->x.n,res->x.n,g);
 873            }
 874            value_clear(g);
 875            return ;
 876        }
 877        else {
 878              if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
 879               /* Product of an expression (polynomial or peririodic) and a rational number */
 880
 881                 emul_rev(e1, res);
 882                  return ;
 883              }
 884              else {
 885                /* Product of a rationel number and an expression (polynomial or peririodic) */
 886
 887                    i = res->x.p->type == modulo ? 1 : 0;
 888                    for (; i<res->x.p->size; i++)
 889                       emul(e1, &res->x.p->arr[i]);
 890
 891                  return ;
 892              }
 893        }
 894    }
 895
 896    return ;
 897 }
 898
 899 /* Frees mask ! */
 900 void emask(evalue *mask, evalue *res) {
 901     int n, i, j;
 902     Polyhedron *d, *fd;
 903     struct section *s;
 904     evalue mone;
 905
 906     assert(mask->x.p->type == partition);
 907     assert(res->x.p->type == partition);
 908
 909     s = (struct section *)
 910             malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
 911                    sizeof(struct section));
 912     assert(s);
 913
 914     value_init(mone.d);
 915     evalue_set_si(&mone, -1, 1);
 916
 917     n = 0;
 918     for (j = 0; j < res->x.p->size/2; ++j) {
 919         assert(mask->x.p->size >= 2);
 920         fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
 921                               EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
 922         if (!emptyQ(fd))
 923             for (i = 1; i < mask->x.p->size/2; ++i) {
 924                 Polyhedron *t = fd;
 925                 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
 926                 Domain_Free(t);
 927                 if (emptyQ(fd))
 928                     break;
 929             }
 930         if (emptyQ(fd)) {
 931             Domain_Free(fd);
 932             continue;
 933         }
 934         value_init(s[n].E.d);
 935         evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
 936         s[n].D = fd;
 937         ++n;
 938     }
 939     for (i = 0; i < mask->x.p->size/2; ++i) {
 940         fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
 941         eadd(&mone, &mask->x.p->arr[2*i+1]);
 942         emul(&mone, &mask->x.p->arr[2*i+1]);
 943         for (j = 0; j < res->x.p->size/2; ++j) {
 944             Polyhedron *t;
 945             d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
 946                                    EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
 947             if (emptyQ(d)) {
 948                 Domain_Free(d);
 949                 continue;
 950             }
 951             t = fd;
 952             fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
 953             if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
 954                 Domain_Free(t);
 955             value_init(s[n].E.d);
 956             evalue_copy(&s[n].E, &mask->x.p->arr[2*i+1]);
 957             emul(&res->x.p->arr[2*j+1], &s[n].E);
 958             s[n].D = d;
 959             ++n;
 960         }
 961         if (!emptyQ(fd)) {
 962             assert(0);          // We don't allow this.
 963         }
 964     }
 965
 966     free_evalue_refs(&mone);
 967     free_evalue_refs(mask);
 968     free_evalue_refs(res);
 969     value_init(res->d);
 970     res->x.p = new_enode(partition, 2*n, -1);
 971     for (j = 0; j < n; ++j) {
 972         EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
 973         value_clear(res->x.p->arr[2*j+1].d);
 974         res->x.p->arr[2*j+1] = s[j].E;
 975     }
 976
 977     free(s);
 978 }
 979
 980 void evalue_copy(evalue *dst, evalue *src)
 981 {
 982     value_assign(dst->d, src->d);
 983     if(value_notzero_p(src->d)) {
 984          value_init(dst->x.n);
 985          value_assign(dst->x.n, src->x.n);
 986     } else
 987          dst->x.p = ecopy(src->x.p);
 988 }
 989
 990 enode *new_enode(enode_type type,int size,int pos) {
 991
 992   enode *res;
 993   int i;
 994
 995   if(size == 0) {
 996     fprintf(stderr, "Allocating enode of size 0 !\n" );
 997     return NULL;
 998   }
 999   res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1000   res->type = type;
1001   res->size = size;
1002   res->pos = pos;
1003   for(i=0; i<size; i++) {
1004     value_init(res->arr[i].d);
1005     value_set_si(res->arr[i].d,0);
1006     res->arr[i].x.p = 0;
1007   }
1008   return res;
1009 } /* new_enode */
1010
1011 enode *ecopy(enode *e) {
1012
1013   enode *res;
1014   int i;
1015
1016   res = new_enode(e->type,e->size,e->pos);
1017   for(i=0;i<e->size;++i) {
1018     value_assign(res->arr[i].d,e->arr[i].d);
1019     if(value_zero_p(res->arr[i].d))
1020       res->arr[i].x.p = ecopy(e->arr[i].x.p);
1021     else if (EVALUE_IS_DOMAIN(res->arr[i]))
1022       EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1023     else {
1024       value_init(res->arr[i].x.n);
1025       value_assign(res->arr[i].x.n,e->arr[i].x.n);
1026     }
1027   }
1028   return(res);
1029 } /* ecopy */
1030
1031 int eequal(evalue *e1,evalue *e2) {
1032
1033     int i;
1034     enode *p1, *p2;
1035
1036     if (value_ne(e1->d,e2->d))
1037         return 0;
1038
1039     /* e1->d == e2->d */
1040     if (value_notzero_p(e1->d)) {
1041         if (value_ne(e1->x.n,e2->x.n))
1042             return 0;
1043
1044         /* e1->d == e2->d != 0  AND e1->n == e2->n */
1045         return 1;
1046     }
1047
1048     /* e1->d == e2->d == 0 */
1049     p1 = e1->x.p;
1050     p2 = e2->x.p;
1051     if (p1->type != p2->type) return 0;
1052     if (p1->size != p2->size) return 0;
1053     if (p1->pos  != p2->pos) return 0;
1054     for (i=0; i<p1->size; i++)
1055         if (!eequal(&p1->arr[i], &p2->arr[i]) )
1056             return 0;
1057     return 1;
1058 } /* eequal */
1059
1060 void free_evalue_refs(evalue *e) {
1061
1062   enode *p;
1063   int i;
1064
1065   if (EVALUE_IS_DOMAIN(*e)) {
1066     Domain_Free(EVALUE_DOMAIN(*e));
1067     value_clear(e->d);
1068     return;
1069   } else if (value_pos_p(e->d)) {
1070
1071     /* 'e' stores a constant */
1072     value_clear(e->d);
1073     value_clear(e->x.n);
1074     return;
1075   }
1076   assert(value_zero_p(e->d));
1077   value_clear(e->d);
1078   p = e->x.p;
1079   if (!p) return;       /* null pointer */
1080   for (i=0; i<p->size; i++) {
1081     free_evalue_refs(&(p->arr[i]));
1082   }
1083   free(p);
1084   return;
1085 } /* free_evalue_refs */
1086
1087 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
1088                         Vector * val, evalue *res)
1089 {
1090     unsigned nparam = periods->Size;
1091
1092     if (p == nparam) {
1093         double d = compute_evalue(e, val->p);
1094         if (d > 0)
1095             d += .25;
1096         else
1097             d -= .25;
1098         value_set_si(res->d, 1);
1099         value_init(res->x.n);
1100         value_set_double(res->x.n, d);
1101         mpz_fdiv_r(res->x.n, res->x.n, m);
1102         return;
1103     }
1104     if (value_one_p(periods->p[p]))
1105         mod2table_r(e, periods, m, p+1, val, res);
1106     else {
1107         Value tmp;
1108         value_init(tmp);
1109
1110         value_assign(tmp, periods->p[p]);
1111         value_set_si(res->d, 0);
1112         res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
1113         do {
1114             value_decrement(tmp, tmp);
1115             value_assign(val->p[p], tmp);
1116             mod2table_r(e, periods, m, p+1, val,
1117                         &res->x.p->arr[VALUE_TO_INT(tmp)]);
1118         } while (value_pos_p(tmp));
1119
1120         value_clear(tmp);
1121     }
1122 }
1123
1124 void evalue_mod2table(evalue *e, int nparam)
1125 {
1126   enode *p;
1127   int i;
1128
1129   if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
1130     return;
1131   p = e->x.p;
1132   for (i=0; i<p->size; i++) {
1133     evalue_mod2table(&(p->arr[i]), nparam);
1134   }
1135   if (p->type == modulo) {
1136     Vector *periods = Vector_Alloc(nparam);
1137     Vector *val = Vector_Alloc(nparam);
1138     Value tmp;
1139     evalue *ev;
1140     evalue EP, res;
1141
1142     value_init(tmp);
1143     value_set_si(tmp, p->pos);
1144     Vector_Set(periods->p, 1, nparam);
1145     Vector_Set(val->p, 0, nparam);
1146     for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
1147       enode *p = ev->x.p;
1148
1149       assert(p->type == polynomial);
1150       assert(p->size == 2);
1151       assert(value_one_p(p->arr[1].d));
1152       Gcd(tmp, p->arr[1].x.n, &periods->p[p->pos-1]);
1153       value_division(periods->p[p->pos-1], tmp, periods->p[p->pos-1]);
1154     }
1155     value_set_si(tmp, p->pos);
1156     value_init(EP.d);
1157     mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
1158
1159     value_init(res.d);
1160     evalue_set_si(&res, 0, 1);
1161     /* Compute the polynomial using Horner's rule */
1162     for (i=p->size-1;i>1;i--) {
1163       eadd(&p->arr[i], &res);
1164       emul(&EP, &res);
1165     }
1166     eadd(&p->arr[1], &res);
1167
1168     free_evalue_refs(e);
1169     free_evalue_refs(&EP);
1170     *e = res;
1171
1172     value_clear(tmp);
1173     Vector_Free(val);
1174     Vector_Free(periods);
1175   }
1176 } /* evalue_mod2table */
1177
1178 /********************************************************/
1179 /* function in domain                                   */
1180 /*    check if the parameters in list_args              */
1181 /*    verifies the constraints of Domain P              */
1182 /********************************************************/
1183 int in_domain(Polyhedron *P, Value *list_args) {
1184
1185   int col,row;
1186   Value v; /* value of the constraint of a row when
1187                parameters are instanciated*/
1188   Value tmp;
1189
1190   value_init(v);
1191   value_init(tmp);
1192
1193   /*P->Constraint constraint matrice of polyhedron P*/
1194   for(row=0;row<P->NbConstraints;row++) {
1195     value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
1196     for(col=1;col<P->Dimension+1;col++) {
1197       value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
1198       value_addto(v,v,tmp);
1199     }
1200     if (value_notzero_p(P->Constraint[row][0])) {
1201
1202       /*if v is not >=0 then this constraint is not respected */
1203       if (value_neg_p(v)) {
1204 next:
1205         value_clear(v);
1206         value_clear(tmp);
1207         return P->next ? in_domain(P->next, list_args) : 0;
1208       }
1209     }
1210     else {
1211
1212       /*if v is not = 0 then this constraint is not respected */
1213       if (value_notzero_p(v))
1214         goto next;
1215     }
1216   }
1217
1218   /*if not return before this point => all
1219     the constraints are respected */
1220   value_clear(v);
1221   value_clear(tmp);
1222   return 1;
1223 } /* in_domain */
1224
1225 /****************************************************/
1226 /* function compute enode                           */
1227 /*     compute the value of enode p with parameters */
1228 /*     list "list_args                              */
1229 /*     compute the polynomial or the periodic       */
1230 /****************************************************/
1231
1232 static double compute_enode(enode *p, Value *list_args) {
1233
1234   int i;
1235   Value m, param;
1236   double res=0.0;
1237
1238   if (!p)
1239     return(0.);
1240
1241   value_init(m);
1242   value_init(param);
1243
1244   if (p->type == polynomial) {
1245     if (p->size > 1)
1246         value_assign(param,list_args[p->pos-1]);
1247
1248     /* Compute the polynomial using Horner's rule */
1249     for (i=p->size-1;i>0;i--) {
1250       res +=compute_evalue(&p->arr[i],list_args);
1251       res *=VALUE_TO_DOUBLE(param);
1252     }
1253     res +=compute_evalue(&p->arr[0],list_args);
1254   }
1255   else if (p->type == modulo) {
1256     double d = compute_evalue(&p->arr[0], list_args);
1257     if (d > 0)
1258         d += .25;
1259     else
1260         d -= .25;
1261     value_set_double(param, d);
1262     value_set_si(m, p->pos);
1263     mpz_fdiv_r(param, param, m);
1264
1265     /* Compute the polynomial using Horner's rule */
1266     for (i=p->size-1;i>1;i--) {
1267       res +=compute_evalue(&p->arr[i],list_args);
1268       res *=VALUE_TO_DOUBLE(param);
1269     }
1270     res +=compute_evalue(&p->arr[1],list_args);
1271   }
1272   else if (p->type == periodic) {
1273     value_assign(param,list_args[p->pos-1]);
1274
1275     /* Choose the right element of the periodic */
1276     value_absolute(m,param);
1277     value_set_si(param,p->size);
1278     value_modulus(m,m,param);
1279     res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
1280   }
1281   else if (p->type == relation) {
1282     if (fabs(compute_evalue(&p->arr[0], list_args)) < 0.5)
1283       res = compute_evalue(&p->arr[1], list_args);
1284     else if (p->size > 2)
1285       res = compute_evalue(&p->arr[2], list_args);
1286   }
1287   else if (p->type == partition) {
1288     for (i = 0; i < p->size/2; ++i)
1289       if (in_domain(EVALUE_DOMAIN(p->arr[2*i]), list_args)) {
1290         res = compute_evalue(&p->arr[2*i+1], list_args);
1291         break;
1292       }
1293   }
1294   value_clear(m);
1295   value_clear(param);
1296   return res;
1297 } /* compute_enode */
1298
1299 /*************************************************/
1300 /* return the value of Ehrhart Polynomial        */
1301 /* It returns a double, because since it is      */
1302 /* a recursive function, some intermediate value */
1303 /* might not be integral                         */
1304 /*************************************************/
1305
1306 double compute_evalue(evalue *e,Value *list_args) {
1307
1308   double res;
1309
1310   if (value_notzero_p(e->d)) {
1311     if (value_notone_p(e->d))
1312       res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
1313     else
1314       res = VALUE_TO_DOUBLE(e->x.n);
1315   }
1316   else
1317     res = compute_enode(e->x.p,list_args);
1318   return res;
1319 } /* compute_evalue */
1320
1321
1322 /****************************************************/
1323 /* function compute_poly :                          */
1324 /* Check for the good validity domain               */
1325 /* return the number of point in the Polyhedron     */
1326 /* in allocated memory                              */
1327 /* Using the Ehrhart pseudo-polynomial              */
1328 /****************************************************/
1329 Value *compute_poly(Enumeration *en,Value *list_args) {
1330
1331   Value *tmp;
1332   /*    double d; int i; */
1333
1334   tmp = (Value *) malloc (sizeof(Value));
1335   assert(tmp != NULL);
1336   value_init(*tmp);
1337   value_set_si(*tmp,0);
1338
1339   if(!en)
1340     return(tmp);        /* no ehrhart polynomial */
1341   if(en->ValidityDomain) {
1342     if(!en->ValidityDomain->Dimension) { /* no parameters */
1343       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1344       return(tmp);
1345     }
1346   }
1347   else
1348     return(tmp);  /* no Validity Domain */
1349   while(en) {
1350     if(in_domain(en->ValidityDomain,list_args)) {
1351
1352 #ifdef EVAL_EHRHART_DEBUG
1353       Print_Domain(stdout,en->ValidityDomain);
1354       print_evalue(stdout,&en->EP);
1355 #endif
1356
1357       /*                        d = compute_evalue(&en->EP,list_args);
1358                                 i = d;
1359                                 printf("(double)%lf = %d\n", d, i ); */
1360       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1361       return(tmp);
1362     }
1363     else
1364       en=en->next;
1365   }
1366   value_set_si(*tmp,0);
1367   return(tmp); /* no compatible domain with the arguments */
1368 } /* compute_poly */
1369