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>=0;i--) {
 110             if (value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n))
 111
 112                 /* Zero coefficient */
 113                 continue;
 114             else
 115                 break;
 116         }
 117         if (i==-1) p->size = 1;
 118         else if (i+1<p->size) p->size = i+1;
 119
 120         /* Try to reduce its strength */
 121         if (p->size == 1) {
 122             value_clear(e->d);
 123             memcpy(e,&p->arr[0],sizeof(evalue));
 124             free(p);
 125         }
 126     }
 127     else if (p->type==modulo) {
 128
 129         /* Try to reduce the degree */
 130         for (i=p->size-1;i>=1;i--) {
 131             if (value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n))
 132
 133                 /* Zero coefficient */
 134                 continue;
 135             else
 136                 break;
 137         }
 138         if (i==0) p->size = 2;
 139         else if (i+1<p->size) p->size = i+1;
 140
 141         /* Try to reduce its strength */
 142         if (p->size == 2) {
 143             value_clear(e->d);
 144             memcpy(e,&p->arr[1],sizeof(evalue));
 145             free(p);
 146         }
 147     }
 148 } /* reduce_evalue */
 149
 150 void print_evalue(FILE *DST,evalue *e,char **pname) {
 151
 152   if(value_notzero_p(e->d)) {
 153     if(value_notone_p(e->d)) {
 154       value_print(DST,VALUE_FMT,e->x.n);
 155       fprintf(DST,"/");
 156       value_print(DST,VALUE_FMT,e->d);
 157     }
 158     else {
 159       value_print(DST,VALUE_FMT,e->x.n);
 160     }
 161   }
 162   else
 163     print_enode(DST,e->x.p,pname);
 164   return;
 165 } /* print_evalue */
 166
 167 void print_enode(FILE *DST,enode *p,char **pname) {
 168
 169   int i;
 170
 171   if (!p) {
 172     fprintf(DST, "NULL");
 173     return;
 174   }
 175   if (p->type == evector) {
 176     fprintf(DST, "{ ");
 177     for (i=0; i<p->size; i++) {
 178       print_evalue(DST, &p->arr[i], pname);
 179       if (i!=(p->size-1))
 180         fprintf(DST, ", ");
 181     }
 182     fprintf(DST, " }\n");
 183   }
 184   else if (p->type == polynomial) {
 185     fprintf(DST, "( ");
 186     for (i=p->size-1; i>=0; i--) {
 187       print_evalue(DST, &p->arr[i], pname);
 188       if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
 189       else if (i>1)
 190         fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
 191     }
 192     fprintf(DST, " )\n");
 193   }
 194   else if (p->type == periodic) {
 195     fprintf(DST, "[ ");
 196     for (i=0; i<p->size; i++) {
 197       print_evalue(DST, &p->arr[i], pname);
 198       if (i!=(p->size-1)) fprintf(DST, ", ");
 199     }
 200     fprintf(DST," ]_%s", pname[p->pos-1]);
 201   }
 202   else if (p->type == modulo) {
 203     fprintf(DST, "( ");
 204     for (i=p->size-1; i>=1; i--) {
 205       print_evalue(DST, &p->arr[i], pname);
 206       if (i >= 2) {
 207         fprintf(DST, " * ");
 208         if (i > 2)
 209           fprintf(DST, "(");
 210         fprintf(DST, "(");
 211         print_evalue(DST, &p->arr[0], pname);
 212         fprintf(DST, ") mod %d", p->pos);
 213         if (i>2)
 214           fprintf(DST, ")^%d + ", i-1);
 215         else
 216           fprintf(DST, " + ", i-1);
 217       }
 218     }
 219     fprintf(DST, " )\n");
 220   }
 221   return;
 222 } /* print_enode */
 223
 224 static int mod_term_smaller(evalue *e1, evalue *e2)
 225 {
 226     if (value_notzero_p(e1->d)) {
 227         if (value_zero_p(e2->d))
 228             return 1;
 229         return value_lt(e1->x.n, e2->x.n);
 230     }
 231     if (value_notzero_p(e2->d))
 232         return 0;
 233     if (e1->x.p->pos < e2->x.p->pos)
 234         return 1;
 235     else if (e1->x.p->pos > e2->x.p->pos)
 236         return 0;
 237     else
 238         return mod_term_smaller(&e1->x.p->arr[0], &e2->x.p->arr[0]);
 239 }
 240
 241 static void eadd_rev(evalue *e1, evalue *res)
 242 {
 243     evalue ev;
 244     value_init(ev.d);
 245     evalue_copy(&ev, e1);
 246     eadd(res, &ev);
 247     free_evalue_refs(res);
 248     *res = ev;
 249 }
 250
 251 static void eadd_rev_cst (evalue *e1, evalue *res)
 252 {
 253     evalue ev;
 254     value_init(ev.d);
 255     evalue_copy(&ev, e1);
 256     eadd(res, &ev.x.p->arr[ev.x.p->type==modulo]);
 257     free_evalue_refs(res);
 258     *res = ev;
 259 }
 260
 261 void eadd(evalue *e1,evalue *res) {
 262
 263  int i;
 264     if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
 265          /* Add two rational numbers */
 266          Value g,m1,m2;
 267          value_init(g);
 268          value_init(m1);
 269          value_init(m2);
 270
 271          value_multiply(m1,e1->x.n,res->d);
 272          value_multiply(m2,res->x.n,e1->d);
 273          value_addto(res->x.n,m1,m2);
 274          value_multiply(res->d,e1->d,res->d);
 275          Gcd(res->x.n,res->d,&g);
 276          if (value_notone_p(g)) {
 277               value_division(res->d,res->d,g);
 278               value_division(res->x.n,res->x.n,g);
 279          }
 280          value_clear(g); value_clear(m1); value_clear(m2);
 281          return ;
 282      }
 283      else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
 284           switch (res->x.p->type) {
 285           case polynomial:
 286               /* Add the constant to the constant term of a polynomial*/
 287                eadd(e1, &res->x.p->arr[0]);
 288                return ;
 289           case periodic:
 290               /* Add the constant to all elements of a periodic number */
 291               for (i=0; i<res->x.p->size; i++) {
 292                   eadd(e1, &res->x.p->arr[i]);
 293               }
 294               return ;
 295           case evector:
 296               fprintf(stderr, "eadd: cannot add const with vector\n");
 297               return;
 298           case modulo:
 299                eadd(e1, &res->x.p->arr[1]);
 300                return ;
 301           }
 302      }
 303      /* add polynomial or periodic to constant
 304       * you have to exchange e1 and res, before doing addition */
 305
 306      else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
 307           eadd_rev(e1, res);
 308           return;
 309      }
 310      else {   // ((e1->d==0) && (res->d==0))
 311                  if ((e1->x.p->type != res->x.p->type) ) {
 312                       /* adding to evalues of different type. two cases are possible
 313                        * res is periodic and e1 is polynomial, you have to exchange
 314                        * e1 and res then to add e1 to the constant term of res */
 315                      if (e1->x.p->type == polynomial) {
 316                           eadd_rev_cst(e1, res);
 317                      }
 318                      else if (res->x.p->type == polynomial) {
 319                           /* res is polynomial and e1 is periodic,
 320                             add e1 to the constant term of res */
 321
 322                           eadd(e1,&res->x.p->arr[0]);
 323                      } else
 324                         assert(0);
 325
 326                      return;
 327                  }
 328                  else if (e1->x.p->pos != res->x.p->pos ||
 329                             (res->x.p->type == modulo &&
 330                              !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
 331                  /* adding evalues of different position (i.e function of different unknowns
 332                   * to case are possible  */
 333
 334                         switch (res->x.p->type) {
 335                         case modulo:
 336                             if(mod_term_smaller(res, e1))
 337                                 eadd(e1,&res->x.p->arr[1]);
 338                             else
 339                                 eadd_rev_cst(e1, res);
 340                             return;
 341                         case polynomial: //  res and e1 are polynomials
 342                                //  add e1 to the constant term of res
 343
 344                             if(res->x.p->pos < e1->x.p->pos)
 345                                eadd(e1,&res->x.p->arr[0]);
 346                             else
 347                                 eadd_rev_cst(e1, res);
 348                               // value_clear(g); value_clear(m1); value_clear(m2);
 349                                return;
 350                         case periodic:  // res and e1 are pointers to periodic numbers
 351                                   //add e1 to all elements of res
 352
 353                             if(res->x.p->pos < e1->x.p->pos)
 354                                   for (i=0;i<res->x.p->size;i++) {
 355                                        eadd(e1,&res->x.p->arr[i]);
 356                                   }
 357                             else
 358                                 eadd_rev(e1, res);
 359                             return;
 360                         }
 361                  }
 362
 363
 364                  //same type , same pos  and same size
 365                  if (e1->x.p->size == res->x.p->size) {
 366                       // add any element in e1 to the corresponding element in res
 367                       if (res->x.p->type == modulo)
 368                         assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
 369                       i = res->x.p->type == modulo ? 1 : 0;
 370                       for (; i<res->x.p->size; i++) {
 371                             eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
 372                       }
 373                       return ;
 374                 }
 375
 376                 /* Sizes are different */
 377                 if (res->x.p->type==polynomial) {
 378                     /* VIN100: if e1-size > res-size you have to copy e1 in a   */
 379                     /* new enode and add res to that new node. If you do not do */
 380                     /* that, you lose the the upper weight part of e1 !         */
 381
 382                      if(e1->x.p->size > res->x.p->size) {
 383                           enode *tmp;
 384                           tmp = ecopy(e1->x.p);
 385                           for(i=0;i<res->x.p->size;++i) {
 386                               eadd(&res->x.p->arr[i], &tmp->arr[i]);
 387                               //  free_evalue_refs(&res->x.p->arr[i]);
 388                           }
 389                           res->x.p = tmp;
 390
 391                      }
 392                      else {
 393
 394                         for (i=0; i<e1->x.p->size ; i++) {
 395                              eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
 396                         }
 397
 398                         return ;
 399                      }
 400                 }
 401
 402     /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
 403                 else if (res->x.p->type==periodic) {
 404                       /* you have to create a new evalue 'ne' in whitch size equals to the lcm
 405                        of the sizes of e1 and res, then to copy res periodicaly in ne, after
 406                        to add periodicaly elements of e1 to elements of ne, and finaly to
 407                        return ne. */
 408                        int x,y,p;
 409                        Value ex, ey ,ep;
 410                        evalue *ne;
 411                        value_init(ex); value_init(ey);value_init(ep);
 412                        x=e1->x.p->size;
 413                        y= res->x.p->size;
 414                        value_set_si(ex,e1->x.p->size);
 415                        value_set_si(ey,res->x.p->size);
 416                        value_assign (ep,*Lcm(ex,ey));
 417                        p=(int)mpz_get_si(ep);
 418                        ne= (evalue *) malloc (sizeof(evalue));
 419                        value_init(ne->d);
 420                        value_set_si( ne->d,0);
 421
 422                        ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
 423                        for(i=0;i<p;i++)  {
 424                               value_assign(ne->x.p->arr[i].d, res->x.p->arr[i%y].d);
 425                               if (value_notzero_p(ne->x.p->arr[i].d))   {
 426                                 value_init(ne->x.p->arr[i].x.n);
 427                                 value_assign(ne->x.p->arr[i].x.n, res->x.p->arr[i%y].x.n);
 428                               }
 429                               else {
 430                                   ne->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%y].x.p);
 431                               }
 432                        }
 433                        for(i=0;i<p;i++)  {
 434                             eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
 435                        }
 436
 437                       value_assign(res->d, ne->d);
 438                       res->x.p=ne->x.p;
 439
 440                         return ;
 441                 }
 442                 else { /* evector */
 443                      fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
 444                      return ;
 445                 }
 446      }
 447      return ;
 448  }/* eadd  */
 449
 450 static void emul_rev (evalue *e1, evalue *res)
 451 {
 452     evalue ev;
 453     value_init(ev.d);
 454     evalue_copy(&ev, e1);
 455     emul(res, &ev);
 456     free_evalue_refs(res);
 457     *res = ev;
 458 }
 459
 460 static void emul_poly (evalue *e1, evalue *res)
 461 {
 462     int i, j, o = res->x.p->type == modulo;
 463     evalue tmp;
 464     int size=(e1->x.p->size + res->x.p->size - o - 1);
 465     value_init(tmp.d);
 466     value_set_si(tmp.d,0);
 467     tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
 468     if (o)
 469         evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
 470     for (i=o; i < e1->x.p->size; i++) {
 471         evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
 472         emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
 473     }
 474     for (; i<size; i++)
 475         evalue_set_si(&tmp.x.p->arr[i], 0, 1);
 476     for (i=o+1; i<res->x.p->size; i++)
 477         for (j=o; j<e1->x.p->size; j++) {
 478             evalue ev;
 479             value_init(ev.d);
 480             evalue_copy(&ev, &e1->x.p->arr[j]);
 481             emul(&res->x.p->arr[i], &ev);
 482             eadd(&ev, &tmp.x.p->arr[i+j-o]);
 483             free_evalue_refs(&ev);
 484         }
 485     free_evalue_refs(res);
 486     *res = tmp;
 487 }
 488
 489 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
 490  * do a copy of "res" befor calling this function if you nead it after. The vector type of
 491  * evalues is not treated here */
 492
 493 void emul (evalue *e1, evalue *res ){
 494     int i,j;
 495
 496 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
 497     fprintf(stderr, "emul: do not proced on evector type !\n");
 498     return;
 499 }
 500
 501    if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
 502        switch(e1->x.p->type) {
 503        case polynomial:
 504            switch(res->x.p->type) {
 505            case polynomial:
 506                if(e1->x.p->pos == res->x.p->pos) {
 507                /* Product of two polynomials of the same variable */
 508                     emul_poly(e1, res);
 509                     return;
 510                }
 511                else {
 512                   /* Product of two polynomials of different variables */
 513
 514                 if(res->x.p->pos < e1->x.p->pos)
 515                      for( i=0; i<res->x.p->size ; i++)
 516                             emul(e1, &res->x.p->arr[i]);
 517                 else
 518                     emul_rev(e1, res);
 519
 520                  return ;
 521                }
 522            case periodic:
 523            case modulo:
 524                 /* Product of a polynomial and a periodic or modulo */
 525                 emul_rev(e1, res);
 526                 return;
 527            }
 528        case periodic:
 529            switch(res->x.p->type) {
 530            case periodic:
 531                 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
 532                  /* Product of two periodics of the same parameter and period */
 533
 534                      for(i=0; i<res->x.p->size;i++)
 535                            emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
 536
 537                      return;
 538                 }
 539                 else{
 540                   if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
 541                    /* Product of two periodics of the same parameter and different periods */
 542                     evalue *newp;
 543                     Value x,y,z;
 544                     int ix,iy,lcm;
 545                     value_init(x); value_init(y);value_init(z);
 546                     ix=e1->x.p->size;
 547                     iy=res->x.p->size;
 548                     value_set_si(x,e1->x.p->size);
 549                     value_set_si(y,res->x.p->size);
 550                     value_assign (z,*Lcm(x,y));
 551                     lcm=(int)mpz_get_si(z);
 552                     newp= (evalue *) malloc (sizeof(evalue));
 553                     value_init(newp->d);
 554                     value_set_si( newp->d,0);
 555                     newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
 556                         for(i=0;i<lcm;i++)  {
 557                            value_assign(newp->x.p->arr[i].d, res->x.p->arr[i%iy].d);
 558                            if (value_notzero_p(newp->x.p->arr[i].d))   {
 559                            value_assign(newp->x.p->arr[i].x.n, res->x.p->arr[i%iy].x.n);
 560                            }
 561                            else {
 562                            newp->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%iy].x.p);
 563                            }
 564
 565                         }
 566                         for(i=0;i<lcm;i++)
 567                             emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
 568
 569                         value_assign(res->d,newp->d);
 570                         res->x.p=newp->x.p;
 571
 572                           value_clear(x); value_clear(y);value_clear(z);
 573                           return ;
 574                   }
 575                   else {
 576                      /* Product of two periodics of different parameters */
 577
 578                         for(i=0; i<res->x.p->size; i++)
 579                             emul(e1, &(res->x.p->arr[i]));
 580
 581                         return;
 582                   }
 583                 }
 584            case polynomial:
 585                   /* Product of a periodic and a polynomial */
 586
 587                        for(i=0; i<res->x.p->size ; i++)
 588                             emul(e1, &(res->x.p->arr[i]));
 589
 590                        return;
 591
 592            }
 593        case modulo:
 594             switch(res->x.p->type) {
 595             case polynomial:
 596                 for(i=0; i<res->x.p->size ; i++)
 597                     emul(e1, &(res->x.p->arr[i]));
 598                 return;
 599             case periodic:
 600                 assert(0);
 601             case modulo:
 602                 if (e1->x.p->pos == res->x.p->pos &&
 603                             eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))
 604                     emul_poly(e1, res);
 605                 else {
 606                     if(mod_term_smaller(res, e1))
 607                         for(i=1; i<res->x.p->size ; i++)
 608                             emul(e1, &(res->x.p->arr[i]));
 609                     else
 610                         emul_rev(e1, res);
 611                     return;
 612                 }
 613             }
 614        }
 615    }
 616    else {
 617        if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
 618            /* Product of two rational numbers */
 619
 620             Value g;
 621             value_init(g);
 622             value_multiply(res->d,e1->d,res->d);
 623             value_multiply(res->x.n,e1->x.n,res->x.n );
 624             Gcd(res->x.n, res->d,&g);
 625            if (value_notone_p(g)) {
 626                value_division(res->d,res->d,g);
 627                value_division(res->x.n,res->x.n,g);
 628            }
 629            value_clear(g);
 630            return ;
 631        }
 632        else {
 633              if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
 634               /* Product of an expression (polynomial or peririodic) and a rational number */
 635
 636                 emul_rev(e1, res);
 637                  return ;
 638              }
 639              else {
 640                /* Product of a rationel number and an expression (polynomial or peririodic) */
 641
 642                    i = res->x.p->type == modulo ? 1 : 0;
 643                    for (; i<res->x.p->size; i++)
 644                       emul(e1, &res->x.p->arr[i]);
 645
 646                  return ;
 647              }
 648        }
 649    }
 650
 651    return ;
 652 }
 653
 654 void evalue_copy(evalue *dst, evalue *src)
 655 {
 656     value_assign(dst->d, src->d);
 657     if(value_notzero_p(src->d)) {
 658          value_init(dst->x.n);
 659          value_assign(dst->x.n, src->x.n);
 660     } else
 661          dst->x.p = ecopy(src->x.p);
 662 }
 663
 664 enode *new_enode(enode_type type,int size,int pos) {
 665
 666   enode *res;
 667   int i;
 668
 669   if(size == 0) {
 670     fprintf(stderr, "Allocating enode of size 0 !\n" );
 671     return NULL;
 672   }
 673   res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
 674   res->type = type;
 675   res->size = size;
 676   res->pos = pos;
 677   for(i=0; i<size; i++) {
 678     value_init(res->arr[i].d);
 679     value_set_si(res->arr[i].d,0);
 680     res->arr[i].x.p = 0;
 681   }
 682   return res;
 683 } /* new_enode */
 684
 685 enode *ecopy(enode *e) {
 686
 687   enode *res;
 688   int i;
 689
 690   res = new_enode(e->type,e->size,e->pos);
 691   for(i=0;i<e->size;++i) {
 692     value_assign(res->arr[i].d,e->arr[i].d);
 693     if(value_zero_p(res->arr[i].d))
 694       res->arr[i].x.p = ecopy(e->arr[i].x.p);
 695     else {
 696       value_init(res->arr[i].x.n);
 697       value_assign(res->arr[i].x.n,e->arr[i].x.n);
 698     }
 699   }
 700   return(res);
 701 } /* ecopy */
 702
 703 int eequal(evalue *e1,evalue *e2) {
 704
 705     int i;
 706     enode *p1, *p2;
 707
 708     if (value_ne(e1->d,e2->d))
 709         return 0;
 710
 711     /* e1->d == e2->d */
 712     if (value_notzero_p(e1->d)) {
 713         if (value_ne(e1->x.n,e2->x.n))
 714             return 0;
 715
 716         /* e1->d == e2->d != 0  AND e1->n == e2->n */
 717         return 1;
 718     }
 719
 720     /* e1->d == e2->d == 0 */
 721     p1 = e1->x.p;
 722     p2 = e2->x.p;
 723     if (p1->type != p2->type) return 0;
 724     if (p1->size != p2->size) return 0;
 725     if (p1->pos  != p2->pos) return 0;
 726     for (i=0; i<p1->size; i++)
 727         if (!eequal(&p1->arr[i], &p2->arr[i]) )
 728             return 0;
 729     return 1;
 730 } /* eequal */
 731
 732 void free_evalue_refs(evalue *e) {
 733
 734   enode *p;
 735   int i;
 736
 737   if (value_notzero_p(e->d)) {
 738
 739     /* 'e' stores a constant */
 740     value_clear(e->d);
 741     value_clear(e->x.n);
 742     return;
 743   }
 744   value_clear(e->d);
 745   p = e->x.p;
 746   if (!p) return;       /* null pointer */
 747   for (i=0; i<p->size; i++) {
 748     free_evalue_refs(&(p->arr[i]));
 749   }
 750   free(p);
 751   return;
 752 } /* free_evalue_refs */
 753
 754 /****************************************************/
 755 /* function compute enode                           */
 756 /*     compute the value of enode p with parameters */
 757 /*     list "list_args                              */
 758 /*     compute the polynomial or the periodic       */
 759 /****************************************************/
 760
 761 static double compute_enode(enode *p, Value *list_args) {
 762
 763   int i;
 764   Value m, param;
 765   double res=0.0;
 766
 767   if (!p)
 768     return(0.);
 769
 770   value_init(m);
 771   value_init(param);
 772
 773   if (p->type == polynomial) {
 774     if (p->size > 1)
 775         value_assign(param,list_args[p->pos-1]);
 776
 777     /* Compute the polynomial using Horner's rule */
 778     for (i=p->size-1;i>0;i--) {
 779       res +=compute_evalue(&p->arr[i],list_args);
 780       res *=VALUE_TO_DOUBLE(param);
 781     }
 782     res +=compute_evalue(&p->arr[0],list_args);
 783   }
 784   else if (p->type == modulo) {
 785     double d = compute_evalue(&p->arr[0], list_args);
 786     if (d > 0)
 787         d += .25;
 788     else
 789         d -= .25;
 790     value_set_double(param, d);
 791     value_set_si(m, p->pos);
 792     mpz_fdiv_r(param, param, m);
 793
 794     /* Compute the polynomial using Horner's rule */
 795     for (i=p->size-1;i>1;i--) {
 796       res +=compute_evalue(&p->arr[i],list_args);
 797       res *=VALUE_TO_DOUBLE(param);
 798     }
 799     res +=compute_evalue(&p->arr[1],list_args);
 800   }
 801   else if (p->type == periodic) {
 802     value_assign(param,list_args[p->pos-1]);
 803
 804     /* Choose the right element of the periodic */
 805     value_absolute(m,param);
 806     value_set_si(param,p->size);
 807     value_modulus(m,m,param);
 808     res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
 809   }
 810   value_clear(m);
 811   value_clear(param);
 812   return res;
 813 } /* compute_enode */
 814
 815 /*************************************************/
 816 /* return the value of Ehrhart Polynomial        */
 817 /* It returns a double, because since it is      */
 818 /* a recursive function, some intermediate value */
 819 /* might not be integral                         */
 820 /*************************************************/
 821
 822 double compute_evalue(evalue *e,Value *list_args) {
 823
 824   double res;
 825
 826   if (value_notzero_p(e->d)) {
 827     if (value_notone_p(e->d))
 828       res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
 829     else
 830       res = VALUE_TO_DOUBLE(e->x.n);
 831   }
 832   else
 833     res = compute_enode(e->x.p,list_args);
 834   return res;
 835 } /* compute_evalue */
 836
 837
 838 /****************************************************/
 839 /* function compute_poly :                          */
 840 /* Check for the good validity domain               */
 841 /* return the number of point in the Polyhedron     */
 842 /* in allocated memory                              */
 843 /* Using the Ehrhart pseudo-polynomial              */
 844 /****************************************************/
 845 Value *compute_poly(Enumeration *en,Value *list_args) {
 846
 847   Value *tmp;
 848   /*    double d; int i; */
 849
 850   tmp = (Value *) malloc (sizeof(Value));
 851   assert(tmp != NULL);
 852   value_init(*tmp);
 853   value_set_si(*tmp,0);
 854
 855   if(!en)
 856     return(tmp);        /* no ehrhart polynomial */
 857   if(en->ValidityDomain) {
 858     if(!en->ValidityDomain->Dimension) { /* no parameters */
 859       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
 860       return(tmp);
 861     }
 862   }
 863   else
 864     return(tmp);  /* no Validity Domain */
 865   while(en) {
 866     if(in_domain(en->ValidityDomain,list_args)) {
 867
 868 #ifdef EVAL_EHRHART_DEBUG
 869       Print_Domain(stdout,en->ValidityDomain);
 870       print_evalue(stdout,&en->EP);
 871 #endif
 872
 873       /*                        d = compute_evalue(&en->EP,list_args);
 874                                 i = d;
 875                                 printf("(double)%lf = %d\n", d, i ); */
 876       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
 877       return(tmp);
 878     }
 879     else
 880       en=en->next;
 881   }
 882   value_set_si(*tmp,0);
 883   return(tmp); /* no compatible domain with the arguments */
 884 } /* compute_poly */
 885