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