5 #include "ev_operations.h"
9 #define ALLOC(p) p = (typeof(p))malloc(sizeof(*p))
10 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
12 void evalue_set_si(evalue
*ev
, int n
, int d
) {
13 value_set_si(ev
->d
, d
);
15 value_set_si(ev
->x
.n
, n
);
18 void evalue_set(evalue
*ev
, Value n
, Value d
) {
19 value_assign(ev
->d
, d
);
21 value_assign(ev
->x
.n
, n
);
24 void aep_evalue(evalue
*e
, int *ref
) {
29 if (value_notzero_p(e
->d
))
30 return; /* a rational number, its already reduced */
32 return; /* hum... an overflow probably occured */
34 /* First check the components of p */
35 for (i
=0;i
<p
->size
;i
++)
36 aep_evalue(&p
->arr
[i
],ref
);
43 p
->pos
= ref
[p
->pos
-1]+1;
49 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
55 if (value_notzero_p(e
->d
))
56 return; /* a rational number, its already reduced */
58 return; /* hum... an overflow probably occured */
61 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
62 for(i
=0;i
<CT
->NbRows
-1;i
++)
63 for(j
=0;j
<CT
->NbColumns
;j
++)
64 if(value_notzero_p(CT
->p
[i
][j
])) {
69 /* Transform the references in e, using ref */
73 } /* addeliminatedparams_evalue */
75 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
83 assert(value_notzero_p(e1
->d
));
84 assert(value_notzero_p(e2
->d
));
85 value_multiply(m
, e1
->x
.n
, e2
->d
);
86 value_multiply(m2
, e2
->x
.n
, e1
->d
);
89 else if (value_gt(m
, m2
))
99 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
101 if (value_notzero_p(e1
->d
)) {
102 if (value_zero_p(e2
->d
))
104 int r
= mod_rational_smaller(e1
, e2
);
105 return r
== -1 ? 0 : r
;
107 if (value_notzero_p(e2
->d
))
109 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
111 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
114 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
116 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
121 static int mod_term_smaller(evalue
*e1
, evalue
*e2
)
123 assert(value_zero_p(e1
->d
));
124 assert(value_zero_p(e2
->d
));
125 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
126 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
127 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
130 /* Negative pos means inequality */
131 /* s is negative of substitution if m is not zero */
140 struct fixed_param
*fixed
;
145 static int relations_depth(evalue
*e
)
150 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
151 e
= &e
->x
.p
->arr
[1], ++d
);
155 static void Lcm3(Value i
, Value j
, Value
*res
)
161 value_multiply(*res
,i
,j
);
162 value_division(*res
, *res
, aux
);
166 static void poly_denom(evalue
*p
, Value
*d
)
170 while (value_zero_p(p
->d
)) {
171 assert(p
->x
.p
->type
== polynomial
);
172 assert(p
->x
.p
->size
== 2);
173 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
174 Lcm3(*d
, p
->x
.p
->arr
[1].d
, d
);
180 #define EVALUE_IS_ONE(ev) (value_one_p((ev).d) && value_one_p((ev).x.n))
182 static void realloc_substitution(struct subst
*s
, int d
)
184 struct fixed_param
*n
;
187 for (i
= 0; i
< s
->n
; ++i
)
194 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
200 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
203 /* May have been reduced already */
204 if (value_notzero_p(m
->d
))
207 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
208 assert(m
->x
.p
->size
== 3);
210 /* fractional was inverted during reduction
211 * invert it back and move constant in
213 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
214 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
215 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
216 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
217 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
218 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
219 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
220 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
221 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
222 value_set_si(m
->x
.p
->arr
[1].d
, 1);
225 /* Oops. Nested identical relations. */
226 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
229 if (s
->n
>= s
->max
) {
230 int d
= relations_depth(r
);
231 realloc_substitution(s
, d
);
235 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
236 assert(p
->x
.p
->size
== 2);
239 assert(value_pos_p(f
->x
.n
));
241 value_init(s
->fixed
[s
->n
].m
);
242 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
243 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
244 value_init(s
->fixed
[s
->n
].d
);
245 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
246 value_init(s
->fixed
[s
->n
].s
.d
);
247 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
253 static int type_offset(enode
*p
)
255 return p
->type
== fractional
? 1 :
256 p
->type
== flooring
? 1 : 0;
259 static void reorder_terms(enode
*p
, evalue
*v
)
262 int offset
= type_offset(p
);
264 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
266 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
267 free_evalue_refs(&(p
->arr
[i
]));
273 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
279 if (value_notzero_p(e
->d
)) {
281 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
282 return; /* a rational number, its already reduced */
286 return; /* hum... an overflow probably occured */
288 /* First reduce the components of p */
289 add
= p
->type
== relation
;
290 for (i
=0; i
<p
->size
; i
++) {
292 add
= add_modulo_substitution(s
, e
);
294 if (i
== 0 && p
->type
==fractional
)
295 _reduce_evalue(&p
->arr
[i
], s
, 1);
297 _reduce_evalue(&p
->arr
[i
], s
, fract
);
299 if (add
&& i
== p
->size
-1) {
301 value_clear(s
->fixed
[s
->n
].m
);
302 value_clear(s
->fixed
[s
->n
].d
);
303 free_evalue_refs(&s
->fixed
[s
->n
].s
);
304 } else if (add
&& i
== 1)
305 s
->fixed
[s
->n
-1].pos
*= -1;
308 if (p
->type
==periodic
) {
310 /* Try to reduce the period */
311 for (i
=1; i
<=(p
->size
)/2; i
++) {
312 if ((p
->size
% i
)==0) {
314 /* Can we reduce the size to i ? */
316 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
317 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
320 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
324 you_lose
: /* OK, lets not do it */
329 /* Try to reduce its strength */
332 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
336 else if (p
->type
==polynomial
) {
337 for (k
= 0; s
&& k
< s
->n
; ++k
) {
338 if (s
->fixed
[k
].pos
== p
->pos
) {
339 int divide
= value_notone_p(s
->fixed
[k
].d
);
342 if (value_notzero_p(s
->fixed
[k
].m
)) {
345 assert(p
->size
== 2);
346 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
348 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
355 value_assign(d
.d
, s
->fixed
[k
].d
);
357 if (value_notzero_p(s
->fixed
[k
].m
))
358 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
360 value_set_si(d
.x
.n
, 1);
363 for (i
=p
->size
-1;i
>=1;i
--) {
364 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
366 emul(&d
, &p
->arr
[i
]);
367 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
368 free_evalue_refs(&(p
->arr
[i
]));
371 _reduce_evalue(&p
->arr
[0], s
, fract
);
374 free_evalue_refs(&d
);
380 /* Try to reduce the degree */
381 for (i
=p
->size
-1;i
>=1;i
--) {
382 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
384 /* Zero coefficient */
385 free_evalue_refs(&(p
->arr
[i
]));
390 /* Try to reduce its strength */
393 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
397 else if (p
->type
==fractional
) {
401 if (value_notzero_p(p
->arr
[0].d
)) {
403 value_assign(v
.d
, p
->arr
[0].d
);
405 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
410 evalue
*pp
= &p
->arr
[0];
411 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
412 assert(pp
->x
.p
->size
== 2);
414 /* search for exact duplicate among the modulo inequalities */
416 f
= &pp
->x
.p
->arr
[1];
417 for (k
= 0; s
&& k
< s
->n
; ++k
) {
418 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
419 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
420 value_eq(s
->fixed
[k
].m
, f
->d
) &&
421 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
428 /* replace { E/m } by { (E-1)/m } + 1/m */
433 evalue_set_si(&extra
, 1, 1);
434 value_assign(extra
.d
, g
);
435 eadd(&extra
, &v
.x
.p
->arr
[1]);
436 free_evalue_refs(&extra
);
438 /* We've been going in circles; stop now */
439 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
440 free_evalue_refs(&v
);
442 evalue_set_si(&v
, 0, 1);
447 value_set_si(v
.d
, 0);
448 v
.x
.p
= new_enode(fractional
, 3, -1);
449 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
450 value_assign(v
.x
.p
->arr
[1].d
, g
);
451 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
452 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
455 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
458 value_division(f
->d
, g
, f
->d
);
459 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
460 value_assign(f
->d
, g
);
461 value_decrement(f
->x
.n
, f
->x
.n
);
462 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
464 Gcd(f
->d
, f
->x
.n
, &g
);
465 value_division(f
->d
, f
->d
, g
);
466 value_division(f
->x
.n
, f
->x
.n
, g
);
475 /* reduction may have made this fractional arg smaller */
476 i
= reorder
? p
->size
: 1;
477 for ( ; i
< p
->size
; ++i
)
478 if (value_zero_p(p
->arr
[i
].d
) &&
479 p
->arr
[i
].x
.p
->type
== fractional
&&
480 !mod_term_smaller(e
, &p
->arr
[i
]))
484 value_set_si(v
.d
, 0);
485 v
.x
.p
= new_enode(fractional
, 3, -1);
486 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
487 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
488 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
498 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
499 pp
= &pp
->x
.p
->arr
[0]) {
500 f
= &pp
->x
.p
->arr
[1];
501 assert(value_pos_p(f
->d
));
502 mpz_mul_ui(twice
, f
->x
.n
, 2);
503 if (value_lt(twice
, f
->d
))
505 if (value_eq(twice
, f
->d
))
513 value_set_si(v
.d
, 0);
514 v
.x
.p
= new_enode(fractional
, 3, -1);
515 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
516 poly_denom(&p
->arr
[0], &twice
);
517 value_assign(v
.x
.p
->arr
[1].d
, twice
);
518 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
519 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
520 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
522 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
523 pp
= &pp
->x
.p
->arr
[0]) {
524 f
= &pp
->x
.p
->arr
[1];
525 value_oppose(f
->x
.n
, f
->x
.n
);
526 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
528 value_division(pp
->d
, twice
, pp
->d
);
529 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
530 value_assign(pp
->d
, twice
);
531 value_oppose(pp
->x
.n
, pp
->x
.n
);
532 value_decrement(pp
->x
.n
, pp
->x
.n
);
533 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
543 reorder_terms(p
, &v
);
544 _reduce_evalue(&p
->arr
[1], s
, fract
);
547 /* Try to reduce the degree */
548 for (i
=p
->size
-1;i
>=2;i
--) {
549 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
551 /* Zero coefficient */
552 free_evalue_refs(&(p
->arr
[i
]));
557 /* Try to reduce its strength */
560 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
561 free_evalue_refs(&(p
->arr
[0]));
565 else if (p
->type
== relation
) {
566 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
567 free_evalue_refs(&(p
->arr
[2]));
568 free_evalue_refs(&(p
->arr
[0]));
575 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
576 free_evalue_refs(&(p
->arr
[2]));
579 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
580 free_evalue_refs(&(p
->arr
[1]));
581 free_evalue_refs(&(p
->arr
[0]));
582 evalue_set_si(e
, 0, 1);
589 /* Relation was reduced by means of an identical
590 * inequality => remove
592 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
595 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
596 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
598 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
600 free_evalue_refs(&(p
->arr
[2]));
604 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
606 evalue_set_si(e
, 0, 1);
607 free_evalue_refs(&(p
->arr
[1]));
609 free_evalue_refs(&(p
->arr
[0]));
615 } /* reduce_evalue */
617 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
622 for (k
= 0; k
< dim
; ++k
)
623 if (value_notzero_p(row
[k
+1]))
626 Vector_Normalize_Positive(row
+1, dim
+1, k
);
627 assert(s
->n
< s
->max
);
628 value_init(s
->fixed
[s
->n
].d
);
629 value_init(s
->fixed
[s
->n
].m
);
630 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
631 s
->fixed
[s
->n
].pos
= k
+1;
632 value_set_si(s
->fixed
[s
->n
].m
, 0);
633 r
= &s
->fixed
[s
->n
].s
;
635 for (l
= k
+1; l
< dim
; ++l
)
636 if (value_notzero_p(row
[l
+1])) {
637 value_set_si(r
->d
, 0);
638 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
639 value_init(r
->x
.p
->arr
[1].x
.n
);
640 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
641 value_set_si(r
->x
.p
->arr
[1].d
, 1);
645 value_oppose(r
->x
.n
, row
[dim
+1]);
646 value_set_si(r
->d
, 1);
650 void reduce_evalue (evalue
*e
) {
651 struct subst s
= { NULL
, 0, 0 };
653 if (value_notzero_p(e
->d
))
654 return; /* a rational number, its already reduced */
656 if (e
->x
.p
->type
== partition
) {
659 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
661 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
662 /* This shouldn't really happen;
663 * Empty domains should not be added.
670 D
= DomainConvex(D
, 0);
671 if (!D
->next
&& D
->NbEq
) {
675 realloc_substitution(&s
, dim
);
677 int d
= relations_depth(&e
->x
.p
->arr
[2*i
+1]);
679 NALLOC(s
.fixed
, s
.max
);
682 for (j
= 0; j
< D
->NbEq
; ++j
)
683 add_substitution(&s
, D
->Constraint
[j
], dim
);
685 if (D
!= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))
687 _reduce_evalue(&e
->x
.p
->arr
[2*i
+1], &s
, 0);
688 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
690 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
691 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
692 value_clear(e
->x
.p
->arr
[2*i
].d
);
694 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
695 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
700 for (j
= 0; j
< s
.n
; ++j
) {
701 value_clear(s
.fixed
[j
].d
);
702 value_clear(s
.fixed
[j
].m
);
703 free_evalue_refs(&s
.fixed
[j
].s
);
707 if (e
->x
.p
->size
== 0) {
709 evalue_set_si(e
, 0, 1);
712 _reduce_evalue(e
, &s
, 0);
717 void print_evalue(FILE *DST
,evalue
*e
,char **pname
) {
719 if(value_notzero_p(e
->d
)) {
720 if(value_notone_p(e
->d
)) {
721 value_print(DST
,VALUE_FMT
,e
->x
.n
);
723 value_print(DST
,VALUE_FMT
,e
->d
);
726 value_print(DST
,VALUE_FMT
,e
->x
.n
);
730 print_enode(DST
,e
->x
.p
,pname
);
734 void print_enode(FILE *DST
,enode
*p
,char **pname
) {
739 fprintf(DST
, "NULL");
745 for (i
=0; i
<p
->size
; i
++) {
746 print_evalue(DST
, &p
->arr
[i
], pname
);
750 fprintf(DST
, " }\n");
754 for (i
=p
->size
-1; i
>=0; i
--) {
755 print_evalue(DST
, &p
->arr
[i
], pname
);
756 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
758 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
760 fprintf(DST
, " )\n");
764 for (i
=0; i
<p
->size
; i
++) {
765 print_evalue(DST
, &p
->arr
[i
], pname
);
766 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
768 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
773 for (i
=p
->size
-1; i
>=1; i
--) {
774 print_evalue(DST
, &p
->arr
[i
], pname
);
777 fprintf(DST
, p
->type
== flooring
? "[" : "{");
778 print_evalue(DST
, &p
->arr
[0], pname
);
779 fprintf(DST
, p
->type
== flooring
? "]" : "}");
781 fprintf(DST
, "^%d + ", i
-1);
786 fprintf(DST
, " )\n");
790 print_evalue(DST
, &p
->arr
[0], pname
);
791 fprintf(DST
, "= 0 ] * \n");
792 print_evalue(DST
, &p
->arr
[1], pname
);
794 fprintf(DST
, " +\n [ ");
795 print_evalue(DST
, &p
->arr
[0], pname
);
796 fprintf(DST
, "!= 0 ] * \n");
797 print_evalue(DST
, &p
->arr
[2], pname
);
801 for (i
=0; i
<p
->size
/2; i
++) {
802 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), pname
);
803 print_evalue(DST
, &p
->arr
[2*i
+1], pname
);
812 static void eadd_rev(evalue
*e1
, evalue
*res
)
816 evalue_copy(&ev
, e1
);
818 free_evalue_refs(res
);
822 static void eadd_rev_cst (evalue
*e1
, evalue
*res
)
826 evalue_copy(&ev
, e1
);
827 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
828 free_evalue_refs(res
);
832 struct section
{ Polyhedron
* D
; evalue E
; };
834 void eadd_partitions (evalue
*e1
,evalue
*res
)
839 s
= (struct section
*)
840 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
841 sizeof(struct section
));
845 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
846 assert(res
->x
.p
->size
>= 2);
847 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
848 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
850 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
852 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
861 value_init(s
[n
].E
.d
);
862 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
866 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
867 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
868 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
870 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
871 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
877 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
878 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
880 value_init(s
[n
].E
.d
);
881 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
882 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
887 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
891 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
894 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
895 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
896 value_clear(res
->x
.p
->arr
[2*i
].d
);
900 res
->x
.p
= new_enode(partition
, 2*n
, -1);
901 for (j
= 0; j
< n
; ++j
) {
902 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
903 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
904 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
905 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
911 static void explicit_complement(evalue
*res
)
913 enode
*rel
= new_enode(relation
, 3, 0);
915 value_clear(rel
->arr
[0].d
);
916 rel
->arr
[0] = res
->x
.p
->arr
[0];
917 value_clear(rel
->arr
[1].d
);
918 rel
->arr
[1] = res
->x
.p
->arr
[1];
919 value_set_si(rel
->arr
[2].d
, 1);
920 value_init(rel
->arr
[2].x
.n
);
921 value_set_si(rel
->arr
[2].x
.n
, 0);
926 void eadd(evalue
*e1
,evalue
*res
) {
929 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
930 /* Add two rational numbers */
936 value_multiply(m1
,e1
->x
.n
,res
->d
);
937 value_multiply(m2
,res
->x
.n
,e1
->d
);
938 value_addto(res
->x
.n
,m1
,m2
);
939 value_multiply(res
->d
,e1
->d
,res
->d
);
940 Gcd(res
->x
.n
,res
->d
,&g
);
941 if (value_notone_p(g
)) {
942 value_division(res
->d
,res
->d
,g
);
943 value_division(res
->x
.n
,res
->x
.n
,g
);
945 value_clear(g
); value_clear(m1
); value_clear(m2
);
948 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
949 switch (res
->x
.p
->type
) {
951 /* Add the constant to the constant term of a polynomial*/
952 eadd(e1
, &res
->x
.p
->arr
[0]);
955 /* Add the constant to all elements of a periodic number */
956 for (i
=0; i
<res
->x
.p
->size
; i
++) {
957 eadd(e1
, &res
->x
.p
->arr
[i
]);
961 fprintf(stderr
, "eadd: cannot add const with vector\n");
965 eadd(e1
, &res
->x
.p
->arr
[1]);
968 assert(EVALUE_IS_ZERO(*e1
));
969 break; /* Do nothing */
971 /* Create (zero) complement if needed */
972 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
973 explicit_complement(res
);
974 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
975 eadd(e1
, &res
->x
.p
->arr
[i
]);
981 /* add polynomial or periodic to constant
982 * you have to exchange e1 and res, before doing addition */
984 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
988 else { // ((e1->d==0) && (res->d==0))
989 assert(!((e1
->x
.p
->type
== partition
) ^
990 (res
->x
.p
->type
== partition
)));
991 if (e1
->x
.p
->type
== partition
) {
992 eadd_partitions(e1
, res
);
995 if (e1
->x
.p
->type
== relation
&&
996 (res
->x
.p
->type
!= relation
||
997 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1001 if (res
->x
.p
->type
== relation
) {
1002 if (e1
->x
.p
->type
== relation
&&
1003 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1004 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1005 explicit_complement(res
);
1006 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1007 explicit_complement(e1
);
1008 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1009 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1012 if (res
->x
.p
->size
< 3)
1013 explicit_complement(res
);
1014 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1015 eadd(e1
, &res
->x
.p
->arr
[i
]);
1018 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
1019 /* adding to evalues of different type. two cases are possible
1020 * res is periodic and e1 is polynomial, you have to exchange
1021 * e1 and res then to add e1 to the constant term of res */
1022 if (e1
->x
.p
->type
== polynomial
) {
1023 eadd_rev_cst(e1
, res
);
1025 else if (res
->x
.p
->type
== polynomial
) {
1026 /* res is polynomial and e1 is periodic,
1027 add e1 to the constant term of res */
1029 eadd(e1
,&res
->x
.p
->arr
[0]);
1035 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
1036 ((res
->x
.p
->type
== fractional
||
1037 res
->x
.p
->type
== flooring
) &&
1038 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1039 /* adding evalues of different position (i.e function of different unknowns
1040 * to case are possible */
1042 switch (res
->x
.p
->type
) {
1045 if(mod_term_smaller(res
, e1
))
1046 eadd(e1
,&res
->x
.p
->arr
[1]);
1048 eadd_rev_cst(e1
, res
);
1050 case polynomial
: // res and e1 are polynomials
1051 // add e1 to the constant term of res
1053 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1054 eadd(e1
,&res
->x
.p
->arr
[0]);
1056 eadd_rev_cst(e1
, res
);
1057 // value_clear(g); value_clear(m1); value_clear(m2);
1059 case periodic
: // res and e1 are pointers to periodic numbers
1060 //add e1 to all elements of res
1062 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1063 for (i
=0;i
<res
->x
.p
->size
;i
++) {
1064 eadd(e1
,&res
->x
.p
->arr
[i
]);
1075 //same type , same pos and same size
1076 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1077 // add any element in e1 to the corresponding element in res
1078 i
= type_offset(res
->x
.p
);
1080 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1081 for (; i
<res
->x
.p
->size
; i
++) {
1082 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1087 /* Sizes are different */
1088 switch(res
->x
.p
->type
) {
1092 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1093 /* new enode and add res to that new node. If you do not do */
1094 /* that, you lose the the upper weight part of e1 ! */
1096 if(e1
->x
.p
->size
> res
->x
.p
->size
)
1099 i
= type_offset(res
->x
.p
);
1101 assert(eequal(&e1
->x
.p
->arr
[0],
1102 &res
->x
.p
->arr
[0]));
1103 for (; i
<e1
->x
.p
->size
; i
++) {
1104 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1111 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1114 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1115 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1116 to add periodicaly elements of e1 to elements of ne, and finaly to
1121 value_init(ex
); value_init(ey
);value_init(ep
);
1124 value_set_si(ex
,e1
->x
.p
->size
);
1125 value_set_si(ey
,res
->x
.p
->size
);
1126 value_assign (ep
,*Lcm(ex
,ey
));
1127 p
=(int)mpz_get_si(ep
);
1128 ne
= (evalue
*) malloc (sizeof(evalue
));
1130 value_set_si( ne
->d
,0);
1132 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1134 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1137 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1140 value_assign(res
->d
, ne
->d
);
1146 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
1155 static void emul_rev (evalue
*e1
, evalue
*res
)
1159 evalue_copy(&ev
, e1
);
1161 free_evalue_refs(res
);
1165 static void emul_poly (evalue
*e1
, evalue
*res
)
1167 int i
, j
, o
= type_offset(res
->x
.p
);
1169 int size
=(e1
->x
.p
->size
+ res
->x
.p
->size
- o
- 1);
1171 value_set_si(tmp
.d
,0);
1172 tmp
.x
.p
=new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1174 evalue_copy(&tmp
.x
.p
->arr
[0], &e1
->x
.p
->arr
[0]);
1175 for (i
=o
; i
< e1
->x
.p
->size
; i
++) {
1176 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1177 emul(&res
->x
.p
->arr
[o
], &tmp
.x
.p
->arr
[i
]);
1180 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1181 for (i
=o
+1; i
<res
->x
.p
->size
; i
++)
1182 for (j
=o
; j
<e1
->x
.p
->size
; j
++) {
1185 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1186 emul(&res
->x
.p
->arr
[i
], &ev
);
1187 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-o
]);
1188 free_evalue_refs(&ev
);
1190 free_evalue_refs(res
);
1194 void emul_partitions (evalue
*e1
,evalue
*res
)
1199 s
= (struct section
*)
1200 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1201 sizeof(struct section
));
1205 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1206 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1207 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1208 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1214 /* This code is only needed because the partitions
1215 are not true partitions.
1217 for (k
= 0; k
< n
; ++k
) {
1218 if (DomainIncludes(s
[k
].D
, d
))
1220 if (DomainIncludes(d
, s
[k
].D
)) {
1221 Domain_Free(s
[k
].D
);
1222 free_evalue_refs(&s
[k
].E
);
1233 value_init(s
[n
].E
.d
);
1234 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1235 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1239 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1240 value_clear(res
->x
.p
->arr
[2*i
].d
);
1241 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1245 res
->x
.p
= new_enode(partition
, 2*n
, -1);
1246 for (j
= 0; j
< n
; ++j
) {
1247 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1248 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1249 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1250 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1256 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1258 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1259 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1260 * evalues is not treated here */
1262 void emul (evalue
*e1
, evalue
*res
){
1265 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1266 fprintf(stderr
, "emul: do not proced on evector type !\n");
1270 if (EVALUE_IS_ZERO(*res
))
1273 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1274 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1275 emul_partitions(e1
, res
);
1278 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1279 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1280 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1282 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1283 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1284 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1285 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1286 explicit_complement(res
);
1287 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1288 explicit_complement(e1
);
1289 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1290 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1293 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1294 emul(e1
, &res
->x
.p
->arr
[i
]);
1296 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1297 switch(e1
->x
.p
->type
) {
1299 switch(res
->x
.p
->type
) {
1301 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1302 /* Product of two polynomials of the same variable */
1307 /* Product of two polynomials of different variables */
1309 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1310 for( i
=0; i
<res
->x
.p
->size
; i
++)
1311 emul(e1
, &res
->x
.p
->arr
[i
]);
1320 /* Product of a polynomial and a periodic or fractional */
1327 switch(res
->x
.p
->type
) {
1329 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1330 /* Product of two periodics of the same parameter and period */
1332 for(i
=0; i
<res
->x
.p
->size
;i
++)
1333 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1338 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1339 /* Product of two periodics of the same parameter and different periods */
1343 value_init(x
); value_init(y
);value_init(z
);
1346 value_set_si(x
,e1
->x
.p
->size
);
1347 value_set_si(y
,res
->x
.p
->size
);
1348 value_assign (z
,*Lcm(x
,y
));
1349 lcm
=(int)mpz_get_si(z
);
1350 newp
= (evalue
*) malloc (sizeof(evalue
));
1351 value_init(newp
->d
);
1352 value_set_si( newp
->d
,0);
1353 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1354 for(i
=0;i
<lcm
;i
++) {
1355 evalue_copy(&newp
->x
.p
->arr
[i
],
1356 &res
->x
.p
->arr
[i
%iy
]);
1359 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1361 value_assign(res
->d
,newp
->d
);
1364 value_clear(x
); value_clear(y
);value_clear(z
);
1368 /* Product of two periodics of different parameters */
1370 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1371 for(i
=0; i
<res
->x
.p
->size
; i
++)
1372 emul(e1
, &(res
->x
.p
->arr
[i
]));
1380 /* Product of a periodic and a polynomial */
1382 for(i
=0; i
<res
->x
.p
->size
; i
++)
1383 emul(e1
, &(res
->x
.p
->arr
[i
]));
1390 switch(res
->x
.p
->type
) {
1392 for(i
=0; i
<res
->x
.p
->size
; i
++)
1393 emul(e1
, &(res
->x
.p
->arr
[i
]));
1400 assert(e1
->x
.p
->type
== res
->x
.p
->type
);
1401 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1402 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1405 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1406 if (e1
->x
.p
->type
!= fractional
|| !value_two_p(d
.d
))
1410 value_set_si(d
.x
.n
, 1);
1412 /* { x }^2 == { x }/2 */
1413 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1414 assert(e1
->x
.p
->size
== 3);
1415 assert(res
->x
.p
->size
== 3);
1417 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1419 eadd(&res
->x
.p
->arr
[1], &tmp
);
1420 emul(&e1
->x
.p
->arr
[2], &tmp
);
1421 emul(&e1
->x
.p
->arr
[1], res
);
1422 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1423 free_evalue_refs(&tmp
);
1428 if(mod_term_smaller(res
, e1
))
1429 for(i
=1; i
<res
->x
.p
->size
; i
++)
1430 emul(e1
, &(res
->x
.p
->arr
[i
]));
1445 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1446 /* Product of two rational numbers */
1450 value_multiply(res
->d
,e1
->d
,res
->d
);
1451 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1452 Gcd(res
->x
.n
, res
->d
,&g
);
1453 if (value_notone_p(g
)) {
1454 value_division(res
->d
,res
->d
,g
);
1455 value_division(res
->x
.n
,res
->x
.n
,g
);
1461 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1462 /* Product of an expression (polynomial or peririodic) and a rational number */
1468 /* Product of a rationel number and an expression (polynomial or peririodic) */
1470 i
= type_offset(res
->x
.p
);
1471 for (; i
<res
->x
.p
->size
; i
++)
1472 emul(e1
, &res
->x
.p
->arr
[i
]);
1482 /* Frees mask content ! */
1483 void emask(evalue
*mask
, evalue
*res
) {
1489 if (EVALUE_IS_ZERO(*res
)) {
1490 free_evalue_refs(mask
);
1494 assert(value_zero_p(mask
->d
));
1495 assert(mask
->x
.p
->type
== partition
);
1496 assert(value_zero_p(res
->d
));
1497 assert(res
->x
.p
->type
== partition
);
1499 s
= (struct section
*)
1500 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1501 sizeof(struct section
));
1505 evalue_set_si(&mone
, -1, 1);
1508 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1509 assert(mask
->x
.p
->size
>= 2);
1510 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1511 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1513 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1515 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1524 value_init(s
[n
].E
.d
);
1525 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1529 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1530 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1533 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1534 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1535 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1536 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1538 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1539 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1545 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1546 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1548 value_init(s
[n
].E
.d
);
1549 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1550 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1556 /* Just ignore; this may have been previously masked off */
1558 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1562 free_evalue_refs(&mone
);
1563 free_evalue_refs(mask
);
1564 free_evalue_refs(res
);
1567 evalue_set_si(res
, 0, 1);
1569 res
->x
.p
= new_enode(partition
, 2*n
, -1);
1570 for (j
= 0; j
< n
; ++j
) {
1571 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1572 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1573 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1580 void evalue_copy(evalue
*dst
, evalue
*src
)
1582 value_assign(dst
->d
, src
->d
);
1583 if(value_notzero_p(src
->d
)) {
1584 value_init(dst
->x
.n
);
1585 value_assign(dst
->x
.n
, src
->x
.n
);
1587 dst
->x
.p
= ecopy(src
->x
.p
);
1590 enode
*new_enode(enode_type type
,int size
,int pos
) {
1596 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1599 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1603 for(i
=0; i
<size
; i
++) {
1604 value_init(res
->arr
[i
].d
);
1605 value_set_si(res
->arr
[i
].d
,0);
1606 res
->arr
[i
].x
.p
= 0;
1611 enode
*ecopy(enode
*e
) {
1616 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1617 for(i
=0;i
<e
->size
;++i
) {
1618 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1619 if(value_zero_p(res
->arr
[i
].d
))
1620 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1621 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1622 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1624 value_init(res
->arr
[i
].x
.n
);
1625 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1631 int ecmp(const evalue
*e1
, const evalue
*e2
)
1637 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
1641 value_multiply(m
, e1
->x
.n
, e2
->d
);
1642 value_multiply(m2
, e2
->x
.n
, e1
->d
);
1644 if (value_lt(m
, m2
))
1646 else if (value_gt(m
, m2
))
1656 if (value_notzero_p(e1
->d
))
1658 if (value_notzero_p(e2
->d
))
1664 if (p1
->type
!= p2
->type
)
1665 return p1
->type
- p2
->type
;
1666 if (p1
->pos
!= p2
->pos
)
1667 return p1
->pos
- p2
->pos
;
1668 if (p1
->size
!= p2
->size
)
1669 return p1
->size
- p2
->size
;
1671 for (i
= p1
->size
-1; i
>= 0; --i
)
1672 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
1678 int eequal(evalue
*e1
,evalue
*e2
) {
1683 if (value_ne(e1
->d
,e2
->d
))
1686 /* e1->d == e2->d */
1687 if (value_notzero_p(e1
->d
)) {
1688 if (value_ne(e1
->x
.n
,e2
->x
.n
))
1691 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1695 /* e1->d == e2->d == 0 */
1698 if (p1
->type
!= p2
->type
) return 0;
1699 if (p1
->size
!= p2
->size
) return 0;
1700 if (p1
->pos
!= p2
->pos
) return 0;
1701 for (i
=0; i
<p1
->size
; i
++)
1702 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
1707 void free_evalue_refs(evalue
*e
) {
1712 if (EVALUE_IS_DOMAIN(*e
)) {
1713 Domain_Free(EVALUE_DOMAIN(*e
));
1716 } else if (value_pos_p(e
->d
)) {
1718 /* 'e' stores a constant */
1720 value_clear(e
->x
.n
);
1723 assert(value_zero_p(e
->d
));
1726 if (!p
) return; /* null pointer */
1727 for (i
=0; i
<p
->size
; i
++) {
1728 free_evalue_refs(&(p
->arr
[i
]));
1732 } /* free_evalue_refs */
1734 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
1735 Vector
* val
, evalue
*res
)
1737 unsigned nparam
= periods
->Size
;
1740 double d
= compute_evalue(e
, val
->p
);
1741 d
*= VALUE_TO_DOUBLE(m
);
1746 value_assign(res
->d
, m
);
1747 value_init(res
->x
.n
);
1748 value_set_double(res
->x
.n
, d
);
1749 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
1752 if (value_one_p(periods
->p
[p
]))
1753 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
1758 value_assign(tmp
, periods
->p
[p
]);
1759 value_set_si(res
->d
, 0);
1760 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
1762 value_decrement(tmp
, tmp
);
1763 value_assign(val
->p
[p
], tmp
);
1764 mod2table_r(e
, periods
, m
, p
+1, val
,
1765 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
1766 } while (value_pos_p(tmp
));
1772 static void rel2table(evalue
*e
, int zero
)
1774 if (value_pos_p(e
->d
)) {
1775 if (value_zero_p(e
->x
.n
) == zero
)
1776 value_set_si(e
->x
.n
, 1);
1778 value_set_si(e
->x
.n
, 0);
1779 value_set_si(e
->d
, 1);
1782 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
1783 rel2table(&e
->x
.p
->arr
[i
], zero
);
1787 void evalue_mod2table(evalue
*e
, int nparam
)
1792 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
1795 for (i
=0; i
<p
->size
; i
++) {
1796 evalue_mod2table(&(p
->arr
[i
]), nparam
);
1798 if (p
->type
== relation
) {
1803 evalue_copy(©
, &p
->arr
[0]);
1805 rel2table(&p
->arr
[0], 1);
1806 emul(&p
->arr
[0], &p
->arr
[1]);
1808 rel2table(©
, 0);
1809 emul(©
, &p
->arr
[2]);
1810 eadd(&p
->arr
[2], &p
->arr
[1]);
1811 free_evalue_refs(&p
->arr
[2]);
1812 free_evalue_refs(©
);
1814 free_evalue_refs(&p
->arr
[0]);
1818 } else if (p
->type
== fractional
) {
1819 Vector
*periods
= Vector_Alloc(nparam
);
1820 Vector
*val
= Vector_Alloc(nparam
);
1826 value_set_si(tmp
, 1);
1827 Vector_Set(periods
->p
, 1, nparam
);
1828 Vector_Set(val
->p
, 0, nparam
);
1829 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
1832 assert(p
->type
== polynomial
);
1833 assert(p
->size
== 2);
1834 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
1835 Lcm3(tmp
, p
->arr
[1].d
, &tmp
);
1838 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
1841 evalue_set_si(&res
, 0, 1);
1842 /* Compute the polynomial using Horner's rule */
1843 for (i
=p
->size
-1;i
>1;i
--) {
1844 eadd(&p
->arr
[i
], &res
);
1847 eadd(&p
->arr
[1], &res
);
1849 free_evalue_refs(e
);
1850 free_evalue_refs(&EP
);
1855 Vector_Free(periods
);
1857 } /* evalue_mod2table */
1859 /********************************************************/
1860 /* function in domain */
1861 /* check if the parameters in list_args */
1862 /* verifies the constraints of Domain P */
1863 /********************************************************/
1864 int in_domain(Polyhedron
*P
, Value
*list_args
) {
1867 Value v
; /* value of the constraint of a row when
1868 parameters are instanciated*/
1874 /*P->Constraint constraint matrice of polyhedron P*/
1875 for(row
=0;row
<P
->NbConstraints
;row
++) {
1876 value_assign(v
,P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
1877 for(col
=1;col
<P
->Dimension
+1;col
++) {
1878 value_multiply(tmp
,P
->Constraint
[row
][col
],list_args
[col
-1]);
1879 value_addto(v
,v
,tmp
);
1881 if (value_notzero_p(P
->Constraint
[row
][0])) {
1883 /*if v is not >=0 then this constraint is not respected */
1884 if (value_neg_p(v
)) {
1888 return P
->next
? in_domain(P
->next
, list_args
) : 0;
1893 /*if v is not = 0 then this constraint is not respected */
1894 if (value_notzero_p(v
))
1899 /*if not return before this point => all
1900 the constraints are respected */
1906 /****************************************************/
1907 /* function compute enode */
1908 /* compute the value of enode p with parameters */
1909 /* list "list_args */
1910 /* compute the polynomial or the periodic */
1911 /****************************************************/
1913 static double compute_enode(enode
*p
, Value
*list_args
) {
1925 if (p
->type
== polynomial
) {
1927 value_assign(param
,list_args
[p
->pos
-1]);
1929 /* Compute the polynomial using Horner's rule */
1930 for (i
=p
->size
-1;i
>0;i
--) {
1931 res
+=compute_evalue(&p
->arr
[i
],list_args
);
1932 res
*=VALUE_TO_DOUBLE(param
);
1934 res
+=compute_evalue(&p
->arr
[0],list_args
);
1936 else if (p
->type
== fractional
) {
1937 double d
= compute_evalue(&p
->arr
[0], list_args
);
1938 d
-= floor(d
+1e-10);
1940 /* Compute the polynomial using Horner's rule */
1941 for (i
=p
->size
-1;i
>1;i
--) {
1942 res
+=compute_evalue(&p
->arr
[i
],list_args
);
1945 res
+=compute_evalue(&p
->arr
[1],list_args
);
1947 else if (p
->type
== periodic
) {
1948 value_assign(param
,list_args
[p
->pos
-1]);
1950 /* Choose the right element of the periodic */
1951 value_absolute(m
,param
);
1952 value_set_si(param
,p
->size
);
1953 value_modulus(m
,m
,param
);
1954 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
1956 else if (p
->type
== relation
) {
1957 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
1958 res
= compute_evalue(&p
->arr
[1], list_args
);
1959 else if (p
->size
> 2)
1960 res
= compute_evalue(&p
->arr
[2], list_args
);
1962 else if (p
->type
== partition
) {
1963 for (i
= 0; i
< p
->size
/2; ++i
)
1964 if (in_domain(EVALUE_DOMAIN(p
->arr
[2*i
]), list_args
)) {
1965 res
= compute_evalue(&p
->arr
[2*i
+1], list_args
);
1974 } /* compute_enode */
1976 /*************************************************/
1977 /* return the value of Ehrhart Polynomial */
1978 /* It returns a double, because since it is */
1979 /* a recursive function, some intermediate value */
1980 /* might not be integral */
1981 /*************************************************/
1983 double compute_evalue(evalue
*e
,Value
*list_args
) {
1987 if (value_notzero_p(e
->d
)) {
1988 if (value_notone_p(e
->d
))
1989 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
1991 res
= VALUE_TO_DOUBLE(e
->x
.n
);
1994 res
= compute_enode(e
->x
.p
,list_args
);
1996 } /* compute_evalue */
1999 /****************************************************/
2000 /* function compute_poly : */
2001 /* Check for the good validity domain */
2002 /* return the number of point in the Polyhedron */
2003 /* in allocated memory */
2004 /* Using the Ehrhart pseudo-polynomial */
2005 /****************************************************/
2006 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2009 /* double d; int i; */
2011 tmp
= (Value
*) malloc (sizeof(Value
));
2012 assert(tmp
!= NULL
);
2014 value_set_si(*tmp
,0);
2017 return(tmp
); /* no ehrhart polynomial */
2018 if(en
->ValidityDomain
) {
2019 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2020 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2025 return(tmp
); /* no Validity Domain */
2027 if(in_domain(en
->ValidityDomain
,list_args
)) {
2029 #ifdef EVAL_EHRHART_DEBUG
2030 Print_Domain(stdout
,en
->ValidityDomain
);
2031 print_evalue(stdout
,&en
->EP
);
2034 /* d = compute_evalue(&en->EP,list_args);
2036 printf("(double)%lf = %d\n", d, i ); */
2037 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2043 value_set_si(*tmp
,0);
2044 return(tmp
); /* no compatible domain with the arguments */
2045 } /* compute_poly */
2047 size_t value_size(Value v
) {
2048 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2049 * sizeof(v
[0]._mp_d
[0]);
2052 size_t domain_size(Polyhedron
*D
)
2055 size_t s
= sizeof(*D
);
2057 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2058 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2059 s
+= value_size(D
->Constraint
[i
][j
]);
2062 for (i = 0; i < D->NbRays; ++i)
2063 for (j = 0; j < D->Dimension+2; ++j)
2064 s += value_size(D->Ray[i][j]);
2067 return D
->next
? s
+domain_size(D
->next
) : s
;
2070 size_t enode_size(enode
*p
) {
2071 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2074 if (p
->type
== partition
)
2075 for (i
= 0; i
< p
->size
/2; ++i
) {
2076 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2077 s
+= evalue_size(&p
->arr
[2*i
+1]);
2080 for (i
= 0; i
< p
->size
; ++i
) {
2081 s
+= evalue_size(&p
->arr
[i
]);
2086 size_t evalue_size(evalue
*e
)
2088 size_t s
= sizeof(*e
);
2089 s
+= value_size(e
->d
);
2090 if (value_notzero_p(e
->d
))
2091 s
+= value_size(e
->x
.n
);
2093 s
+= enode_size(e
->x
.p
);
2097 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2099 evalue
*found
= NULL
;
2104 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2107 value_init(offset
.d
);
2108 value_init(offset
.x
.n
);
2109 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2110 Lcm3(m
, offset
.d
, &offset
.d
);
2111 value_set_si(offset
.x
.n
, 1);
2114 evalue_copy(©
, cst
);
2117 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2119 if (eequal(base
, &e
->x
.p
->arr
[0]))
2120 found
= &e
->x
.p
->arr
[0];
2122 value_set_si(offset
.x
.n
, -2);
2125 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2127 if (eequal(base
, &e
->x
.p
->arr
[0]))
2130 free_evalue_refs(cst
);
2131 free_evalue_refs(&offset
);
2134 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2135 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2140 static evalue
*find_relation_pair(evalue
*e
)
2143 evalue
*found
= NULL
;
2145 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2148 if (e
->x
.p
->type
== fractional
) {
2153 poly_denom(&e
->x
.p
->arr
[0], &m
);
2155 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2156 cst
= &cst
->x
.p
->arr
[0])
2159 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2160 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2165 i
= e
->x
.p
->type
== relation
;
2166 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2167 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2172 void evalue_mod2relation(evalue
*e
) {
2175 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2178 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2179 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2180 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2181 value_clear(e
->x
.p
->arr
[2*i
].d
);
2182 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2184 if (2*i
< e
->x
.p
->size
) {
2185 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2186 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2191 if (e
->x
.p
->size
== 0) {
2193 evalue_set_si(e
, 0, 1);
2199 while ((d
= find_relation_pair(e
)) != NULL
) {
2203 value_init(split
.d
);
2204 value_set_si(split
.d
, 0);
2205 split
.x
.p
= new_enode(relation
, 3, 0);
2206 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2207 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2209 ev
= &split
.x
.p
->arr
[0];
2210 value_set_si(ev
->d
, 0);
2211 ev
->x
.p
= new_enode(fractional
, 3, -1);
2212 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2213 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2214 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2220 free_evalue_refs(&split
);
2224 static int evalue_comp(const void * a
, const void * b
)
2226 const evalue
*e1
= *(const evalue
**)a
;
2227 const evalue
*e2
= *(const evalue
**)b
;
2228 return ecmp(e1
, e2
);
2231 void evalue_combine(evalue
*e
)
2238 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2241 NALLOC(evs
, e
->x
.p
->size
/2);
2242 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2243 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2244 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2245 p
= new_enode(partition
, e
->x
.p
->size
, -1);
2246 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2247 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2248 value_clear(p
->arr
[2*k
].d
);
2249 value_clear(p
->arr
[2*k
+1].d
);
2250 p
->arr
[2*k
] = *(evs
[i
]-1);
2251 p
->arr
[2*k
+1] = *(evs
[i
]);
2254 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2257 value_clear((evs
[i
]-1)->d
);
2261 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2262 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2263 free_evalue_refs(evs
[i
]);
2267 for (i
= 2*k
; i
< p
->size
; ++i
)
2268 value_clear(p
->arr
[i
].d
);
2275 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2277 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2279 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2282 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2283 Polyhedron
*D
, *N
, **P
;
2286 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2293 if (D
->NbEq
<= H
->NbEq
) {
2299 tmp
.x
.p
= new_enode(partition
, 2, -1);
2300 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2301 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2302 reduce_evalue(&tmp
);
2303 if (value_notzero_p(tmp
.d
) ||
2304 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2307 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2308 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2311 free_evalue_refs(&tmp
);
2317 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2319 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2321 value_clear(e
->x
.p
->arr
[2*i
].d
);
2322 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2324 if (2*i
< e
->x
.p
->size
) {
2325 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2326 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2333 H
= DomainConvex(D
, 0);
2334 E
= DomainDifference(H
, D
, 0);
2336 D
= DomainDifference(H
, E
, 0);
2339 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2343 /* May change coefficients to become non-standard if fiddle is set
2344 * => reduce p afterwards to correct
2346 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2347 Matrix
**R
, int fiddle
)
2351 unsigned dim
= D
->Dimension
;
2352 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2357 assert(p
->type
== fractional
);
2359 value_set_si(T
->p
[1][dim
], 1);
2361 while (value_zero_p(pp
->d
)) {
2362 assert(pp
->x
.p
->type
== polynomial
);
2363 assert(pp
->x
.p
->size
== 2);
2364 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2365 value_division(T
->p
[0][pp
->x
.p
->pos
-1], *d
, pp
->x
.p
->arr
[1].d
);
2366 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2367 if (fiddle
&& value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2368 value_substract(pp
->x
.p
->arr
[1].x
.n
,
2369 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2370 value_multiply(T
->p
[0][pp
->x
.p
->pos
-1],
2371 T
->p
[0][pp
->x
.p
->pos
-1], pp
->x
.p
->arr
[1].x
.n
);
2372 pp
= &pp
->x
.p
->arr
[0];
2374 value_division(T
->p
[0][dim
], *d
, pp
->d
);
2375 value_multiply(T
->p
[0][dim
], T
->p
[0][dim
], pp
->x
.n
);
2376 I
= DomainImage(D
, T
, 0);
2377 H
= DomainConvex(I
, 0);
2386 static int reduce_in_domain(evalue
*e
, Polyhedron
*D
)
2395 if (value_notzero_p(e
->d
))
2400 if (p
->type
== relation
) {
2406 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
, 1);
2407 line_minmax(I
, &min
, &max
); /* frees I */
2408 equal
= value_eq(min
, max
);
2409 mpz_cdiv_q(min
, min
, d
);
2410 mpz_fdiv_q(max
, max
, d
);
2412 if (value_gt(min
, max
)) {
2418 evalue_set_si(e
, 0, 1);
2421 free_evalue_refs(&(p
->arr
[1]));
2422 free_evalue_refs(&(p
->arr
[0]));
2428 return r
? r
: reduce_in_domain(e
, D
);
2432 free_evalue_refs(&(p
->arr
[2]));
2435 free_evalue_refs(&(p
->arr
[0]));
2441 return reduce_in_domain(e
, D
);
2442 } else if (value_eq(min
, max
)) {
2443 /* zero for a single value */
2445 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2446 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2447 value_multiply(min
, min
, d
);
2448 value_substract(M
->p
[0][D
->Dimension
+1],
2449 M
->p
[0][D
->Dimension
+1], min
);
2450 E
= DomainAddConstraints(D
, M
, 0);
2456 r
= reduce_in_domain(&p
->arr
[1], E
);
2458 r
|= reduce_in_domain(&p
->arr
[2], D
);
2460 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2468 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2471 i
= p
->type
== relation
? 1 :
2472 p
->type
== fractional
? 1 : 0;
2473 for (; i
<p
->size
; i
++)
2474 r
|= reduce_in_domain(&p
->arr
[i
], D
);
2476 if (p
->type
!= fractional
) {
2477 if (r
&& p
->type
== polynomial
) {
2480 value_set_si(f
.d
, 0);
2481 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2482 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2483 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2484 reorder_terms(p
, &f
);
2495 I
= polynomial_projection(p
, D
, &d
, &T
, 1);
2497 line_minmax(I
, &min
, &max
); /* frees I */
2498 mpz_fdiv_q(min
, min
, d
);
2499 mpz_fdiv_q(max
, max
, d
);
2501 if (value_eq(min
, max
)) {
2504 value_init(inc
.x
.n
);
2505 value_set_si(inc
.d
, 1);
2506 value_oppose(inc
.x
.n
, min
);
2507 eadd(&inc
, &p
->arr
[0]);
2508 reorder_terms(p
, &p
->arr
[0]); /* frees arr[0] */
2512 free_evalue_refs(&inc
);
2515 _reduce_evalue(&p
->arr
[0], 0, 1);
2519 value_set_si(f
.d
, 0);
2520 f
.x
.p
= new_enode(fractional
, 3, -1);
2521 value_clear(f
.x
.p
->arr
[0].d
);
2522 f
.x
.p
->arr
[0] = p
->arr
[0];
2523 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2524 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
2525 reorder_terms(p
, &f
);
2539 void evalue_range_reduction(evalue
*e
)
2542 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2545 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2546 if (reduce_in_domain(&e
->x
.p
->arr
[2*i
+1],
2547 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
2548 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
2550 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2551 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2552 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
2553 value_clear(e
->x
.p
->arr
[2*i
].d
);
2555 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2556 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2564 Enumeration
* partition2enumeration(evalue
*EP
)
2567 Enumeration
*en
, *res
= NULL
;
2569 if (EVALUE_IS_ZERO(*EP
)) {
2574 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
2575 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
2578 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
2579 value_clear(EP
->x
.p
->arr
[2*i
].d
);
2580 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
2588 static int frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
2598 if (value_notzero_p(e
->d
))
2603 i
= p
->type
== relation
? 1 :
2604 p
->type
== fractional
? 1 : 0;
2605 for (; i
<p
->size
; i
++)
2606 r
|= frac2floor_in_domain(&p
->arr
[i
], D
);
2608 if (p
->type
!= fractional
) {
2609 if (r
&& p
->type
== polynomial
) {
2612 value_set_si(f
.d
, 0);
2613 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2614 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2615 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2616 reorder_terms(p
, &f
);
2625 I
= polynomial_projection(p
, D
, &d
, &T
, 0);
2628 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2631 assert(I
->NbEq
== 0); /* Should have been reduced */
2634 for (i
= 0; i
< I
->NbConstraints
; ++i
)
2635 if (value_pos_p(I
->Constraint
[i
][1]))
2638 assert(i
< I
->NbConstraints
);
2640 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
2641 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
2642 if (value_neg_p(min
)) {
2644 mpz_fdiv_q(min
, min
, d
);
2645 value_init(offset
.d
);
2646 value_set_si(offset
.d
, 1);
2647 value_init(offset
.x
.n
);
2648 value_oppose(offset
.x
.n
, min
);
2649 eadd(&offset
, &p
->arr
[0]);
2650 free_evalue_refs(&offset
);
2659 value_set_si(fl
.d
, 0);
2660 fl
.x
.p
= new_enode(flooring
, 3, -1);
2661 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
2662 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
2663 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
2665 eadd(&fl
, &p
->arr
[0]);
2666 reorder_terms(p
, &p
->arr
[0]);
2669 free_evalue_refs(&fl
);
2674 void evalue_frac2floor(evalue
*e
)
2677 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2680 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2681 if (frac2floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
2682 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
])))
2683 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
2686 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
2691 int nparam
= D
->Dimension
- nvar
;
2694 nr
= D
->NbConstraints
+ 2;
2695 nc
= D
->Dimension
+ 2 + 1;
2696 C
= Matrix_Alloc(nr
, nc
);
2697 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
2698 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
2699 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
2700 D
->Dimension
+ 1 - nvar
);
2705 nc
= C
->NbColumns
+ 1;
2706 C
= Matrix_Alloc(nr
, nc
);
2707 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
2708 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
2709 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
2710 oldC
->NbColumns
- 1 - nvar
);
2713 value_set_si(C
->p
[nr
-2][0], 1);
2714 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
2715 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
2717 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
2718 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
2724 evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
2731 evalue
*factor
= NULL
;
2734 if (EVALUE_IS_ZERO(*e
))
2737 if (value_notzero_p(e
->d
)) {
2739 int nparam
= D
->Dimension
- nvar
;
2743 D
= Constraints2Polyhedron(C
, 0);
2747 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
2752 if (!EVALUE_IS_ONE(*e
))
2758 switch (e
->x
.p
->type
) {
2760 evalue
*pp
= &e
->x
.p
->arr
[0];
2761 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
2762 poly_denom(pp
, &row
->p
[1 + nvar
]);
2763 value_set_si(row
->p
[0], 1);
2764 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
2765 pp
= &pp
->x
.p
->arr
[0]) {
2767 assert(pp
->x
.p
->type
== polynomial
);
2769 if (pos
>= 1 + nvar
)
2771 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
2772 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
2773 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
2775 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
2776 value_division(row
->p
[1 + D
->Dimension
+ 1],
2777 row
->p
[1 + D
->Dimension
+ 1],
2779 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
2780 row
->p
[1 + D
->Dimension
+ 1],
2782 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
2786 int pos
= e
->x
.p
->pos
;
2790 value_init(factor
->d
);
2791 value_set_si(factor
->d
, 0);
2792 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
2793 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
2794 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
2798 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
2799 for (i
= 0; i
< D
->NbRays
; ++i
)
2800 if (value_notzero_p(D
->Ray
[i
][pos
]))
2802 assert(i
< D
->NbRays
);
2803 if (value_neg_p(D
->Ray
[i
][pos
])) {
2805 value_init(factor
->d
);
2806 evalue_set_si(factor
, -1, 1);
2808 value_set_si(row
->p
[0], 1);
2809 value_set_si(row
->p
[pos
], 1);
2810 value_set_si(row
->p
[1 + nvar
], -1);
2817 i
= type_offset(e
->x
.p
);
2819 res
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
2824 evalue_copy(&cum
, factor
);
2828 for (; i
< e
->x
.p
->size
; ++i
) {
2832 C
= esum_add_constraint(nvar
, D
, C
, row
);
2838 Vector_Print(stderr, P_VALUE_FMT, row);
2840 Matrix_Print(stderr, P_VALUE_FMT, C);
2842 t
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
2851 free_evalue_refs(t
);
2854 if (factor
&& i
+1 < e
->x
.p
->size
)
2861 free_evalue_refs(factor
);
2862 free_evalue_refs(&cum
);
2872 evalue
*esum(evalue
*e
, int nvar
)
2880 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
2881 evalue_copy(res
, e
);
2885 evalue_set_si(res
, 0, 1);
2887 assert(value_zero_p(e
->d
));
2888 assert(e
->x
.p
->type
== partition
);
2890 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2891 char *test
[] = { "A", "B", "C", "D", "E", "F", "G" };
2893 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
2894 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2895 Polyhedron_Print(stderr
, P_VALUE_FMT
, EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
2896 print_evalue(stderr
, t
, test
);
2898 free_evalue_refs(t
);