ev_operations.c

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