3 #include "ev_operations.h"
5 void evalue_set_si(evalue
*ev
, int n
, int d
) {
6 value_set_si(ev
->d
, d
);
8 value_set_si(ev
->x
.n
, n
);
11 void evalue_set(evalue
*ev
, Value n
, Value d
) {
12 value_assign(ev
->d
, d
);
14 value_assign(ev
->x
.n
, n
);
17 void aep_evalue(evalue
*e
, int *ref
) {
22 if (value_notzero_p(e
->d
))
23 return; /* a rational number, its already reduced */
25 return; /* hum... an overflow probably occured */
27 /* First check the components of p */
28 for (i
=0;i
<p
->size
;i
++)
29 aep_evalue(&p
->arr
[i
],ref
);
32 if (p
->type
!= modulo
)
33 p
->pos
= ref
[p
->pos
-1]+1;
38 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
44 if (value_notzero_p(e
->d
))
45 return; /* a rational number, its already reduced */
47 return; /* hum... an overflow probably occured */
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
])) {
58 /* Transform the references in e, using ref */
62 } /* addeliminatedparams_evalue */
64 void reduce_evalue (evalue
*e
) {
69 if (value_notzero_p(e
->d
))
70 return; /* a rational number, its already reduced */
72 return; /* hum... an overflow probably occured */
74 /* First reduce the components of p */
75 for (i
=0; i
<p
->size
; i
++)
76 reduce_evalue(&p
->arr
[i
]);
78 if (p
->type
==periodic
) {
80 /* Try to reduce the period */
81 for (i
=1; i
<=(p
->size
)/2; i
++) {
82 if ((p
->size
% i
)==0) {
84 /* Can we reduce the size to i ? */
86 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
87 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
90 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
94 you_lose
: /* OK, lets not do it */
99 /* Try to reduce its strength */
102 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
106 else if (p
->type
==polynomial
) {
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
)))
112 /* Zero coefficient */
113 free_evalue_refs(&(p
->arr
[i
]));
118 /* Try to reduce its strength */
121 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
125 else if (p
->type
==modulo
) {
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
)))
131 /* Zero coefficient */
132 free_evalue_refs(&(p
->arr
[i
]));
137 /* Try to reduce its strength */
140 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
141 free_evalue_refs(&(p
->arr
[0]));
145 } /* reduce_evalue */
147 void print_evalue(FILE *DST
,evalue
*e
,char **pname
) {
149 if(value_notzero_p(e
->d
)) {
150 if(value_notone_p(e
->d
)) {
151 value_print(DST
,VALUE_FMT
,e
->x
.n
);
153 value_print(DST
,VALUE_FMT
,e
->d
);
156 value_print(DST
,VALUE_FMT
,e
->x
.n
);
160 print_enode(DST
,e
->x
.p
,pname
);
164 void print_enode(FILE *DST
,enode
*p
,char **pname
) {
169 fprintf(DST
, "NULL");
175 for (i
=0; i
<p
->size
; i
++) {
176 print_evalue(DST
, &p
->arr
[i
], pname
);
180 fprintf(DST
, " }\n");
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]);
188 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
190 fprintf(DST
, " )\n");
194 for (i
=0; i
<p
->size
; i
++) {
195 print_evalue(DST
, &p
->arr
[i
], pname
);
196 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
198 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
202 for (i
=p
->size
-1; i
>=1; i
--) {
203 print_evalue(DST
, &p
->arr
[i
], pname
);
209 print_evalue(DST
, &p
->arr
[0], pname
);
210 fprintf(DST
, ") mod %d", p
->pos
);
212 fprintf(DST
, ")^%d + ", i
-1);
214 fprintf(DST
, " + ", i
-1);
217 fprintf(DST
, " )\n");
221 print_evalue(DST
, &p
->arr
[0], pname
);
222 fprintf(DST
, "= 0 ] * \n");
223 print_evalue(DST
, &p
->arr
[1], pname
);
226 for (i
=0; i
<p
->size
/2; i
++) {
227 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), pname
);
228 print_evalue(DST
, &p
->arr
[2*i
+1], pname
);
237 static int mod_term_smaller(evalue
*e1
, evalue
*e2
)
239 if (value_notzero_p(e1
->d
)) {
240 if (value_zero_p(e2
->d
))
242 return value_lt(e1
->x
.n
, e2
->x
.n
);
244 if (value_notzero_p(e2
->d
))
246 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
248 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
251 return mod_term_smaller(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
254 static void eadd_rev(evalue
*e1
, evalue
*res
)
258 evalue_copy(&ev
, e1
);
260 free_evalue_refs(res
);
264 static void eadd_rev_cst (evalue
*e1
, evalue
*res
)
268 evalue_copy(&ev
, e1
);
269 eadd(res
, &ev
.x
.p
->arr
[ev
.x
.p
->type
==modulo
]);
270 free_evalue_refs(res
);
274 struct section
{ Polyhedron
* D
; evalue E
; };
276 void eadd_partitions (evalue
*e1
,evalue
*res
)
281 s
= (struct section
*)
282 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
283 sizeof(struct section
));
287 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
288 assert(res
->x
.p
->size
>= 2);
289 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
290 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
292 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
294 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
303 value_init(s
[n
].E
.d
);
304 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
308 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
309 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
310 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
312 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
313 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
319 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
320 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
322 value_init(s
[n
].E
.d
);
323 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
324 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
329 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
330 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
331 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
338 res
->x
.p
= new_enode(partition
, 2*n
, -1);
339 for (j
= 0; j
< n
; ++j
) {
340 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
341 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
347 void eadd(evalue
*e1
,evalue
*res
) {
350 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
351 /* Add two rational numbers */
357 value_multiply(m1
,e1
->x
.n
,res
->d
);
358 value_multiply(m2
,res
->x
.n
,e1
->d
);
359 value_addto(res
->x
.n
,m1
,m2
);
360 value_multiply(res
->d
,e1
->d
,res
->d
);
361 Gcd(res
->x
.n
,res
->d
,&g
);
362 if (value_notone_p(g
)) {
363 value_division(res
->d
,res
->d
,g
);
364 value_division(res
->x
.n
,res
->x
.n
,g
);
366 value_clear(g
); value_clear(m1
); value_clear(m2
);
369 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
370 switch (res
->x
.p
->type
) {
372 /* Add the constant to the constant term of a polynomial*/
373 eadd(e1
, &res
->x
.p
->arr
[0]);
376 /* Add the constant to all elements of a periodic number */
377 for (i
=0; i
<res
->x
.p
->size
; i
++) {
378 eadd(e1
, &res
->x
.p
->arr
[i
]);
382 fprintf(stderr
, "eadd: cannot add const with vector\n");
385 eadd(e1
, &res
->x
.p
->arr
[1]);
388 assert(EVALUE_IS_ZERO(*e1
));
389 break; /* Do nothing */
394 /* add polynomial or periodic to constant
395 * you have to exchange e1 and res, before doing addition */
397 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
401 else { // ((e1->d==0) && (res->d==0))
402 assert(!((e1
->x
.p
->type
== partition
) ^
403 (res
->x
.p
->type
== partition
)));
404 if (e1
->x
.p
->type
== partition
) {
405 eadd_partitions(e1
, res
);
408 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
409 /* adding to evalues of different type. two cases are possible
410 * res is periodic and e1 is polynomial, you have to exchange
411 * e1 and res then to add e1 to the constant term of res */
412 if (e1
->x
.p
->type
== polynomial
) {
413 eadd_rev_cst(e1
, res
);
415 else if (res
->x
.p
->type
== polynomial
) {
416 /* res is polynomial and e1 is periodic,
417 add e1 to the constant term of res */
419 eadd(e1
,&res
->x
.p
->arr
[0]);
425 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
426 (res
->x
.p
->type
== modulo
&&
427 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
428 /* adding evalues of different position (i.e function of different unknowns
429 * to case are possible */
431 switch (res
->x
.p
->type
) {
433 if(mod_term_smaller(res
, e1
))
434 eadd(e1
,&res
->x
.p
->arr
[1]);
436 eadd_rev_cst(e1
, res
);
438 case polynomial
: // res and e1 are polynomials
439 // add e1 to the constant term of res
441 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
442 eadd(e1
,&res
->x
.p
->arr
[0]);
444 eadd_rev_cst(e1
, res
);
445 // value_clear(g); value_clear(m1); value_clear(m2);
447 case periodic
: // res and e1 are pointers to periodic numbers
448 //add e1 to all elements of res
450 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
451 for (i
=0;i
<res
->x
.p
->size
;i
++) {
452 eadd(e1
,&res
->x
.p
->arr
[i
]);
461 //same type , same pos and same size
462 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
463 // add any element in e1 to the corresponding element in res
464 if (res
->x
.p
->type
== modulo
)
465 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
466 i
= res
->x
.p
->type
== modulo
? 1 : 0;
467 for (; i
<res
->x
.p
->size
; i
++) {
468 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
473 /* Sizes are different */
474 switch(res
->x
.p
->type
) {
477 /* VIN100: if e1-size > res-size you have to copy e1 in a */
478 /* new enode and add res to that new node. If you do not do */
479 /* that, you lose the the upper weight part of e1 ! */
481 if(e1
->x
.p
->size
> res
->x
.p
->size
)
485 if (res
->x
.p
->type
== modulo
)
486 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
487 i
= res
->x
.p
->type
== modulo
? 1 : 0;
488 for (; i
<e1
->x
.p
->size
; i
++) {
489 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
496 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
499 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
500 of the sizes of e1 and res, then to copy res periodicaly in ne, after
501 to add periodicaly elements of e1 to elements of ne, and finaly to
506 value_init(ex
); value_init(ey
);value_init(ep
);
509 value_set_si(ex
,e1
->x
.p
->size
);
510 value_set_si(ey
,res
->x
.p
->size
);
511 value_assign (ep
,*Lcm(ex
,ey
));
512 p
=(int)mpz_get_si(ep
);
513 ne
= (evalue
*) malloc (sizeof(evalue
));
515 value_set_si( ne
->d
,0);
517 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
519 value_assign(ne
->x
.p
->arr
[i
].d
, res
->x
.p
->arr
[i
%y
].d
);
520 if (value_notzero_p(ne
->x
.p
->arr
[i
].d
)) {
521 value_init(ne
->x
.p
->arr
[i
].x
.n
);
522 value_assign(ne
->x
.p
->arr
[i
].x
.n
, res
->x
.p
->arr
[i
%y
].x
.n
);
525 ne
->x
.p
->arr
[i
].x
.p
=ecopy(res
->x
.p
->arr
[i
%y
].x
.p
);
529 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
532 value_assign(res
->d
, ne
->d
);
538 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
545 static void emul_rev (evalue
*e1
, evalue
*res
)
549 evalue_copy(&ev
, e1
);
551 free_evalue_refs(res
);
555 static void emul_poly (evalue
*e1
, evalue
*res
)
557 int i
, j
, o
= res
->x
.p
->type
== modulo
;
559 int size
=(e1
->x
.p
->size
+ res
->x
.p
->size
- o
- 1);
561 value_set_si(tmp
.d
,0);
562 tmp
.x
.p
=new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
564 evalue_copy(&tmp
.x
.p
->arr
[0], &e1
->x
.p
->arr
[0]);
565 for (i
=o
; i
< e1
->x
.p
->size
; i
++) {
566 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
567 emul(&res
->x
.p
->arr
[o
], &tmp
.x
.p
->arr
[i
]);
570 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
571 for (i
=o
+1; i
<res
->x
.p
->size
; i
++)
572 for (j
=o
; j
<e1
->x
.p
->size
; j
++) {
575 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
576 emul(&res
->x
.p
->arr
[i
], &ev
);
577 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-o
]);
578 free_evalue_refs(&ev
);
580 free_evalue_refs(res
);
584 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
585 * do a copy of "res" befor calling this function if you nead it after. The vector type of
586 * evalues is not treated here */
588 void emul (evalue
*e1
, evalue
*res
){
591 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
592 fprintf(stderr
, "emul: do not proced on evector type !\n");
596 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
)
597 emul(e1
, &res
->x
.p
->arr
[1]);
599 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
600 switch(e1
->x
.p
->type
) {
602 switch(res
->x
.p
->type
) {
604 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
605 /* Product of two polynomials of the same variable */
610 /* Product of two polynomials of different variables */
612 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
613 for( i
=0; i
<res
->x
.p
->size
; i
++)
614 emul(e1
, &res
->x
.p
->arr
[i
]);
622 /* Product of a polynomial and a periodic or modulo */
627 switch(res
->x
.p
->type
) {
629 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
630 /* Product of two periodics of the same parameter and period */
632 for(i
=0; i
<res
->x
.p
->size
;i
++)
633 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
638 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
639 /* Product of two periodics of the same parameter and different periods */
643 value_init(x
); value_init(y
);value_init(z
);
646 value_set_si(x
,e1
->x
.p
->size
);
647 value_set_si(y
,res
->x
.p
->size
);
648 value_assign (z
,*Lcm(x
,y
));
649 lcm
=(int)mpz_get_si(z
);
650 newp
= (evalue
*) malloc (sizeof(evalue
));
652 value_set_si( newp
->d
,0);
653 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
655 value_assign(newp
->x
.p
->arr
[i
].d
, res
->x
.p
->arr
[i
%iy
].d
);
656 if (value_notzero_p(newp
->x
.p
->arr
[i
].d
)) {
657 value_assign(newp
->x
.p
->arr
[i
].x
.n
, res
->x
.p
->arr
[i
%iy
].x
.n
);
660 newp
->x
.p
->arr
[i
].x
.p
=ecopy(res
->x
.p
->arr
[i
%iy
].x
.p
);
665 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
667 value_assign(res
->d
,newp
->d
);
670 value_clear(x
); value_clear(y
);value_clear(z
);
674 /* Product of two periodics of different parameters */
676 for(i
=0; i
<res
->x
.p
->size
; i
++)
677 emul(e1
, &(res
->x
.p
->arr
[i
]));
683 /* Product of a periodic and a polynomial */
685 for(i
=0; i
<res
->x
.p
->size
; i
++)
686 emul(e1
, &(res
->x
.p
->arr
[i
]));
692 switch(res
->x
.p
->type
) {
694 for(i
=0; i
<res
->x
.p
->size
; i
++)
695 emul(e1
, &(res
->x
.p
->arr
[i
]));
700 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
701 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
702 if (e1
->x
.p
->pos
!= 2)
705 /* x mod 2 == (x mod 2)^2 */
706 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1) (x mod 2) */
707 assert(e1
->x
.p
->size
== 3);
708 assert(res
->x
.p
->size
== 3);
711 evalue_copy(&tmp
, &res
->x
.p
->arr
[1]);
712 eadd(&res
->x
.p
->arr
[2], &tmp
);
713 emul(&e1
->x
.p
->arr
[2], &tmp
);
714 emul(&e1
->x
.p
->arr
[1], res
);
715 eadd(&tmp
, &res
->x
.p
->arr
[2]);
716 free_evalue_refs(&tmp
);
719 if(mod_term_smaller(res
, e1
))
720 for(i
=1; i
<res
->x
.p
->size
; i
++)
721 emul(e1
, &(res
->x
.p
->arr
[i
]));
736 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
737 /* Product of two rational numbers */
741 value_multiply(res
->d
,e1
->d
,res
->d
);
742 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
743 Gcd(res
->x
.n
, res
->d
,&g
);
744 if (value_notone_p(g
)) {
745 value_division(res
->d
,res
->d
,g
);
746 value_division(res
->x
.n
,res
->x
.n
,g
);
752 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
753 /* Product of an expression (polynomial or peririodic) and a rational number */
759 /* Product of a rationel number and an expression (polynomial or peririodic) */
761 i
= res
->x
.p
->type
== modulo
? 1 : 0;
762 for (; i
<res
->x
.p
->size
; i
++)
763 emul(e1
, &res
->x
.p
->arr
[i
]);
773 void evalue_copy(evalue
*dst
, evalue
*src
)
775 value_assign(dst
->d
, src
->d
);
776 if(value_notzero_p(src
->d
)) {
777 value_init(dst
->x
.n
);
778 value_assign(dst
->x
.n
, src
->x
.n
);
780 dst
->x
.p
= ecopy(src
->x
.p
);
783 enode
*new_enode(enode_type type
,int size
,int pos
) {
789 fprintf(stderr
, "Allocating enode of size 0 !\n" );
792 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
796 for(i
=0; i
<size
; i
++) {
797 value_init(res
->arr
[i
].d
);
798 value_set_si(res
->arr
[i
].d
,0);
804 enode
*ecopy(enode
*e
) {
809 res
= new_enode(e
->type
,e
->size
,e
->pos
);
810 for(i
=0;i
<e
->size
;++i
) {
811 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
812 if(value_zero_p(res
->arr
[i
].d
))
813 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
814 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
815 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
817 value_init(res
->arr
[i
].x
.n
);
818 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
824 int eequal(evalue
*e1
,evalue
*e2
) {
829 if (value_ne(e1
->d
,e2
->d
))
833 if (value_notzero_p(e1
->d
)) {
834 if (value_ne(e1
->x
.n
,e2
->x
.n
))
837 /* e1->d == e2->d != 0 AND e1->n == e2->n */
841 /* e1->d == e2->d == 0 */
844 if (p1
->type
!= p2
->type
) return 0;
845 if (p1
->size
!= p2
->size
) return 0;
846 if (p1
->pos
!= p2
->pos
) return 0;
847 for (i
=0; i
<p1
->size
; i
++)
848 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
853 void free_evalue_refs(evalue
*e
) {
858 if (EVALUE_IS_DOMAIN(*e
)) {
859 Domain_Free(EVALUE_DOMAIN(*e
));
861 } else if (value_pos_p(e
->d
)) {
863 /* 'e' stores a constant */
868 assert(value_zero_p(e
->d
));
871 if (!p
) return; /* null pointer */
872 for (i
=0; i
<p
->size
; i
++) {
873 free_evalue_refs(&(p
->arr
[i
]));
877 } /* free_evalue_refs */
879 /********************************************************/
880 /* function in domain */
881 /* check if the parameters in list_args */
882 /* verifies the constraints of Domain P */
883 /********************************************************/
884 int in_domain(Polyhedron
*P
, Value
*list_args
) {
887 Value v
; /* value of the constraint of a row when
888 parameters are instanciated*/
894 /*P->Constraint constraint matrice of polyhedron P*/
895 for(row
=0;row
<P
->NbConstraints
;row
++) {
896 value_assign(v
,P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
897 for(col
=1;col
<P
->Dimension
+1;col
++) {
898 value_multiply(tmp
,P
->Constraint
[row
][col
],list_args
[col
-1]);
899 value_addto(v
,v
,tmp
);
901 if (value_notzero_p(P
->Constraint
[row
][0])) {
903 /*if v is not >=0 then this constraint is not respected */
904 if (value_neg_p(v
)) {
908 return P
->next
? in_domain(P
->next
, list_args
) : 0;
913 /*if v is not = 0 then this constraint is not respected */
914 if (value_notzero_p(v
))
919 /*if not return before this point => all
920 the constraints are respected */
926 /****************************************************/
927 /* function compute enode */
928 /* compute the value of enode p with parameters */
929 /* list "list_args */
930 /* compute the polynomial or the periodic */
931 /****************************************************/
933 static double compute_enode(enode
*p
, Value
*list_args
) {
945 if (p
->type
== polynomial
) {
947 value_assign(param
,list_args
[p
->pos
-1]);
949 /* Compute the polynomial using Horner's rule */
950 for (i
=p
->size
-1;i
>0;i
--) {
951 res
+=compute_evalue(&p
->arr
[i
],list_args
);
952 res
*=VALUE_TO_DOUBLE(param
);
954 res
+=compute_evalue(&p
->arr
[0],list_args
);
956 else if (p
->type
== modulo
) {
957 double d
= compute_evalue(&p
->arr
[0], list_args
);
962 value_set_double(param
, d
);
963 value_set_si(m
, p
->pos
);
964 mpz_fdiv_r(param
, param
, m
);
966 /* Compute the polynomial using Horner's rule */
967 for (i
=p
->size
-1;i
>1;i
--) {
968 res
+=compute_evalue(&p
->arr
[i
],list_args
);
969 res
*=VALUE_TO_DOUBLE(param
);
971 res
+=compute_evalue(&p
->arr
[1],list_args
);
973 else if (p
->type
== periodic
) {
974 value_assign(param
,list_args
[p
->pos
-1]);
976 /* Choose the right element of the periodic */
977 value_absolute(m
,param
);
978 value_set_si(param
,p
->size
);
979 value_modulus(m
,m
,param
);
980 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
982 else if (p
->type
== relation
) {
983 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 0.5)
984 res
= compute_evalue(&p
->arr
[1], list_args
);
989 } /* compute_enode */
991 /*************************************************/
992 /* return the value of Ehrhart Polynomial */
993 /* It returns a double, because since it is */
994 /* a recursive function, some intermediate value */
995 /* might not be integral */
996 /*************************************************/
998 double compute_evalue(evalue
*e
,Value
*list_args
) {
1002 if (value_notzero_p(e
->d
)) {
1003 if (value_notone_p(e
->d
))
1004 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
1006 res
= VALUE_TO_DOUBLE(e
->x
.n
);
1009 res
= compute_enode(e
->x
.p
,list_args
);
1011 } /* compute_evalue */
1014 /****************************************************/
1015 /* function compute_poly : */
1016 /* Check for the good validity domain */
1017 /* return the number of point in the Polyhedron */
1018 /* in allocated memory */
1019 /* Using the Ehrhart pseudo-polynomial */
1020 /****************************************************/
1021 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
1024 /* double d; int i; */
1026 tmp
= (Value
*) malloc (sizeof(Value
));
1027 assert(tmp
!= NULL
);
1029 value_set_si(*tmp
,0);
1032 return(tmp
); /* no ehrhart polynomial */
1033 if(en
->ValidityDomain
) {
1034 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
1035 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
1040 return(tmp
); /* no Validity Domain */
1042 if(in_domain(en
->ValidityDomain
,list_args
)) {
1044 #ifdef EVAL_EHRHART_DEBUG
1045 Print_Domain(stdout
,en
->ValidityDomain
);
1046 print_evalue(stdout
,&en
->EP
);
1049 /* d = compute_evalue(&en->EP,list_args);
1051 printf("(double)%lf = %d\n", d, i ); */
1052 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
1058 value_set_si(*tmp
,0);
1059 return(tmp
); /* no compatible domain with the arguments */
1060 } /* compute_poly */