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