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 static void mod2table_r(evalue *e, Vector *periods, int p, Vector * val,
 880                         evalue *res)
 881 {
 882     unsigned nparam = periods->Size;
 883
 884     if (p == nparam) {
 885         double d = compute_evalue(e, val->p);
 886         if (d > 0)
 887             d += .25;
 888         else
 889             d -= .25;
 890         value_set_si(res->d, 1);
 891         value_init(res->x.n);
 892         value_set_double(res->x.n, d);
 893         return;
 894     }
 895     if (value_one_p(periods->p[p]))
 896         mod2table_r(e, periods, p+1, val, res);
 897     else {
 898         Value tmp;
 899         value_init(tmp);
 900
 901         value_assign(tmp, periods->p[p]);
 902         value_set_si(res->d, 0);
 903         res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
 904         do {
 905             value_decrement(tmp, tmp);
 906             value_assign(val->p[p], tmp);
 907             mod2table_r(e, periods, p+1, val,
 908                         &res->x.p->arr[VALUE_TO_INT(tmp)]);
 909         } while (value_pos_p(tmp));
 910
 911         value_clear(tmp);
 912     }
 913 }
 914
 915 void evalue_mod2table(evalue *e, int nparam)
 916 {
 917   enode *p;
 918   int i;
 919
 920   if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
 921     return;
 922   p = e->x.p;
 923   for (i=0; i<p->size; i++) {
 924     evalue_mod2table(&(p->arr[i]), nparam);
 925   }
 926   if (p->type == modulo) {
 927     Vector *periods = Vector_Alloc(nparam);
 928     Vector *val = Vector_Alloc(nparam);
 929     Value tmp;
 930     evalue *ev;
 931     evalue EP, res;
 932
 933     value_init(tmp);
 934     value_set_si(tmp, p->pos);
 935     Vector_Set(periods->p, 1, nparam);
 936     Vector_Set(val->p, 0, nparam);
 937     for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
 938       enode *p = ev->x.p;
 939
 940       assert(p->type == polynomial);
 941       assert(p->size == 2);
 942       assert(value_one_p(p->arr[1].d));
 943       Gcd(tmp, p->arr[1].x.n, &periods->p[p->pos-1]);
 944       value_division(periods->p[p->pos-1], tmp, periods->p[p->pos-1]);
 945     }
 946     value_init(EP.d);
 947     mod2table_r(&p->arr[0], periods, 0, val, &EP);
 948
 949     value_init(res.d);
 950     evalue_set_si(&res, 0, 1);
 951     /* Compute the polynomial using Horner's rule */
 952     for (i=p->size-1;i>1;i--) {
 953       eadd(&p->arr[i], &res);
 954       emul(&EP, &res);
 955     }
 956     eadd(&p->arr[1], &res);
 957
 958     free_evalue_refs(e);
 959     free_evalue_refs(&EP);
 960     *e = res;
 961
 962     value_clear(tmp);
 963     Vector_Free(val);
 964     Vector_Free(periods);
 965   }
 966 } /* evalue_mod2table */
 967
 968 /********************************************************/
 969 /* function in domain                                   */
 970 /*    check if the parameters in list_args              */
 971 /*    verifies the constraints of Domain P              */
 972 /********************************************************/
 973 int in_domain(Polyhedron *P, Value *list_args) {
 974
 975   int col,row;
 976   Value v; /* value of the constraint of a row when
 977                parameters are instanciated*/
 978   Value tmp;
 979
 980   value_init(v);
 981   value_init(tmp);
 982
 983   /*P->Constraint constraint matrice of polyhedron P*/
 984   for(row=0;row<P->NbConstraints;row++) {
 985     value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
 986     for(col=1;col<P->Dimension+1;col++) {
 987       value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
 988       value_addto(v,v,tmp);
 989     }
 990     if (value_notzero_p(P->Constraint[row][0])) {
 991
 992       /*if v is not >=0 then this constraint is not respected */
 993       if (value_neg_p(v)) {
 994 next:
 995         value_clear(v);
 996         value_clear(tmp);
 997         return P->next ? in_domain(P->next, list_args) : 0;
 998       }
 999     }
1000     else {
1001
1002       /*if v is not = 0 then this constraint is not respected */
1003       if (value_notzero_p(v))
1004         goto next;
1005     }
1006   }
1007
1008   /*if not return before this point => all
1009     the constraints are respected */
1010   value_clear(v);
1011   value_clear(tmp);
1012   return 1;
1013 } /* in_domain */
1014
1015 /****************************************************/
1016 /* function compute enode                           */
1017 /*     compute the value of enode p with parameters */
1018 /*     list "list_args                              */
1019 /*     compute the polynomial or the periodic       */
1020 /****************************************************/
1021
1022 static double compute_enode(enode *p, Value *list_args) {
1023
1024   int i;
1025   Value m, param;
1026   double res=0.0;
1027
1028   if (!p)
1029     return(0.);
1030
1031   value_init(m);
1032   value_init(param);
1033
1034   if (p->type == polynomial) {
1035     if (p->size > 1)
1036         value_assign(param,list_args[p->pos-1]);
1037
1038     /* Compute the polynomial using Horner's rule */
1039     for (i=p->size-1;i>0;i--) {
1040       res +=compute_evalue(&p->arr[i],list_args);
1041       res *=VALUE_TO_DOUBLE(param);
1042     }
1043     res +=compute_evalue(&p->arr[0],list_args);
1044   }
1045   else if (p->type == modulo) {
1046     double d = compute_evalue(&p->arr[0], list_args);
1047     if (d > 0)
1048         d += .25;
1049     else
1050         d -= .25;
1051     value_set_double(param, d);
1052     value_set_si(m, p->pos);
1053     mpz_fdiv_r(param, param, m);
1054
1055     /* Compute the polynomial using Horner's rule */
1056     for (i=p->size-1;i>1;i--) {
1057       res +=compute_evalue(&p->arr[i],list_args);
1058       res *=VALUE_TO_DOUBLE(param);
1059     }
1060     res +=compute_evalue(&p->arr[1],list_args);
1061   }
1062   else if (p->type == periodic) {
1063     value_assign(param,list_args[p->pos-1]);
1064
1065     /* Choose the right element of the periodic */
1066     value_absolute(m,param);
1067     value_set_si(param,p->size);
1068     value_modulus(m,m,param);
1069     res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
1070   }
1071   else if (p->type == relation) {
1072     if (fabs(compute_evalue(&p->arr[0], list_args)) < 0.5)
1073       res = compute_evalue(&p->arr[1], list_args);
1074   }
1075   value_clear(m);
1076   value_clear(param);
1077   return res;
1078 } /* compute_enode */
1079
1080 /*************************************************/
1081 /* return the value of Ehrhart Polynomial        */
1082 /* It returns a double, because since it is      */
1083 /* a recursive function, some intermediate value */
1084 /* might not be integral                         */
1085 /*************************************************/
1086
1087 double compute_evalue(evalue *e,Value *list_args) {
1088
1089   double res;
1090
1091   if (value_notzero_p(e->d)) {
1092     if (value_notone_p(e->d))
1093       res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
1094     else
1095       res = VALUE_TO_DOUBLE(e->x.n);
1096   }
1097   else
1098     res = compute_enode(e->x.p,list_args);
1099   return res;
1100 } /* compute_evalue */
1101
1102
1103 /****************************************************/
1104 /* function compute_poly :                          */
1105 /* Check for the good validity domain               */
1106 /* return the number of point in the Polyhedron     */
1107 /* in allocated memory                              */
1108 /* Using the Ehrhart pseudo-polynomial              */
1109 /****************************************************/
1110 Value *compute_poly(Enumeration *en,Value *list_args) {
1111
1112   Value *tmp;
1113   /*    double d; int i; */
1114
1115   tmp = (Value *) malloc (sizeof(Value));
1116   assert(tmp != NULL);
1117   value_init(*tmp);
1118   value_set_si(*tmp,0);
1119
1120   if(!en)
1121     return(tmp);        /* no ehrhart polynomial */
1122   if(en->ValidityDomain) {
1123     if(!en->ValidityDomain->Dimension) { /* no parameters */
1124       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1125       return(tmp);
1126     }
1127   }
1128   else
1129     return(tmp);  /* no Validity Domain */
1130   while(en) {
1131     if(in_domain(en->ValidityDomain,list_args)) {
1132
1133 #ifdef EVAL_EHRHART_DEBUG
1134       Print_Domain(stdout,en->ValidityDomain);
1135       print_evalue(stdout,&en->EP);
1136 #endif
1137
1138       /*                        d = compute_evalue(&en->EP,list_args);
1139                                 i = d;
1140                                 printf("(double)%lf = %d\n", d, i ); */
1141       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1142       return(tmp);
1143     }
1144     else
1145       en=en->next;
1146   }
1147   value_set_si(*tmp,0);
1148   return(tmp); /* no compatible domain with the arguments */
1149 } /* compute_poly */
1150