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 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
333 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
334 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
335 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
336 value_clear(res
->x
.p
->arr
[2*i
].d
);
340 res
->x
.p
= new_enode(partition
, 2*n
, -1);
341 for (j
= 0; j
< n
; ++j
) {
342 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
343 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
344 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
350 void eadd(evalue
*e1
,evalue
*res
) {
353 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
354 /* Add two rational numbers */
360 value_multiply(m1
,e1
->x
.n
,res
->d
);
361 value_multiply(m2
,res
->x
.n
,e1
->d
);
362 value_addto(res
->x
.n
,m1
,m2
);
363 value_multiply(res
->d
,e1
->d
,res
->d
);
364 Gcd(res
->x
.n
,res
->d
,&g
);
365 if (value_notone_p(g
)) {
366 value_division(res
->d
,res
->d
,g
);
367 value_division(res
->x
.n
,res
->x
.n
,g
);
369 value_clear(g
); value_clear(m1
); value_clear(m2
);
372 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
373 switch (res
->x
.p
->type
) {
375 /* Add the constant to the constant term of a polynomial*/
376 eadd(e1
, &res
->x
.p
->arr
[0]);
379 /* Add the constant to all elements of a periodic number */
380 for (i
=0; i
<res
->x
.p
->size
; i
++) {
381 eadd(e1
, &res
->x
.p
->arr
[i
]);
385 fprintf(stderr
, "eadd: cannot add const with vector\n");
388 eadd(e1
, &res
->x
.p
->arr
[1]);
391 assert(EVALUE_IS_ZERO(*e1
));
392 break; /* Do nothing */
397 /* add polynomial or periodic to constant
398 * you have to exchange e1 and res, before doing addition */
400 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
404 else { // ((e1->d==0) && (res->d==0))
405 assert(!((e1
->x
.p
->type
== partition
) ^
406 (res
->x
.p
->type
== partition
)));
407 if (e1
->x
.p
->type
== partition
) {
408 eadd_partitions(e1
, res
);
411 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
412 /* adding to evalues of different type. two cases are possible
413 * res is periodic and e1 is polynomial, you have to exchange
414 * e1 and res then to add e1 to the constant term of res */
415 if (e1
->x
.p
->type
== polynomial
) {
416 eadd_rev_cst(e1
, res
);
418 else if (res
->x
.p
->type
== polynomial
) {
419 /* res is polynomial and e1 is periodic,
420 add e1 to the constant term of res */
422 eadd(e1
,&res
->x
.p
->arr
[0]);
428 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
429 (res
->x
.p
->type
== modulo
&&
430 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
431 /* adding evalues of different position (i.e function of different unknowns
432 * to case are possible */
434 switch (res
->x
.p
->type
) {
436 if(mod_term_smaller(res
, e1
))
437 eadd(e1
,&res
->x
.p
->arr
[1]);
439 eadd_rev_cst(e1
, res
);
441 case polynomial
: // res and e1 are polynomials
442 // add e1 to the constant term of res
444 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
445 eadd(e1
,&res
->x
.p
->arr
[0]);
447 eadd_rev_cst(e1
, res
);
448 // value_clear(g); value_clear(m1); value_clear(m2);
450 case periodic
: // res and e1 are pointers to periodic numbers
451 //add e1 to all elements of res
453 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
454 for (i
=0;i
<res
->x
.p
->size
;i
++) {
455 eadd(e1
,&res
->x
.p
->arr
[i
]);
464 //same type , same pos and same size
465 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
466 // add any element in e1 to the corresponding element in res
467 if (res
->x
.p
->type
== modulo
)
468 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
469 i
= res
->x
.p
->type
== modulo
? 1 : 0;
470 for (; i
<res
->x
.p
->size
; i
++) {
471 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
476 /* Sizes are different */
477 switch(res
->x
.p
->type
) {
480 /* VIN100: if e1-size > res-size you have to copy e1 in a */
481 /* new enode and add res to that new node. If you do not do */
482 /* that, you lose the the upper weight part of e1 ! */
484 if(e1
->x
.p
->size
> res
->x
.p
->size
)
488 if (res
->x
.p
->type
== modulo
)
489 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
490 i
= res
->x
.p
->type
== modulo
? 1 : 0;
491 for (; i
<e1
->x
.p
->size
; i
++) {
492 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
499 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
502 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
503 of the sizes of e1 and res, then to copy res periodicaly in ne, after
504 to add periodicaly elements of e1 to elements of ne, and finaly to
509 value_init(ex
); value_init(ey
);value_init(ep
);
512 value_set_si(ex
,e1
->x
.p
->size
);
513 value_set_si(ey
,res
->x
.p
->size
);
514 value_assign (ep
,*Lcm(ex
,ey
));
515 p
=(int)mpz_get_si(ep
);
516 ne
= (evalue
*) malloc (sizeof(evalue
));
518 value_set_si( ne
->d
,0);
520 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
522 value_assign(ne
->x
.p
->arr
[i
].d
, res
->x
.p
->arr
[i
%y
].d
);
523 if (value_notzero_p(ne
->x
.p
->arr
[i
].d
)) {
524 value_init(ne
->x
.p
->arr
[i
].x
.n
);
525 value_assign(ne
->x
.p
->arr
[i
].x
.n
, res
->x
.p
->arr
[i
%y
].x
.n
);
528 ne
->x
.p
->arr
[i
].x
.p
=ecopy(res
->x
.p
->arr
[i
%y
].x
.p
);
532 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
535 value_assign(res
->d
, ne
->d
);
541 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
548 static void emul_rev (evalue
*e1
, evalue
*res
)
552 evalue_copy(&ev
, e1
);
554 free_evalue_refs(res
);
558 static void emul_poly (evalue
*e1
, evalue
*res
)
560 int i
, j
, o
= res
->x
.p
->type
== modulo
;
562 int size
=(e1
->x
.p
->size
+ res
->x
.p
->size
- o
- 1);
564 value_set_si(tmp
.d
,0);
565 tmp
.x
.p
=new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
567 evalue_copy(&tmp
.x
.p
->arr
[0], &e1
->x
.p
->arr
[0]);
568 for (i
=o
; i
< e1
->x
.p
->size
; i
++) {
569 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
570 emul(&res
->x
.p
->arr
[o
], &tmp
.x
.p
->arr
[i
]);
573 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
574 for (i
=o
+1; i
<res
->x
.p
->size
; i
++)
575 for (j
=o
; j
<e1
->x
.p
->size
; j
++) {
578 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
579 emul(&res
->x
.p
->arr
[i
], &ev
);
580 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-o
]);
581 free_evalue_refs(&ev
);
583 free_evalue_refs(res
);
587 void emul_partitions (evalue
*e1
,evalue
*res
)
592 s
= (struct section
*)
593 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
594 sizeof(struct section
));
598 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
599 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
600 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
601 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
607 /* This code is only needed because the partitions
608 are not true partitions.
610 for (k
= 0; k
< n
; ++k
) {
611 if (DomainIncludes(s
[k
].D
, d
))
613 if (DomainIncludes(d
, s
[k
].D
)) {
615 free_evalue_refs(&s
[k
].E
);
626 value_init(s
[n
].E
.d
);
627 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
628 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
632 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
633 value_clear(res
->x
.p
->arr
[2*i
].d
);
634 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
638 res
->x
.p
= new_enode(partition
, 2*n
, -1);
639 for (j
= 0; j
< n
; ++j
) {
640 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
641 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
642 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
648 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
649 * do a copy of "res" befor calling this function if you nead it after. The vector type of
650 * evalues is not treated here */
652 void emul (evalue
*e1
, evalue
*res
){
655 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
656 fprintf(stderr
, "emul: do not proced on evector type !\n");
660 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
661 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
662 emul_partitions(e1
, res
);
665 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
666 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
667 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
669 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
)
670 emul(e1
, &res
->x
.p
->arr
[1]);
672 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
673 switch(e1
->x
.p
->type
) {
675 switch(res
->x
.p
->type
) {
677 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
678 /* Product of two polynomials of the same variable */
683 /* Product of two polynomials of different variables */
685 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
686 for( i
=0; i
<res
->x
.p
->size
; i
++)
687 emul(e1
, &res
->x
.p
->arr
[i
]);
695 /* Product of a polynomial and a periodic or modulo */
700 switch(res
->x
.p
->type
) {
702 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
703 /* Product of two periodics of the same parameter and period */
705 for(i
=0; i
<res
->x
.p
->size
;i
++)
706 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
711 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
712 /* Product of two periodics of the same parameter and different periods */
716 value_init(x
); value_init(y
);value_init(z
);
719 value_set_si(x
,e1
->x
.p
->size
);
720 value_set_si(y
,res
->x
.p
->size
);
721 value_assign (z
,*Lcm(x
,y
));
722 lcm
=(int)mpz_get_si(z
);
723 newp
= (evalue
*) malloc (sizeof(evalue
));
725 value_set_si( newp
->d
,0);
726 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
728 value_assign(newp
->x
.p
->arr
[i
].d
, res
->x
.p
->arr
[i
%iy
].d
);
729 if (value_notzero_p(newp
->x
.p
->arr
[i
].d
)) {
730 value_assign(newp
->x
.p
->arr
[i
].x
.n
, res
->x
.p
->arr
[i
%iy
].x
.n
);
733 newp
->x
.p
->arr
[i
].x
.p
=ecopy(res
->x
.p
->arr
[i
%iy
].x
.p
);
738 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
740 value_assign(res
->d
,newp
->d
);
743 value_clear(x
); value_clear(y
);value_clear(z
);
747 /* Product of two periodics of different parameters */
749 for(i
=0; i
<res
->x
.p
->size
; i
++)
750 emul(e1
, &(res
->x
.p
->arr
[i
]));
756 /* Product of a periodic and a polynomial */
758 for(i
=0; i
<res
->x
.p
->size
; i
++)
759 emul(e1
, &(res
->x
.p
->arr
[i
]));
765 switch(res
->x
.p
->type
) {
767 for(i
=0; i
<res
->x
.p
->size
; i
++)
768 emul(e1
, &(res
->x
.p
->arr
[i
]));
773 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
774 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
775 if (e1
->x
.p
->pos
!= 2)
778 /* x mod 2 == (x mod 2)^2 */
779 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1) (x mod 2) */
780 assert(e1
->x
.p
->size
== 3);
781 assert(res
->x
.p
->size
== 3);
784 evalue_copy(&tmp
, &res
->x
.p
->arr
[1]);
785 eadd(&res
->x
.p
->arr
[2], &tmp
);
786 emul(&e1
->x
.p
->arr
[2], &tmp
);
787 emul(&e1
->x
.p
->arr
[1], res
);
788 eadd(&tmp
, &res
->x
.p
->arr
[2]);
789 free_evalue_refs(&tmp
);
792 if(mod_term_smaller(res
, e1
))
793 for(i
=1; i
<res
->x
.p
->size
; i
++)
794 emul(e1
, &(res
->x
.p
->arr
[i
]));
809 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
810 /* Product of two rational numbers */
814 value_multiply(res
->d
,e1
->d
,res
->d
);
815 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
816 Gcd(res
->x
.n
, res
->d
,&g
);
817 if (value_notone_p(g
)) {
818 value_division(res
->d
,res
->d
,g
);
819 value_division(res
->x
.n
,res
->x
.n
,g
);
825 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
826 /* Product of an expression (polynomial or peririodic) and a rational number */
832 /* Product of a rationel number and an expression (polynomial or peririodic) */
834 i
= res
->x
.p
->type
== modulo
? 1 : 0;
835 for (; i
<res
->x
.p
->size
; i
++)
836 emul(e1
, &res
->x
.p
->arr
[i
]);
846 void evalue_copy(evalue
*dst
, evalue
*src
)
848 value_assign(dst
->d
, src
->d
);
849 if(value_notzero_p(src
->d
)) {
850 value_init(dst
->x
.n
);
851 value_assign(dst
->x
.n
, src
->x
.n
);
853 dst
->x
.p
= ecopy(src
->x
.p
);
856 enode
*new_enode(enode_type type
,int size
,int pos
) {
862 fprintf(stderr
, "Allocating enode of size 0 !\n" );
865 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
869 for(i
=0; i
<size
; i
++) {
870 value_init(res
->arr
[i
].d
);
871 value_set_si(res
->arr
[i
].d
,0);
877 enode
*ecopy(enode
*e
) {
882 res
= new_enode(e
->type
,e
->size
,e
->pos
);
883 for(i
=0;i
<e
->size
;++i
) {
884 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
885 if(value_zero_p(res
->arr
[i
].d
))
886 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
887 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
888 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
890 value_init(res
->arr
[i
].x
.n
);
891 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
897 int eequal(evalue
*e1
,evalue
*e2
) {
902 if (value_ne(e1
->d
,e2
->d
))
906 if (value_notzero_p(e1
->d
)) {
907 if (value_ne(e1
->x
.n
,e2
->x
.n
))
910 /* e1->d == e2->d != 0 AND e1->n == e2->n */
914 /* e1->d == e2->d == 0 */
917 if (p1
->type
!= p2
->type
) return 0;
918 if (p1
->size
!= p2
->size
) return 0;
919 if (p1
->pos
!= p2
->pos
) return 0;
920 for (i
=0; i
<p1
->size
; i
++)
921 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
926 void free_evalue_refs(evalue
*e
) {
931 if (EVALUE_IS_DOMAIN(*e
)) {
932 Domain_Free(EVALUE_DOMAIN(*e
));
935 } else if (value_pos_p(e
->d
)) {
937 /* 'e' stores a constant */
942 assert(value_zero_p(e
->d
));
945 if (!p
) return; /* null pointer */
946 for (i
=0; i
<p
->size
; i
++) {
947 free_evalue_refs(&(p
->arr
[i
]));
951 } /* free_evalue_refs */
953 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
954 Vector
* val
, evalue
*res
)
956 unsigned nparam
= periods
->Size
;
959 double d
= compute_evalue(e
, val
->p
);
964 value_set_si(res
->d
, 1);
965 value_init(res
->x
.n
);
966 value_set_double(res
->x
.n
, d
);
967 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
970 if (value_one_p(periods
->p
[p
]))
971 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
976 value_assign(tmp
, periods
->p
[p
]);
977 value_set_si(res
->d
, 0);
978 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
980 value_decrement(tmp
, tmp
);
981 value_assign(val
->p
[p
], tmp
);
982 mod2table_r(e
, periods
, m
, p
+1, val
,
983 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
984 } while (value_pos_p(tmp
));
990 void evalue_mod2table(evalue
*e
, int nparam
)
995 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
998 for (i
=0; i
<p
->size
; i
++) {
999 evalue_mod2table(&(p
->arr
[i
]), nparam
);
1001 if (p
->type
== modulo
) {
1002 Vector
*periods
= Vector_Alloc(nparam
);
1003 Vector
*val
= Vector_Alloc(nparam
);
1009 value_set_si(tmp
, p
->pos
);
1010 Vector_Set(periods
->p
, 1, nparam
);
1011 Vector_Set(val
->p
, 0, nparam
);
1012 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
1015 assert(p
->type
== polynomial
);
1016 assert(p
->size
== 2);
1017 assert(value_one_p(p
->arr
[1].d
));
1018 Gcd(tmp
, p
->arr
[1].x
.n
, &periods
->p
[p
->pos
-1]);
1019 value_division(periods
->p
[p
->pos
-1], tmp
, periods
->p
[p
->pos
-1]);
1021 value_set_si(tmp
, p
->pos
);
1023 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
1026 evalue_set_si(&res
, 0, 1);
1027 /* Compute the polynomial using Horner's rule */
1028 for (i
=p
->size
-1;i
>1;i
--) {
1029 eadd(&p
->arr
[i
], &res
);
1032 eadd(&p
->arr
[1], &res
);
1034 free_evalue_refs(e
);
1035 free_evalue_refs(&EP
);
1040 Vector_Free(periods
);
1042 } /* evalue_mod2table */
1044 /********************************************************/
1045 /* function in domain */
1046 /* check if the parameters in list_args */
1047 /* verifies the constraints of Domain P */
1048 /********************************************************/
1049 int in_domain(Polyhedron
*P
, Value
*list_args
) {
1052 Value v
; /* value of the constraint of a row when
1053 parameters are instanciated*/
1059 /*P->Constraint constraint matrice of polyhedron P*/
1060 for(row
=0;row
<P
->NbConstraints
;row
++) {
1061 value_assign(v
,P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
1062 for(col
=1;col
<P
->Dimension
+1;col
++) {
1063 value_multiply(tmp
,P
->Constraint
[row
][col
],list_args
[col
-1]);
1064 value_addto(v
,v
,tmp
);
1066 if (value_notzero_p(P
->Constraint
[row
][0])) {
1068 /*if v is not >=0 then this constraint is not respected */
1069 if (value_neg_p(v
)) {
1073 return P
->next
? in_domain(P
->next
, list_args
) : 0;
1078 /*if v is not = 0 then this constraint is not respected */
1079 if (value_notzero_p(v
))
1084 /*if not return before this point => all
1085 the constraints are respected */
1091 /****************************************************/
1092 /* function compute enode */
1093 /* compute the value of enode p with parameters */
1094 /* list "list_args */
1095 /* compute the polynomial or the periodic */
1096 /****************************************************/
1098 static double compute_enode(enode
*p
, Value
*list_args
) {
1110 if (p
->type
== polynomial
) {
1112 value_assign(param
,list_args
[p
->pos
-1]);
1114 /* Compute the polynomial using Horner's rule */
1115 for (i
=p
->size
-1;i
>0;i
--) {
1116 res
+=compute_evalue(&p
->arr
[i
],list_args
);
1117 res
*=VALUE_TO_DOUBLE(param
);
1119 res
+=compute_evalue(&p
->arr
[0],list_args
);
1121 else if (p
->type
== modulo
) {
1122 double d
= compute_evalue(&p
->arr
[0], list_args
);
1127 value_set_double(param
, d
);
1128 value_set_si(m
, p
->pos
);
1129 mpz_fdiv_r(param
, param
, m
);
1131 /* Compute the polynomial using Horner's rule */
1132 for (i
=p
->size
-1;i
>1;i
--) {
1133 res
+=compute_evalue(&p
->arr
[i
],list_args
);
1134 res
*=VALUE_TO_DOUBLE(param
);
1136 res
+=compute_evalue(&p
->arr
[1],list_args
);
1138 else if (p
->type
== periodic
) {
1139 value_assign(param
,list_args
[p
->pos
-1]);
1141 /* Choose the right element of the periodic */
1142 value_absolute(m
,param
);
1143 value_set_si(param
,p
->size
);
1144 value_modulus(m
,m
,param
);
1145 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
1147 else if (p
->type
== relation
) {
1148 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 0.5)
1149 res
= compute_evalue(&p
->arr
[1], list_args
);
1154 } /* compute_enode */
1156 /*************************************************/
1157 /* return the value of Ehrhart Polynomial */
1158 /* It returns a double, because since it is */
1159 /* a recursive function, some intermediate value */
1160 /* might not be integral */
1161 /*************************************************/
1163 double compute_evalue(evalue
*e
,Value
*list_args
) {
1167 if (value_notzero_p(e
->d
)) {
1168 if (value_notone_p(e
->d
))
1169 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
1171 res
= VALUE_TO_DOUBLE(e
->x
.n
);
1174 res
= compute_enode(e
->x
.p
,list_args
);
1176 } /* compute_evalue */
1179 /****************************************************/
1180 /* function compute_poly : */
1181 /* Check for the good validity domain */
1182 /* return the number of point in the Polyhedron */
1183 /* in allocated memory */
1184 /* Using the Ehrhart pseudo-polynomial */
1185 /****************************************************/
1186 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
1189 /* double d; int i; */
1191 tmp
= (Value
*) malloc (sizeof(Value
));
1192 assert(tmp
!= NULL
);
1194 value_set_si(*tmp
,0);
1197 return(tmp
); /* no ehrhart polynomial */
1198 if(en
->ValidityDomain
) {
1199 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
1200 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
1205 return(tmp
); /* no Validity Domain */
1207 if(in_domain(en
->ValidityDomain
,list_args
)) {
1209 #ifdef EVAL_EHRHART_DEBUG
1210 Print_Domain(stdout
,en
->ValidityDomain
);
1211 print_evalue(stdout
,&en
->EP
);
1214 /* d = compute_evalue(&en->EP,list_args);
1216 printf("(double)%lf = %d\n", d, i ); */
1217 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
1223 value_set_si(*tmp
,0);
1224 return(tmp
); /* no compatible domain with the arguments */
1225 } /* compute_poly */