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 void addeliminatedparams_enum(evalue
*e
, Matrix
*CT
, Polyhedron
*CEq
,
76 unsigned MaxRays
, unsigned nparam
)
81 if (CT
->NbRows
== CT
->NbColumns
)
84 if (EVALUE_IS_ZERO(*e
))
87 if (value_notzero_p(e
->d
)) {
90 value_set_si(res
.d
, 0);
91 res
.x
.p
= new_enode(partition
, 2, nparam
);
92 EVALUE_SET_DOMAIN(res
.x
.p
->arr
[0],
93 DomainConstraintSimplify(Polyhedron_Copy(CEq
), MaxRays
));
94 value_clear(res
.x
.p
->arr
[1].d
);
102 assert(p
->type
== partition
);
105 for (i
=0; i
<p
->size
/2; i
++) {
106 Polyhedron
*D
= EVALUE_DOMAIN(p
->arr
[2*i
]);
107 Polyhedron
*T
= DomainPreimage(D
, CT
, MaxRays
);
110 T
= DomainIntersection(D
, CEq
, MaxRays
);
112 EVALUE_SET_DOMAIN(p
->arr
[2*i
], T
);
113 addeliminatedparams_evalue(&p
->arr
[2*i
+1], CT
);
117 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
125 assert(value_notzero_p(e1
->d
));
126 assert(value_notzero_p(e2
->d
));
127 value_multiply(m
, e1
->x
.n
, e2
->d
);
128 value_multiply(m2
, e2
->x
.n
, e1
->d
);
131 else if (value_gt(m
, m2
))
141 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
143 if (value_notzero_p(e1
->d
)) {
144 if (value_zero_p(e2
->d
))
146 int r
= mod_rational_smaller(e1
, e2
);
147 return r
== -1 ? 0 : r
;
149 if (value_notzero_p(e2
->d
))
151 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
153 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
156 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
158 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
163 static int mod_term_smaller(evalue
*e1
, evalue
*e2
)
165 assert(value_zero_p(e1
->d
));
166 assert(value_zero_p(e2
->d
));
167 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
168 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
169 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
172 /* Negative pos means inequality */
173 /* s is negative of substitution if m is not zero */
182 struct fixed_param
*fixed
;
187 static int relations_depth(evalue
*e
)
192 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
193 e
= &e
->x
.p
->arr
[1], ++d
);
197 static void poly_denom(evalue
*p
, Value
*d
)
201 while (value_zero_p(p
->d
)) {
202 assert(p
->x
.p
->type
== polynomial
);
203 assert(p
->x
.p
->size
== 2);
204 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
205 value_lcm(*d
, p
->x
.p
->arr
[1].d
, d
);
208 value_lcm(*d
, p
->d
, d
);
211 #define EVALUE_IS_ONE(ev) (value_one_p((ev).d) && value_one_p((ev).x.n))
213 static void realloc_substitution(struct subst
*s
, int d
)
215 struct fixed_param
*n
;
218 for (i
= 0; i
< s
->n
; ++i
)
225 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
231 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
234 /* May have been reduced already */
235 if (value_notzero_p(m
->d
))
238 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
239 assert(m
->x
.p
->size
== 3);
241 /* fractional was inverted during reduction
242 * invert it back and move constant in
244 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
245 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
246 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
247 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
248 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
249 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
250 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
251 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
252 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
253 value_set_si(m
->x
.p
->arr
[1].d
, 1);
256 /* Oops. Nested identical relations. */
257 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
260 if (s
->n
>= s
->max
) {
261 int d
= relations_depth(r
);
262 realloc_substitution(s
, d
);
266 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
267 assert(p
->x
.p
->size
== 2);
270 assert(value_pos_p(f
->x
.n
));
272 value_init(s
->fixed
[s
->n
].m
);
273 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
274 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
275 value_init(s
->fixed
[s
->n
].d
);
276 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
277 value_init(s
->fixed
[s
->n
].s
.d
);
278 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
284 static int type_offset(enode
*p
)
286 return p
->type
== fractional
? 1 :
287 p
->type
== flooring
? 1 : 0;
290 static void reorder_terms(enode
*p
, evalue
*v
)
293 int offset
= type_offset(p
);
295 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
297 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
298 free_evalue_refs(&(p
->arr
[i
]));
304 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
310 if (value_notzero_p(e
->d
)) {
312 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
313 return; /* a rational number, its already reduced */
317 return; /* hum... an overflow probably occured */
319 /* First reduce the components of p */
320 add
= p
->type
== relation
;
321 for (i
=0; i
<p
->size
; i
++) {
323 add
= add_modulo_substitution(s
, e
);
325 if (i
== 0 && p
->type
==fractional
)
326 _reduce_evalue(&p
->arr
[i
], s
, 1);
328 _reduce_evalue(&p
->arr
[i
], s
, fract
);
330 if (add
&& i
== p
->size
-1) {
332 value_clear(s
->fixed
[s
->n
].m
);
333 value_clear(s
->fixed
[s
->n
].d
);
334 free_evalue_refs(&s
->fixed
[s
->n
].s
);
335 } else if (add
&& i
== 1)
336 s
->fixed
[s
->n
-1].pos
*= -1;
339 if (p
->type
==periodic
) {
341 /* Try to reduce the period */
342 for (i
=1; i
<=(p
->size
)/2; i
++) {
343 if ((p
->size
% i
)==0) {
345 /* Can we reduce the size to i ? */
347 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
348 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
351 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
355 you_lose
: /* OK, lets not do it */
360 /* Try to reduce its strength */
363 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
367 else if (p
->type
==polynomial
) {
368 for (k
= 0; s
&& k
< s
->n
; ++k
) {
369 if (s
->fixed
[k
].pos
== p
->pos
) {
370 int divide
= value_notone_p(s
->fixed
[k
].d
);
373 if (value_notzero_p(s
->fixed
[k
].m
)) {
376 assert(p
->size
== 2);
377 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
379 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
386 value_assign(d
.d
, s
->fixed
[k
].d
);
388 if (value_notzero_p(s
->fixed
[k
].m
))
389 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
391 value_set_si(d
.x
.n
, 1);
394 for (i
=p
->size
-1;i
>=1;i
--) {
395 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
397 emul(&d
, &p
->arr
[i
]);
398 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
399 free_evalue_refs(&(p
->arr
[i
]));
402 _reduce_evalue(&p
->arr
[0], s
, fract
);
405 free_evalue_refs(&d
);
411 /* Try to reduce the degree */
412 for (i
=p
->size
-1;i
>=1;i
--) {
413 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
415 /* Zero coefficient */
416 free_evalue_refs(&(p
->arr
[i
]));
421 /* Try to reduce its strength */
424 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
428 else if (p
->type
==fractional
) {
432 if (value_notzero_p(p
->arr
[0].d
)) {
434 value_assign(v
.d
, p
->arr
[0].d
);
436 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
441 evalue
*pp
= &p
->arr
[0];
442 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
443 assert(pp
->x
.p
->size
== 2);
445 /* search for exact duplicate among the modulo inequalities */
447 f
= &pp
->x
.p
->arr
[1];
448 for (k
= 0; s
&& k
< s
->n
; ++k
) {
449 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
450 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
451 value_eq(s
->fixed
[k
].m
, f
->d
) &&
452 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
459 /* replace { E/m } by { (E-1)/m } + 1/m */
464 evalue_set_si(&extra
, 1, 1);
465 value_assign(extra
.d
, g
);
466 eadd(&extra
, &v
.x
.p
->arr
[1]);
467 free_evalue_refs(&extra
);
469 /* We've been going in circles; stop now */
470 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
471 free_evalue_refs(&v
);
473 evalue_set_si(&v
, 0, 1);
478 value_set_si(v
.d
, 0);
479 v
.x
.p
= new_enode(fractional
, 3, -1);
480 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
481 value_assign(v
.x
.p
->arr
[1].d
, g
);
482 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
483 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
486 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
489 value_division(f
->d
, g
, f
->d
);
490 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
491 value_assign(f
->d
, g
);
492 value_decrement(f
->x
.n
, f
->x
.n
);
493 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
495 Gcd(f
->d
, f
->x
.n
, &g
);
496 value_division(f
->d
, f
->d
, g
);
497 value_division(f
->x
.n
, f
->x
.n
, g
);
506 /* reduction may have made this fractional arg smaller */
507 i
= reorder
? p
->size
: 1;
508 for ( ; i
< p
->size
; ++i
)
509 if (value_zero_p(p
->arr
[i
].d
) &&
510 p
->arr
[i
].x
.p
->type
== fractional
&&
511 !mod_term_smaller(e
, &p
->arr
[i
]))
515 value_set_si(v
.d
, 0);
516 v
.x
.p
= new_enode(fractional
, 3, -1);
517 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
518 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
519 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
529 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
530 pp
= &pp
->x
.p
->arr
[0]) {
531 f
= &pp
->x
.p
->arr
[1];
532 assert(value_pos_p(f
->d
));
533 mpz_mul_ui(twice
, f
->x
.n
, 2);
534 if (value_lt(twice
, f
->d
))
536 if (value_eq(twice
, f
->d
))
544 value_set_si(v
.d
, 0);
545 v
.x
.p
= new_enode(fractional
, 3, -1);
546 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
547 poly_denom(&p
->arr
[0], &twice
);
548 value_assign(v
.x
.p
->arr
[1].d
, twice
);
549 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
550 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
551 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
553 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
554 pp
= &pp
->x
.p
->arr
[0]) {
555 f
= &pp
->x
.p
->arr
[1];
556 value_oppose(f
->x
.n
, f
->x
.n
);
557 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
559 value_division(pp
->d
, twice
, pp
->d
);
560 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
561 value_assign(pp
->d
, twice
);
562 value_oppose(pp
->x
.n
, pp
->x
.n
);
563 value_decrement(pp
->x
.n
, pp
->x
.n
);
564 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
574 reorder_terms(p
, &v
);
575 _reduce_evalue(&p
->arr
[1], s
, fract
);
578 /* Try to reduce the degree */
579 for (i
=p
->size
-1;i
>=2;i
--) {
580 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
582 /* Zero coefficient */
583 free_evalue_refs(&(p
->arr
[i
]));
588 /* Try to reduce its strength */
591 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
592 free_evalue_refs(&(p
->arr
[0]));
596 else if (p
->type
== flooring
) {
597 /* Try to reduce the degree */
598 for (i
=p
->size
-1;i
>=2;i
--) {
599 if (!EVALUE_IS_ZERO(p
->arr
[i
]))
601 /* Zero coefficient */
602 free_evalue_refs(&(p
->arr
[i
]));
607 /* Try to reduce its strength */
610 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
611 free_evalue_refs(&(p
->arr
[0]));
615 else if (p
->type
== relation
) {
616 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
617 free_evalue_refs(&(p
->arr
[2]));
618 free_evalue_refs(&(p
->arr
[0]));
625 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
626 free_evalue_refs(&(p
->arr
[2]));
629 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
630 free_evalue_refs(&(p
->arr
[1]));
631 free_evalue_refs(&(p
->arr
[0]));
632 evalue_set_si(e
, 0, 1);
639 /* Relation was reduced by means of an identical
640 * inequality => remove
642 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
645 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
646 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
648 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
650 free_evalue_refs(&(p
->arr
[2]));
654 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
656 evalue_set_si(e
, 0, 1);
657 free_evalue_refs(&(p
->arr
[1]));
659 free_evalue_refs(&(p
->arr
[0]));
665 } /* reduce_evalue */
667 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
672 for (k
= 0; k
< dim
; ++k
)
673 if (value_notzero_p(row
[k
+1]))
676 Vector_Normalize_Positive(row
+1, dim
+1, k
);
677 assert(s
->n
< s
->max
);
678 value_init(s
->fixed
[s
->n
].d
);
679 value_init(s
->fixed
[s
->n
].m
);
680 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
681 s
->fixed
[s
->n
].pos
= k
+1;
682 value_set_si(s
->fixed
[s
->n
].m
, 0);
683 r
= &s
->fixed
[s
->n
].s
;
685 for (l
= k
+1; l
< dim
; ++l
)
686 if (value_notzero_p(row
[l
+1])) {
687 value_set_si(r
->d
, 0);
688 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
689 value_init(r
->x
.p
->arr
[1].x
.n
);
690 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
691 value_set_si(r
->x
.p
->arr
[1].d
, 1);
695 value_oppose(r
->x
.n
, row
[dim
+1]);
696 value_set_si(r
->d
, 1);
700 void reduce_evalue (evalue
*e
) {
701 struct subst s
= { NULL
, 0, 0 };
703 if (value_notzero_p(e
->d
))
704 return; /* a rational number, its already reduced */
706 if (e
->x
.p
->type
== partition
) {
709 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
711 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
712 /* This shouldn't really happen;
713 * Empty domains should not be added.
720 D
= DomainConvex(D
, 0);
721 if (!D
->next
&& D
->NbEq
) {
725 realloc_substitution(&s
, dim
);
727 int d
= relations_depth(&e
->x
.p
->arr
[2*i
+1]);
729 NALLOC(s
.fixed
, s
.max
);
732 for (j
= 0; j
< D
->NbEq
; ++j
)
733 add_substitution(&s
, D
->Constraint
[j
], dim
);
735 if (D
!= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))
737 _reduce_evalue(&e
->x
.p
->arr
[2*i
+1], &s
, 0);
738 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
740 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
741 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
742 value_clear(e
->x
.p
->arr
[2*i
].d
);
744 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
745 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
750 for (j
= 0; j
< s
.n
; ++j
) {
751 value_clear(s
.fixed
[j
].d
);
752 value_clear(s
.fixed
[j
].m
);
753 free_evalue_refs(&s
.fixed
[j
].s
);
757 if (e
->x
.p
->size
== 0) {
759 evalue_set_si(e
, 0, 1);
762 _reduce_evalue(e
, &s
, 0);
767 void print_evalue(FILE *DST
,evalue
*e
,char **pname
) {
769 if(value_notzero_p(e
->d
)) {
770 if(value_notone_p(e
->d
)) {
771 value_print(DST
,VALUE_FMT
,e
->x
.n
);
773 value_print(DST
,VALUE_FMT
,e
->d
);
776 value_print(DST
,VALUE_FMT
,e
->x
.n
);
780 print_enode(DST
,e
->x
.p
,pname
);
784 void print_enode(FILE *DST
,enode
*p
,char **pname
) {
789 fprintf(DST
, "NULL");
795 for (i
=0; i
<p
->size
; i
++) {
796 print_evalue(DST
, &p
->arr
[i
], pname
);
800 fprintf(DST
, " }\n");
804 for (i
=p
->size
-1; i
>=0; i
--) {
805 print_evalue(DST
, &p
->arr
[i
], pname
);
806 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
808 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
810 fprintf(DST
, " )\n");
814 for (i
=0; i
<p
->size
; i
++) {
815 print_evalue(DST
, &p
->arr
[i
], pname
);
816 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
818 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
823 for (i
=p
->size
-1; i
>=1; i
--) {
824 print_evalue(DST
, &p
->arr
[i
], pname
);
827 fprintf(DST
, p
->type
== flooring
? "[" : "{");
828 print_evalue(DST
, &p
->arr
[0], pname
);
829 fprintf(DST
, p
->type
== flooring
? "]" : "}");
831 fprintf(DST
, "^%d + ", i
-1);
836 fprintf(DST
, " )\n");
840 print_evalue(DST
, &p
->arr
[0], pname
);
841 fprintf(DST
, "= 0 ] * \n");
842 print_evalue(DST
, &p
->arr
[1], pname
);
844 fprintf(DST
, " +\n [ ");
845 print_evalue(DST
, &p
->arr
[0], pname
);
846 fprintf(DST
, "!= 0 ] * \n");
847 print_evalue(DST
, &p
->arr
[2], pname
);
851 char **names
= pname
;
852 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
853 if (p
->pos
< maxdim
) {
854 NALLOC(names
, maxdim
);
855 for (i
= 0; i
< p
->pos
; ++i
)
857 for ( ; i
< maxdim
; ++i
) {
858 NALLOC(names
[i
], 10);
859 snprintf(names
[i
], 10, "_p%d", i
);
863 for (i
=0; i
<p
->size
/2; i
++) {
864 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
865 print_evalue(DST
, &p
->arr
[2*i
+1], names
);
868 if (p
->pos
< maxdim
) {
869 for (i
= p
->pos
; i
< maxdim
; ++i
)
882 static void eadd_rev(evalue
*e1
, evalue
*res
)
886 evalue_copy(&ev
, e1
);
888 free_evalue_refs(res
);
892 static void eadd_rev_cst (evalue
*e1
, evalue
*res
)
896 evalue_copy(&ev
, e1
);
897 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
898 free_evalue_refs(res
);
902 struct section
{ Polyhedron
* D
; evalue E
; };
904 void eadd_partitions (evalue
*e1
,evalue
*res
)
909 s
= (struct section
*)
910 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
911 sizeof(struct section
));
913 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
914 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
917 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
918 assert(res
->x
.p
->size
>= 2);
919 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
920 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
922 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
924 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
933 value_init(s
[n
].E
.d
);
934 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
938 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
939 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
940 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
942 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
943 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
949 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
950 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
952 value_init(s
[n
].E
.d
);
953 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
954 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
959 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
963 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
966 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
967 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
968 value_clear(res
->x
.p
->arr
[2*i
].d
);
973 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
974 for (j
= 0; j
< n
; ++j
) {
975 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
976 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
977 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
978 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
984 static void explicit_complement(evalue
*res
)
986 enode
*rel
= new_enode(relation
, 3, 0);
988 value_clear(rel
->arr
[0].d
);
989 rel
->arr
[0] = res
->x
.p
->arr
[0];
990 value_clear(rel
->arr
[1].d
);
991 rel
->arr
[1] = res
->x
.p
->arr
[1];
992 value_set_si(rel
->arr
[2].d
, 1);
993 value_init(rel
->arr
[2].x
.n
);
994 value_set_si(rel
->arr
[2].x
.n
, 0);
999 void eadd(evalue
*e1
,evalue
*res
) {
1002 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1003 /* Add two rational numbers */
1009 value_multiply(m1
,e1
->x
.n
,res
->d
);
1010 value_multiply(m2
,res
->x
.n
,e1
->d
);
1011 value_addto(res
->x
.n
,m1
,m2
);
1012 value_multiply(res
->d
,e1
->d
,res
->d
);
1013 Gcd(res
->x
.n
,res
->d
,&g
);
1014 if (value_notone_p(g
)) {
1015 value_division(res
->d
,res
->d
,g
);
1016 value_division(res
->x
.n
,res
->x
.n
,g
);
1018 value_clear(g
); value_clear(m1
); value_clear(m2
);
1021 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
1022 switch (res
->x
.p
->type
) {
1024 /* Add the constant to the constant term of a polynomial*/
1025 eadd(e1
, &res
->x
.p
->arr
[0]);
1028 /* Add the constant to all elements of a periodic number */
1029 for (i
=0; i
<res
->x
.p
->size
; i
++) {
1030 eadd(e1
, &res
->x
.p
->arr
[i
]);
1034 fprintf(stderr
, "eadd: cannot add const with vector\n");
1038 eadd(e1
, &res
->x
.p
->arr
[1]);
1041 assert(EVALUE_IS_ZERO(*e1
));
1042 break; /* Do nothing */
1044 /* Create (zero) complement if needed */
1045 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
1046 explicit_complement(res
);
1047 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1048 eadd(e1
, &res
->x
.p
->arr
[i
]);
1054 /* add polynomial or periodic to constant
1055 * you have to exchange e1 and res, before doing addition */
1057 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1061 else { // ((e1->d==0) && (res->d==0))
1062 assert(!((e1
->x
.p
->type
== partition
) ^
1063 (res
->x
.p
->type
== partition
)));
1064 if (e1
->x
.p
->type
== partition
) {
1065 eadd_partitions(e1
, res
);
1068 if (e1
->x
.p
->type
== relation
&&
1069 (res
->x
.p
->type
!= relation
||
1070 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1074 if (res
->x
.p
->type
== relation
) {
1075 if (e1
->x
.p
->type
== relation
&&
1076 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1077 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1078 explicit_complement(res
);
1079 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1080 explicit_complement(e1
);
1081 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1082 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1085 if (res
->x
.p
->size
< 3)
1086 explicit_complement(res
);
1087 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1088 eadd(e1
, &res
->x
.p
->arr
[i
]);
1091 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
1092 /* adding to evalues of different type. two cases are possible
1093 * res is periodic and e1 is polynomial, you have to exchange
1094 * e1 and res then to add e1 to the constant term of res */
1095 if (e1
->x
.p
->type
== polynomial
) {
1096 eadd_rev_cst(e1
, res
);
1098 else if (res
->x
.p
->type
== polynomial
) {
1099 /* res is polynomial and e1 is periodic,
1100 add e1 to the constant term of res */
1102 eadd(e1
,&res
->x
.p
->arr
[0]);
1108 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
1109 ((res
->x
.p
->type
== fractional
||
1110 res
->x
.p
->type
== flooring
) &&
1111 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1112 /* adding evalues of different position (i.e function of different unknowns
1113 * to case are possible */
1115 switch (res
->x
.p
->type
) {
1118 if(mod_term_smaller(res
, e1
))
1119 eadd(e1
,&res
->x
.p
->arr
[1]);
1121 eadd_rev_cst(e1
, res
);
1123 case polynomial
: // res and e1 are polynomials
1124 // add e1 to the constant term of res
1126 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1127 eadd(e1
,&res
->x
.p
->arr
[0]);
1129 eadd_rev_cst(e1
, res
);
1130 // value_clear(g); value_clear(m1); value_clear(m2);
1132 case periodic
: // res and e1 are pointers to periodic numbers
1133 //add e1 to all elements of res
1135 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1136 for (i
=0;i
<res
->x
.p
->size
;i
++) {
1137 eadd(e1
,&res
->x
.p
->arr
[i
]);
1148 //same type , same pos and same size
1149 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1150 // add any element in e1 to the corresponding element in res
1151 i
= type_offset(res
->x
.p
);
1153 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1154 for (; i
<res
->x
.p
->size
; i
++) {
1155 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1160 /* Sizes are different */
1161 switch(res
->x
.p
->type
) {
1165 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1166 /* new enode and add res to that new node. If you do not do */
1167 /* that, you lose the the upper weight part of e1 ! */
1169 if(e1
->x
.p
->size
> res
->x
.p
->size
)
1172 i
= type_offset(res
->x
.p
);
1174 assert(eequal(&e1
->x
.p
->arr
[0],
1175 &res
->x
.p
->arr
[0]));
1176 for (; i
<e1
->x
.p
->size
; i
++) {
1177 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1184 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1187 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1188 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1189 to add periodicaly elements of e1 to elements of ne, and finaly to
1194 value_init(ex
); value_init(ey
);value_init(ep
);
1197 value_set_si(ex
,e1
->x
.p
->size
);
1198 value_set_si(ey
,res
->x
.p
->size
);
1199 value_assign (ep
,*Lcm(ex
,ey
));
1200 p
=(int)mpz_get_si(ep
);
1201 ne
= (evalue
*) malloc (sizeof(evalue
));
1203 value_set_si( ne
->d
,0);
1205 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1207 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1210 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1213 value_assign(res
->d
, ne
->d
);
1219 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
1228 static void emul_rev (evalue
*e1
, evalue
*res
)
1232 evalue_copy(&ev
, e1
);
1234 free_evalue_refs(res
);
1238 static void emul_poly (evalue
*e1
, evalue
*res
)
1240 int i
, j
, o
= type_offset(res
->x
.p
);
1242 int size
=(e1
->x
.p
->size
+ res
->x
.p
->size
- o
- 1);
1244 value_set_si(tmp
.d
,0);
1245 tmp
.x
.p
=new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1247 evalue_copy(&tmp
.x
.p
->arr
[0], &e1
->x
.p
->arr
[0]);
1248 for (i
=o
; i
< e1
->x
.p
->size
; i
++) {
1249 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1250 emul(&res
->x
.p
->arr
[o
], &tmp
.x
.p
->arr
[i
]);
1253 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1254 for (i
=o
+1; i
<res
->x
.p
->size
; i
++)
1255 for (j
=o
; j
<e1
->x
.p
->size
; j
++) {
1258 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1259 emul(&res
->x
.p
->arr
[i
], &ev
);
1260 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-o
]);
1261 free_evalue_refs(&ev
);
1263 free_evalue_refs(res
);
1267 void emul_partitions (evalue
*e1
,evalue
*res
)
1272 s
= (struct section
*)
1273 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1274 sizeof(struct section
));
1276 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1277 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1280 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1281 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1282 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1283 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1289 /* This code is only needed because the partitions
1290 are not true partitions.
1292 for (k
= 0; k
< n
; ++k
) {
1293 if (DomainIncludes(s
[k
].D
, d
))
1295 if (DomainIncludes(d
, s
[k
].D
)) {
1296 Domain_Free(s
[k
].D
);
1297 free_evalue_refs(&s
[k
].E
);
1308 value_init(s
[n
].E
.d
);
1309 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1310 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1314 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1315 value_clear(res
->x
.p
->arr
[2*i
].d
);
1316 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1321 evalue_set_si(res
, 0, 1);
1323 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1324 for (j
= 0; j
< n
; ++j
) {
1325 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1326 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1327 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1328 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1335 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1337 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1338 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1339 * evalues is not treated here */
1341 void emul (evalue
*e1
, evalue
*res
){
1344 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1345 fprintf(stderr
, "emul: do not proced on evector type !\n");
1349 if (EVALUE_IS_ZERO(*res
))
1352 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1353 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1354 emul_partitions(e1
, res
);
1357 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1358 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1359 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1361 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1362 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1363 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1364 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1365 explicit_complement(res
);
1366 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1367 explicit_complement(e1
);
1368 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1369 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1372 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1373 emul(e1
, &res
->x
.p
->arr
[i
]);
1375 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1376 switch(e1
->x
.p
->type
) {
1378 switch(res
->x
.p
->type
) {
1380 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1381 /* Product of two polynomials of the same variable */
1386 /* Product of two polynomials of different variables */
1388 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1389 for( i
=0; i
<res
->x
.p
->size
; i
++)
1390 emul(e1
, &res
->x
.p
->arr
[i
]);
1399 /* Product of a polynomial and a periodic or fractional */
1406 switch(res
->x
.p
->type
) {
1408 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1409 /* Product of two periodics of the same parameter and period */
1411 for(i
=0; i
<res
->x
.p
->size
;i
++)
1412 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1417 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1418 /* Product of two periodics of the same parameter and different periods */
1422 value_init(x
); value_init(y
);value_init(z
);
1425 value_set_si(x
,e1
->x
.p
->size
);
1426 value_set_si(y
,res
->x
.p
->size
);
1427 value_assign (z
,*Lcm(x
,y
));
1428 lcm
=(int)mpz_get_si(z
);
1429 newp
= (evalue
*) malloc (sizeof(evalue
));
1430 value_init(newp
->d
);
1431 value_set_si( newp
->d
,0);
1432 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1433 for(i
=0;i
<lcm
;i
++) {
1434 evalue_copy(&newp
->x
.p
->arr
[i
],
1435 &res
->x
.p
->arr
[i
%iy
]);
1438 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1440 value_assign(res
->d
,newp
->d
);
1443 value_clear(x
); value_clear(y
);value_clear(z
);
1447 /* Product of two periodics of different parameters */
1449 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1450 for(i
=0; i
<res
->x
.p
->size
; i
++)
1451 emul(e1
, &(res
->x
.p
->arr
[i
]));
1459 /* Product of a periodic and a polynomial */
1461 for(i
=0; i
<res
->x
.p
->size
; i
++)
1462 emul(e1
, &(res
->x
.p
->arr
[i
]));
1469 switch(res
->x
.p
->type
) {
1471 for(i
=0; i
<res
->x
.p
->size
; i
++)
1472 emul(e1
, &(res
->x
.p
->arr
[i
]));
1479 assert(e1
->x
.p
->type
== res
->x
.p
->type
);
1480 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1481 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1484 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1485 if (e1
->x
.p
->type
!= fractional
|| !value_two_p(d
.d
))
1489 value_set_si(d
.x
.n
, 1);
1491 /* { x }^2 == { x }/2 */
1492 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1493 assert(e1
->x
.p
->size
== 3);
1494 assert(res
->x
.p
->size
== 3);
1496 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1498 eadd(&res
->x
.p
->arr
[1], &tmp
);
1499 emul(&e1
->x
.p
->arr
[2], &tmp
);
1500 emul(&e1
->x
.p
->arr
[1], res
);
1501 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1502 free_evalue_refs(&tmp
);
1507 if(mod_term_smaller(res
, e1
))
1508 for(i
=1; i
<res
->x
.p
->size
; i
++)
1509 emul(e1
, &(res
->x
.p
->arr
[i
]));
1524 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1525 /* Product of two rational numbers */
1529 value_multiply(res
->d
,e1
->d
,res
->d
);
1530 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1531 Gcd(res
->x
.n
, res
->d
,&g
);
1532 if (value_notone_p(g
)) {
1533 value_division(res
->d
,res
->d
,g
);
1534 value_division(res
->x
.n
,res
->x
.n
,g
);
1540 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1541 /* Product of an expression (polynomial or peririodic) and a rational number */
1547 /* Product of a rationel number and an expression (polynomial or peririodic) */
1549 i
= type_offset(res
->x
.p
);
1550 for (; i
<res
->x
.p
->size
; i
++)
1551 emul(e1
, &res
->x
.p
->arr
[i
]);
1561 /* Frees mask content ! */
1562 void emask(evalue
*mask
, evalue
*res
) {
1569 if (EVALUE_IS_ZERO(*res
)) {
1570 free_evalue_refs(mask
);
1574 assert(value_zero_p(mask
->d
));
1575 assert(mask
->x
.p
->type
== partition
);
1576 assert(value_zero_p(res
->d
));
1577 assert(res
->x
.p
->type
== partition
);
1578 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1579 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1580 pos
= res
->x
.p
->pos
;
1582 s
= (struct section
*)
1583 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1584 sizeof(struct section
));
1588 evalue_set_si(&mone
, -1, 1);
1591 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1592 assert(mask
->x
.p
->size
>= 2);
1593 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1594 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1596 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1598 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1607 value_init(s
[n
].E
.d
);
1608 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1612 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1613 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1616 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1617 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1618 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1619 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1621 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1622 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1628 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1629 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1631 value_init(s
[n
].E
.d
);
1632 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1633 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1639 /* Just ignore; this may have been previously masked off */
1641 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1645 free_evalue_refs(&mone
);
1646 free_evalue_refs(mask
);
1647 free_evalue_refs(res
);
1650 evalue_set_si(res
, 0, 1);
1652 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1653 for (j
= 0; j
< n
; ++j
) {
1654 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1655 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1656 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1663 void evalue_copy(evalue
*dst
, evalue
*src
)
1665 value_assign(dst
->d
, src
->d
);
1666 if(value_notzero_p(src
->d
)) {
1667 value_init(dst
->x
.n
);
1668 value_assign(dst
->x
.n
, src
->x
.n
);
1670 dst
->x
.p
= ecopy(src
->x
.p
);
1673 enode
*new_enode(enode_type type
,int size
,int pos
) {
1679 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1682 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1686 for(i
=0; i
<size
; i
++) {
1687 value_init(res
->arr
[i
].d
);
1688 value_set_si(res
->arr
[i
].d
,0);
1689 res
->arr
[i
].x
.p
= 0;
1694 enode
*ecopy(enode
*e
) {
1699 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1700 for(i
=0;i
<e
->size
;++i
) {
1701 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1702 if(value_zero_p(res
->arr
[i
].d
))
1703 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1704 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1705 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1707 value_init(res
->arr
[i
].x
.n
);
1708 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1714 int ecmp(const evalue
*e1
, const evalue
*e2
)
1720 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
1724 value_multiply(m
, e1
->x
.n
, e2
->d
);
1725 value_multiply(m2
, e2
->x
.n
, e1
->d
);
1727 if (value_lt(m
, m2
))
1729 else if (value_gt(m
, m2
))
1739 if (value_notzero_p(e1
->d
))
1741 if (value_notzero_p(e2
->d
))
1747 if (p1
->type
!= p2
->type
)
1748 return p1
->type
- p2
->type
;
1749 if (p1
->pos
!= p2
->pos
)
1750 return p1
->pos
- p2
->pos
;
1751 if (p1
->size
!= p2
->size
)
1752 return p1
->size
- p2
->size
;
1754 for (i
= p1
->size
-1; i
>= 0; --i
)
1755 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
1761 int eequal(evalue
*e1
,evalue
*e2
) {
1766 if (value_ne(e1
->d
,e2
->d
))
1769 /* e1->d == e2->d */
1770 if (value_notzero_p(e1
->d
)) {
1771 if (value_ne(e1
->x
.n
,e2
->x
.n
))
1774 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1778 /* e1->d == e2->d == 0 */
1781 if (p1
->type
!= p2
->type
) return 0;
1782 if (p1
->size
!= p2
->size
) return 0;
1783 if (p1
->pos
!= p2
->pos
) return 0;
1784 for (i
=0; i
<p1
->size
; i
++)
1785 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
1790 void free_evalue_refs(evalue
*e
) {
1795 if (EVALUE_IS_DOMAIN(*e
)) {
1796 Domain_Free(EVALUE_DOMAIN(*e
));
1799 } else if (value_pos_p(e
->d
)) {
1801 /* 'e' stores a constant */
1803 value_clear(e
->x
.n
);
1806 assert(value_zero_p(e
->d
));
1809 if (!p
) return; /* null pointer */
1810 for (i
=0; i
<p
->size
; i
++) {
1811 free_evalue_refs(&(p
->arr
[i
]));
1815 } /* free_evalue_refs */
1817 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
1818 Vector
* val
, evalue
*res
)
1820 unsigned nparam
= periods
->Size
;
1823 double d
= compute_evalue(e
, val
->p
);
1824 d
*= VALUE_TO_DOUBLE(m
);
1829 value_assign(res
->d
, m
);
1830 value_init(res
->x
.n
);
1831 value_set_double(res
->x
.n
, d
);
1832 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
1835 if (value_one_p(periods
->p
[p
]))
1836 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
1841 value_assign(tmp
, periods
->p
[p
]);
1842 value_set_si(res
->d
, 0);
1843 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
1845 value_decrement(tmp
, tmp
);
1846 value_assign(val
->p
[p
], tmp
);
1847 mod2table_r(e
, periods
, m
, p
+1, val
,
1848 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
1849 } while (value_pos_p(tmp
));
1855 static void rel2table(evalue
*e
, int zero
)
1857 if (value_pos_p(e
->d
)) {
1858 if (value_zero_p(e
->x
.n
) == zero
)
1859 value_set_si(e
->x
.n
, 1);
1861 value_set_si(e
->x
.n
, 0);
1862 value_set_si(e
->d
, 1);
1865 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
1866 rel2table(&e
->x
.p
->arr
[i
], zero
);
1870 void evalue_mod2table(evalue
*e
, int nparam
)
1875 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
1878 for (i
=0; i
<p
->size
; i
++) {
1879 evalue_mod2table(&(p
->arr
[i
]), nparam
);
1881 if (p
->type
== relation
) {
1886 evalue_copy(©
, &p
->arr
[0]);
1888 rel2table(&p
->arr
[0], 1);
1889 emul(&p
->arr
[0], &p
->arr
[1]);
1891 rel2table(©
, 0);
1892 emul(©
, &p
->arr
[2]);
1893 eadd(&p
->arr
[2], &p
->arr
[1]);
1894 free_evalue_refs(&p
->arr
[2]);
1895 free_evalue_refs(©
);
1897 free_evalue_refs(&p
->arr
[0]);
1901 } else if (p
->type
== fractional
) {
1902 Vector
*periods
= Vector_Alloc(nparam
);
1903 Vector
*val
= Vector_Alloc(nparam
);
1909 value_set_si(tmp
, 1);
1910 Vector_Set(periods
->p
, 1, nparam
);
1911 Vector_Set(val
->p
, 0, nparam
);
1912 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
1915 assert(p
->type
== polynomial
);
1916 assert(p
->size
== 2);
1917 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
1918 value_lcm(tmp
, p
->arr
[1].d
, &tmp
);
1921 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
1924 evalue_set_si(&res
, 0, 1);
1925 /* Compute the polynomial using Horner's rule */
1926 for (i
=p
->size
-1;i
>1;i
--) {
1927 eadd(&p
->arr
[i
], &res
);
1930 eadd(&p
->arr
[1], &res
);
1932 free_evalue_refs(e
);
1933 free_evalue_refs(&EP
);
1938 Vector_Free(periods
);
1940 } /* evalue_mod2table */
1942 /********************************************************/
1943 /* function in domain */
1944 /* check if the parameters in list_args */
1945 /* verifies the constraints of Domain P */
1946 /********************************************************/
1947 int in_domain(Polyhedron
*P
, Value
*list_args
) {
1950 Value v
; /* value of the constraint of a row when
1951 parameters are instanciated*/
1957 /*P->Constraint constraint matrice of polyhedron P*/
1958 for(row
=0;row
<P
->NbConstraints
;row
++) {
1959 value_assign(v
,P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
1960 for(col
=1;col
<P
->Dimension
+1;col
++) {
1961 value_multiply(tmp
,P
->Constraint
[row
][col
],list_args
[col
-1]);
1962 value_addto(v
,v
,tmp
);
1964 if (value_notzero_p(P
->Constraint
[row
][0])) {
1966 /*if v is not >=0 then this constraint is not respected */
1967 if (value_neg_p(v
)) {
1971 return P
->next
? in_domain(P
->next
, list_args
) : 0;
1976 /*if v is not = 0 then this constraint is not respected */
1977 if (value_notzero_p(v
))
1982 /*if not return before this point => all
1983 the constraints are respected */
1989 /****************************************************/
1990 /* function compute enode */
1991 /* compute the value of enode p with parameters */
1992 /* list "list_args */
1993 /* compute the polynomial or the periodic */
1994 /****************************************************/
1996 static double compute_enode(enode
*p
, Value
*list_args
) {
2008 if (p
->type
== polynomial
) {
2010 value_assign(param
,list_args
[p
->pos
-1]);
2012 /* Compute the polynomial using Horner's rule */
2013 for (i
=p
->size
-1;i
>0;i
--) {
2014 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2015 res
*=VALUE_TO_DOUBLE(param
);
2017 res
+=compute_evalue(&p
->arr
[0],list_args
);
2019 else if (p
->type
== fractional
) {
2020 double d
= compute_evalue(&p
->arr
[0], list_args
);
2021 d
-= floor(d
+1e-10);
2023 /* Compute the polynomial using Horner's rule */
2024 for (i
=p
->size
-1;i
>1;i
--) {
2025 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2028 res
+=compute_evalue(&p
->arr
[1],list_args
);
2030 else if (p
->type
== flooring
) {
2031 double d
= compute_evalue(&p
->arr
[0], list_args
);
2034 /* Compute the polynomial using Horner's rule */
2035 for (i
=p
->size
-1;i
>1;i
--) {
2036 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2039 res
+=compute_evalue(&p
->arr
[1],list_args
);
2041 else if (p
->type
== periodic
) {
2042 value_assign(param
,list_args
[p
->pos
-1]);
2044 /* Choose the right element of the periodic */
2045 value_absolute(m
,param
);
2046 value_set_si(param
,p
->size
);
2047 value_modulus(m
,m
,param
);
2048 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2050 else if (p
->type
== relation
) {
2051 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2052 res
= compute_evalue(&p
->arr
[1], list_args
);
2053 else if (p
->size
> 2)
2054 res
= compute_evalue(&p
->arr
[2], list_args
);
2056 else if (p
->type
== partition
) {
2057 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2058 Value
*vals
= list_args
;
2061 for (i
= 0; i
< dim
; ++i
) {
2062 value_init(vals
[i
]);
2064 value_assign(vals
[i
], list_args
[i
]);
2067 for (i
= 0; i
< p
->size
/2; ++i
)
2068 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2069 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2073 for (i
= 0; i
< dim
; ++i
)
2074 value_clear(vals
[i
]);
2083 } /* compute_enode */
2085 /*************************************************/
2086 /* return the value of Ehrhart Polynomial */
2087 /* It returns a double, because since it is */
2088 /* a recursive function, some intermediate value */
2089 /* might not be integral */
2090 /*************************************************/
2092 double compute_evalue(evalue
*e
,Value
*list_args
) {
2096 if (value_notzero_p(e
->d
)) {
2097 if (value_notone_p(e
->d
))
2098 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2100 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2103 res
= compute_enode(e
->x
.p
,list_args
);
2105 } /* compute_evalue */
2108 /****************************************************/
2109 /* function compute_poly : */
2110 /* Check for the good validity domain */
2111 /* return the number of point in the Polyhedron */
2112 /* in allocated memory */
2113 /* Using the Ehrhart pseudo-polynomial */
2114 /****************************************************/
2115 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2118 /* double d; int i; */
2120 tmp
= (Value
*) malloc (sizeof(Value
));
2121 assert(tmp
!= NULL
);
2123 value_set_si(*tmp
,0);
2126 return(tmp
); /* no ehrhart polynomial */
2127 if(en
->ValidityDomain
) {
2128 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2129 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2134 return(tmp
); /* no Validity Domain */
2136 if(in_domain(en
->ValidityDomain
,list_args
)) {
2138 #ifdef EVAL_EHRHART_DEBUG
2139 Print_Domain(stdout
,en
->ValidityDomain
);
2140 print_evalue(stdout
,&en
->EP
);
2143 /* d = compute_evalue(&en->EP,list_args);
2145 printf("(double)%lf = %d\n", d, i ); */
2146 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2152 value_set_si(*tmp
,0);
2153 return(tmp
); /* no compatible domain with the arguments */
2154 } /* compute_poly */
2156 size_t value_size(Value v
) {
2157 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2158 * sizeof(v
[0]._mp_d
[0]);
2161 size_t domain_size(Polyhedron
*D
)
2164 size_t s
= sizeof(*D
);
2166 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2167 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2168 s
+= value_size(D
->Constraint
[i
][j
]);
2171 for (i = 0; i < D->NbRays; ++i)
2172 for (j = 0; j < D->Dimension+2; ++j)
2173 s += value_size(D->Ray[i][j]);
2176 return D
->next
? s
+domain_size(D
->next
) : s
;
2179 size_t enode_size(enode
*p
) {
2180 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2183 if (p
->type
== partition
)
2184 for (i
= 0; i
< p
->size
/2; ++i
) {
2185 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2186 s
+= evalue_size(&p
->arr
[2*i
+1]);
2189 for (i
= 0; i
< p
->size
; ++i
) {
2190 s
+= evalue_size(&p
->arr
[i
]);
2195 size_t evalue_size(evalue
*e
)
2197 size_t s
= sizeof(*e
);
2198 s
+= value_size(e
->d
);
2199 if (value_notzero_p(e
->d
))
2200 s
+= value_size(e
->x
.n
);
2202 s
+= enode_size(e
->x
.p
);
2206 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2208 evalue
*found
= NULL
;
2213 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2216 value_init(offset
.d
);
2217 value_init(offset
.x
.n
);
2218 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2219 value_lcm(m
, offset
.d
, &offset
.d
);
2220 value_set_si(offset
.x
.n
, 1);
2223 evalue_copy(©
, cst
);
2226 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2228 if (eequal(base
, &e
->x
.p
->arr
[0]))
2229 found
= &e
->x
.p
->arr
[0];
2231 value_set_si(offset
.x
.n
, -2);
2234 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2236 if (eequal(base
, &e
->x
.p
->arr
[0]))
2239 free_evalue_refs(cst
);
2240 free_evalue_refs(&offset
);
2243 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2244 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2249 static evalue
*find_relation_pair(evalue
*e
)
2252 evalue
*found
= NULL
;
2254 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2257 if (e
->x
.p
->type
== fractional
) {
2262 poly_denom(&e
->x
.p
->arr
[0], &m
);
2264 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2265 cst
= &cst
->x
.p
->arr
[0])
2268 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2269 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2274 i
= e
->x
.p
->type
== relation
;
2275 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2276 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2281 void evalue_mod2relation(evalue
*e
) {
2284 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2287 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2288 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2289 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2290 value_clear(e
->x
.p
->arr
[2*i
].d
);
2291 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2293 if (2*i
< e
->x
.p
->size
) {
2294 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2295 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2300 if (e
->x
.p
->size
== 0) {
2302 evalue_set_si(e
, 0, 1);
2308 while ((d
= find_relation_pair(e
)) != NULL
) {
2312 value_init(split
.d
);
2313 value_set_si(split
.d
, 0);
2314 split
.x
.p
= new_enode(relation
, 3, 0);
2315 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2316 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2318 ev
= &split
.x
.p
->arr
[0];
2319 value_set_si(ev
->d
, 0);
2320 ev
->x
.p
= new_enode(fractional
, 3, -1);
2321 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2322 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2323 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2329 free_evalue_refs(&split
);
2333 static int evalue_comp(const void * a
, const void * b
)
2335 const evalue
*e1
= *(const evalue
**)a
;
2336 const evalue
*e2
= *(const evalue
**)b
;
2337 return ecmp(e1
, e2
);
2340 void evalue_combine(evalue
*e
)
2347 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2350 NALLOC(evs
, e
->x
.p
->size
/2);
2351 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2352 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2353 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2354 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2355 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2356 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2357 value_clear(p
->arr
[2*k
].d
);
2358 value_clear(p
->arr
[2*k
+1].d
);
2359 p
->arr
[2*k
] = *(evs
[i
]-1);
2360 p
->arr
[2*k
+1] = *(evs
[i
]);
2363 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2366 value_clear((evs
[i
]-1)->d
);
2370 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2371 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2372 free_evalue_refs(evs
[i
]);
2376 for (i
= 2*k
; i
< p
->size
; ++i
)
2377 value_clear(p
->arr
[i
].d
);
2384 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2386 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2388 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2391 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2392 Polyhedron
*D
, *N
, **P
;
2395 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2402 if (D
->NbEq
<= H
->NbEq
) {
2408 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2409 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2410 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2411 reduce_evalue(&tmp
);
2412 if (value_notzero_p(tmp
.d
) ||
2413 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2416 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2417 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2420 free_evalue_refs(&tmp
);
2426 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2428 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2430 value_clear(e
->x
.p
->arr
[2*i
].d
);
2431 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2433 if (2*i
< e
->x
.p
->size
) {
2434 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2435 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2442 H
= DomainConvex(D
, 0);
2443 E
= DomainDifference(H
, D
, 0);
2445 D
= DomainDifference(H
, E
, 0);
2448 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2452 /* May change coefficients to become non-standard if fiddle is set
2453 * => reduce p afterwards to correct
2455 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2456 Matrix
**R
, int fiddle
)
2460 unsigned dim
= D
->Dimension
;
2461 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2466 assert(p
->type
== fractional
);
2468 value_set_si(T
->p
[1][dim
], 1);
2470 while (value_zero_p(pp
->d
)) {
2471 assert(pp
->x
.p
->type
== polynomial
);
2472 assert(pp
->x
.p
->size
== 2);
2473 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2474 value_division(T
->p
[0][pp
->x
.p
->pos
-1], *d
, pp
->x
.p
->arr
[1].d
);
2475 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2476 if (fiddle
&& value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2477 value_substract(pp
->x
.p
->arr
[1].x
.n
,
2478 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2479 value_multiply(T
->p
[0][pp
->x
.p
->pos
-1],
2480 T
->p
[0][pp
->x
.p
->pos
-1], pp
->x
.p
->arr
[1].x
.n
);
2481 pp
= &pp
->x
.p
->arr
[0];
2483 value_division(T
->p
[0][dim
], *d
, pp
->d
);
2484 value_multiply(T
->p
[0][dim
], T
->p
[0][dim
], pp
->x
.n
);
2485 I
= DomainImage(D
, T
, 0);
2486 H
= DomainConvex(I
, 0);
2495 static int reduce_in_domain(evalue
*e
, Polyhedron
*D
)
2505 if (value_notzero_p(e
->d
))
2510 if (p
->type
== relation
) {
2516 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
, 1);
2517 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2518 equal
= value_eq(min
, max
);
2519 mpz_cdiv_q(min
, min
, d
);
2520 mpz_fdiv_q(max
, max
, d
);
2522 if (bounded
&& value_gt(min
, max
)) {
2528 evalue_set_si(e
, 0, 1);
2531 free_evalue_refs(&(p
->arr
[1]));
2532 free_evalue_refs(&(p
->arr
[0]));
2538 return r
? r
: reduce_in_domain(e
, D
);
2539 } else if (bounded
&& equal
) {
2542 free_evalue_refs(&(p
->arr
[2]));
2545 free_evalue_refs(&(p
->arr
[0]));
2551 return reduce_in_domain(e
, D
);
2552 } else if (bounded
&& value_eq(min
, max
)) {
2553 /* zero for a single value */
2555 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2556 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2557 value_multiply(min
, min
, d
);
2558 value_substract(M
->p
[0][D
->Dimension
+1],
2559 M
->p
[0][D
->Dimension
+1], min
);
2560 E
= DomainAddConstraints(D
, M
, 0);
2566 r
= reduce_in_domain(&p
->arr
[1], E
);
2568 r
|= reduce_in_domain(&p
->arr
[2], D
);
2570 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2578 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2581 i
= p
->type
== relation
? 1 :
2582 p
->type
== fractional
? 1 : 0;
2583 for (; i
<p
->size
; i
++)
2584 r
|= reduce_in_domain(&p
->arr
[i
], D
);
2586 if (p
->type
!= fractional
) {
2587 if (r
&& p
->type
== polynomial
) {
2590 value_set_si(f
.d
, 0);
2591 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2592 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2593 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2594 reorder_terms(p
, &f
);
2605 I
= polynomial_projection(p
, D
, &d
, &T
, 1);
2607 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2608 mpz_fdiv_q(min
, min
, d
);
2609 mpz_fdiv_q(max
, max
, d
);
2610 value_substract(d
, max
, min
);
2612 if (bounded
&& value_eq(min
, max
)) {
2615 value_init(inc
.x
.n
);
2616 value_set_si(inc
.d
, 1);
2617 value_oppose(inc
.x
.n
, min
);
2618 eadd(&inc
, &p
->arr
[0]);
2619 reorder_terms(p
, &p
->arr
[0]); /* frees arr[0] */
2623 free_evalue_refs(&inc
);
2625 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
2628 value_set_si(rem
.d
, 0);
2629 rem
.x
.p
= new_enode(fractional
, 3, -1);
2630 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
2631 rem
.x
.p
->arr
[1] = p
->arr
[1];
2632 rem
.x
.p
->arr
[2] = p
->arr
[2];
2633 for (i
= 3; i
< p
->size
; ++i
)
2634 p
->arr
[i
-2] = p
->arr
[i
];
2639 value_init(inc
.x
.n
);
2640 value_set_si(inc
.d
, 1);
2641 value_oppose(inc
.x
.n
, min
);
2645 evalue_copy(&t
, &p
->arr
[0]);
2650 value_set_si(f
.d
, 0);
2651 f
.x
.p
= new_enode(fractional
, 3, -1);
2652 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
2653 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2654 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
2657 value_init(factor
.d
);
2658 evalue_set_si(&factor
, -1, 1);
2664 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2665 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
2671 free_evalue_refs(&inc
);
2672 free_evalue_refs(&t
);
2673 free_evalue_refs(&f
);
2674 free_evalue_refs(&factor
);
2675 free_evalue_refs(&rem
);
2677 reduce_in_domain(e
, D
);
2681 _reduce_evalue(&p
->arr
[0], 0, 1);
2685 value_set_si(f
.d
, 0);
2686 f
.x
.p
= new_enode(fractional
, 3, -1);
2687 value_clear(f
.x
.p
->arr
[0].d
);
2688 f
.x
.p
->arr
[0] = p
->arr
[0];
2689 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2690 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
2691 reorder_terms(p
, &f
);
2705 void evalue_range_reduction(evalue
*e
)
2708 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2711 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2712 if (reduce_in_domain(&e
->x
.p
->arr
[2*i
+1],
2713 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
2714 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
2716 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2717 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2718 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
2719 value_clear(e
->x
.p
->arr
[2*i
].d
);
2721 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2722 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2730 Enumeration
* partition2enumeration(evalue
*EP
)
2733 Enumeration
*en
, *res
= NULL
;
2735 if (EVALUE_IS_ZERO(*EP
)) {
2740 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
2741 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
2742 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
2745 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
2746 value_clear(EP
->x
.p
->arr
[2*i
].d
);
2747 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
2755 static int frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
2765 if (value_notzero_p(e
->d
))
2770 i
= p
->type
== relation
? 1 :
2771 p
->type
== fractional
? 1 : 0;
2772 for (; i
<p
->size
; i
++)
2773 r
|= frac2floor_in_domain(&p
->arr
[i
], D
);
2775 if (p
->type
!= fractional
) {
2776 if (r
&& p
->type
== polynomial
) {
2779 value_set_si(f
.d
, 0);
2780 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2781 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2782 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2783 reorder_terms(p
, &f
);
2792 I
= polynomial_projection(p
, D
, &d
, &T
, 0);
2795 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2798 assert(I
->NbEq
== 0); /* Should have been reduced */
2801 for (i
= 0; i
< I
->NbConstraints
; ++i
)
2802 if (value_pos_p(I
->Constraint
[i
][1]))
2805 assert(i
< I
->NbConstraints
);
2807 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
2808 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
2809 if (value_neg_p(min
)) {
2811 mpz_fdiv_q(min
, min
, d
);
2812 value_init(offset
.d
);
2813 value_set_si(offset
.d
, 1);
2814 value_init(offset
.x
.n
);
2815 value_oppose(offset
.x
.n
, min
);
2816 eadd(&offset
, &p
->arr
[0]);
2817 free_evalue_refs(&offset
);
2826 value_set_si(fl
.d
, 0);
2827 fl
.x
.p
= new_enode(flooring
, 3, -1);
2828 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
2829 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
2830 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
2832 eadd(&fl
, &p
->arr
[0]);
2833 reorder_terms(p
, &p
->arr
[0]);
2836 free_evalue_refs(&fl
);
2841 void evalue_frac2floor(evalue
*e
)
2844 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2847 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2848 if (frac2floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
2849 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
])))
2850 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
2853 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
2858 int nparam
= D
->Dimension
- nvar
;
2861 nr
= D
->NbConstraints
+ 2;
2862 nc
= D
->Dimension
+ 2 + 1;
2863 C
= Matrix_Alloc(nr
, nc
);
2864 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
2865 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
2866 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
2867 D
->Dimension
+ 1 - nvar
);
2872 nc
= C
->NbColumns
+ 1;
2873 C
= Matrix_Alloc(nr
, nc
);
2874 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
2875 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
2876 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
2877 oldC
->NbColumns
- 1 - nvar
);
2880 value_set_si(C
->p
[nr
-2][0], 1);
2881 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
2882 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
2884 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
2885 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
2891 static void floor2frac_r(evalue
*e
, int nvar
)
2898 if (value_notzero_p(e
->d
))
2903 assert(p
->type
== flooring
);
2904 for (i
= 1; i
< p
->size
; i
++)
2905 floor2frac_r(&p
->arr
[i
], nvar
);
2907 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
2908 assert(pp
->x
.p
->type
== polynomial
);
2909 pp
->x
.p
->pos
-= nvar
;
2913 value_set_si(f
.d
, 0);
2914 f
.x
.p
= new_enode(fractional
, 3, -1);
2915 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2916 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
2917 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
2919 eadd(&f
, &p
->arr
[0]);
2920 reorder_terms(p
, &p
->arr
[0]);
2923 free_evalue_refs(&f
);
2926 /* Convert flooring back to fractional and shift position
2927 * of the parameters by nvar
2929 static void floor2frac(evalue
*e
, int nvar
)
2931 floor2frac_r(e
, nvar
);
2935 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
2938 int nparam
= D
->Dimension
- nvar
;
2942 D
= Constraints2Polyhedron(C
, 0);
2946 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
2948 /* Double check that D was not unbounded. */
2949 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
2957 evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
2964 evalue
*factor
= NULL
;
2967 if (EVALUE_IS_ZERO(*e
))
2971 Polyhedron
*DD
= Disjoint_Domain(D
, 0, 0);
2978 res
= esum_over_domain(e
, nvar
, Q
, C
);
2981 for (Q
= DD
; Q
; Q
= DD
) {
2987 t
= esum_over_domain(e
, nvar
, Q
, C
);
2994 free_evalue_refs(t
);
3001 if (value_notzero_p(e
->d
)) {
3004 t
= esum_over_domain_cst(nvar
, D
, C
);
3006 if (!EVALUE_IS_ONE(*e
))
3012 switch (e
->x
.p
->type
) {
3014 evalue
*pp
= &e
->x
.p
->arr
[0];
3016 if (pp
->x
.p
->pos
> nvar
) {
3017 /* remainder is independent of the summated vars */
3023 floor2frac(&f
, nvar
);
3025 t
= esum_over_domain_cst(nvar
, D
, C
);
3029 free_evalue_refs(&f
);
3034 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3035 poly_denom(pp
, &row
->p
[1 + nvar
]);
3036 value_set_si(row
->p
[0], 1);
3037 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3038 pp
= &pp
->x
.p
->arr
[0]) {
3040 assert(pp
->x
.p
->type
== polynomial
);
3042 if (pos
>= 1 + nvar
)
3044 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3045 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3046 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3048 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3049 value_division(row
->p
[1 + D
->Dimension
+ 1],
3050 row
->p
[1 + D
->Dimension
+ 1],
3052 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3053 row
->p
[1 + D
->Dimension
+ 1],
3055 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3059 int pos
= e
->x
.p
->pos
;
3063 value_init(factor
->d
);
3064 value_set_si(factor
->d
, 0);
3065 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3066 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3067 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3071 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3072 for (i
= 0; i
< D
->NbRays
; ++i
)
3073 if (value_notzero_p(D
->Ray
[i
][pos
]))
3075 assert(i
< D
->NbRays
);
3076 if (value_neg_p(D
->Ray
[i
][pos
])) {
3078 value_init(factor
->d
);
3079 evalue_set_si(factor
, -1, 1);
3081 value_set_si(row
->p
[0], 1);
3082 value_set_si(row
->p
[pos
], 1);
3083 value_set_si(row
->p
[1 + nvar
], -1);
3090 i
= type_offset(e
->x
.p
);
3092 res
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
3097 evalue_copy(&cum
, factor
);
3101 for (; i
< e
->x
.p
->size
; ++i
) {
3105 C
= esum_add_constraint(nvar
, D
, C
, row
);
3111 Vector_Print(stderr, P_VALUE_FMT, row);
3113 Matrix_Print(stderr, P_VALUE_FMT, C);
3115 t
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
3124 free_evalue_refs(t
);
3127 if (factor
&& i
+1 < e
->x
.p
->size
)
3134 free_evalue_refs(factor
);
3135 free_evalue_refs(&cum
);
3147 evalue
*esum(evalue
*e
, int nvar
)
3155 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3156 evalue_copy(res
, e
);
3160 evalue_set_si(res
, 0, 1);
3162 assert(value_zero_p(e
->d
));
3163 assert(e
->x
.p
->type
== partition
);
3165 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3167 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3168 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
3170 free_evalue_refs(t
);
3179 /* Initial silly implementation */
3180 void eor(evalue
*e1
, evalue
*res
)
3186 evalue_set_si(&mone
, -1, 1);
3188 evalue_copy(&E
, res
);
3194 free_evalue_refs(&E
);
3195 free_evalue_refs(&mone
);