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             s[n].E = res->x.p->arr[2*i+1];
 330             s[n].D = fd;
 331             ++n;
 332         } else
 333             free_evalue_refs(&res->x.p->arr[2*i+1]);
 334         if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
 335             Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
 336         value_clear(res->x.p->arr[2*i].d);
 337     }
 338
 339     free(res->x.p);
 340     res->x.p = new_enode(partition, 2*n, -1);
 341     for (j = 0; j < n; ++j) {
 342         EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
 343         value_clear(res->x.p->arr[2*j+1].d);
 344         res->x.p->arr[2*j+1] = s[j].E;
 345     }
 346
 347     free(s);
 348 }
 349
 350 void eadd(evalue *e1,evalue *res) {
 351
 352  int i;
 353     if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
 354          /* Add two rational numbers */
 355          Value g,m1,m2;
 356          value_init(g);
 357          value_init(m1);
 358          value_init(m2);
 359
 360          value_multiply(m1,e1->x.n,res->d);
 361          value_multiply(m2,res->x.n,e1->d);
 362          value_addto(res->x.n,m1,m2);
 363          value_multiply(res->d,e1->d,res->d);
 364          Gcd(res->x.n,res->d,&g);
 365          if (value_notone_p(g)) {
 366               value_division(res->d,res->d,g);
 367               value_division(res->x.n,res->x.n,g);
 368          }
 369          value_clear(g); value_clear(m1); value_clear(m2);
 370          return ;
 371      }
 372      else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
 373           switch (res->x.p->type) {
 374           case polynomial:
 375               /* Add the constant to the constant term of a polynomial*/
 376                eadd(e1, &res->x.p->arr[0]);
 377                return ;
 378           case periodic:
 379               /* Add the constant to all elements of a periodic number */
 380               for (i=0; i<res->x.p->size; i++) {
 381                   eadd(e1, &res->x.p->arr[i]);
 382               }
 383               return ;
 384           case evector:
 385               fprintf(stderr, "eadd: cannot add const with vector\n");
 386               return;
 387           case modulo:
 388                eadd(e1, &res->x.p->arr[1]);
 389                return ;
 390           case partition:
 391                 assert(EVALUE_IS_ZERO(*e1));
 392                 break;                          /* Do nothing */
 393           default:
 394                 assert(0);
 395           }
 396      }
 397      /* add polynomial or periodic to constant
 398       * you have to exchange e1 and res, before doing addition */
 399
 400      else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
 401           eadd_rev(e1, res);
 402           return;
 403      }
 404      else {   // ((e1->d==0) && (res->d==0))
 405         assert(!((e1->x.p->type == partition) ^
 406                  (res->x.p->type == partition)));
 407         if (e1->x.p->type == partition) {
 408             eadd_partitions(e1, res);
 409             return;
 410         }
 411                  if ((e1->x.p->type != res->x.p->type) ) {
 412                       /* adding to evalues of different type. two cases are possible
 413                        * res is periodic and e1 is polynomial, you have to exchange
 414                        * e1 and res then to add e1 to the constant term of res */
 415                      if (e1->x.p->type == polynomial) {
 416                           eadd_rev_cst(e1, res);
 417                      }
 418                      else if (res->x.p->type == polynomial) {
 419                           /* res is polynomial and e1 is periodic,
 420                             add e1 to the constant term of res */
 421
 422                           eadd(e1,&res->x.p->arr[0]);
 423                      } else
 424                         assert(0);
 425
 426                      return;
 427                  }
 428                  else if (e1->x.p->pos != res->x.p->pos ||
 429                             (res->x.p->type == modulo &&
 430                              !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
 431                  /* adding evalues of different position (i.e function of different unknowns
 432                   * to case are possible  */
 433
 434                         switch (res->x.p->type) {
 435                         case modulo:
 436                             if(mod_term_smaller(res, e1))
 437                                 eadd(e1,&res->x.p->arr[1]);
 438                             else
 439                                 eadd_rev_cst(e1, res);
 440                             return;
 441                         case polynomial: //  res and e1 are polynomials
 442                                //  add e1 to the constant term of res
 443
 444                             if(res->x.p->pos < e1->x.p->pos)
 445                                eadd(e1,&res->x.p->arr[0]);
 446                             else
 447                                 eadd_rev_cst(e1, res);
 448                               // value_clear(g); value_clear(m1); value_clear(m2);
 449                                return;
 450                         case periodic:  // res and e1 are pointers to periodic numbers
 451                                   //add e1 to all elements of res
 452
 453                             if(res->x.p->pos < e1->x.p->pos)
 454                                   for (i=0;i<res->x.p->size;i++) {
 455                                        eadd(e1,&res->x.p->arr[i]);
 456                                   }
 457                             else
 458                                 eadd_rev(e1, res);
 459                             return;
 460                         }
 461                  }
 462
 463
 464                  //same type , same pos  and same size
 465                  if (e1->x.p->size == res->x.p->size) {
 466                       // add any element in e1 to the corresponding element in res
 467                       if (res->x.p->type == modulo)
 468                         assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
 469                       i = res->x.p->type == modulo ? 1 : 0;
 470                       for (; i<res->x.p->size; i++) {
 471                             eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
 472                       }
 473                       return ;
 474                 }
 475
 476                 /* Sizes are different */
 477                 switch(res->x.p->type) {
 478                 case polynomial:
 479                 case modulo:
 480                     /* VIN100: if e1-size > res-size you have to copy e1 in a   */
 481                     /* new enode and add res to that new node. If you do not do */
 482                     /* that, you lose the the upper weight part of e1 !         */
 483
 484                      if(e1->x.p->size > res->x.p->size)
 485                           eadd_rev(e1, res);
 486                      else {
 487
 488                         if (res->x.p->type == modulo)
 489                             assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
 490                         i = res->x.p->type == modulo ? 1 : 0;
 491                         for (; i<e1->x.p->size ; i++) {
 492                              eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
 493                         }
 494
 495                         return ;
 496                      }
 497                      break;
 498
 499     /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
 500                 case periodic:
 501                 {
 502                       /* you have to create a new evalue 'ne' in whitch size equals to the lcm
 503                        of the sizes of e1 and res, then to copy res periodicaly in ne, after
 504                        to add periodicaly elements of e1 to elements of ne, and finaly to
 505                        return ne. */
 506                        int x,y,p;
 507                        Value ex, ey ,ep;
 508                        evalue *ne;
 509                        value_init(ex); value_init(ey);value_init(ep);
 510                        x=e1->x.p->size;
 511                        y= res->x.p->size;
 512                        value_set_si(ex,e1->x.p->size);
 513                        value_set_si(ey,res->x.p->size);
 514                        value_assign (ep,*Lcm(ex,ey));
 515                        p=(int)mpz_get_si(ep);
 516                        ne= (evalue *) malloc (sizeof(evalue));
 517                        value_init(ne->d);
 518                        value_set_si( ne->d,0);
 519
 520                        ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
 521                        for(i=0;i<p;i++)  {
 522                               value_assign(ne->x.p->arr[i].d, res->x.p->arr[i%y].d);
 523                               if (value_notzero_p(ne->x.p->arr[i].d))   {
 524                                 value_init(ne->x.p->arr[i].x.n);
 525                                 value_assign(ne->x.p->arr[i].x.n, res->x.p->arr[i%y].x.n);
 526                               }
 527                               else {
 528                                   ne->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%y].x.p);
 529                               }
 530                        }
 531                        for(i=0;i<p;i++)  {
 532                             eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
 533                        }
 534
 535                       value_assign(res->d, ne->d);
 536                       res->x.p=ne->x.p;
 537
 538                         return ;
 539                 }
 540                 case evector:
 541                      fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
 542                      return ;
 543                 }
 544      }
 545      return ;
 546  }/* eadd  */
 547
 548 static void emul_rev (evalue *e1, evalue *res)
 549 {
 550     evalue ev;
 551     value_init(ev.d);
 552     evalue_copy(&ev, e1);
 553     emul(res, &ev);
 554     free_evalue_refs(res);
 555     *res = ev;
 556 }
 557
 558 static void emul_poly (evalue *e1, evalue *res)
 559 {
 560     int i, j, o = res->x.p->type == modulo;
 561     evalue tmp;
 562     int size=(e1->x.p->size + res->x.p->size - o - 1);
 563     value_init(tmp.d);
 564     value_set_si(tmp.d,0);
 565     tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
 566     if (o)
 567         evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
 568     for (i=o; i < e1->x.p->size; i++) {
 569         evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
 570         emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
 571     }
 572     for (; i<size; i++)
 573         evalue_set_si(&tmp.x.p->arr[i], 0, 1);
 574     for (i=o+1; i<res->x.p->size; i++)
 575         for (j=o; j<e1->x.p->size; j++) {
 576             evalue ev;
 577             value_init(ev.d);
 578             evalue_copy(&ev, &e1->x.p->arr[j]);
 579             emul(&res->x.p->arr[i], &ev);
 580             eadd(&ev, &tmp.x.p->arr[i+j-o]);
 581             free_evalue_refs(&ev);
 582         }
 583     free_evalue_refs(res);
 584     *res = tmp;
 585 }
 586
 587 void emul_partitions (evalue *e1,evalue *res)
 588 {
 589     int n, i, j, k;
 590     Polyhedron *d;
 591     struct section *s;
 592     s = (struct section *)
 593             malloc((e1->x.p->size/2) * (res->x.p->size/2) *
 594                    sizeof(struct section));
 595     assert(s);
 596
 597     n = 0;
 598     for (i = 0; i < res->x.p->size/2; ++i) {
 599         for (j = 0; j < e1->x.p->size/2; ++j) {
 600             d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
 601                                    EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
 602             if (emptyQ(d)) {
 603                 Domain_Free(d);
 604                 continue;
 605             }
 606
 607             /* This code is only needed because the partitions
 608                are not true partitions.
 609              */
 610             for (k = 0; k < n; ++k) {
 611                 if (DomainIncludes(s[k].D, d))
 612                     break;
 613                 if (DomainIncludes(d, s[k].D)) {
 614                     Domain_Free(s[k].D);
 615                     free_evalue_refs(&s[k].E);
 616                     if (n > k)
 617                         s[k] = s[--n];
 618                     --k;
 619                 }
 620             }
 621             if (k < n) {
 622                 Domain_Free(d);
 623                 continue;
 624             }
 625
 626             value_init(s[n].E.d);
 627             evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
 628             emul(&e1->x.p->arr[2*j+1], &s[n].E);
 629             s[n].D = d;
 630             ++n;
 631         }
 632         Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
 633         value_clear(res->x.p->arr[2*i].d);
 634         free_evalue_refs(&res->x.p->arr[2*i+1]);
 635     }
 636
 637     free(res->x.p);
 638     res->x.p = new_enode(partition, 2*n, -1);
 639     for (j = 0; j < n; ++j) {
 640         EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
 641         value_clear(res->x.p->arr[2*j+1].d);
 642         res->x.p->arr[2*j+1] = s[j].E;
 643     }
 644
 645     free(s);
 646 }
 647
 648 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
 649  * do a copy of "res" befor calling this function if you nead it after. The vector type of
 650  * evalues is not treated here */
 651
 652 void emul (evalue *e1, evalue *res ){
 653     int i,j;
 654
 655 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
 656     fprintf(stderr, "emul: do not proced on evector type !\n");
 657     return;
 658 }
 659
 660     if (value_zero_p(e1->d) && e1->x.p->type == partition) {
 661         if (value_zero_p(res->d) && res->x.p->type == partition)
 662             emul_partitions(e1, res);
 663         else
 664             emul_rev(e1, res);
 665     } else if (value_zero_p(res->d) && res->x.p->type == partition) {
 666         for (i = 0; i < res->x.p->size/2; ++i)
 667             emul(e1, &res->x.p->arr[2*i+1]);
 668     } else
 669    if (value_zero_p(res->d) && res->x.p->type == relation)
 670         emul(e1, &res->x.p->arr[1]);
 671    else
 672    if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
 673        switch(e1->x.p->type) {
 674        case polynomial:
 675            switch(res->x.p->type) {
 676            case polynomial:
 677                if(e1->x.p->pos == res->x.p->pos) {
 678                /* Product of two polynomials of the same variable */
 679                     emul_poly(e1, res);
 680                     return;
 681                }
 682                else {
 683                   /* Product of two polynomials of different variables */
 684
 685                 if(res->x.p->pos < e1->x.p->pos)
 686                      for( i=0; i<res->x.p->size ; i++)
 687                             emul(e1, &res->x.p->arr[i]);
 688                 else
 689                     emul_rev(e1, res);
 690
 691                  return ;
 692                }
 693            case periodic:
 694            case modulo:
 695                 /* Product of a polynomial and a periodic or modulo */
 696                 emul_rev(e1, res);
 697                 return;
 698            }
 699        case periodic:
 700            switch(res->x.p->type) {
 701            case periodic:
 702                 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
 703                  /* Product of two periodics of the same parameter and period */
 704
 705                      for(i=0; i<res->x.p->size;i++)
 706                            emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
 707
 708                      return;
 709                 }
 710                 else{
 711                   if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
 712                    /* Product of two periodics of the same parameter and different periods */
 713                     evalue *newp;
 714                     Value x,y,z;
 715                     int ix,iy,lcm;
 716                     value_init(x); value_init(y);value_init(z);
 717                     ix=e1->x.p->size;
 718                     iy=res->x.p->size;
 719                     value_set_si(x,e1->x.p->size);
 720                     value_set_si(y,res->x.p->size);
 721                     value_assign (z,*Lcm(x,y));
 722                     lcm=(int)mpz_get_si(z);
 723                     newp= (evalue *) malloc (sizeof(evalue));
 724                     value_init(newp->d);
 725                     value_set_si( newp->d,0);
 726                     newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
 727                         for(i=0;i<lcm;i++)  {
 728                            value_assign(newp->x.p->arr[i].d, res->x.p->arr[i%iy].d);
 729                            if (value_notzero_p(newp->x.p->arr[i].d))   {
 730                            value_assign(newp->x.p->arr[i].x.n, res->x.p->arr[i%iy].x.n);
 731                            }
 732                            else {
 733                            newp->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%iy].x.p);
 734                            }
 735
 736                         }
 737                         for(i=0;i<lcm;i++)
 738                             emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
 739
 740                         value_assign(res->d,newp->d);
 741                         res->x.p=newp->x.p;
 742
 743                           value_clear(x); value_clear(y);value_clear(z);
 744                           return ;
 745                   }
 746                   else {
 747                      /* Product of two periodics of different parameters */
 748
 749                         for(i=0; i<res->x.p->size; i++)
 750                             emul(e1, &(res->x.p->arr[i]));
 751
 752                         return;
 753                   }
 754                 }
 755            case polynomial:
 756                   /* Product of a periodic and a polynomial */
 757
 758                        for(i=0; i<res->x.p->size ; i++)
 759                             emul(e1, &(res->x.p->arr[i]));
 760
 761                        return;
 762
 763            }
 764        case modulo:
 765             switch(res->x.p->type) {
 766             case polynomial:
 767                 for(i=0; i<res->x.p->size ; i++)
 768                     emul(e1, &(res->x.p->arr[i]));
 769                 return;
 770             case periodic:
 771                 assert(0);
 772             case modulo:
 773                 if (e1->x.p->pos == res->x.p->pos &&
 774                             eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
 775                     if (e1->x.p->pos != 2)
 776                         emul_poly(e1, res);
 777                     else {
 778                         evalue tmp;
 779                         /* x mod 2 == (x mod 2)^2 */
 780                         /* a0 b0 + (a0 b1 + a1 b0 + a1 b1) (x mod 2) */
 781                         assert(e1->x.p->size == 3);
 782                         assert(res->x.p->size == 3);
 783                         value_init(tmp.d);
 784                         evalue_copy(&tmp, &res->x.p->arr[1]);
 785                         eadd(&res->x.p->arr[2], &tmp);
 786                         emul(&e1->x.p->arr[2], &tmp);
 787                         emul(&e1->x.p->arr[1], res);
 788                         eadd(&tmp, &res->x.p->arr[2]);
 789                         free_evalue_refs(&tmp);
 790                     }
 791                 } else {
 792                     if(mod_term_smaller(res, e1))
 793                         for(i=1; i<res->x.p->size ; i++)
 794                             emul(e1, &(res->x.p->arr[i]));
 795                     else
 796                         emul_rev(e1, res);
 797                     return;
 798                 }
 799             }
 800             break;
 801        case relation:
 802             emul_rev(e1, res);
 803             break;
 804        default:
 805             assert(0);
 806        }
 807    }
 808    else {
 809        if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
 810            /* Product of two rational numbers */
 811
 812             Value g;
 813             value_init(g);
 814             value_multiply(res->d,e1->d,res->d);
 815             value_multiply(res->x.n,e1->x.n,res->x.n );
 816             Gcd(res->x.n, res->d,&g);
 817            if (value_notone_p(g)) {
 818                value_division(res->d,res->d,g);
 819                value_division(res->x.n,res->x.n,g);
 820            }
 821            value_clear(g);
 822            return ;
 823        }
 824        else {
 825              if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
 826               /* Product of an expression (polynomial or peririodic) and a rational number */
 827
 828                 emul_rev(e1, res);
 829                  return ;
 830              }
 831              else {
 832                /* Product of a rationel number and an expression (polynomial or peririodic) */
 833
 834                    i = res->x.p->type == modulo ? 1 : 0;
 835                    for (; i<res->x.p->size; i++)
 836                       emul(e1, &res->x.p->arr[i]);
 837
 838                  return ;
 839              }
 840        }
 841    }
 842
 843    return ;
 844 }
 845
 846 void evalue_copy(evalue *dst, evalue *src)
 847 {
 848     value_assign(dst->d, src->d);
 849     if(value_notzero_p(src->d)) {
 850          value_init(dst->x.n);
 851          value_assign(dst->x.n, src->x.n);
 852     } else
 853          dst->x.p = ecopy(src->x.p);
 854 }
 855
 856 enode *new_enode(enode_type type,int size,int pos) {
 857
 858   enode *res;
 859   int i;
 860
 861   if(size == 0) {
 862     fprintf(stderr, "Allocating enode of size 0 !\n" );
 863     return NULL;
 864   }
 865   res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
 866   res->type = type;
 867   res->size = size;
 868   res->pos = pos;
 869   for(i=0; i<size; i++) {
 870     value_init(res->arr[i].d);
 871     value_set_si(res->arr[i].d,0);
 872     res->arr[i].x.p = 0;
 873   }
 874   return res;
 875 } /* new_enode */
 876
 877 enode *ecopy(enode *e) {
 878
 879   enode *res;
 880   int i;
 881
 882   res = new_enode(e->type,e->size,e->pos);
 883   for(i=0;i<e->size;++i) {
 884     value_assign(res->arr[i].d,e->arr[i].d);
 885     if(value_zero_p(res->arr[i].d))
 886       res->arr[i].x.p = ecopy(e->arr[i].x.p);
 887     else if (EVALUE_IS_DOMAIN(res->arr[i]))
 888       EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
 889     else {
 890       value_init(res->arr[i].x.n);
 891       value_assign(res->arr[i].x.n,e->arr[i].x.n);
 892     }
 893   }
 894   return(res);
 895 } /* ecopy */
 896
 897 int eequal(evalue *e1,evalue *e2) {
 898
 899     int i;
 900     enode *p1, *p2;
 901
 902     if (value_ne(e1->d,e2->d))
 903         return 0;
 904
 905     /* e1->d == e2->d */
 906     if (value_notzero_p(e1->d)) {
 907         if (value_ne(e1->x.n,e2->x.n))
 908             return 0;
 909
 910         /* e1->d == e2->d != 0  AND e1->n == e2->n */
 911         return 1;
 912     }
 913
 914     /* e1->d == e2->d == 0 */
 915     p1 = e1->x.p;
 916     p2 = e2->x.p;
 917     if (p1->type != p2->type) return 0;
 918     if (p1->size != p2->size) return 0;
 919     if (p1->pos  != p2->pos) return 0;
 920     for (i=0; i<p1->size; i++)
 921         if (!eequal(&p1->arr[i], &p2->arr[i]) )
 922             return 0;
 923     return 1;
 924 } /* eequal */
 925
 926 void free_evalue_refs(evalue *e) {
 927
 928   enode *p;
 929   int i;
 930
 931   if (EVALUE_IS_DOMAIN(*e)) {
 932     Domain_Free(EVALUE_DOMAIN(*e));
 933     value_clear(e->d);
 934     return;
 935   } else if (value_pos_p(e->d)) {
 936
 937     /* 'e' stores a constant */
 938     value_clear(e->d);
 939     value_clear(e->x.n);
 940     return;
 941   }
 942   assert(value_zero_p(e->d));
 943   value_clear(e->d);
 944   p = e->x.p;
 945   if (!p) return;       /* null pointer */
 946   for (i=0; i<p->size; i++) {
 947     free_evalue_refs(&(p->arr[i]));
 948   }
 949   free(p);
 950   return;
 951 } /* free_evalue_refs */
 952
 953 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
 954                         Vector * val, evalue *res)
 955 {
 956     unsigned nparam = periods->Size;
 957
 958     if (p == nparam) {
 959         double d = compute_evalue(e, val->p);
 960         if (d > 0)
 961             d += .25;
 962         else
 963             d -= .25;
 964         value_set_si(res->d, 1);
 965         value_init(res->x.n);
 966         value_set_double(res->x.n, d);
 967         mpz_fdiv_r(res->x.n, res->x.n, m);
 968         return;
 969     }
 970     if (value_one_p(periods->p[p]))
 971         mod2table_r(e, periods, m, p+1, val, res);
 972     else {
 973         Value tmp;
 974         value_init(tmp);
 975
 976         value_assign(tmp, periods->p[p]);
 977         value_set_si(res->d, 0);
 978         res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
 979         do {
 980             value_decrement(tmp, tmp);
 981             value_assign(val->p[p], tmp);
 982             mod2table_r(e, periods, m, p+1, val,
 983                         &res->x.p->arr[VALUE_TO_INT(tmp)]);
 984         } while (value_pos_p(tmp));
 985
 986         value_clear(tmp);
 987     }
 988 }
 989
 990 void evalue_mod2table(evalue *e, int nparam)
 991 {
 992   enode *p;
 993   int i;
 994
 995   if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
 996     return;
 997   p = e->x.p;
 998   for (i=0; i<p->size; i++) {
 999     evalue_mod2table(&(p->arr[i]), nparam);
1000   }
1001   if (p->type == modulo) {
1002     Vector *periods = Vector_Alloc(nparam);
1003     Vector *val = Vector_Alloc(nparam);
1004     Value tmp;
1005     evalue *ev;
1006     evalue EP, res;
1007
1008     value_init(tmp);
1009     value_set_si(tmp, p->pos);
1010     Vector_Set(periods->p, 1, nparam);
1011     Vector_Set(val->p, 0, nparam);
1012     for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
1013       enode *p = ev->x.p;
1014
1015       assert(p->type == polynomial);
1016       assert(p->size == 2);
1017       assert(value_one_p(p->arr[1].d));
1018       Gcd(tmp, p->arr[1].x.n, &periods->p[p->pos-1]);
1019       value_division(periods->p[p->pos-1], tmp, periods->p[p->pos-1]);
1020     }
1021     value_set_si(tmp, p->pos);
1022     value_init(EP.d);
1023     mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
1024
1025     value_init(res.d);
1026     evalue_set_si(&res, 0, 1);
1027     /* Compute the polynomial using Horner's rule */
1028     for (i=p->size-1;i>1;i--) {
1029       eadd(&p->arr[i], &res);
1030       emul(&EP, &res);
1031     }
1032     eadd(&p->arr[1], &res);
1033
1034     free_evalue_refs(e);
1035     free_evalue_refs(&EP);
1036     *e = res;
1037
1038     value_clear(tmp);
1039     Vector_Free(val);
1040     Vector_Free(periods);
1041   }
1042 } /* evalue_mod2table */
1043
1044 /********************************************************/
1045 /* function in domain                                   */
1046 /*    check if the parameters in list_args              */
1047 /*    verifies the constraints of Domain P              */
1048 /********************************************************/
1049 int in_domain(Polyhedron *P, Value *list_args) {
1050
1051   int col,row;
1052   Value v; /* value of the constraint of a row when
1053                parameters are instanciated*/
1054   Value tmp;
1055
1056   value_init(v);
1057   value_init(tmp);
1058
1059   /*P->Constraint constraint matrice of polyhedron P*/
1060   for(row=0;row<P->NbConstraints;row++) {
1061     value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
1062     for(col=1;col<P->Dimension+1;col++) {
1063       value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
1064       value_addto(v,v,tmp);
1065     }
1066     if (value_notzero_p(P->Constraint[row][0])) {
1067
1068       /*if v is not >=0 then this constraint is not respected */
1069       if (value_neg_p(v)) {
1070 next:
1071         value_clear(v);
1072         value_clear(tmp);
1073         return P->next ? in_domain(P->next, list_args) : 0;
1074       }
1075     }
1076     else {
1077
1078       /*if v is not = 0 then this constraint is not respected */
1079       if (value_notzero_p(v))
1080         goto next;
1081     }
1082   }
1083
1084   /*if not return before this point => all
1085     the constraints are respected */
1086   value_clear(v);
1087   value_clear(tmp);
1088   return 1;
1089 } /* in_domain */
1090
1091 /****************************************************/
1092 /* function compute enode                           */
1093 /*     compute the value of enode p with parameters */
1094 /*     list "list_args                              */
1095 /*     compute the polynomial or the periodic       */
1096 /****************************************************/
1097
1098 static double compute_enode(enode *p, Value *list_args) {
1099
1100   int i;
1101   Value m, param;
1102   double res=0.0;
1103
1104   if (!p)
1105     return(0.);
1106
1107   value_init(m);
1108   value_init(param);
1109
1110   if (p->type == polynomial) {
1111     if (p->size > 1)
1112         value_assign(param,list_args[p->pos-1]);
1113
1114     /* Compute the polynomial using Horner's rule */
1115     for (i=p->size-1;i>0;i--) {
1116       res +=compute_evalue(&p->arr[i],list_args);
1117       res *=VALUE_TO_DOUBLE(param);
1118     }
1119     res +=compute_evalue(&p->arr[0],list_args);
1120   }
1121   else if (p->type == modulo) {
1122     double d = compute_evalue(&p->arr[0], list_args);
1123     if (d > 0)
1124         d += .25;
1125     else
1126         d -= .25;
1127     value_set_double(param, d);
1128     value_set_si(m, p->pos);
1129     mpz_fdiv_r(param, param, m);
1130
1131     /* Compute the polynomial using Horner's rule */
1132     for (i=p->size-1;i>1;i--) {
1133       res +=compute_evalue(&p->arr[i],list_args);
1134       res *=VALUE_TO_DOUBLE(param);
1135     }
1136     res +=compute_evalue(&p->arr[1],list_args);
1137   }
1138   else if (p->type == periodic) {
1139     value_assign(param,list_args[p->pos-1]);
1140
1141     /* Choose the right element of the periodic */
1142     value_absolute(m,param);
1143     value_set_si(param,p->size);
1144     value_modulus(m,m,param);
1145     res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
1146   }
1147   else if (p->type == relation) {
1148     if (fabs(compute_evalue(&p->arr[0], list_args)) < 0.5)
1149       res = compute_evalue(&p->arr[1], list_args);
1150   }
1151   value_clear(m);
1152   value_clear(param);
1153   return res;
1154 } /* compute_enode */
1155
1156 /*************************************************/
1157 /* return the value of Ehrhart Polynomial        */
1158 /* It returns a double, because since it is      */
1159 /* a recursive function, some intermediate value */
1160 /* might not be integral                         */
1161 /*************************************************/
1162
1163 double compute_evalue(evalue *e,Value *list_args) {
1164
1165   double res;
1166
1167   if (value_notzero_p(e->d)) {
1168     if (value_notone_p(e->d))
1169       res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
1170     else
1171       res = VALUE_TO_DOUBLE(e->x.n);
1172   }
1173   else
1174     res = compute_enode(e->x.p,list_args);
1175   return res;
1176 } /* compute_evalue */
1177
1178
1179 /****************************************************/
1180 /* function compute_poly :                          */
1181 /* Check for the good validity domain               */
1182 /* return the number of point in the Polyhedron     */
1183 /* in allocated memory                              */
1184 /* Using the Ehrhart pseudo-polynomial              */
1185 /****************************************************/
1186 Value *compute_poly(Enumeration *en,Value *list_args) {
1187
1188   Value *tmp;
1189   /*    double d; int i; */
1190
1191   tmp = (Value *) malloc (sizeof(Value));
1192   assert(tmp != NULL);
1193   value_init(*tmp);
1194   value_set_si(*tmp,0);
1195
1196   if(!en)
1197     return(tmp);        /* no ehrhart polynomial */
1198   if(en->ValidityDomain) {
1199     if(!en->ValidityDomain->Dimension) { /* no parameters */
1200       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1201       return(tmp);
1202     }
1203   }
1204   else
1205     return(tmp);  /* no Validity Domain */
1206   while(en) {
1207     if(in_domain(en->ValidityDomain,list_args)) {
1208
1209 #ifdef EVAL_EHRHART_DEBUG
1210       Print_Domain(stdout,en->ValidityDomain);
1211       print_evalue(stdout,&en->EP);
1212 #endif
1213
1214       /*                        d = compute_evalue(&en->EP,list_args);
1215                                 i = d;
1216                                 printf("(double)%lf = %d\n", d, i ); */
1217       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1218       return(tmp);
1219     }
1220     else
1221       en=en->next;
1222   }
1223   value_set_si(*tmp,0);
1224   return(tmp); /* no compatible domain with the arguments */
1225 } /* compute_poly */
1226