5 #include "ev_operations.h"
10 #define value_pmodulus(ref,val1,val2) (mpz_fdiv_r((ref),(val1),(val2)))
13 #define ALLOC(p) p = (typeof(p))malloc(sizeof(*p))
14 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
16 void evalue_set_si(evalue
*ev
, int n
, int d
) {
17 value_set_si(ev
->d
, d
);
19 value_set_si(ev
->x
.n
, n
);
22 void evalue_set(evalue
*ev
, Value n
, Value d
) {
23 value_assign(ev
->d
, d
);
25 value_assign(ev
->x
.n
, n
);
28 void aep_evalue(evalue
*e
, int *ref
) {
33 if (value_notzero_p(e
->d
))
34 return; /* a rational number, its already reduced */
36 return; /* hum... an overflow probably occured */
38 /* First check the components of p */
39 for (i
=0;i
<p
->size
;i
++)
40 aep_evalue(&p
->arr
[i
],ref
);
47 p
->pos
= ref
[p
->pos
-1]+1;
53 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
59 if (value_notzero_p(e
->d
))
60 return; /* a rational number, its already reduced */
62 return; /* hum... an overflow probably occured */
65 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
66 for(i
=0;i
<CT
->NbRows
-1;i
++)
67 for(j
=0;j
<CT
->NbColumns
;j
++)
68 if(value_notzero_p(CT
->p
[i
][j
])) {
73 /* Transform the references in e, using ref */
77 } /* addeliminatedparams_evalue */
79 void addeliminatedparams_enum(evalue
*e
, Matrix
*CT
, Polyhedron
*CEq
,
80 unsigned MaxRays
, unsigned nparam
)
85 if (CT
->NbRows
== CT
->NbColumns
)
88 if (EVALUE_IS_ZERO(*e
))
91 if (value_notzero_p(e
->d
)) {
94 value_set_si(res
.d
, 0);
95 res
.x
.p
= new_enode(partition
, 2, nparam
);
96 EVALUE_SET_DOMAIN(res
.x
.p
->arr
[0],
97 DomainConstraintSimplify(Polyhedron_Copy(CEq
), MaxRays
));
98 value_clear(res
.x
.p
->arr
[1].d
);
106 assert(p
->type
== partition
);
109 for (i
=0; i
<p
->size
/2; i
++) {
110 Polyhedron
*D
= EVALUE_DOMAIN(p
->arr
[2*i
]);
111 Polyhedron
*T
= DomainPreimage(D
, CT
, MaxRays
);
114 T
= DomainIntersection(D
, CEq
, MaxRays
);
116 EVALUE_SET_DOMAIN(p
->arr
[2*i
], T
);
117 addeliminatedparams_evalue(&p
->arr
[2*i
+1], CT
);
121 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
129 assert(value_notzero_p(e1
->d
));
130 assert(value_notzero_p(e2
->d
));
131 value_multiply(m
, e1
->x
.n
, e2
->d
);
132 value_multiply(m2
, e2
->x
.n
, e1
->d
);
135 else if (value_gt(m
, m2
))
145 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
147 if (value_notzero_p(e1
->d
)) {
149 if (value_zero_p(e2
->d
))
151 r
= mod_rational_smaller(e1
, e2
);
152 return r
== -1 ? 0 : r
;
154 if (value_notzero_p(e2
->d
))
156 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
158 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
161 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
163 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
168 static int mod_term_smaller(evalue
*e1
, evalue
*e2
)
170 assert(value_zero_p(e1
->d
));
171 assert(value_zero_p(e2
->d
));
172 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
173 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
174 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
177 /* Negative pos means inequality */
178 /* s is negative of substitution if m is not zero */
187 struct fixed_param
*fixed
;
192 static int relations_depth(evalue
*e
)
197 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
198 e
= &e
->x
.p
->arr
[1], ++d
);
202 static void poly_denom(evalue
*p
, Value
*d
)
206 while (value_zero_p(p
->d
)) {
207 assert(p
->x
.p
->type
== polynomial
);
208 assert(p
->x
.p
->size
== 2);
209 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
210 value_lcm(*d
, p
->x
.p
->arr
[1].d
, d
);
213 value_lcm(*d
, p
->d
, d
);
216 #define EVALUE_IS_ONE(ev) (value_one_p((ev).d) && value_one_p((ev).x.n))
218 static void realloc_substitution(struct subst
*s
, int d
)
220 struct fixed_param
*n
;
223 for (i
= 0; i
< s
->n
; ++i
)
230 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
236 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
239 /* May have been reduced already */
240 if (value_notzero_p(m
->d
))
243 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
244 assert(m
->x
.p
->size
== 3);
246 /* fractional was inverted during reduction
247 * invert it back and move constant in
249 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
250 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
251 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
252 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
253 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
254 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
255 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
256 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
257 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
258 value_set_si(m
->x
.p
->arr
[1].d
, 1);
261 /* Oops. Nested identical relations. */
262 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
265 if (s
->n
>= s
->max
) {
266 int d
= relations_depth(r
);
267 realloc_substitution(s
, d
);
271 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
272 assert(p
->x
.p
->size
== 2);
275 assert(value_pos_p(f
->x
.n
));
277 value_init(s
->fixed
[s
->n
].m
);
278 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
279 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
280 value_init(s
->fixed
[s
->n
].d
);
281 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
282 value_init(s
->fixed
[s
->n
].s
.d
);
283 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
289 static int type_offset(enode
*p
)
291 return p
->type
== fractional
? 1 :
292 p
->type
== flooring
? 1 : 0;
295 static void reorder_terms(enode
*p
, evalue
*v
)
298 int offset
= type_offset(p
);
300 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
302 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
303 free_evalue_refs(&(p
->arr
[i
]));
309 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
315 if (value_notzero_p(e
->d
)) {
317 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
318 return; /* a rational number, its already reduced */
322 return; /* hum... an overflow probably occured */
324 /* First reduce the components of p */
325 add
= p
->type
== relation
;
326 for (i
=0; i
<p
->size
; i
++) {
328 add
= add_modulo_substitution(s
, e
);
330 if (i
== 0 && p
->type
==fractional
)
331 _reduce_evalue(&p
->arr
[i
], s
, 1);
333 _reduce_evalue(&p
->arr
[i
], s
, fract
);
335 if (add
&& i
== p
->size
-1) {
337 value_clear(s
->fixed
[s
->n
].m
);
338 value_clear(s
->fixed
[s
->n
].d
);
339 free_evalue_refs(&s
->fixed
[s
->n
].s
);
340 } else if (add
&& i
== 1)
341 s
->fixed
[s
->n
-1].pos
*= -1;
344 if (p
->type
==periodic
) {
346 /* Try to reduce the period */
347 for (i
=1; i
<=(p
->size
)/2; i
++) {
348 if ((p
->size
% i
)==0) {
350 /* Can we reduce the size to i ? */
352 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
353 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
356 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
360 you_lose
: /* OK, lets not do it */
365 /* Try to reduce its strength */
368 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
372 else if (p
->type
==polynomial
) {
373 for (k
= 0; s
&& k
< s
->n
; ++k
) {
374 if (s
->fixed
[k
].pos
== p
->pos
) {
375 int divide
= value_notone_p(s
->fixed
[k
].d
);
378 if (value_notzero_p(s
->fixed
[k
].m
)) {
381 assert(p
->size
== 2);
382 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
384 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
391 value_assign(d
.d
, s
->fixed
[k
].d
);
393 if (value_notzero_p(s
->fixed
[k
].m
))
394 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
396 value_set_si(d
.x
.n
, 1);
399 for (i
=p
->size
-1;i
>=1;i
--) {
400 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
402 emul(&d
, &p
->arr
[i
]);
403 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
404 free_evalue_refs(&(p
->arr
[i
]));
407 _reduce_evalue(&p
->arr
[0], s
, fract
);
410 free_evalue_refs(&d
);
416 /* Try to reduce the degree */
417 for (i
=p
->size
-1;i
>=1;i
--) {
418 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
420 /* Zero coefficient */
421 free_evalue_refs(&(p
->arr
[i
]));
426 /* Try to reduce its strength */
429 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
433 else if (p
->type
==fractional
) {
437 if (value_notzero_p(p
->arr
[0].d
)) {
439 value_assign(v
.d
, p
->arr
[0].d
);
441 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
446 evalue
*pp
= &p
->arr
[0];
447 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
448 assert(pp
->x
.p
->size
== 2);
450 /* search for exact duplicate among the modulo inequalities */
452 f
= &pp
->x
.p
->arr
[1];
453 for (k
= 0; s
&& k
< s
->n
; ++k
) {
454 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
455 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
456 value_eq(s
->fixed
[k
].m
, f
->d
) &&
457 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
464 /* replace { E/m } by { (E-1)/m } + 1/m */
469 evalue_set_si(&extra
, 1, 1);
470 value_assign(extra
.d
, g
);
471 eadd(&extra
, &v
.x
.p
->arr
[1]);
472 free_evalue_refs(&extra
);
474 /* We've been going in circles; stop now */
475 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
476 free_evalue_refs(&v
);
478 evalue_set_si(&v
, 0, 1);
483 value_set_si(v
.d
, 0);
484 v
.x
.p
= new_enode(fractional
, 3, -1);
485 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
486 value_assign(v
.x
.p
->arr
[1].d
, g
);
487 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
488 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
491 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
494 value_division(f
->d
, g
, f
->d
);
495 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
496 value_assign(f
->d
, g
);
497 value_decrement(f
->x
.n
, f
->x
.n
);
498 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
500 Gcd(f
->d
, f
->x
.n
, &g
);
501 value_division(f
->d
, f
->d
, g
);
502 value_division(f
->x
.n
, f
->x
.n
, g
);
511 /* reduction may have made this fractional arg smaller */
512 i
= reorder
? p
->size
: 1;
513 for ( ; i
< p
->size
; ++i
)
514 if (value_zero_p(p
->arr
[i
].d
) &&
515 p
->arr
[i
].x
.p
->type
== fractional
&&
516 !mod_term_smaller(e
, &p
->arr
[i
]))
520 value_set_si(v
.d
, 0);
521 v
.x
.p
= new_enode(fractional
, 3, -1);
522 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
523 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
524 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
534 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
535 pp
= &pp
->x
.p
->arr
[0]) {
536 f
= &pp
->x
.p
->arr
[1];
537 assert(value_pos_p(f
->d
));
538 mpz_mul_ui(twice
, f
->x
.n
, 2);
539 if (value_lt(twice
, f
->d
))
541 if (value_eq(twice
, f
->d
))
549 value_set_si(v
.d
, 0);
550 v
.x
.p
= new_enode(fractional
, 3, -1);
551 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
552 poly_denom(&p
->arr
[0], &twice
);
553 value_assign(v
.x
.p
->arr
[1].d
, twice
);
554 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
555 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
556 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
558 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
559 pp
= &pp
->x
.p
->arr
[0]) {
560 f
= &pp
->x
.p
->arr
[1];
561 value_oppose(f
->x
.n
, f
->x
.n
);
562 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
564 value_division(pp
->d
, twice
, pp
->d
);
565 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
566 value_assign(pp
->d
, twice
);
567 value_oppose(pp
->x
.n
, pp
->x
.n
);
568 value_decrement(pp
->x
.n
, pp
->x
.n
);
569 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
579 reorder_terms(p
, &v
);
580 _reduce_evalue(&p
->arr
[1], s
, fract
);
583 /* Try to reduce the degree */
584 for (i
=p
->size
-1;i
>=2;i
--) {
585 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
587 /* Zero coefficient */
588 free_evalue_refs(&(p
->arr
[i
]));
593 /* Try to reduce its strength */
596 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
597 free_evalue_refs(&(p
->arr
[0]));
601 else if (p
->type
== flooring
) {
602 /* Try to reduce the degree */
603 for (i
=p
->size
-1;i
>=2;i
--) {
604 if (!EVALUE_IS_ZERO(p
->arr
[i
]))
606 /* Zero coefficient */
607 free_evalue_refs(&(p
->arr
[i
]));
612 /* Try to reduce its strength */
615 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
616 free_evalue_refs(&(p
->arr
[0]));
620 else if (p
->type
== relation
) {
621 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
622 free_evalue_refs(&(p
->arr
[2]));
623 free_evalue_refs(&(p
->arr
[0]));
630 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
631 free_evalue_refs(&(p
->arr
[2]));
634 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
635 free_evalue_refs(&(p
->arr
[1]));
636 free_evalue_refs(&(p
->arr
[0]));
637 evalue_set_si(e
, 0, 1);
644 /* Relation was reduced by means of an identical
645 * inequality => remove
647 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
650 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
651 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
653 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
655 free_evalue_refs(&(p
->arr
[2]));
659 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
661 evalue_set_si(e
, 0, 1);
662 free_evalue_refs(&(p
->arr
[1]));
664 free_evalue_refs(&(p
->arr
[0]));
670 } /* reduce_evalue */
672 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
677 for (k
= 0; k
< dim
; ++k
)
678 if (value_notzero_p(row
[k
+1]))
681 Vector_Normalize_Positive(row
+1, dim
+1, k
);
682 assert(s
->n
< s
->max
);
683 value_init(s
->fixed
[s
->n
].d
);
684 value_init(s
->fixed
[s
->n
].m
);
685 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
686 s
->fixed
[s
->n
].pos
= k
+1;
687 value_set_si(s
->fixed
[s
->n
].m
, 0);
688 r
= &s
->fixed
[s
->n
].s
;
690 for (l
= k
+1; l
< dim
; ++l
)
691 if (value_notzero_p(row
[l
+1])) {
692 value_set_si(r
->d
, 0);
693 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
694 value_init(r
->x
.p
->arr
[1].x
.n
);
695 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
696 value_set_si(r
->x
.p
->arr
[1].d
, 1);
700 value_oppose(r
->x
.n
, row
[dim
+1]);
701 value_set_si(r
->d
, 1);
705 void reduce_evalue (evalue
*e
) {
706 struct subst s
= { NULL
, 0, 0 };
708 if (value_notzero_p(e
->d
))
709 return; /* a rational number, its already reduced */
711 if (e
->x
.p
->type
== partition
) {
714 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
715 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
717 /* This shouldn't really happen;
718 * Empty domains should not be added.
725 D
= DomainConvex(D
, 0);
726 if (!D
->next
&& D
->NbEq
) {
730 realloc_substitution(&s
, dim
);
732 int d
= relations_depth(&e
->x
.p
->arr
[2*i
+1]);
734 NALLOC(s
.fixed
, s
.max
);
737 for (j
= 0; j
< D
->NbEq
; ++j
)
738 add_substitution(&s
, D
->Constraint
[j
], dim
);
740 if (D
!= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))
742 _reduce_evalue(&e
->x
.p
->arr
[2*i
+1], &s
, 0);
743 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
745 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
746 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
747 value_clear(e
->x
.p
->arr
[2*i
].d
);
749 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
750 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
755 for (j
= 0; j
< s
.n
; ++j
) {
756 value_clear(s
.fixed
[j
].d
);
757 value_clear(s
.fixed
[j
].m
);
758 free_evalue_refs(&s
.fixed
[j
].s
);
762 if (e
->x
.p
->size
== 0) {
764 evalue_set_si(e
, 0, 1);
767 _reduce_evalue(e
, &s
, 0);
772 void print_evalue(FILE *DST
,evalue
*e
,char **pname
) {
774 if(value_notzero_p(e
->d
)) {
775 if(value_notone_p(e
->d
)) {
776 value_print(DST
,VALUE_FMT
,e
->x
.n
);
778 value_print(DST
,VALUE_FMT
,e
->d
);
781 value_print(DST
,VALUE_FMT
,e
->x
.n
);
785 print_enode(DST
,e
->x
.p
,pname
);
789 void print_enode(FILE *DST
,enode
*p
,char **pname
) {
794 fprintf(DST
, "NULL");
800 for (i
=0; i
<p
->size
; i
++) {
801 print_evalue(DST
, &p
->arr
[i
], pname
);
805 fprintf(DST
, " }\n");
809 for (i
=p
->size
-1; i
>=0; i
--) {
810 print_evalue(DST
, &p
->arr
[i
], pname
);
811 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
813 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
815 fprintf(DST
, " )\n");
819 for (i
=0; i
<p
->size
; i
++) {
820 print_evalue(DST
, &p
->arr
[i
], pname
);
821 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
823 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
828 for (i
=p
->size
-1; i
>=1; i
--) {
829 print_evalue(DST
, &p
->arr
[i
], pname
);
832 fprintf(DST
, p
->type
== flooring
? "[" : "{");
833 print_evalue(DST
, &p
->arr
[0], pname
);
834 fprintf(DST
, p
->type
== flooring
? "]" : "}");
836 fprintf(DST
, "^%d + ", i
-1);
841 fprintf(DST
, " )\n");
845 print_evalue(DST
, &p
->arr
[0], pname
);
846 fprintf(DST
, "= 0 ] * \n");
847 print_evalue(DST
, &p
->arr
[1], pname
);
849 fprintf(DST
, " +\n [ ");
850 print_evalue(DST
, &p
->arr
[0], pname
);
851 fprintf(DST
, "!= 0 ] * \n");
852 print_evalue(DST
, &p
->arr
[2], pname
);
856 char **names
= pname
;
857 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
858 if (p
->pos
< maxdim
) {
859 NALLOC(names
, maxdim
);
860 for (i
= 0; i
< p
->pos
; ++i
)
862 for ( ; i
< maxdim
; ++i
) {
863 NALLOC(names
[i
], 10);
864 snprintf(names
[i
], 10, "_p%d", i
);
868 for (i
=0; i
<p
->size
/2; i
++) {
869 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
870 print_evalue(DST
, &p
->arr
[2*i
+1], names
);
873 if (p
->pos
< maxdim
) {
874 for (i
= p
->pos
; i
< maxdim
; ++i
)
887 static void eadd_rev(evalue
*e1
, evalue
*res
)
891 evalue_copy(&ev
, e1
);
893 free_evalue_refs(res
);
897 static void eadd_rev_cst (evalue
*e1
, evalue
*res
)
901 evalue_copy(&ev
, e1
);
902 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
903 free_evalue_refs(res
);
907 struct section
{ Polyhedron
* D
; evalue E
; };
909 void eadd_partitions (evalue
*e1
,evalue
*res
)
914 s
= (struct section
*)
915 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
916 sizeof(struct section
));
918 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
919 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
920 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
923 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
924 assert(res
->x
.p
->size
>= 2);
925 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
926 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
928 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
930 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
939 value_init(s
[n
].E
.d
);
940 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
944 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
945 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
946 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
948 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
949 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
955 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
956 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
958 value_init(s
[n
].E
.d
);
959 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
960 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
965 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
969 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
972 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
973 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
974 value_clear(res
->x
.p
->arr
[2*i
].d
);
979 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
980 for (j
= 0; j
< n
; ++j
) {
981 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
982 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
983 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
984 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
990 static void explicit_complement(evalue
*res
)
992 enode
*rel
= new_enode(relation
, 3, 0);
994 value_clear(rel
->arr
[0].d
);
995 rel
->arr
[0] = res
->x
.p
->arr
[0];
996 value_clear(rel
->arr
[1].d
);
997 rel
->arr
[1] = res
->x
.p
->arr
[1];
998 value_set_si(rel
->arr
[2].d
, 1);
999 value_init(rel
->arr
[2].x
.n
);
1000 value_set_si(rel
->arr
[2].x
.n
, 0);
1005 void eadd(evalue
*e1
,evalue
*res
) {
1008 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1009 /* Add two rational numbers */
1015 value_multiply(m1
,e1
->x
.n
,res
->d
);
1016 value_multiply(m2
,res
->x
.n
,e1
->d
);
1017 value_addto(res
->x
.n
,m1
,m2
);
1018 value_multiply(res
->d
,e1
->d
,res
->d
);
1019 Gcd(res
->x
.n
,res
->d
,&g
);
1020 if (value_notone_p(g
)) {
1021 value_division(res
->d
,res
->d
,g
);
1022 value_division(res
->x
.n
,res
->x
.n
,g
);
1024 value_clear(g
); value_clear(m1
); value_clear(m2
);
1027 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
1028 switch (res
->x
.p
->type
) {
1030 /* Add the constant to the constant term of a polynomial*/
1031 eadd(e1
, &res
->x
.p
->arr
[0]);
1034 /* Add the constant to all elements of a periodic number */
1035 for (i
=0; i
<res
->x
.p
->size
; i
++) {
1036 eadd(e1
, &res
->x
.p
->arr
[i
]);
1040 fprintf(stderr
, "eadd: cannot add const with vector\n");
1044 eadd(e1
, &res
->x
.p
->arr
[1]);
1047 assert(EVALUE_IS_ZERO(*e1
));
1048 break; /* Do nothing */
1050 /* Create (zero) complement if needed */
1051 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
1052 explicit_complement(res
);
1053 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1054 eadd(e1
, &res
->x
.p
->arr
[i
]);
1060 /* add polynomial or periodic to constant
1061 * you have to exchange e1 and res, before doing addition */
1063 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1067 else { // ((e1->d==0) && (res->d==0))
1068 assert(!((e1
->x
.p
->type
== partition
) ^
1069 (res
->x
.p
->type
== partition
)));
1070 if (e1
->x
.p
->type
== partition
) {
1071 eadd_partitions(e1
, res
);
1074 if (e1
->x
.p
->type
== relation
&&
1075 (res
->x
.p
->type
!= relation
||
1076 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1080 if (res
->x
.p
->type
== relation
) {
1081 if (e1
->x
.p
->type
== relation
&&
1082 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1083 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1084 explicit_complement(res
);
1085 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1086 explicit_complement(e1
);
1087 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1088 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1091 if (res
->x
.p
->size
< 3)
1092 explicit_complement(res
);
1093 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1094 eadd(e1
, &res
->x
.p
->arr
[i
]);
1097 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
1098 /* adding to evalues of different type. two cases are possible
1099 * res is periodic and e1 is polynomial, you have to exchange
1100 * e1 and res then to add e1 to the constant term of res */
1101 if (e1
->x
.p
->type
== polynomial
) {
1102 eadd_rev_cst(e1
, res
);
1104 else if (res
->x
.p
->type
== polynomial
) {
1105 /* res is polynomial and e1 is periodic,
1106 add e1 to the constant term of res */
1108 eadd(e1
,&res
->x
.p
->arr
[0]);
1114 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
1115 ((res
->x
.p
->type
== fractional
||
1116 res
->x
.p
->type
== flooring
) &&
1117 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1118 /* adding evalues of different position (i.e function of different unknowns
1119 * to case are possible */
1121 switch (res
->x
.p
->type
) {
1124 if(mod_term_smaller(res
, e1
))
1125 eadd(e1
,&res
->x
.p
->arr
[1]);
1127 eadd_rev_cst(e1
, res
);
1129 case polynomial
: // res and e1 are polynomials
1130 // add e1 to the constant term of res
1132 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1133 eadd(e1
,&res
->x
.p
->arr
[0]);
1135 eadd_rev_cst(e1
, res
);
1136 // value_clear(g); value_clear(m1); value_clear(m2);
1138 case periodic
: // res and e1 are pointers to periodic numbers
1139 //add e1 to all elements of res
1141 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1142 for (i
=0;i
<res
->x
.p
->size
;i
++) {
1143 eadd(e1
,&res
->x
.p
->arr
[i
]);
1154 //same type , same pos and same size
1155 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1156 // add any element in e1 to the corresponding element in res
1157 i
= type_offset(res
->x
.p
);
1159 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1160 for (; i
<res
->x
.p
->size
; i
++) {
1161 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1166 /* Sizes are different */
1167 switch(res
->x
.p
->type
) {
1171 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1172 /* new enode and add res to that new node. If you do not do */
1173 /* that, you lose the the upper weight part of e1 ! */
1175 if(e1
->x
.p
->size
> res
->x
.p
->size
)
1178 i
= type_offset(res
->x
.p
);
1180 assert(eequal(&e1
->x
.p
->arr
[0],
1181 &res
->x
.p
->arr
[0]));
1182 for (; i
<e1
->x
.p
->size
; i
++) {
1183 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1190 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1193 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1194 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1195 to add periodicaly elements of e1 to elements of ne, and finaly to
1200 value_init(ex
); value_init(ey
);value_init(ep
);
1203 value_set_si(ex
,e1
->x
.p
->size
);
1204 value_set_si(ey
,res
->x
.p
->size
);
1205 value_assign (ep
,*Lcm(ex
,ey
));
1206 p
=(int)mpz_get_si(ep
);
1207 ne
= (evalue
*) malloc (sizeof(evalue
));
1209 value_set_si( ne
->d
,0);
1211 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1213 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1216 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1219 value_assign(res
->d
, ne
->d
);
1225 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
1234 static void emul_rev (evalue
*e1
, evalue
*res
)
1238 evalue_copy(&ev
, e1
);
1240 free_evalue_refs(res
);
1244 static void emul_poly (evalue
*e1
, evalue
*res
)
1246 int i
, j
, o
= type_offset(res
->x
.p
);
1248 int size
=(e1
->x
.p
->size
+ res
->x
.p
->size
- o
- 1);
1250 value_set_si(tmp
.d
,0);
1251 tmp
.x
.p
=new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1253 evalue_copy(&tmp
.x
.p
->arr
[0], &e1
->x
.p
->arr
[0]);
1254 for (i
=o
; i
< e1
->x
.p
->size
; i
++) {
1255 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1256 emul(&res
->x
.p
->arr
[o
], &tmp
.x
.p
->arr
[i
]);
1259 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1260 for (i
=o
+1; i
<res
->x
.p
->size
; i
++)
1261 for (j
=o
; j
<e1
->x
.p
->size
; j
++) {
1264 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1265 emul(&res
->x
.p
->arr
[i
], &ev
);
1266 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-o
]);
1267 free_evalue_refs(&ev
);
1269 free_evalue_refs(res
);
1273 void emul_partitions (evalue
*e1
,evalue
*res
)
1278 s
= (struct section
*)
1279 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1280 sizeof(struct section
));
1282 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1283 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1284 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1287 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1288 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1289 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1290 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1296 /* This code is only needed because the partitions
1297 are not true partitions.
1299 for (k
= 0; k
< n
; ++k
) {
1300 if (DomainIncludes(s
[k
].D
, d
))
1302 if (DomainIncludes(d
, s
[k
].D
)) {
1303 Domain_Free(s
[k
].D
);
1304 free_evalue_refs(&s
[k
].E
);
1315 value_init(s
[n
].E
.d
);
1316 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1317 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1321 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1322 value_clear(res
->x
.p
->arr
[2*i
].d
);
1323 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1328 evalue_set_si(res
, 0, 1);
1330 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1331 for (j
= 0; j
< n
; ++j
) {
1332 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1333 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1334 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1335 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1342 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1344 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1345 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1346 * evalues is not treated here */
1348 void emul (evalue
*e1
, evalue
*res
){
1351 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1352 fprintf(stderr
, "emul: do not proced on evector type !\n");
1356 if (EVALUE_IS_ZERO(*res
))
1359 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1360 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1361 emul_partitions(e1
, res
);
1364 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1365 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1366 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1368 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1369 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1370 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1371 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1372 explicit_complement(res
);
1373 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1374 explicit_complement(e1
);
1375 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1376 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1379 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1380 emul(e1
, &res
->x
.p
->arr
[i
]);
1382 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1383 switch(e1
->x
.p
->type
) {
1385 switch(res
->x
.p
->type
) {
1387 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1388 /* Product of two polynomials of the same variable */
1393 /* Product of two polynomials of different variables */
1395 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1396 for( i
=0; i
<res
->x
.p
->size
; i
++)
1397 emul(e1
, &res
->x
.p
->arr
[i
]);
1406 /* Product of a polynomial and a periodic or fractional */
1413 switch(res
->x
.p
->type
) {
1415 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1416 /* Product of two periodics of the same parameter and period */
1418 for(i
=0; i
<res
->x
.p
->size
;i
++)
1419 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1424 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1425 /* Product of two periodics of the same parameter and different periods */
1429 value_init(x
); value_init(y
);value_init(z
);
1432 value_set_si(x
,e1
->x
.p
->size
);
1433 value_set_si(y
,res
->x
.p
->size
);
1434 value_assign (z
,*Lcm(x
,y
));
1435 lcm
=(int)mpz_get_si(z
);
1436 newp
= (evalue
*) malloc (sizeof(evalue
));
1437 value_init(newp
->d
);
1438 value_set_si( newp
->d
,0);
1439 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1440 for(i
=0;i
<lcm
;i
++) {
1441 evalue_copy(&newp
->x
.p
->arr
[i
],
1442 &res
->x
.p
->arr
[i
%iy
]);
1445 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1447 value_assign(res
->d
,newp
->d
);
1450 value_clear(x
); value_clear(y
);value_clear(z
);
1454 /* Product of two periodics of different parameters */
1456 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1457 for(i
=0; i
<res
->x
.p
->size
; i
++)
1458 emul(e1
, &(res
->x
.p
->arr
[i
]));
1466 /* Product of a periodic and a polynomial */
1468 for(i
=0; i
<res
->x
.p
->size
; i
++)
1469 emul(e1
, &(res
->x
.p
->arr
[i
]));
1476 switch(res
->x
.p
->type
) {
1478 for(i
=0; i
<res
->x
.p
->size
; i
++)
1479 emul(e1
, &(res
->x
.p
->arr
[i
]));
1486 assert(e1
->x
.p
->type
== res
->x
.p
->type
);
1487 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1488 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1491 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1492 if (e1
->x
.p
->type
!= fractional
|| !value_two_p(d
.d
))
1497 value_set_si(d
.x
.n
, 1);
1498 /* { x }^2 == { x }/2 */
1499 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1500 assert(e1
->x
.p
->size
== 3);
1501 assert(res
->x
.p
->size
== 3);
1503 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1505 eadd(&res
->x
.p
->arr
[1], &tmp
);
1506 emul(&e1
->x
.p
->arr
[2], &tmp
);
1507 emul(&e1
->x
.p
->arr
[1], res
);
1508 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1509 free_evalue_refs(&tmp
);
1514 if(mod_term_smaller(res
, e1
))
1515 for(i
=1; i
<res
->x
.p
->size
; i
++)
1516 emul(e1
, &(res
->x
.p
->arr
[i
]));
1531 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1532 /* Product of two rational numbers */
1536 value_multiply(res
->d
,e1
->d
,res
->d
);
1537 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1538 Gcd(res
->x
.n
, res
->d
,&g
);
1539 if (value_notone_p(g
)) {
1540 value_division(res
->d
,res
->d
,g
);
1541 value_division(res
->x
.n
,res
->x
.n
,g
);
1547 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1548 /* Product of an expression (polynomial or peririodic) and a rational number */
1554 /* Product of a rationel number and an expression (polynomial or peririodic) */
1556 i
= type_offset(res
->x
.p
);
1557 for (; i
<res
->x
.p
->size
; i
++)
1558 emul(e1
, &res
->x
.p
->arr
[i
]);
1568 /* Frees mask content ! */
1569 void emask(evalue
*mask
, evalue
*res
) {
1576 if (EVALUE_IS_ZERO(*res
)) {
1577 free_evalue_refs(mask
);
1581 assert(value_zero_p(mask
->d
));
1582 assert(mask
->x
.p
->type
== partition
);
1583 assert(value_zero_p(res
->d
));
1584 assert(res
->x
.p
->type
== partition
);
1585 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1586 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1587 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1588 pos
= res
->x
.p
->pos
;
1590 s
= (struct section
*)
1591 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1592 sizeof(struct section
));
1596 evalue_set_si(&mone
, -1, 1);
1599 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1600 assert(mask
->x
.p
->size
>= 2);
1601 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1602 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1604 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1606 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1615 value_init(s
[n
].E
.d
);
1616 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1620 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1621 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1624 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1625 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1626 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1627 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1629 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1630 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1636 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1637 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1639 value_init(s
[n
].E
.d
);
1640 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1641 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1647 /* Just ignore; this may have been previously masked off */
1649 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1653 free_evalue_refs(&mone
);
1654 free_evalue_refs(mask
);
1655 free_evalue_refs(res
);
1658 evalue_set_si(res
, 0, 1);
1660 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1661 for (j
= 0; j
< n
; ++j
) {
1662 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1663 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1664 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1671 void evalue_copy(evalue
*dst
, evalue
*src
)
1673 value_assign(dst
->d
, src
->d
);
1674 if(value_notzero_p(src
->d
)) {
1675 value_init(dst
->x
.n
);
1676 value_assign(dst
->x
.n
, src
->x
.n
);
1678 dst
->x
.p
= ecopy(src
->x
.p
);
1681 enode
*new_enode(enode_type type
,int size
,int pos
) {
1687 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1690 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1694 for(i
=0; i
<size
; i
++) {
1695 value_init(res
->arr
[i
].d
);
1696 value_set_si(res
->arr
[i
].d
,0);
1697 res
->arr
[i
].x
.p
= 0;
1702 enode
*ecopy(enode
*e
) {
1707 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1708 for(i
=0;i
<e
->size
;++i
) {
1709 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1710 if(value_zero_p(res
->arr
[i
].d
))
1711 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1712 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1713 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1715 value_init(res
->arr
[i
].x
.n
);
1716 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1722 int ecmp(const evalue
*e1
, const evalue
*e2
)
1728 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
1732 value_multiply(m
, e1
->x
.n
, e2
->d
);
1733 value_multiply(m2
, e2
->x
.n
, e1
->d
);
1735 if (value_lt(m
, m2
))
1737 else if (value_gt(m
, m2
))
1747 if (value_notzero_p(e1
->d
))
1749 if (value_notzero_p(e2
->d
))
1755 if (p1
->type
!= p2
->type
)
1756 return p1
->type
- p2
->type
;
1757 if (p1
->pos
!= p2
->pos
)
1758 return p1
->pos
- p2
->pos
;
1759 if (p1
->size
!= p2
->size
)
1760 return p1
->size
- p2
->size
;
1762 for (i
= p1
->size
-1; i
>= 0; --i
)
1763 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
1769 int eequal(evalue
*e1
,evalue
*e2
) {
1774 if (value_ne(e1
->d
,e2
->d
))
1777 /* e1->d == e2->d */
1778 if (value_notzero_p(e1
->d
)) {
1779 if (value_ne(e1
->x
.n
,e2
->x
.n
))
1782 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1786 /* e1->d == e2->d == 0 */
1789 if (p1
->type
!= p2
->type
) return 0;
1790 if (p1
->size
!= p2
->size
) return 0;
1791 if (p1
->pos
!= p2
->pos
) return 0;
1792 for (i
=0; i
<p1
->size
; i
++)
1793 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
1798 void free_evalue_refs(evalue
*e
) {
1803 if (EVALUE_IS_DOMAIN(*e
)) {
1804 Domain_Free(EVALUE_DOMAIN(*e
));
1807 } else if (value_pos_p(e
->d
)) {
1809 /* 'e' stores a constant */
1811 value_clear(e
->x
.n
);
1814 assert(value_zero_p(e
->d
));
1817 if (!p
) return; /* null pointer */
1818 for (i
=0; i
<p
->size
; i
++) {
1819 free_evalue_refs(&(p
->arr
[i
]));
1823 } /* free_evalue_refs */
1825 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
1826 Vector
* val
, evalue
*res
)
1828 unsigned nparam
= periods
->Size
;
1831 double d
= compute_evalue(e
, val
->p
);
1832 d
*= VALUE_TO_DOUBLE(m
);
1837 value_assign(res
->d
, m
);
1838 value_init(res
->x
.n
);
1839 value_set_double(res
->x
.n
, d
);
1840 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
1843 if (value_one_p(periods
->p
[p
]))
1844 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
1849 value_assign(tmp
, periods
->p
[p
]);
1850 value_set_si(res
->d
, 0);
1851 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
1853 value_decrement(tmp
, tmp
);
1854 value_assign(val
->p
[p
], tmp
);
1855 mod2table_r(e
, periods
, m
, p
+1, val
,
1856 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
1857 } while (value_pos_p(tmp
));
1863 static void rel2table(evalue
*e
, int zero
)
1865 if (value_pos_p(e
->d
)) {
1866 if (value_zero_p(e
->x
.n
) == zero
)
1867 value_set_si(e
->x
.n
, 1);
1869 value_set_si(e
->x
.n
, 0);
1870 value_set_si(e
->d
, 1);
1873 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
1874 rel2table(&e
->x
.p
->arr
[i
], zero
);
1878 void evalue_mod2table(evalue
*e
, int nparam
)
1883 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
1886 for (i
=0; i
<p
->size
; i
++) {
1887 evalue_mod2table(&(p
->arr
[i
]), nparam
);
1889 if (p
->type
== relation
) {
1894 evalue_copy(©
, &p
->arr
[0]);
1896 rel2table(&p
->arr
[0], 1);
1897 emul(&p
->arr
[0], &p
->arr
[1]);
1899 rel2table(©
, 0);
1900 emul(©
, &p
->arr
[2]);
1901 eadd(&p
->arr
[2], &p
->arr
[1]);
1902 free_evalue_refs(&p
->arr
[2]);
1903 free_evalue_refs(©
);
1905 free_evalue_refs(&p
->arr
[0]);
1909 } else if (p
->type
== fractional
) {
1910 Vector
*periods
= Vector_Alloc(nparam
);
1911 Vector
*val
= Vector_Alloc(nparam
);
1917 value_set_si(tmp
, 1);
1918 Vector_Set(periods
->p
, 1, nparam
);
1919 Vector_Set(val
->p
, 0, nparam
);
1920 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
1923 assert(p
->type
== polynomial
);
1924 assert(p
->size
== 2);
1925 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
1926 value_lcm(tmp
, p
->arr
[1].d
, &tmp
);
1928 value_lcm(tmp
, ev
->d
, &tmp
);
1930 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
1933 evalue_set_si(&res
, 0, 1);
1934 /* Compute the polynomial using Horner's rule */
1935 for (i
=p
->size
-1;i
>1;i
--) {
1936 eadd(&p
->arr
[i
], &res
);
1939 eadd(&p
->arr
[1], &res
);
1941 free_evalue_refs(e
);
1942 free_evalue_refs(&EP
);
1947 Vector_Free(periods
);
1949 } /* evalue_mod2table */
1951 /********************************************************/
1952 /* function in domain */
1953 /* check if the parameters in list_args */
1954 /* verifies the constraints of Domain P */
1955 /********************************************************/
1956 int in_domain(Polyhedron
*P
, Value
*list_args
) {
1959 Value v
; /* value of the constraint of a row when
1960 parameters are instanciated*/
1966 /*P->Constraint constraint matrice of polyhedron P*/
1967 for(row
=0;row
<P
->NbConstraints
;row
++) {
1968 value_assign(v
,P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
1969 for(col
=1;col
<P
->Dimension
+1;col
++) {
1970 value_multiply(tmp
,P
->Constraint
[row
][col
],list_args
[col
-1]);
1971 value_addto(v
,v
,tmp
);
1973 if (value_notzero_p(P
->Constraint
[row
][0])) {
1975 /*if v is not >=0 then this constraint is not respected */
1976 if (value_neg_p(v
)) {
1980 return P
->next
? in_domain(P
->next
, list_args
) : 0;
1985 /*if v is not = 0 then this constraint is not respected */
1986 if (value_notzero_p(v
))
1991 /*if not return before this point => all
1992 the constraints are respected */
1998 /****************************************************/
1999 /* function compute enode */
2000 /* compute the value of enode p with parameters */
2001 /* list "list_args */
2002 /* compute the polynomial or the periodic */
2003 /****************************************************/
2005 static double compute_enode(enode
*p
, Value
*list_args
) {
2017 if (p
->type
== polynomial
) {
2019 value_assign(param
,list_args
[p
->pos
-1]);
2021 /* Compute the polynomial using Horner's rule */
2022 for (i
=p
->size
-1;i
>0;i
--) {
2023 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2024 res
*=VALUE_TO_DOUBLE(param
);
2026 res
+=compute_evalue(&p
->arr
[0],list_args
);
2028 else if (p
->type
== fractional
) {
2029 double d
= compute_evalue(&p
->arr
[0], list_args
);
2030 d
-= floor(d
+1e-10);
2032 /* Compute the polynomial using Horner's rule */
2033 for (i
=p
->size
-1;i
>1;i
--) {
2034 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2037 res
+=compute_evalue(&p
->arr
[1],list_args
);
2039 else if (p
->type
== flooring
) {
2040 double d
= compute_evalue(&p
->arr
[0], list_args
);
2043 /* Compute the polynomial using Horner's rule */
2044 for (i
=p
->size
-1;i
>1;i
--) {
2045 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2048 res
+=compute_evalue(&p
->arr
[1],list_args
);
2050 else if (p
->type
== periodic
) {
2051 value_assign(m
,list_args
[p
->pos
-1]);
2053 /* Choose the right element of the periodic */
2054 value_set_si(param
,p
->size
);
2055 value_pmodulus(m
,m
,param
);
2056 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2058 else if (p
->type
== relation
) {
2059 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2060 res
= compute_evalue(&p
->arr
[1], list_args
);
2061 else if (p
->size
> 2)
2062 res
= compute_evalue(&p
->arr
[2], list_args
);
2064 else if (p
->type
== partition
) {
2065 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2066 Value
*vals
= list_args
;
2069 for (i
= 0; i
< dim
; ++i
) {
2070 value_init(vals
[i
]);
2072 value_assign(vals
[i
], list_args
[i
]);
2075 for (i
= 0; i
< p
->size
/2; ++i
)
2076 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2077 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2081 for (i
= 0; i
< dim
; ++i
)
2082 value_clear(vals
[i
]);
2091 } /* compute_enode */
2093 /*************************************************/
2094 /* return the value of Ehrhart Polynomial */
2095 /* It returns a double, because since it is */
2096 /* a recursive function, some intermediate value */
2097 /* might not be integral */
2098 /*************************************************/
2100 double compute_evalue(evalue
*e
,Value
*list_args
) {
2104 if (value_notzero_p(e
->d
)) {
2105 if (value_notone_p(e
->d
))
2106 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2108 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2111 res
= compute_enode(e
->x
.p
,list_args
);
2113 } /* compute_evalue */
2116 /****************************************************/
2117 /* function compute_poly : */
2118 /* Check for the good validity domain */
2119 /* return the number of point in the Polyhedron */
2120 /* in allocated memory */
2121 /* Using the Ehrhart pseudo-polynomial */
2122 /****************************************************/
2123 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2126 /* double d; int i; */
2128 tmp
= (Value
*) malloc (sizeof(Value
));
2129 assert(tmp
!= NULL
);
2131 value_set_si(*tmp
,0);
2134 return(tmp
); /* no ehrhart polynomial */
2135 if(en
->ValidityDomain
) {
2136 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2137 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2142 return(tmp
); /* no Validity Domain */
2144 if(in_domain(en
->ValidityDomain
,list_args
)) {
2146 #ifdef EVAL_EHRHART_DEBUG
2147 Print_Domain(stdout
,en
->ValidityDomain
);
2148 print_evalue(stdout
,&en
->EP
);
2151 /* d = compute_evalue(&en->EP,list_args);
2153 printf("(double)%lf = %d\n", d, i ); */
2154 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2160 value_set_si(*tmp
,0);
2161 return(tmp
); /* no compatible domain with the arguments */
2162 } /* compute_poly */
2164 size_t value_size(Value v
) {
2165 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2166 * sizeof(v
[0]._mp_d
[0]);
2169 size_t domain_size(Polyhedron
*D
)
2172 size_t s
= sizeof(*D
);
2174 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2175 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2176 s
+= value_size(D
->Constraint
[i
][j
]);
2179 for (i = 0; i < D->NbRays; ++i)
2180 for (j = 0; j < D->Dimension+2; ++j)
2181 s += value_size(D->Ray[i][j]);
2184 return D
->next
? s
+domain_size(D
->next
) : s
;
2187 size_t enode_size(enode
*p
) {
2188 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2191 if (p
->type
== partition
)
2192 for (i
= 0; i
< p
->size
/2; ++i
) {
2193 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2194 s
+= evalue_size(&p
->arr
[2*i
+1]);
2197 for (i
= 0; i
< p
->size
; ++i
) {
2198 s
+= evalue_size(&p
->arr
[i
]);
2203 size_t evalue_size(evalue
*e
)
2205 size_t s
= sizeof(*e
);
2206 s
+= value_size(e
->d
);
2207 if (value_notzero_p(e
->d
))
2208 s
+= value_size(e
->x
.n
);
2210 s
+= enode_size(e
->x
.p
);
2214 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2216 evalue
*found
= NULL
;
2221 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2224 value_init(offset
.d
);
2225 value_init(offset
.x
.n
);
2226 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2227 value_lcm(m
, offset
.d
, &offset
.d
);
2228 value_set_si(offset
.x
.n
, 1);
2231 evalue_copy(©
, cst
);
2234 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2236 if (eequal(base
, &e
->x
.p
->arr
[0]))
2237 found
= &e
->x
.p
->arr
[0];
2239 value_set_si(offset
.x
.n
, -2);
2242 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2244 if (eequal(base
, &e
->x
.p
->arr
[0]))
2247 free_evalue_refs(cst
);
2248 free_evalue_refs(&offset
);
2251 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2252 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2257 static evalue
*find_relation_pair(evalue
*e
)
2260 evalue
*found
= NULL
;
2262 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2265 if (e
->x
.p
->type
== fractional
) {
2270 poly_denom(&e
->x
.p
->arr
[0], &m
);
2272 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2273 cst
= &cst
->x
.p
->arr
[0])
2276 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2277 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2282 i
= e
->x
.p
->type
== relation
;
2283 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2284 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2289 void evalue_mod2relation(evalue
*e
) {
2292 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2295 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2296 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2297 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2298 value_clear(e
->x
.p
->arr
[2*i
].d
);
2299 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2301 if (2*i
< e
->x
.p
->size
) {
2302 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2303 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2308 if (e
->x
.p
->size
== 0) {
2310 evalue_set_si(e
, 0, 1);
2316 while ((d
= find_relation_pair(e
)) != NULL
) {
2320 value_init(split
.d
);
2321 value_set_si(split
.d
, 0);
2322 split
.x
.p
= new_enode(relation
, 3, 0);
2323 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2324 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2326 ev
= &split
.x
.p
->arr
[0];
2327 value_set_si(ev
->d
, 0);
2328 ev
->x
.p
= new_enode(fractional
, 3, -1);
2329 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2330 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2331 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2337 free_evalue_refs(&split
);
2341 static int evalue_comp(const void * a
, const void * b
)
2343 const evalue
*e1
= *(const evalue
**)a
;
2344 const evalue
*e2
= *(const evalue
**)b
;
2345 return ecmp(e1
, e2
);
2348 void evalue_combine(evalue
*e
)
2355 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2358 NALLOC(evs
, e
->x
.p
->size
/2);
2359 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2360 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2361 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2362 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2363 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2364 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2365 value_clear(p
->arr
[2*k
].d
);
2366 value_clear(p
->arr
[2*k
+1].d
);
2367 p
->arr
[2*k
] = *(evs
[i
]-1);
2368 p
->arr
[2*k
+1] = *(evs
[i
]);
2371 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2374 value_clear((evs
[i
]-1)->d
);
2378 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2379 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2380 free_evalue_refs(evs
[i
]);
2384 for (i
= 2*k
; i
< p
->size
; ++i
)
2385 value_clear(p
->arr
[i
].d
);
2392 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2394 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2396 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2399 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2400 Polyhedron
*D
, *N
, **P
;
2403 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2410 if (D
->NbEq
<= H
->NbEq
) {
2416 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2417 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2418 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2419 reduce_evalue(&tmp
);
2420 if (value_notzero_p(tmp
.d
) ||
2421 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2424 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2425 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2428 free_evalue_refs(&tmp
);
2434 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2436 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2438 value_clear(e
->x
.p
->arr
[2*i
].d
);
2439 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2441 if (2*i
< e
->x
.p
->size
) {
2442 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2443 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2450 H
= DomainConvex(D
, 0);
2451 E
= DomainDifference(H
, D
, 0);
2453 D
= DomainDifference(H
, E
, 0);
2456 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2460 /* May change coefficients to become non-standard if fiddle is set
2461 * => reduce p afterwards to correct
2463 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2464 Matrix
**R
, int fiddle
)
2468 unsigned dim
= D
->Dimension
;
2469 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2474 assert(p
->type
== fractional
);
2476 value_set_si(T
->p
[1][dim
], 1);
2478 while (value_zero_p(pp
->d
)) {
2479 assert(pp
->x
.p
->type
== polynomial
);
2480 assert(pp
->x
.p
->size
== 2);
2481 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2482 value_division(T
->p
[0][pp
->x
.p
->pos
-1], *d
, pp
->x
.p
->arr
[1].d
);
2483 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2484 if (fiddle
&& value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2485 value_substract(pp
->x
.p
->arr
[1].x
.n
,
2486 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2487 value_multiply(T
->p
[0][pp
->x
.p
->pos
-1],
2488 T
->p
[0][pp
->x
.p
->pos
-1], pp
->x
.p
->arr
[1].x
.n
);
2489 pp
= &pp
->x
.p
->arr
[0];
2491 value_division(T
->p
[0][dim
], *d
, pp
->d
);
2492 value_multiply(T
->p
[0][dim
], T
->p
[0][dim
], pp
->x
.n
);
2493 I
= DomainImage(D
, T
, 0);
2494 H
= DomainConvex(I
, 0);
2503 static int reduce_in_domain(evalue
*e
, Polyhedron
*D
)
2513 if (value_notzero_p(e
->d
))
2518 if (p
->type
== relation
) {
2524 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
, 1);
2525 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2526 equal
= value_eq(min
, max
);
2527 mpz_cdiv_q(min
, min
, d
);
2528 mpz_fdiv_q(max
, max
, d
);
2530 if (bounded
&& value_gt(min
, max
)) {
2536 evalue_set_si(e
, 0, 1);
2539 free_evalue_refs(&(p
->arr
[1]));
2540 free_evalue_refs(&(p
->arr
[0]));
2546 return r
? r
: reduce_in_domain(e
, D
);
2547 } else if (bounded
&& equal
) {
2550 free_evalue_refs(&(p
->arr
[2]));
2553 free_evalue_refs(&(p
->arr
[0]));
2559 return reduce_in_domain(e
, D
);
2560 } else if (bounded
&& value_eq(min
, max
)) {
2561 /* zero for a single value */
2563 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2564 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2565 value_multiply(min
, min
, d
);
2566 value_substract(M
->p
[0][D
->Dimension
+1],
2567 M
->p
[0][D
->Dimension
+1], min
);
2568 E
= DomainAddConstraints(D
, M
, 0);
2574 r
= reduce_in_domain(&p
->arr
[1], E
);
2576 r
|= reduce_in_domain(&p
->arr
[2], D
);
2578 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2586 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2589 i
= p
->type
== relation
? 1 :
2590 p
->type
== fractional
? 1 : 0;
2591 for (; i
<p
->size
; i
++)
2592 r
|= reduce_in_domain(&p
->arr
[i
], D
);
2594 if (p
->type
!= fractional
) {
2595 if (r
&& p
->type
== polynomial
) {
2598 value_set_si(f
.d
, 0);
2599 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2600 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2601 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2602 reorder_terms(p
, &f
);
2613 I
= polynomial_projection(p
, D
, &d
, &T
, 1);
2615 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2616 mpz_fdiv_q(min
, min
, d
);
2617 mpz_fdiv_q(max
, max
, d
);
2618 value_substract(d
, max
, min
);
2620 if (bounded
&& value_eq(min
, max
)) {
2623 value_init(inc
.x
.n
);
2624 value_set_si(inc
.d
, 1);
2625 value_oppose(inc
.x
.n
, min
);
2626 eadd(&inc
, &p
->arr
[0]);
2627 reorder_terms(p
, &p
->arr
[0]); /* frees arr[0] */
2631 free_evalue_refs(&inc
);
2633 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
2640 value_set_si(rem
.d
, 0);
2641 rem
.x
.p
= new_enode(fractional
, 3, -1);
2642 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
2643 rem
.x
.p
->arr
[1] = p
->arr
[1];
2644 rem
.x
.p
->arr
[2] = p
->arr
[2];
2645 for (i
= 3; i
< p
->size
; ++i
)
2646 p
->arr
[i
-2] = p
->arr
[i
];
2650 value_init(inc
.x
.n
);
2651 value_set_si(inc
.d
, 1);
2652 value_oppose(inc
.x
.n
, min
);
2655 evalue_copy(&t
, &p
->arr
[0]);
2659 value_set_si(f
.d
, 0);
2660 f
.x
.p
= new_enode(fractional
, 3, -1);
2661 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
2662 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2663 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
2665 value_init(factor
.d
);
2666 evalue_set_si(&factor
, -1, 1);
2672 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2673 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
2679 free_evalue_refs(&inc
);
2680 free_evalue_refs(&t
);
2681 free_evalue_refs(&f
);
2682 free_evalue_refs(&factor
);
2683 free_evalue_refs(&rem
);
2685 reduce_in_domain(e
, D
);
2689 _reduce_evalue(&p
->arr
[0], 0, 1);
2693 value_set_si(f
.d
, 0);
2694 f
.x
.p
= new_enode(fractional
, 3, -1);
2695 value_clear(f
.x
.p
->arr
[0].d
);
2696 f
.x
.p
->arr
[0] = p
->arr
[0];
2697 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2698 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
2699 reorder_terms(p
, &f
);
2713 void evalue_range_reduction(evalue
*e
)
2716 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2719 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2720 if (reduce_in_domain(&e
->x
.p
->arr
[2*i
+1],
2721 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
2722 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
2724 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2725 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2726 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
2727 value_clear(e
->x
.p
->arr
[2*i
].d
);
2729 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2730 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2738 Enumeration
* partition2enumeration(evalue
*EP
)
2741 Enumeration
*en
, *res
= NULL
;
2743 if (EVALUE_IS_ZERO(*EP
)) {
2748 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
2749 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
2750 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
2753 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
2754 value_clear(EP
->x
.p
->arr
[2*i
].d
);
2755 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
2763 static int frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
2773 if (value_notzero_p(e
->d
))
2778 i
= p
->type
== relation
? 1 :
2779 p
->type
== fractional
? 1 : 0;
2780 for (; i
<p
->size
; i
++)
2781 r
|= frac2floor_in_domain(&p
->arr
[i
], D
);
2783 if (p
->type
!= fractional
) {
2784 if (r
&& p
->type
== polynomial
) {
2787 value_set_si(f
.d
, 0);
2788 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2789 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2790 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2791 reorder_terms(p
, &f
);
2800 I
= polynomial_projection(p
, D
, &d
, &T
, 0);
2803 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2806 assert(I
->NbEq
== 0); /* Should have been reduced */
2809 for (i
= 0; i
< I
->NbConstraints
; ++i
)
2810 if (value_pos_p(I
->Constraint
[i
][1]))
2813 assert(i
< I
->NbConstraints
);
2815 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
2816 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
2817 if (value_neg_p(min
)) {
2819 mpz_fdiv_q(min
, min
, d
);
2820 value_init(offset
.d
);
2821 value_set_si(offset
.d
, 1);
2822 value_init(offset
.x
.n
);
2823 value_oppose(offset
.x
.n
, min
);
2824 eadd(&offset
, &p
->arr
[0]);
2825 free_evalue_refs(&offset
);
2834 value_set_si(fl
.d
, 0);
2835 fl
.x
.p
= new_enode(flooring
, 3, -1);
2836 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
2837 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
2838 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
2840 eadd(&fl
, &p
->arr
[0]);
2841 reorder_terms(p
, &p
->arr
[0]);
2844 free_evalue_refs(&fl
);
2849 void evalue_frac2floor(evalue
*e
)
2852 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2855 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2856 if (frac2floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
2857 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
])))
2858 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
2861 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
2866 int nparam
= D
->Dimension
- nvar
;
2869 nr
= D
->NbConstraints
+ 2;
2870 nc
= D
->Dimension
+ 2 + 1;
2871 C
= Matrix_Alloc(nr
, nc
);
2872 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
2873 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
2874 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
2875 D
->Dimension
+ 1 - nvar
);
2880 nc
= C
->NbColumns
+ 1;
2881 C
= Matrix_Alloc(nr
, nc
);
2882 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
2883 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
2884 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
2885 oldC
->NbColumns
- 1 - nvar
);
2888 value_set_si(C
->p
[nr
-2][0], 1);
2889 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
2890 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
2892 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
2893 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
2899 static void floor2frac_r(evalue
*e
, int nvar
)
2906 if (value_notzero_p(e
->d
))
2911 assert(p
->type
== flooring
);
2912 for (i
= 1; i
< p
->size
; i
++)
2913 floor2frac_r(&p
->arr
[i
], nvar
);
2915 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
2916 assert(pp
->x
.p
->type
== polynomial
);
2917 pp
->x
.p
->pos
-= nvar
;
2921 value_set_si(f
.d
, 0);
2922 f
.x
.p
= new_enode(fractional
, 3, -1);
2923 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2924 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
2925 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
2927 eadd(&f
, &p
->arr
[0]);
2928 reorder_terms(p
, &p
->arr
[0]);
2931 free_evalue_refs(&f
);
2934 /* Convert flooring back to fractional and shift position
2935 * of the parameters by nvar
2937 static void floor2frac(evalue
*e
, int nvar
)
2939 floor2frac_r(e
, nvar
);
2943 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
2946 int nparam
= D
->Dimension
- nvar
;
2950 D
= Constraints2Polyhedron(C
, 0);
2954 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
2956 /* Double check that D was not unbounded. */
2957 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
2965 evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
2972 evalue
*factor
= NULL
;
2975 if (EVALUE_IS_ZERO(*e
))
2979 Polyhedron
*DD
= Disjoint_Domain(D
, 0, 0);
2986 res
= esum_over_domain(e
, nvar
, Q
, C
);
2989 for (Q
= DD
; Q
; Q
= DD
) {
2995 t
= esum_over_domain(e
, nvar
, Q
, C
);
3002 free_evalue_refs(t
);
3009 if (value_notzero_p(e
->d
)) {
3012 t
= esum_over_domain_cst(nvar
, D
, C
);
3014 if (!EVALUE_IS_ONE(*e
))
3020 switch (e
->x
.p
->type
) {
3022 evalue
*pp
= &e
->x
.p
->arr
[0];
3024 if (pp
->x
.p
->pos
> nvar
) {
3025 /* remainder is independent of the summated vars */
3031 floor2frac(&f
, nvar
);
3033 t
= esum_over_domain_cst(nvar
, D
, C
);
3037 free_evalue_refs(&f
);
3042 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3043 poly_denom(pp
, &row
->p
[1 + nvar
]);
3044 value_set_si(row
->p
[0], 1);
3045 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3046 pp
= &pp
->x
.p
->arr
[0]) {
3048 assert(pp
->x
.p
->type
== polynomial
);
3050 if (pos
>= 1 + nvar
)
3052 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3053 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3054 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3056 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3057 value_division(row
->p
[1 + D
->Dimension
+ 1],
3058 row
->p
[1 + D
->Dimension
+ 1],
3060 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3061 row
->p
[1 + D
->Dimension
+ 1],
3063 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3067 int pos
= e
->x
.p
->pos
;
3071 value_init(factor
->d
);
3072 value_set_si(factor
->d
, 0);
3073 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3074 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3075 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3079 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3080 for (i
= 0; i
< D
->NbRays
; ++i
)
3081 if (value_notzero_p(D
->Ray
[i
][pos
]))
3083 assert(i
< D
->NbRays
);
3084 if (value_neg_p(D
->Ray
[i
][pos
])) {
3086 value_init(factor
->d
);
3087 evalue_set_si(factor
, -1, 1);
3089 value_set_si(row
->p
[0], 1);
3090 value_set_si(row
->p
[pos
], 1);
3091 value_set_si(row
->p
[1 + nvar
], -1);
3098 i
= type_offset(e
->x
.p
);
3100 res
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
3105 evalue_copy(&cum
, factor
);
3109 for (; i
< e
->x
.p
->size
; ++i
) {
3113 C
= esum_add_constraint(nvar
, D
, C
, row
);
3119 Vector_Print(stderr, P_VALUE_FMT, row);
3121 Matrix_Print(stderr, P_VALUE_FMT, C);
3123 t
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
3132 free_evalue_refs(t
);
3135 if (factor
&& i
+1 < e
->x
.p
->size
)
3142 free_evalue_refs(factor
);
3143 free_evalue_refs(&cum
);
3155 evalue
*esum(evalue
*e
, int nvar
)
3163 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3164 evalue_copy(res
, e
);
3168 evalue_set_si(res
, 0, 1);
3170 assert(value_zero_p(e
->d
));
3171 assert(e
->x
.p
->type
== partition
);
3173 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3175 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3176 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
3178 free_evalue_refs(t
);
3187 /* Initial silly implementation */
3188 void eor(evalue
*e1
, evalue
*res
)
3194 evalue_set_si(&mone
, -1, 1);
3196 evalue_copy(&E
, res
);
3202 free_evalue_refs(&E
);
3203 free_evalue_refs(&mone
);