1 /***********************************************************************/
2 /* copyright 1997, Doran Wilde */
3 /* copyright 1997-2000, Vincent Loechner */
4 /* copyright 2003-2006, Sven Verdoolaege */
5 /* Permission is granted to copy, use, and distribute */
6 /* for any commercial or noncommercial purpose under the terms */
7 /* of the GNU General Public license, version 2, June 1991 */
8 /* (see file : LICENSE). */
9 /***********************************************************************/
16 #include <barvinok/evalue.h>
17 #include <barvinok/barvinok.h>
18 #include <barvinok/util.h>
20 #ifndef value_pmodulus
21 #define value_pmodulus(ref,val1,val2) (mpz_fdiv_r((ref),(val1),(val2)))
24 #define ALLOC(type) (type*)malloc(sizeof(type))
25 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
28 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
30 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
33 void evalue_set_si(evalue
*ev
, int n
, int d
) {
34 value_set_si(ev
->d
, d
);
36 value_set_si(ev
->x
.n
, n
);
39 void evalue_set(evalue
*ev
, Value n
, Value d
) {
40 value_assign(ev
->d
, d
);
42 value_assign(ev
->x
.n
, n
);
47 evalue
*EP
= ALLOC(evalue
);
49 evalue_set_si(EP
, 0, 1);
53 /* returns an evalue that corresponds to
57 evalue
*evalue_var(int var
)
59 evalue
*EP
= ALLOC(evalue
);
61 value_set_si(EP
->d
,0);
62 EP
->x
.p
= new_enode(polynomial
, 2, var
+ 1);
63 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
64 evalue_set_si(&EP
->x
.p
->arr
[1], 1, 1);
68 void aep_evalue(evalue
*e
, int *ref
) {
73 if (value_notzero_p(e
->d
))
74 return; /* a rational number, its already reduced */
76 return; /* hum... an overflow probably occured */
78 /* First check the components of p */
79 for (i
=0;i
<p
->size
;i
++)
80 aep_evalue(&p
->arr
[i
],ref
);
87 p
->pos
= ref
[p
->pos
-1]+1;
93 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
99 if (value_notzero_p(e
->d
))
100 return; /* a rational number, its already reduced */
102 return; /* hum... an overflow probably occured */
105 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
106 for(i
=0;i
<CT
->NbRows
-1;i
++)
107 for(j
=0;j
<CT
->NbColumns
;j
++)
108 if(value_notzero_p(CT
->p
[i
][j
])) {
113 /* Transform the references in e, using ref */
117 } /* addeliminatedparams_evalue */
119 static void addeliminatedparams_partition(enode
*p
, Matrix
*CT
, Polyhedron
*CEq
,
120 unsigned nparam
, unsigned MaxRays
)
123 assert(p
->type
== partition
);
126 for (i
= 0; i
< p
->size
/2; i
++) {
127 Polyhedron
*D
= EVALUE_DOMAIN(p
->arr
[2*i
]);
128 Polyhedron
*T
= DomainPreimage(D
, CT
, MaxRays
);
132 T
= DomainIntersection(D
, CEq
, MaxRays
);
135 EVALUE_SET_DOMAIN(p
->arr
[2*i
], T
);
139 void addeliminatedparams_enum(evalue
*e
, Matrix
*CT
, Polyhedron
*CEq
,
140 unsigned MaxRays
, unsigned nparam
)
145 if (CT
->NbRows
== CT
->NbColumns
)
148 if (EVALUE_IS_ZERO(*e
))
151 if (value_notzero_p(e
->d
)) {
154 value_set_si(res
.d
, 0);
155 res
.x
.p
= new_enode(partition
, 2, nparam
);
156 EVALUE_SET_DOMAIN(res
.x
.p
->arr
[0],
157 DomainConstraintSimplify(Polyhedron_Copy(CEq
), MaxRays
));
158 value_clear(res
.x
.p
->arr
[1].d
);
159 res
.x
.p
->arr
[1] = *e
;
167 addeliminatedparams_partition(p
, CT
, CEq
, nparam
, MaxRays
);
168 for (i
= 0; i
< p
->size
/2; i
++)
169 addeliminatedparams_evalue(&p
->arr
[2*i
+1], CT
);
172 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
180 assert(value_notzero_p(e1
->d
));
181 assert(value_notzero_p(e2
->d
));
182 value_multiply(m
, e1
->x
.n
, e2
->d
);
183 value_multiply(m2
, e2
->x
.n
, e1
->d
);
186 else if (value_gt(m
, m2
))
196 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
198 if (value_notzero_p(e1
->d
)) {
200 if (value_zero_p(e2
->d
))
202 r
= mod_rational_smaller(e1
, e2
);
203 return r
== -1 ? 0 : r
;
205 if (value_notzero_p(e2
->d
))
207 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
209 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
212 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
214 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
219 static int mod_term_smaller(const evalue
*e1
, const evalue
*e2
)
221 assert(value_zero_p(e1
->d
));
222 assert(value_zero_p(e2
->d
));
223 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
224 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
225 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
228 static void check_order(const evalue
*e
)
233 if (value_notzero_p(e
->d
))
236 switch (e
->x
.p
->type
) {
238 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
239 check_order(&e
->x
.p
->arr
[2*i
+1]);
242 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
244 if (value_notzero_p(a
->d
))
246 switch (a
->x
.p
->type
) {
248 assert(mod_term_smaller(&e
->x
.p
->arr
[0], &a
->x
.p
->arr
[0]));
257 for (i
= 0; i
< e
->x
.p
->size
; ++i
) {
259 if (value_notzero_p(a
->d
))
261 switch (a
->x
.p
->type
) {
263 assert(e
->x
.p
->pos
< a
->x
.p
->pos
);
274 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
276 if (value_notzero_p(a
->d
))
278 switch (a
->x
.p
->type
) {
289 /* Negative pos means inequality */
290 /* s is negative of substitution if m is not zero */
299 struct fixed_param
*fixed
;
304 static int relations_depth(evalue
*e
)
309 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
310 e
= &e
->x
.p
->arr
[1], ++d
);
314 static void poly_denom_not_constant(evalue
**pp
, Value
*d
)
319 while (value_zero_p(p
->d
)) {
320 assert(p
->x
.p
->type
== polynomial
);
321 assert(p
->x
.p
->size
== 2);
322 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
323 value_lcm(*d
, *d
, p
->x
.p
->arr
[1].d
);
329 static void poly_denom(evalue
*p
, Value
*d
)
331 poly_denom_not_constant(&p
, d
);
332 value_lcm(*d
, *d
, p
->d
);
335 static void realloc_substitution(struct subst
*s
, int d
)
337 struct fixed_param
*n
;
340 for (i
= 0; i
< s
->n
; ++i
)
347 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
353 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
356 /* May have been reduced already */
357 if (value_notzero_p(m
->d
))
360 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
361 assert(m
->x
.p
->size
== 3);
363 /* fractional was inverted during reduction
364 * invert it back and move constant in
366 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
367 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
368 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
369 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
370 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
371 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
372 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
373 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
374 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
375 value_set_si(m
->x
.p
->arr
[1].d
, 1);
378 /* Oops. Nested identical relations. */
379 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
382 if (s
->n
>= s
->max
) {
383 int d
= relations_depth(r
);
384 realloc_substitution(s
, d
);
388 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
389 assert(p
->x
.p
->size
== 2);
392 assert(value_pos_p(f
->x
.n
));
394 value_init(s
->fixed
[s
->n
].m
);
395 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
396 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
397 value_init(s
->fixed
[s
->n
].d
);
398 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
399 value_init(s
->fixed
[s
->n
].s
.d
);
400 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
406 static int type_offset(enode
*p
)
408 return p
->type
== fractional
? 1 :
409 p
->type
== flooring
? 1 :
410 p
->type
== relation
? 1 : 0;
413 static void reorder_terms_about(enode
*p
, evalue
*v
)
416 int offset
= type_offset(p
);
418 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
420 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
421 free_evalue_refs(&(p
->arr
[i
]));
427 static void reorder_terms(evalue
*e
)
432 assert(value_zero_p(e
->d
));
434 assert(p
->type
== fractional
); /* for now */
437 value_set_si(f
.d
, 0);
438 f
.x
.p
= new_enode(fractional
, 3, -1);
439 value_clear(f
.x
.p
->arr
[0].d
);
440 f
.x
.p
->arr
[0] = p
->arr
[0];
441 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
442 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
443 reorder_terms_about(p
, &f
);
449 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
455 if (value_notzero_p(e
->d
)) {
457 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
458 return; /* a rational number, its already reduced */
462 return; /* hum... an overflow probably occured */
464 /* First reduce the components of p */
465 add
= p
->type
== relation
;
466 for (i
=0; i
<p
->size
; i
++) {
468 add
= add_modulo_substitution(s
, e
);
470 if (i
== 0 && p
->type
==fractional
)
471 _reduce_evalue(&p
->arr
[i
], s
, 1);
473 _reduce_evalue(&p
->arr
[i
], s
, fract
);
475 if (add
&& i
== p
->size
-1) {
477 value_clear(s
->fixed
[s
->n
].m
);
478 value_clear(s
->fixed
[s
->n
].d
);
479 free_evalue_refs(&s
->fixed
[s
->n
].s
);
480 } else if (add
&& i
== 1)
481 s
->fixed
[s
->n
-1].pos
*= -1;
484 if (p
->type
==periodic
) {
486 /* Try to reduce the period */
487 for (i
=1; i
<=(p
->size
)/2; i
++) {
488 if ((p
->size
% i
)==0) {
490 /* Can we reduce the size to i ? */
492 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
493 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
496 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
500 you_lose
: /* OK, lets not do it */
505 /* Try to reduce its strength */
508 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
512 else if (p
->type
==polynomial
) {
513 for (k
= 0; s
&& k
< s
->n
; ++k
) {
514 if (s
->fixed
[k
].pos
== p
->pos
) {
515 int divide
= value_notone_p(s
->fixed
[k
].d
);
518 if (value_notzero_p(s
->fixed
[k
].m
)) {
521 assert(p
->size
== 2);
522 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
524 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
531 value_assign(d
.d
, s
->fixed
[k
].d
);
533 if (value_notzero_p(s
->fixed
[k
].m
))
534 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
536 value_set_si(d
.x
.n
, 1);
539 for (i
=p
->size
-1;i
>=1;i
--) {
540 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
542 emul(&d
, &p
->arr
[i
]);
543 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
544 free_evalue_refs(&(p
->arr
[i
]));
547 _reduce_evalue(&p
->arr
[0], s
, fract
);
550 free_evalue_refs(&d
);
556 /* Try to reduce the degree */
557 for (i
=p
->size
-1;i
>=1;i
--) {
558 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
560 /* Zero coefficient */
561 free_evalue_refs(&(p
->arr
[i
]));
566 /* Try to reduce its strength */
569 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
573 else if (p
->type
==fractional
) {
577 if (value_notzero_p(p
->arr
[0].d
)) {
579 value_assign(v
.d
, p
->arr
[0].d
);
581 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
586 evalue
*pp
= &p
->arr
[0];
587 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
588 assert(pp
->x
.p
->size
== 2);
590 /* search for exact duplicate among the modulo inequalities */
592 f
= &pp
->x
.p
->arr
[1];
593 for (k
= 0; s
&& k
< s
->n
; ++k
) {
594 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
595 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
596 value_eq(s
->fixed
[k
].m
, f
->d
) &&
597 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
604 /* replace { E/m } by { (E-1)/m } + 1/m */
609 evalue_set_si(&extra
, 1, 1);
610 value_assign(extra
.d
, g
);
611 eadd(&extra
, &v
.x
.p
->arr
[1]);
612 free_evalue_refs(&extra
);
614 /* We've been going in circles; stop now */
615 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
616 free_evalue_refs(&v
);
618 evalue_set_si(&v
, 0, 1);
623 value_set_si(v
.d
, 0);
624 v
.x
.p
= new_enode(fractional
, 3, -1);
625 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
626 value_assign(v
.x
.p
->arr
[1].d
, g
);
627 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
628 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
631 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
634 value_division(f
->d
, g
, f
->d
);
635 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
636 value_assign(f
->d
, g
);
637 value_decrement(f
->x
.n
, f
->x
.n
);
638 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
640 value_gcd(g
, f
->d
, f
->x
.n
);
641 value_division(f
->d
, f
->d
, g
);
642 value_division(f
->x
.n
, f
->x
.n
, g
);
651 /* reduction may have made this fractional arg smaller */
652 i
= reorder
? p
->size
: 1;
653 for ( ; i
< p
->size
; ++i
)
654 if (value_zero_p(p
->arr
[i
].d
) &&
655 p
->arr
[i
].x
.p
->type
== fractional
&&
656 !mod_term_smaller(e
, &p
->arr
[i
]))
660 value_set_si(v
.d
, 0);
661 v
.x
.p
= new_enode(fractional
, 3, -1);
662 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
663 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
664 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
672 evalue
*pp
= &p
->arr
[0];
675 poly_denom_not_constant(&pp
, &m
);
676 mpz_fdiv_r(r
, m
, pp
->d
);
677 if (value_notzero_p(r
)) {
679 value_set_si(v
.d
, 0);
680 v
.x
.p
= new_enode(fractional
, 3, -1);
682 value_multiply(r
, m
, pp
->x
.n
);
683 value_multiply(v
.x
.p
->arr
[1].d
, m
, pp
->d
);
684 value_init(v
.x
.p
->arr
[1].x
.n
);
685 mpz_fdiv_r(v
.x
.p
->arr
[1].x
.n
, r
, pp
->d
);
686 mpz_fdiv_q(r
, r
, pp
->d
);
688 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
689 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
691 while (value_zero_p(pp
->d
))
692 pp
= &pp
->x
.p
->arr
[0];
694 value_assign(pp
->d
, m
);
695 value_assign(pp
->x
.n
, r
);
697 value_gcd(r
, pp
->d
, pp
->x
.n
);
698 value_division(pp
->d
, pp
->d
, r
);
699 value_division(pp
->x
.n
, pp
->x
.n
, r
);
712 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
713 pp
= &pp
->x
.p
->arr
[0]) {
714 f
= &pp
->x
.p
->arr
[1];
715 assert(value_pos_p(f
->d
));
716 mpz_mul_ui(twice
, f
->x
.n
, 2);
717 if (value_lt(twice
, f
->d
))
719 if (value_eq(twice
, f
->d
))
727 value_set_si(v
.d
, 0);
728 v
.x
.p
= new_enode(fractional
, 3, -1);
729 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
730 poly_denom(&p
->arr
[0], &twice
);
731 value_assign(v
.x
.p
->arr
[1].d
, twice
);
732 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
733 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
734 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
736 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
737 pp
= &pp
->x
.p
->arr
[0]) {
738 f
= &pp
->x
.p
->arr
[1];
739 value_oppose(f
->x
.n
, f
->x
.n
);
740 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
742 value_division(pp
->d
, twice
, pp
->d
);
743 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
744 value_assign(pp
->d
, twice
);
745 value_oppose(pp
->x
.n
, pp
->x
.n
);
746 value_decrement(pp
->x
.n
, pp
->x
.n
);
747 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
749 /* Maybe we should do this during reduction of
752 value_gcd(twice
, pp
->d
, pp
->x
.n
);
753 value_division(pp
->d
, pp
->d
, twice
);
754 value_division(pp
->x
.n
, pp
->x
.n
, twice
);
764 reorder_terms_about(p
, &v
);
765 _reduce_evalue(&p
->arr
[1], s
, fract
);
768 /* Try to reduce the degree */
769 for (i
=p
->size
-1;i
>=2;i
--) {
770 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
772 /* Zero coefficient */
773 free_evalue_refs(&(p
->arr
[i
]));
778 /* Try to reduce its strength */
781 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
782 free_evalue_refs(&(p
->arr
[0]));
786 else if (p
->type
== flooring
) {
787 /* Try to reduce the degree */
788 for (i
=p
->size
-1;i
>=2;i
--) {
789 if (!EVALUE_IS_ZERO(p
->arr
[i
]))
791 /* Zero coefficient */
792 free_evalue_refs(&(p
->arr
[i
]));
797 /* Try to reduce its strength */
800 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
801 free_evalue_refs(&(p
->arr
[0]));
805 else if (p
->type
== relation
) {
806 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
807 free_evalue_refs(&(p
->arr
[2]));
808 free_evalue_refs(&(p
->arr
[0]));
815 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
816 free_evalue_refs(&(p
->arr
[2]));
819 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
820 free_evalue_refs(&(p
->arr
[1]));
821 free_evalue_refs(&(p
->arr
[0]));
822 evalue_set_si(e
, 0, 1);
829 /* Relation was reduced by means of an identical
830 * inequality => remove
832 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
835 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
836 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
838 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
840 free_evalue_refs(&(p
->arr
[2]));
844 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
846 evalue_set_si(e
, 0, 1);
847 free_evalue_refs(&(p
->arr
[1]));
849 free_evalue_refs(&(p
->arr
[0]));
855 } /* reduce_evalue */
857 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
862 for (k
= 0; k
< dim
; ++k
)
863 if (value_notzero_p(row
[k
+1]))
866 Vector_Normalize_Positive(row
+1, dim
+1, k
);
867 assert(s
->n
< s
->max
);
868 value_init(s
->fixed
[s
->n
].d
);
869 value_init(s
->fixed
[s
->n
].m
);
870 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
871 s
->fixed
[s
->n
].pos
= k
+1;
872 value_set_si(s
->fixed
[s
->n
].m
, 0);
873 r
= &s
->fixed
[s
->n
].s
;
875 for (l
= k
+1; l
< dim
; ++l
)
876 if (value_notzero_p(row
[l
+1])) {
877 value_set_si(r
->d
, 0);
878 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
879 value_init(r
->x
.p
->arr
[1].x
.n
);
880 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
881 value_set_si(r
->x
.p
->arr
[1].d
, 1);
885 value_oppose(r
->x
.n
, row
[dim
+1]);
886 value_set_si(r
->d
, 1);
890 static void _reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
, struct subst
*s
)
893 Polyhedron
*orig
= D
;
898 D
= DomainConvex(D
, 0);
899 /* We don't perform any substitutions if the domain is a union.
900 * We may therefore miss out on some possible simplifications,
901 * e.g., if a variable is always even in the whole union,
902 * while there is a relation in the evalue that evaluates
903 * to zero for even values of the variable.
905 if (!D
->next
&& D
->NbEq
) {
909 realloc_substitution(s
, dim
);
911 int d
= relations_depth(e
);
913 NALLOC(s
->fixed
, s
->max
);
916 for (j
= 0; j
< D
->NbEq
; ++j
)
917 add_substitution(s
, D
->Constraint
[j
], dim
);
921 _reduce_evalue(e
, s
, 0);
924 for (j
= 0; j
< s
->n
; ++j
) {
925 value_clear(s
->fixed
[j
].d
);
926 value_clear(s
->fixed
[j
].m
);
927 free_evalue_refs(&s
->fixed
[j
].s
);
932 void reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
)
934 struct subst s
= { NULL
, 0, 0 };
936 if (EVALUE_IS_ZERO(*e
))
940 evalue_set_si(e
, 0, 1);
943 _reduce_evalue_in_domain(e
, D
, &s
);
948 void reduce_evalue (evalue
*e
) {
949 struct subst s
= { NULL
, 0, 0 };
951 if (value_notzero_p(e
->d
))
952 return; /* a rational number, its already reduced */
954 if (e
->x
.p
->type
== partition
) {
957 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
958 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
960 /* This shouldn't really happen;
961 * Empty domains should not be added.
963 POL_ENSURE_VERTICES(D
);
965 _reduce_evalue_in_domain(&e
->x
.p
->arr
[2*i
+1], D
, &s
);
967 if (emptyQ(D
) || EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
968 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
969 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
970 value_clear(e
->x
.p
->arr
[2*i
].d
);
972 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
973 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
977 if (e
->x
.p
->size
== 0) {
979 evalue_set_si(e
, 0, 1);
982 _reduce_evalue(e
, &s
, 0);
987 static void print_evalue_r(FILE *DST
, const evalue
*e
, const char *const *pname
)
989 if(value_notzero_p(e
->d
)) {
990 if(value_notone_p(e
->d
)) {
991 value_print(DST
,VALUE_FMT
,e
->x
.n
);
993 value_print(DST
,VALUE_FMT
,e
->d
);
996 value_print(DST
,VALUE_FMT
,e
->x
.n
);
1000 print_enode(DST
,e
->x
.p
,pname
);
1002 } /* print_evalue */
1004 void print_evalue(FILE *DST
, const evalue
*e
, const char * const *pname
)
1006 print_evalue_r(DST
, e
, pname
);
1007 if (value_notzero_p(e
->d
))
1011 void print_enode(FILE *DST
, enode
*p
, const char *const *pname
)
1016 fprintf(DST
, "NULL");
1022 for (i
=0; i
<p
->size
; i
++) {
1023 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1027 fprintf(DST
, " }\n");
1031 for (i
=p
->size
-1; i
>=0; i
--) {
1032 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1033 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
1035 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
1037 fprintf(DST
, " )\n");
1041 for (i
=0; i
<p
->size
; i
++) {
1042 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1043 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
1045 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
1050 for (i
=p
->size
-1; i
>=1; i
--) {
1051 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1053 fprintf(DST
, " * ");
1054 fprintf(DST
, p
->type
== flooring
? "[" : "{");
1055 print_evalue_r(DST
, &p
->arr
[0], pname
);
1056 fprintf(DST
, p
->type
== flooring
? "]" : "}");
1058 fprintf(DST
, "^%d + ", i
-1);
1060 fprintf(DST
, " + ");
1063 fprintf(DST
, " )\n");
1067 print_evalue_r(DST
, &p
->arr
[0], pname
);
1068 fprintf(DST
, "= 0 ] * \n");
1069 print_evalue_r(DST
, &p
->arr
[1], pname
);
1071 fprintf(DST
, " +\n [ ");
1072 print_evalue_r(DST
, &p
->arr
[0], pname
);
1073 fprintf(DST
, "!= 0 ] * \n");
1074 print_evalue_r(DST
, &p
->arr
[2], pname
);
1078 char **new_names
= NULL
;
1079 const char *const *names
= pname
;
1080 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
1081 if (!pname
|| p
->pos
< maxdim
) {
1082 new_names
= ALLOCN(char *, maxdim
);
1083 for (i
= 0; i
< p
->pos
; ++i
) {
1085 new_names
[i
] = (char *)pname
[i
];
1087 new_names
[i
] = ALLOCN(char, 10);
1088 snprintf(new_names
[i
], 10, "%c", 'P'+i
);
1091 for ( ; i
< maxdim
; ++i
) {
1092 new_names
[i
] = ALLOCN(char, 10);
1093 snprintf(new_names
[i
], 10, "_p%d", i
);
1095 names
= (const char**)new_names
;
1098 for (i
=0; i
<p
->size
/2; i
++) {
1099 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
1100 print_evalue_r(DST
, &p
->arr
[2*i
+1], names
);
1101 if (value_notzero_p(p
->arr
[2*i
+1].d
))
1105 if (!pname
|| p
->pos
< maxdim
) {
1106 for (i
= pname
? p
->pos
: 0; i
< maxdim
; ++i
)
1120 * 0 if toplevels of e1 and e2 are at the same level
1121 * <0 if toplevel of e1 should be outside of toplevel of e2
1122 * >0 if toplevel of e2 should be outside of toplevel of e1
1124 static int evalue_level_cmp(const evalue
*e1
, const evalue
*e2
)
1126 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
))
1128 if (value_notzero_p(e1
->d
))
1130 if (value_notzero_p(e2
->d
))
1132 if (e1
->x
.p
->type
== partition
&& e2
->x
.p
->type
== partition
)
1134 if (e1
->x
.p
->type
== partition
)
1136 if (e2
->x
.p
->type
== partition
)
1138 if (e1
->x
.p
->type
== relation
&& e2
->x
.p
->type
== relation
) {
1139 if (eequal(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]))
1141 if (mod_term_smaller(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]))
1146 if (e1
->x
.p
->type
== relation
)
1148 if (e2
->x
.p
->type
== relation
)
1150 if (e1
->x
.p
->type
== polynomial
&& e2
->x
.p
->type
== polynomial
)
1151 return e1
->x
.p
->pos
- e2
->x
.p
->pos
;
1152 if (e1
->x
.p
->type
== polynomial
)
1154 if (e2
->x
.p
->type
== polynomial
)
1156 if (e1
->x
.p
->type
== periodic
&& e2
->x
.p
->type
== periodic
)
1157 return e1
->x
.p
->pos
- e2
->x
.p
->pos
;
1158 assert(e1
->x
.p
->type
!= periodic
);
1159 assert(e2
->x
.p
->type
!= periodic
);
1160 assert(e1
->x
.p
->type
== e2
->x
.p
->type
);
1161 if (eequal(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]))
1163 if (mod_term_smaller(e1
, e2
))
1169 static void eadd_rev(const evalue
*e1
, evalue
*res
)
1173 evalue_copy(&ev
, e1
);
1175 free_evalue_refs(res
);
1179 static void eadd_rev_cst(const evalue
*e1
, evalue
*res
)
1183 evalue_copy(&ev
, e1
);
1184 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
1185 free_evalue_refs(res
);
1189 static int is_zero_on(evalue
*e
, Polyhedron
*D
)
1194 tmp
.x
.p
= new_enode(partition
, 2, D
->Dimension
);
1195 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Domain_Copy(D
));
1196 evalue_copy(&tmp
.x
.p
->arr
[1], e
);
1197 reduce_evalue(&tmp
);
1198 is_zero
= EVALUE_IS_ZERO(tmp
);
1199 free_evalue_refs(&tmp
);
1203 struct section
{ Polyhedron
* D
; evalue E
; };
1205 void eadd_partitions(const evalue
*e1
, evalue
*res
)
1210 s
= (struct section
*)
1211 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
1212 sizeof(struct section
));
1214 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1215 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1216 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1219 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1220 assert(res
->x
.p
->size
>= 2);
1221 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1222 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
1224 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
1226 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1231 fd
= DomainConstraintSimplify(fd
, 0);
1236 /* See if we can extend one of the domains in res to cover fd */
1237 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1238 if (is_zero_on(&res
->x
.p
->arr
[2*i
+1], fd
))
1240 if (i
< res
->x
.p
->size
/2) {
1241 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
1242 DomainConcat(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
])));
1245 value_init(s
[n
].E
.d
);
1246 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
1250 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1251 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
1252 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1254 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1255 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1256 d
= DomainConstraintSimplify(d
, 0);
1262 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
1263 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1265 value_init(s
[n
].E
.d
);
1266 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1267 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1268 fd
= DomainConstraintSimplify(fd
, 0);
1269 if (!emptyQ(fd
) && is_zero_on(&e1
->x
.p
->arr
[2*j
+1], fd
)) {
1270 d
= DomainConcat(fd
, d
);
1271 fd
= Empty_Polyhedron(fd
->Dimension
);
1277 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
1281 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1284 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1285 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1286 value_clear(res
->x
.p
->arr
[2*i
].d
);
1291 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1292 for (j
= 0; j
< n
; ++j
) {
1293 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1294 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1295 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1301 static void explicit_complement(evalue
*res
)
1303 enode
*rel
= new_enode(relation
, 3, 0);
1305 value_clear(rel
->arr
[0].d
);
1306 rel
->arr
[0] = res
->x
.p
->arr
[0];
1307 value_clear(rel
->arr
[1].d
);
1308 rel
->arr
[1] = res
->x
.p
->arr
[1];
1309 value_set_si(rel
->arr
[2].d
, 1);
1310 value_init(rel
->arr
[2].x
.n
);
1311 value_set_si(rel
->arr
[2].x
.n
, 0);
1316 static void reduce_constant(evalue
*e
)
1321 value_gcd(g
, e
->x
.n
, e
->d
);
1322 if (value_notone_p(g
)) {
1323 value_division(e
->d
, e
->d
,g
);
1324 value_division(e
->x
.n
, e
->x
.n
,g
);
1329 /* Add two rational numbers */
1330 static void eadd_rationals(const evalue
*e1
, evalue
*res
)
1332 if (value_eq(e1
->d
, res
->d
))
1333 value_addto(res
->x
.n
, res
->x
.n
, e1
->x
.n
);
1335 value_multiply(res
->x
.n
, res
->x
.n
, e1
->d
);
1336 value_addmul(res
->x
.n
, e1
->x
.n
, res
->d
);
1337 value_multiply(res
->d
,e1
->d
,res
->d
);
1339 reduce_constant(res
);
1342 static void eadd_relations(const evalue
*e1
, evalue
*res
)
1346 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1347 explicit_complement(res
);
1348 for (i
= 1; i
< e1
->x
.p
->size
; ++i
)
1349 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1352 static void eadd_arrays(const evalue
*e1
, evalue
*res
, int n
)
1356 // add any element in e1 to the corresponding element in res
1357 i
= type_offset(res
->x
.p
);
1359 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1361 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1364 static void eadd_poly(const evalue
*e1
, evalue
*res
)
1366 if (e1
->x
.p
->size
> res
->x
.p
->size
)
1369 eadd_arrays(e1
, res
, e1
->x
.p
->size
);
1372 static void eadd_periodics(const evalue
*e1
, evalue
*res
)
1379 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1380 eadd_arrays(e1
, res
, e1
->x
.p
->size
);
1383 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1384 * of the sizes of e1 and res, then to copy res periodicaly in ne, after
1385 * to add periodicaly elements of e1 to elements of ne, and finaly to
1393 value_set_si(ex
, e1
->x
.p
->size
);
1394 value_set_si(ey
, res
->x
.p
->size
);
1395 value_lcm(ep
, ex
, ey
);
1396 p
= (int)mpz_get_si(ep
);
1397 ne
= (evalue
*) malloc(sizeof(evalue
));
1399 value_set_si(ne
->d
, 0);
1401 ne
->x
.p
= new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1402 for (i
= 0; i
< p
; i
++)
1403 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1404 for (i
= 0; i
< p
; i
++)
1405 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1407 value_assign(res
->d
, ne
->d
);
1414 void eadd(const evalue
*e1
, evalue
*res
)
1418 if (EVALUE_IS_ZERO(*e1
))
1421 if (EVALUE_IS_ZERO(*res
)) {
1422 if (value_notzero_p(e1
->d
)) {
1423 value_assign(res
->d
, e1
->d
);
1424 value_assign(res
->x
.n
, e1
->x
.n
);
1426 value_clear(res
->x
.n
);
1427 value_set_si(res
->d
, 0);
1428 res
->x
.p
= ecopy(e1
->x
.p
);
1433 cmp
= evalue_level_cmp(res
, e1
);
1435 switch (e1
->x
.p
->type
) {
1439 eadd_rev_cst(e1
, res
);
1444 } else if (cmp
== 0) {
1445 if (value_notzero_p(e1
->d
)) {
1446 eadd_rationals(e1
, res
);
1448 switch (e1
->x
.p
->type
) {
1450 eadd_partitions(e1
, res
);
1453 eadd_relations(e1
, res
);
1456 assert(e1
->x
.p
->size
== res
->x
.p
->size
);
1463 eadd_periodics(e1
, res
);
1471 switch (res
->x
.p
->type
) {
1475 /* Add to the constant term of a polynomial */
1476 eadd(e1
, &res
->x
.p
->arr
[type_offset(res
->x
.p
)]);
1479 /* Add to all elements of a periodic number */
1480 for (i
= 0; i
< res
->x
.p
->size
; i
++)
1481 eadd(e1
, &res
->x
.p
->arr
[i
]);
1484 fprintf(stderr
, "eadd: cannot add const with vector\n");
1489 /* Create (zero) complement if needed */
1490 if (res
->x
.p
->size
< 3)
1491 explicit_complement(res
);
1492 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1493 eadd(e1
, &res
->x
.p
->arr
[i
]);
1501 static void emul_rev(const evalue
*e1
, evalue
*res
)
1505 evalue_copy(&ev
, e1
);
1507 free_evalue_refs(res
);
1511 static void emul_poly(const evalue
*e1
, evalue
*res
)
1513 int i
, j
, offset
= type_offset(res
->x
.p
);
1516 int size
= (e1
->x
.p
->size
+ res
->x
.p
->size
- offset
- 1);
1518 p
= new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1520 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1521 if (!EVALUE_IS_ZERO(e1
->x
.p
->arr
[i
]))
1524 /* special case pure power */
1525 if (i
== e1
->x
.p
->size
-1) {
1527 value_clear(p
->arr
[0].d
);
1528 p
->arr
[0] = res
->x
.p
->arr
[0];
1530 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1531 evalue_set_si(&p
->arr
[i
], 0, 1);
1532 for (i
= offset
; i
< res
->x
.p
->size
; ++i
) {
1533 value_clear(p
->arr
[i
+e1
->x
.p
->size
-offset
-1].d
);
1534 p
->arr
[i
+e1
->x
.p
->size
-offset
-1] = res
->x
.p
->arr
[i
];
1535 emul(&e1
->x
.p
->arr
[e1
->x
.p
->size
-1],
1536 &p
->arr
[i
+e1
->x
.p
->size
-offset
-1]);
1544 value_set_si(tmp
.d
,0);
1547 evalue_copy(&p
->arr
[0], &e1
->x
.p
->arr
[0]);
1548 for (i
= offset
; i
< e1
->x
.p
->size
; i
++) {
1549 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1550 emul(&res
->x
.p
->arr
[offset
], &tmp
.x
.p
->arr
[i
]);
1553 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1554 for (i
= offset
+1; i
<res
->x
.p
->size
; i
++)
1555 for (j
= offset
; j
<e1
->x
.p
->size
; j
++) {
1558 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1559 emul(&res
->x
.p
->arr
[i
], &ev
);
1560 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-offset
]);
1561 free_evalue_refs(&ev
);
1563 free_evalue_refs(res
);
1567 void emul_partitions(const evalue
*e1
, evalue
*res
)
1572 s
= (struct section
*)
1573 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1574 sizeof(struct section
));
1576 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1577 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1578 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1581 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1582 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1583 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1584 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1585 d
= DomainConstraintSimplify(d
, 0);
1591 /* This code is only needed because the partitions
1592 are not true partitions.
1594 for (k
= 0; k
< n
; ++k
) {
1595 if (DomainIncludes(s
[k
].D
, d
))
1597 if (DomainIncludes(d
, s
[k
].D
)) {
1598 Domain_Free(s
[k
].D
);
1599 free_evalue_refs(&s
[k
].E
);
1610 value_init(s
[n
].E
.d
);
1611 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1612 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1616 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1617 value_clear(res
->x
.p
->arr
[2*i
].d
);
1618 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1623 evalue_set_si(res
, 0, 1);
1625 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1626 for (j
= 0; j
< n
; ++j
) {
1627 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1628 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1629 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1636 /* Product of two rational numbers */
1637 static void emul_rationals(const evalue
*e1
, evalue
*res
)
1639 value_multiply(res
->d
, e1
->d
, res
->d
);
1640 value_multiply(res
->x
.n
, e1
->x
.n
, res
->x
.n
);
1641 reduce_constant(res
);
1644 static void emul_relations(const evalue
*e1
, evalue
*res
)
1648 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3) {
1649 free_evalue_refs(&res
->x
.p
->arr
[2]);
1652 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1653 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1656 static void emul_periodics(const evalue
*e1
, evalue
*res
)
1663 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1664 /* Product of two periodics of the same parameter and period */
1665 for (i
= 0; i
< res
->x
.p
->size
; i
++)
1666 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1670 /* Product of two periodics of the same parameter and different periods */
1675 iy
= res
->x
.p
->size
;
1676 value_set_si(x
, e1
->x
.p
->size
);
1677 value_set_si(y
, res
->x
.p
->size
);
1679 lcm
= (int)mpz_get_si(z
);
1680 newp
= (evalue
*) malloc(sizeof(evalue
));
1681 value_init(newp
->d
);
1682 value_set_si(newp
->d
, 0);
1683 newp
->x
.p
= new_enode(periodic
, lcm
, e1
->x
.p
->pos
);
1684 for (i
= 0; i
< lcm
; i
++)
1685 evalue_copy(&newp
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%iy
]);
1686 for (i
= 0; i
< lcm
; i
++)
1687 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1688 value_assign(res
->d
, newp
->d
);
1689 res
->x
.p
= newp
->x
.p
;
1695 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1697 static void emul_fractionals(const evalue
*e1
, evalue
*res
)
1701 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1702 if (!value_two_p(d
.d
))
1707 value_set_si(d
.x
.n
, 1);
1708 /* { x }^2 == { x }/2 */
1709 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1710 assert(e1
->x
.p
->size
== 3);
1711 assert(res
->x
.p
->size
== 3);
1713 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1715 eadd(&res
->x
.p
->arr
[1], &tmp
);
1716 emul(&e1
->x
.p
->arr
[2], &tmp
);
1717 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[1]);
1718 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[2]);
1719 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1720 free_evalue_refs(&tmp
);
1726 /* Computes the product of two evalues "e1" and "res" and puts
1727 * the result in "res". You need to make a copy of "res"
1728 * before calling this function if you still need it afterward.
1729 * The vector type of evalues is not treated here
1731 void emul(const evalue
*e1
, evalue
*res
)
1735 assert(!(value_zero_p(e1
->d
) && e1
->x
.p
->type
== evector
));
1736 assert(!(value_zero_p(res
->d
) && res
->x
.p
->type
== evector
));
1738 if (EVALUE_IS_ZERO(*res
))
1741 if (EVALUE_IS_ONE(*e1
))
1744 if (EVALUE_IS_ZERO(*e1
)) {
1745 if (value_notzero_p(res
->d
)) {
1746 value_assign(res
->d
, e1
->d
);
1747 value_assign(res
->x
.n
, e1
->x
.n
);
1749 free_evalue_refs(res
);
1751 evalue_set_si(res
, 0, 1);
1756 cmp
= evalue_level_cmp(res
, e1
);
1759 } else if (cmp
== 0) {
1760 if (value_notzero_p(e1
->d
)) {
1761 emul_rationals(e1
, res
);
1763 switch (e1
->x
.p
->type
) {
1765 emul_partitions(e1
, res
);
1768 emul_relations(e1
, res
);
1775 emul_periodics(e1
, res
);
1778 emul_fractionals(e1
, res
);
1784 switch (res
->x
.p
->type
) {
1786 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1787 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1794 for (i
= type_offset(res
->x
.p
); i
< res
->x
.p
->size
; ++i
)
1795 emul(e1
, &res
->x
.p
->arr
[i
]);
1801 /* Frees mask content ! */
1802 void emask(evalue
*mask
, evalue
*res
) {
1809 if (EVALUE_IS_ZERO(*res
)) {
1810 free_evalue_refs(mask
);
1814 assert(value_zero_p(mask
->d
));
1815 assert(mask
->x
.p
->type
== partition
);
1816 assert(value_zero_p(res
->d
));
1817 assert(res
->x
.p
->type
== partition
);
1818 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1819 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1820 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1821 pos
= res
->x
.p
->pos
;
1823 s
= (struct section
*)
1824 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1825 sizeof(struct section
));
1829 evalue_set_si(&mone
, -1, 1);
1832 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1833 assert(mask
->x
.p
->size
>= 2);
1834 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1835 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1837 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1839 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1848 value_init(s
[n
].E
.d
);
1849 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1853 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1854 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1857 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1858 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1859 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1860 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1862 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1863 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1869 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1870 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1872 value_init(s
[n
].E
.d
);
1873 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1874 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1880 /* Just ignore; this may have been previously masked off */
1882 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1886 free_evalue_refs(&mone
);
1887 free_evalue_refs(mask
);
1888 free_evalue_refs(res
);
1891 evalue_set_si(res
, 0, 1);
1893 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1894 for (j
= 0; j
< n
; ++j
) {
1895 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1896 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1897 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1904 void evalue_copy(evalue
*dst
, const evalue
*src
)
1906 value_assign(dst
->d
, src
->d
);
1907 if(value_notzero_p(src
->d
)) {
1908 value_init(dst
->x
.n
);
1909 value_assign(dst
->x
.n
, src
->x
.n
);
1911 dst
->x
.p
= ecopy(src
->x
.p
);
1914 evalue
*evalue_dup(const evalue
*e
)
1916 evalue
*res
= ALLOC(evalue
);
1918 evalue_copy(res
, e
);
1922 enode
*new_enode(enode_type type
,int size
,int pos
) {
1928 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1931 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1935 for(i
=0; i
<size
; i
++) {
1936 value_init(res
->arr
[i
].d
);
1937 value_set_si(res
->arr
[i
].d
,0);
1938 res
->arr
[i
].x
.p
= 0;
1943 enode
*ecopy(enode
*e
) {
1948 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1949 for(i
=0;i
<e
->size
;++i
) {
1950 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1951 if(value_zero_p(res
->arr
[i
].d
))
1952 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1953 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1954 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1956 value_init(res
->arr
[i
].x
.n
);
1957 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1963 int ecmp(const evalue
*e1
, const evalue
*e2
)
1969 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
1973 value_multiply(m
, e1
->x
.n
, e2
->d
);
1974 value_multiply(m2
, e2
->x
.n
, e1
->d
);
1976 if (value_lt(m
, m2
))
1978 else if (value_gt(m
, m2
))
1988 if (value_notzero_p(e1
->d
))
1990 if (value_notzero_p(e2
->d
))
1996 if (p1
->type
!= p2
->type
)
1997 return p1
->type
- p2
->type
;
1998 if (p1
->pos
!= p2
->pos
)
1999 return p1
->pos
- p2
->pos
;
2000 if (p1
->size
!= p2
->size
)
2001 return p1
->size
- p2
->size
;
2003 for (i
= p1
->size
-1; i
>= 0; --i
)
2004 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
2010 int eequal(const evalue
*e1
, const evalue
*e2
)
2015 if (value_ne(e1
->d
,e2
->d
))
2018 /* e1->d == e2->d */
2019 if (value_notzero_p(e1
->d
)) {
2020 if (value_ne(e1
->x
.n
,e2
->x
.n
))
2023 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2027 /* e1->d == e2->d == 0 */
2030 if (p1
->type
!= p2
->type
) return 0;
2031 if (p1
->size
!= p2
->size
) return 0;
2032 if (p1
->pos
!= p2
->pos
) return 0;
2033 for (i
=0; i
<p1
->size
; i
++)
2034 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
2039 void free_evalue_refs(evalue
*e
) {
2044 if (EVALUE_IS_DOMAIN(*e
)) {
2045 Domain_Free(EVALUE_DOMAIN(*e
));
2048 } else if (value_pos_p(e
->d
)) {
2050 /* 'e' stores a constant */
2052 value_clear(e
->x
.n
);
2055 assert(value_zero_p(e
->d
));
2058 if (!p
) return; /* null pointer */
2059 for (i
=0; i
<p
->size
; i
++) {
2060 free_evalue_refs(&(p
->arr
[i
]));
2064 } /* free_evalue_refs */
2066 void evalue_free(evalue
*e
)
2068 free_evalue_refs(e
);
2072 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
2073 Vector
* val
, evalue
*res
)
2075 unsigned nparam
= periods
->Size
;
2078 double d
= compute_evalue(e
, val
->p
);
2079 d
*= VALUE_TO_DOUBLE(m
);
2084 value_assign(res
->d
, m
);
2085 value_init(res
->x
.n
);
2086 value_set_double(res
->x
.n
, d
);
2087 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
2090 if (value_one_p(periods
->p
[p
]))
2091 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
2096 value_assign(tmp
, periods
->p
[p
]);
2097 value_set_si(res
->d
, 0);
2098 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
2100 value_decrement(tmp
, tmp
);
2101 value_assign(val
->p
[p
], tmp
);
2102 mod2table_r(e
, periods
, m
, p
+1, val
,
2103 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
2104 } while (value_pos_p(tmp
));
2110 static void rel2table(evalue
*e
, int zero
)
2112 if (value_pos_p(e
->d
)) {
2113 if (value_zero_p(e
->x
.n
) == zero
)
2114 value_set_si(e
->x
.n
, 1);
2116 value_set_si(e
->x
.n
, 0);
2117 value_set_si(e
->d
, 1);
2120 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
2121 rel2table(&e
->x
.p
->arr
[i
], zero
);
2125 void evalue_mod2table(evalue
*e
, int nparam
)
2130 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2133 for (i
=0; i
<p
->size
; i
++) {
2134 evalue_mod2table(&(p
->arr
[i
]), nparam
);
2136 if (p
->type
== relation
) {
2141 evalue_copy(©
, &p
->arr
[0]);
2143 rel2table(&p
->arr
[0], 1);
2144 emul(&p
->arr
[0], &p
->arr
[1]);
2146 rel2table(©
, 0);
2147 emul(©
, &p
->arr
[2]);
2148 eadd(&p
->arr
[2], &p
->arr
[1]);
2149 free_evalue_refs(&p
->arr
[2]);
2150 free_evalue_refs(©
);
2152 free_evalue_refs(&p
->arr
[0]);
2156 } else if (p
->type
== fractional
) {
2157 Vector
*periods
= Vector_Alloc(nparam
);
2158 Vector
*val
= Vector_Alloc(nparam
);
2164 value_set_si(tmp
, 1);
2165 Vector_Set(periods
->p
, 1, nparam
);
2166 Vector_Set(val
->p
, 0, nparam
);
2167 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
2170 assert(p
->type
== polynomial
);
2171 assert(p
->size
== 2);
2172 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
2173 value_lcm(tmp
, tmp
, p
->arr
[1].d
);
2175 value_lcm(tmp
, tmp
, ev
->d
);
2177 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
2180 evalue_set_si(&res
, 0, 1);
2181 /* Compute the polynomial using Horner's rule */
2182 for (i
=p
->size
-1;i
>1;i
--) {
2183 eadd(&p
->arr
[i
], &res
);
2186 eadd(&p
->arr
[1], &res
);
2188 free_evalue_refs(e
);
2189 free_evalue_refs(&EP
);
2194 Vector_Free(periods
);
2196 } /* evalue_mod2table */
2198 /********************************************************/
2199 /* function in domain */
2200 /* check if the parameters in list_args */
2201 /* verifies the constraints of Domain P */
2202 /********************************************************/
2203 int in_domain(Polyhedron
*P
, Value
*list_args
)
2206 Value v
; /* value of the constraint of a row when
2207 parameters are instantiated*/
2211 for (row
= 0; row
< P
->NbConstraints
; row
++) {
2212 Inner_Product(P
->Constraint
[row
]+1, list_args
, P
->Dimension
, &v
);
2213 value_addto(v
, v
, P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
2214 if (value_neg_p(v
) ||
2215 value_zero_p(P
->Constraint
[row
][0]) && value_notzero_p(v
)) {
2222 return in
|| (P
->next
&& in_domain(P
->next
, list_args
));
2225 /****************************************************/
2226 /* function compute enode */
2227 /* compute the value of enode p with parameters */
2228 /* list "list_args */
2229 /* compute the polynomial or the periodic */
2230 /****************************************************/
2232 static double compute_enode(enode
*p
, Value
*list_args
) {
2244 if (p
->type
== polynomial
) {
2246 value_assign(param
,list_args
[p
->pos
-1]);
2248 /* Compute the polynomial using Horner's rule */
2249 for (i
=p
->size
-1;i
>0;i
--) {
2250 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2251 res
*=VALUE_TO_DOUBLE(param
);
2253 res
+=compute_evalue(&p
->arr
[0],list_args
);
2255 else if (p
->type
== fractional
) {
2256 double d
= compute_evalue(&p
->arr
[0], list_args
);
2257 d
-= floor(d
+1e-10);
2259 /* Compute the polynomial using Horner's rule */
2260 for (i
=p
->size
-1;i
>1;i
--) {
2261 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2264 res
+=compute_evalue(&p
->arr
[1],list_args
);
2266 else if (p
->type
== flooring
) {
2267 double d
= compute_evalue(&p
->arr
[0], list_args
);
2270 /* Compute the polynomial using Horner's rule */
2271 for (i
=p
->size
-1;i
>1;i
--) {
2272 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2275 res
+=compute_evalue(&p
->arr
[1],list_args
);
2277 else if (p
->type
== periodic
) {
2278 value_assign(m
,list_args
[p
->pos
-1]);
2280 /* Choose the right element of the periodic */
2281 value_set_si(param
,p
->size
);
2282 value_pmodulus(m
,m
,param
);
2283 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2285 else if (p
->type
== relation
) {
2286 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2287 res
= compute_evalue(&p
->arr
[1], list_args
);
2288 else if (p
->size
> 2)
2289 res
= compute_evalue(&p
->arr
[2], list_args
);
2291 else if (p
->type
== partition
) {
2292 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2293 Value
*vals
= list_args
;
2296 for (i
= 0; i
< dim
; ++i
) {
2297 value_init(vals
[i
]);
2299 value_assign(vals
[i
], list_args
[i
]);
2302 for (i
= 0; i
< p
->size
/2; ++i
)
2303 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2304 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2308 for (i
= 0; i
< dim
; ++i
)
2309 value_clear(vals
[i
]);
2318 } /* compute_enode */
2320 /*************************************************/
2321 /* return the value of Ehrhart Polynomial */
2322 /* It returns a double, because since it is */
2323 /* a recursive function, some intermediate value */
2324 /* might not be integral */
2325 /*************************************************/
2327 double compute_evalue(const evalue
*e
, Value
*list_args
)
2331 if (value_notzero_p(e
->d
)) {
2332 if (value_notone_p(e
->d
))
2333 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2335 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2338 res
= compute_enode(e
->x
.p
,list_args
);
2340 } /* compute_evalue */
2343 /****************************************************/
2344 /* function compute_poly : */
2345 /* Check for the good validity domain */
2346 /* return the number of point in the Polyhedron */
2347 /* in allocated memory */
2348 /* Using the Ehrhart pseudo-polynomial */
2349 /****************************************************/
2350 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2353 /* double d; int i; */
2355 tmp
= (Value
*) malloc (sizeof(Value
));
2356 assert(tmp
!= NULL
);
2358 value_set_si(*tmp
,0);
2361 return(tmp
); /* no ehrhart polynomial */
2362 if(en
->ValidityDomain
) {
2363 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2364 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2369 return(tmp
); /* no Validity Domain */
2371 if(in_domain(en
->ValidityDomain
,list_args
)) {
2373 #ifdef EVAL_EHRHART_DEBUG
2374 Print_Domain(stdout
,en
->ValidityDomain
);
2375 print_evalue(stdout
,&en
->EP
);
2378 /* d = compute_evalue(&en->EP,list_args);
2380 printf("(double)%lf = %d\n", d, i ); */
2381 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2387 value_set_si(*tmp
,0);
2388 return(tmp
); /* no compatible domain with the arguments */
2389 } /* compute_poly */
2391 static evalue
*eval_polynomial(const enode
*p
, int offset
,
2392 evalue
*base
, Value
*values
)
2397 res
= evalue_zero();
2398 for (i
= p
->size
-1; i
> offset
; --i
) {
2399 c
= evalue_eval(&p
->arr
[i
], values
);
2404 c
= evalue_eval(&p
->arr
[offset
], values
);
2411 evalue
*evalue_eval(const evalue
*e
, Value
*values
)
2418 if (value_notzero_p(e
->d
)) {
2419 res
= ALLOC(evalue
);
2421 evalue_copy(res
, e
);
2424 switch (e
->x
.p
->type
) {
2426 value_init(param
.x
.n
);
2427 value_assign(param
.x
.n
, values
[e
->x
.p
->pos
-1]);
2428 value_init(param
.d
);
2429 value_set_si(param
.d
, 1);
2431 res
= eval_polynomial(e
->x
.p
, 0, ¶m
, values
);
2432 free_evalue_refs(¶m
);
2435 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2436 mpz_fdiv_r(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2438 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2439 evalue_free(param2
);
2442 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2443 mpz_fdiv_q(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2444 value_set_si(param2
->d
, 1);
2446 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2447 evalue_free(param2
);
2450 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2451 if (value_zero_p(param2
->x
.n
))
2452 res
= evalue_eval(&e
->x
.p
->arr
[1], values
);
2453 else if (e
->x
.p
->size
> 2)
2454 res
= evalue_eval(&e
->x
.p
->arr
[2], values
);
2456 res
= evalue_zero();
2457 evalue_free(param2
);
2460 assert(e
->x
.p
->pos
== EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
);
2461 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2462 if (in_domain(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), values
)) {
2463 res
= evalue_eval(&e
->x
.p
->arr
[2*i
+1], values
);
2467 res
= evalue_zero();
2475 size_t value_size(Value v
) {
2476 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2477 * sizeof(v
[0]._mp_d
[0]);
2480 size_t domain_size(Polyhedron
*D
)
2483 size_t s
= sizeof(*D
);
2485 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2486 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2487 s
+= value_size(D
->Constraint
[i
][j
]);
2490 for (i = 0; i < D->NbRays; ++i)
2491 for (j = 0; j < D->Dimension+2; ++j)
2492 s += value_size(D->Ray[i][j]);
2495 return D
->next
? s
+domain_size(D
->next
) : s
;
2498 size_t enode_size(enode
*p
) {
2499 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2502 if (p
->type
== partition
)
2503 for (i
= 0; i
< p
->size
/2; ++i
) {
2504 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2505 s
+= evalue_size(&p
->arr
[2*i
+1]);
2508 for (i
= 0; i
< p
->size
; ++i
) {
2509 s
+= evalue_size(&p
->arr
[i
]);
2514 size_t evalue_size(evalue
*e
)
2516 size_t s
= sizeof(*e
);
2517 s
+= value_size(e
->d
);
2518 if (value_notzero_p(e
->d
))
2519 s
+= value_size(e
->x
.n
);
2521 s
+= enode_size(e
->x
.p
);
2525 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2527 evalue
*found
= NULL
;
2532 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2535 value_init(offset
.d
);
2536 value_init(offset
.x
.n
);
2537 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2538 value_lcm(offset
.d
, m
, offset
.d
);
2539 value_set_si(offset
.x
.n
, 1);
2542 evalue_copy(©
, cst
);
2545 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2547 if (eequal(base
, &e
->x
.p
->arr
[0]))
2548 found
= &e
->x
.p
->arr
[0];
2550 value_set_si(offset
.x
.n
, -2);
2553 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2555 if (eequal(base
, &e
->x
.p
->arr
[0]))
2558 free_evalue_refs(cst
);
2559 free_evalue_refs(&offset
);
2562 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2563 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2568 static evalue
*find_relation_pair(evalue
*e
)
2571 evalue
*found
= NULL
;
2573 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2576 if (e
->x
.p
->type
== fractional
) {
2581 poly_denom(&e
->x
.p
->arr
[0], &m
);
2583 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2584 cst
= &cst
->x
.p
->arr
[0])
2587 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2588 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2593 i
= e
->x
.p
->type
== relation
;
2594 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2595 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2600 void evalue_mod2relation(evalue
*e
) {
2603 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2606 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2607 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2608 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2609 value_clear(e
->x
.p
->arr
[2*i
].d
);
2610 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2612 if (2*i
< e
->x
.p
->size
) {
2613 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2614 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2619 if (e
->x
.p
->size
== 0) {
2621 evalue_set_si(e
, 0, 1);
2627 while ((d
= find_relation_pair(e
)) != NULL
) {
2631 value_init(split
.d
);
2632 value_set_si(split
.d
, 0);
2633 split
.x
.p
= new_enode(relation
, 3, 0);
2634 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2635 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2637 ev
= &split
.x
.p
->arr
[0];
2638 value_set_si(ev
->d
, 0);
2639 ev
->x
.p
= new_enode(fractional
, 3, -1);
2640 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2641 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2642 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2648 free_evalue_refs(&split
);
2652 static int evalue_comp(const void * a
, const void * b
)
2654 const evalue
*e1
= *(const evalue
**)a
;
2655 const evalue
*e2
= *(const evalue
**)b
;
2656 return ecmp(e1
, e2
);
2659 void evalue_combine(evalue
*e
)
2666 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2669 NALLOC(evs
, e
->x
.p
->size
/2);
2670 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2671 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2672 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2673 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2674 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2675 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2676 value_clear(p
->arr
[2*k
].d
);
2677 value_clear(p
->arr
[2*k
+1].d
);
2678 p
->arr
[2*k
] = *(evs
[i
]-1);
2679 p
->arr
[2*k
+1] = *(evs
[i
]);
2682 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2685 value_clear((evs
[i
]-1)->d
);
2689 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2690 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2691 free_evalue_refs(evs
[i
]);
2695 for (i
= 2*k
; i
< p
->size
; ++i
)
2696 value_clear(p
->arr
[i
].d
);
2703 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2705 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2707 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2710 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2711 Polyhedron
*D
, *N
, **P
;
2714 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2721 if (D
->NbEq
<= H
->NbEq
) {
2727 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2728 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2729 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2730 reduce_evalue(&tmp
);
2731 if (value_notzero_p(tmp
.d
) ||
2732 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2735 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2736 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2739 free_evalue_refs(&tmp
);
2745 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2747 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2749 value_clear(e
->x
.p
->arr
[2*i
].d
);
2750 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2752 if (2*i
< e
->x
.p
->size
) {
2753 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2754 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2761 H
= DomainConvex(D
, 0);
2762 E
= DomainDifference(H
, D
, 0);
2764 D
= DomainDifference(H
, E
, 0);
2767 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2771 /* Use smallest representative for coefficients in affine form in
2772 * argument of fractional.
2773 * Since any change will make the argument non-standard,
2774 * the containing evalue will have to be reduced again afterward.
2776 static void fractional_minimal_coefficients(enode
*p
)
2782 assert(p
->type
== fractional
);
2784 while (value_zero_p(pp
->d
)) {
2785 assert(pp
->x
.p
->type
== polynomial
);
2786 assert(pp
->x
.p
->size
== 2);
2787 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2788 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2789 if (value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2790 value_subtract(pp
->x
.p
->arr
[1].x
.n
,
2791 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2792 pp
= &pp
->x
.p
->arr
[0];
2798 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2803 unsigned dim
= D
->Dimension
;
2804 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2807 assert(p
->type
== fractional
|| p
->type
== flooring
);
2808 value_set_si(T
->p
[1][dim
], 1);
2809 evalue_extract_affine(&p
->arr
[0], T
->p
[0], &T
->p
[0][dim
], d
);
2810 I
= DomainImage(D
, T
, 0);
2811 H
= DomainConvex(I
, 0);
2821 static void replace_by_affine(evalue
*e
, Value offset
)
2828 value_init(inc
.x
.n
);
2829 value_set_si(inc
.d
, 1);
2830 value_oppose(inc
.x
.n
, offset
);
2831 eadd(&inc
, &p
->arr
[0]);
2832 reorder_terms_about(p
, &p
->arr
[0]); /* frees arr[0] */
2836 free_evalue_refs(&inc
);
2839 int evalue_range_reduction_in_domain(evalue
*e
, Polyhedron
*D
)
2848 if (value_notzero_p(e
->d
))
2853 if (p
->type
== relation
) {
2860 fractional_minimal_coefficients(p
->arr
[0].x
.p
);
2861 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
);
2862 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2863 equal
= value_eq(min
, max
);
2864 mpz_cdiv_q(min
, min
, d
);
2865 mpz_fdiv_q(max
, max
, d
);
2867 if (bounded
&& value_gt(min
, max
)) {
2873 evalue_set_si(e
, 0, 1);
2876 free_evalue_refs(&(p
->arr
[1]));
2877 free_evalue_refs(&(p
->arr
[0]));
2883 return r
? r
: evalue_range_reduction_in_domain(e
, D
);
2884 } else if (bounded
&& equal
) {
2887 free_evalue_refs(&(p
->arr
[2]));
2890 free_evalue_refs(&(p
->arr
[0]));
2896 return evalue_range_reduction_in_domain(e
, D
);
2897 } else if (bounded
&& value_eq(min
, max
)) {
2898 /* zero for a single value */
2900 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2901 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2902 value_multiply(min
, min
, d
);
2903 value_subtract(M
->p
[0][D
->Dimension
+1],
2904 M
->p
[0][D
->Dimension
+1], min
);
2905 E
= DomainAddConstraints(D
, M
, 0);
2911 r
= evalue_range_reduction_in_domain(&p
->arr
[1], E
);
2913 r
|= evalue_range_reduction_in_domain(&p
->arr
[2], D
);
2915 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2923 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2926 i
= p
->type
== relation
? 1 :
2927 p
->type
== fractional
? 1 : 0;
2928 for (; i
<p
->size
; i
++)
2929 r
|= evalue_range_reduction_in_domain(&p
->arr
[i
], D
);
2931 if (p
->type
!= fractional
) {
2932 if (r
&& p
->type
== polynomial
) {
2935 value_set_si(f
.d
, 0);
2936 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2937 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2938 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2939 reorder_terms_about(p
, &f
);
2950 fractional_minimal_coefficients(p
);
2951 I
= polynomial_projection(p
, D
, &d
, NULL
);
2952 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2953 mpz_fdiv_q(min
, min
, d
);
2954 mpz_fdiv_q(max
, max
, d
);
2955 value_subtract(d
, max
, min
);
2957 if (bounded
&& value_eq(min
, max
)) {
2958 replace_by_affine(e
, min
);
2960 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
2961 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
2962 * See pages 199-200 of PhD thesis.
2970 value_set_si(rem
.d
, 0);
2971 rem
.x
.p
= new_enode(fractional
, 3, -1);
2972 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
2973 value_clear(rem
.x
.p
->arr
[1].d
);
2974 value_clear(rem
.x
.p
->arr
[2].d
);
2975 rem
.x
.p
->arr
[1] = p
->arr
[1];
2976 rem
.x
.p
->arr
[2] = p
->arr
[2];
2977 for (i
= 3; i
< p
->size
; ++i
)
2978 p
->arr
[i
-2] = p
->arr
[i
];
2982 value_init(inc
.x
.n
);
2983 value_set_si(inc
.d
, 1);
2984 value_oppose(inc
.x
.n
, min
);
2987 evalue_copy(&t
, &p
->arr
[0]);
2991 value_set_si(f
.d
, 0);
2992 f
.x
.p
= new_enode(fractional
, 3, -1);
2993 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
2994 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2995 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
2997 value_init(factor
.d
);
2998 evalue_set_si(&factor
, -1, 1);
3004 value_clear(f
.x
.p
->arr
[1].x
.n
);
3005 value_clear(f
.x
.p
->arr
[2].x
.n
);
3006 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3007 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3011 reorder_terms(&rem
);
3018 free_evalue_refs(&inc
);
3019 free_evalue_refs(&t
);
3020 free_evalue_refs(&f
);
3021 free_evalue_refs(&factor
);
3022 free_evalue_refs(&rem
);
3024 evalue_range_reduction_in_domain(e
, D
);
3028 _reduce_evalue(&p
->arr
[0], 0, 1);
3040 void evalue_range_reduction(evalue
*e
)
3043 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3046 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3047 if (evalue_range_reduction_in_domain(&e
->x
.p
->arr
[2*i
+1],
3048 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
3049 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3051 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
3052 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
3053 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3054 value_clear(e
->x
.p
->arr
[2*i
].d
);
3056 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
3057 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
3065 Enumeration
* partition2enumeration(evalue
*EP
)
3068 Enumeration
*en
, *res
= NULL
;
3070 if (EVALUE_IS_ZERO(*EP
)) {
3075 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
3076 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
3077 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
3080 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
3081 value_clear(EP
->x
.p
->arr
[2*i
].d
);
3082 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
3090 int evalue_frac2floor_in_domain3(evalue
*e
, Polyhedron
*D
, int shift
)
3099 if (value_notzero_p(e
->d
))
3104 i
= p
->type
== relation
? 1 :
3105 p
->type
== fractional
? 1 : 0;
3106 for (; i
<p
->size
; i
++)
3107 r
|= evalue_frac2floor_in_domain3(&p
->arr
[i
], D
, shift
);
3109 if (p
->type
!= fractional
) {
3110 if (r
&& p
->type
== polynomial
) {
3113 value_set_si(f
.d
, 0);
3114 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3115 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3116 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3117 reorder_terms_about(p
, &f
);
3127 I
= polynomial_projection(p
, D
, &d
, NULL
);
3130 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3133 assert(I
->NbEq
== 0); /* Should have been reduced */
3136 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3137 if (value_pos_p(I
->Constraint
[i
][1]))
3140 if (i
< I
->NbConstraints
) {
3142 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3143 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3144 if (value_neg_p(min
)) {
3146 mpz_fdiv_q(min
, min
, d
);
3147 value_init(offset
.d
);
3148 value_set_si(offset
.d
, 1);
3149 value_init(offset
.x
.n
);
3150 value_oppose(offset
.x
.n
, min
);
3151 eadd(&offset
, &p
->arr
[0]);
3152 free_evalue_refs(&offset
);
3162 value_set_si(fl
.d
, 0);
3163 fl
.x
.p
= new_enode(flooring
, 3, -1);
3164 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
3165 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
3166 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
3168 eadd(&fl
, &p
->arr
[0]);
3169 reorder_terms_about(p
, &p
->arr
[0]);
3173 free_evalue_refs(&fl
);
3178 int evalue_frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
3180 return evalue_frac2floor_in_domain3(e
, D
, 1);
3183 void evalue_frac2floor2(evalue
*e
, int shift
)
3186 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
3188 if (evalue_frac2floor_in_domain3(e
, NULL
, 0))
3194 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3195 if (evalue_frac2floor_in_domain3(&e
->x
.p
->arr
[2*i
+1],
3196 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), shift
))
3197 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3200 void evalue_frac2floor(evalue
*e
)
3202 evalue_frac2floor2(e
, 1);
3205 /* Add a new variable with lower bound 1 and upper bound specified
3206 * by row. If negative is true, then the new variable has upper
3207 * bound -1 and lower bound specified by row.
3209 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
3210 Vector
*row
, int negative
)
3214 int nparam
= D
->Dimension
- nvar
;
3217 nr
= D
->NbConstraints
+ 2;
3218 nc
= D
->Dimension
+ 2 + 1;
3219 C
= Matrix_Alloc(nr
, nc
);
3220 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
3221 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
3222 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3223 D
->Dimension
+ 1 - nvar
);
3228 nc
= C
->NbColumns
+ 1;
3229 C
= Matrix_Alloc(nr
, nc
);
3230 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
3231 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
3232 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3233 oldC
->NbColumns
- 1 - nvar
);
3236 value_set_si(C
->p
[nr
-2][0], 1);
3238 value_set_si(C
->p
[nr
-2][1 + nvar
], -1);
3240 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
3241 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
3243 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
3244 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
3250 static void floor2frac_r(evalue
*e
, int nvar
)
3257 if (value_notzero_p(e
->d
))
3262 assert(p
->type
== flooring
);
3263 for (i
= 1; i
< p
->size
; i
++)
3264 floor2frac_r(&p
->arr
[i
], nvar
);
3266 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
3267 assert(pp
->x
.p
->type
== polynomial
);
3268 pp
->x
.p
->pos
-= nvar
;
3272 value_set_si(f
.d
, 0);
3273 f
.x
.p
= new_enode(fractional
, 3, -1);
3274 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3275 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3276 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3278 eadd(&f
, &p
->arr
[0]);
3279 reorder_terms_about(p
, &p
->arr
[0]);
3283 free_evalue_refs(&f
);
3286 /* Convert flooring back to fractional and shift position
3287 * of the parameters by nvar
3289 static void floor2frac(evalue
*e
, int nvar
)
3291 floor2frac_r(e
, nvar
);
3295 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
3298 int nparam
= D
->Dimension
- nvar
;
3302 D
= Constraints2Polyhedron(C
, 0);
3306 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
3308 /* Double check that D was not unbounded. */
3309 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
3317 static evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
3318 int *signs
, Matrix
*C
, unsigned MaxRays
)
3324 evalue
*factor
= NULL
;
3328 if (EVALUE_IS_ZERO(*e
))
3332 Polyhedron
*DD
= Disjoint_Domain(D
, 0, MaxRays
);
3339 res
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3342 for (Q
= DD
; Q
; Q
= DD
) {
3348 t
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3361 if (value_notzero_p(e
->d
)) {
3364 t
= esum_over_domain_cst(nvar
, D
, C
);
3366 if (!EVALUE_IS_ONE(*e
))
3372 switch (e
->x
.p
->type
) {
3374 evalue
*pp
= &e
->x
.p
->arr
[0];
3376 if (pp
->x
.p
->pos
> nvar
) {
3377 /* remainder is independent of the summated vars */
3383 floor2frac(&f
, nvar
);
3385 t
= esum_over_domain_cst(nvar
, D
, C
);
3389 free_evalue_refs(&f
);
3394 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3395 poly_denom(pp
, &row
->p
[1 + nvar
]);
3396 value_set_si(row
->p
[0], 1);
3397 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3398 pp
= &pp
->x
.p
->arr
[0]) {
3400 assert(pp
->x
.p
->type
== polynomial
);
3402 if (pos
>= 1 + nvar
)
3404 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3405 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3406 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3408 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3409 value_division(row
->p
[1 + D
->Dimension
+ 1],
3410 row
->p
[1 + D
->Dimension
+ 1],
3412 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3413 row
->p
[1 + D
->Dimension
+ 1],
3415 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3419 int pos
= e
->x
.p
->pos
;
3422 factor
= ALLOC(evalue
);
3423 value_init(factor
->d
);
3424 value_set_si(factor
->d
, 0);
3425 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3426 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3427 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3431 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3432 negative
= signs
[pos
-1] < 0;
3433 value_set_si(row
->p
[0], 1);
3435 value_set_si(row
->p
[pos
], -1);
3436 value_set_si(row
->p
[1 + nvar
], 1);
3438 value_set_si(row
->p
[pos
], 1);
3439 value_set_si(row
->p
[1 + nvar
], -1);
3447 offset
= type_offset(e
->x
.p
);
3449 res
= esum_over_domain(&e
->x
.p
->arr
[offset
], nvar
, D
, signs
, C
, MaxRays
);
3453 evalue_copy(&cum
, factor
);
3457 for (i
= 1; offset
+i
< e
->x
.p
->size
; ++i
) {
3461 C
= esum_add_constraint(nvar
, D
, C
, row
, negative
);
3467 Vector_Print(stderr, P_VALUE_FMT, row);
3469 Matrix_Print(stderr, P_VALUE_FMT, C);
3471 t
= esum_over_domain(&e
->x
.p
->arr
[offset
+i
], nvar
, D
, signs
, C
, MaxRays
);
3476 if (negative
&& (i
% 2))
3486 if (factor
&& offset
+i
+1 < e
->x
.p
->size
)
3493 free_evalue_refs(&cum
);
3494 evalue_free(factor
);
3505 static void domain_signs(Polyhedron
*D
, int *signs
)
3509 POL_ENSURE_VERTICES(D
);
3510 for (j
= 0; j
< D
->Dimension
; ++j
) {
3512 for (k
= 0; k
< D
->NbRays
; ++k
) {
3513 signs
[j
] = value_sign(D
->Ray
[k
][1+j
]);
3520 static void shift_floor_in_domain(evalue
*e
, Polyhedron
*D
)
3527 if (value_notzero_p(e
->d
))
3532 for (i
= type_offset(p
); i
< p
->size
; ++i
)
3533 shift_floor_in_domain(&p
->arr
[i
], D
);
3535 if (p
->type
!= flooring
)
3541 I
= polynomial_projection(p
, D
, &d
, NULL
);
3542 assert(I
->NbEq
== 0); /* Should have been reduced */
3544 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3545 if (value_pos_p(I
->Constraint
[i
][1]))
3547 assert(i
< I
->NbConstraints
);
3548 if (i
< I
->NbConstraints
) {
3549 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3550 mpz_fdiv_q(m
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3551 if (value_neg_p(m
)) {
3552 /* replace [e] by [e-m]+m such that e-m >= 0 */
3557 value_set_si(f
.d
, 1);
3558 value_oppose(f
.x
.n
, m
);
3559 eadd(&f
, &p
->arr
[0]);
3562 value_set_si(f
.d
, 0);
3563 f
.x
.p
= new_enode(flooring
, 3, -1);
3564 value_clear(f
.x
.p
->arr
[0].d
);
3565 f
.x
.p
->arr
[0] = p
->arr
[0];
3566 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
3567 value_set_si(f
.x
.p
->arr
[1].d
, 1);
3568 value_init(f
.x
.p
->arr
[1].x
.n
);
3569 value_assign(f
.x
.p
->arr
[1].x
.n
, m
);
3570 reorder_terms_about(p
, &f
);
3581 /* Make arguments of all floors non-negative */
3582 static void shift_floor_arguments(evalue
*e
)
3586 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3589 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3590 shift_floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
3591 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3594 evalue
*evalue_sum(evalue
*e
, int nvar
, unsigned MaxRays
)
3598 evalue
*res
= ALLOC(evalue
);
3602 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3603 evalue_copy(res
, e
);
3607 evalue_split_domains_into_orthants(e
, MaxRays
);
3608 evalue_frac2floor2(e
, 0);
3609 evalue_set_si(res
, 0, 1);
3611 assert(value_zero_p(e
->d
));
3612 assert(e
->x
.p
->type
== partition
);
3613 shift_floor_arguments(e
);
3615 assert(e
->x
.p
->size
>= 2);
3616 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3618 signs
= alloca(sizeof(int) * dim
);
3620 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3622 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3623 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3624 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
, 0,
3635 evalue
*esum(evalue
*e
, int nvar
)
3637 return evalue_sum(e
, nvar
, 0);
3640 /* Initial silly implementation */
3641 void eor(evalue
*e1
, evalue
*res
)
3647 evalue_set_si(&mone
, -1, 1);
3649 evalue_copy(&E
, res
);
3655 free_evalue_refs(&E
);
3656 free_evalue_refs(&mone
);
3659 /* computes denominator of polynomial evalue
3660 * d should point to a value initialized to 1
3662 void evalue_denom(const evalue
*e
, Value
*d
)
3666 if (value_notzero_p(e
->d
)) {
3667 value_lcm(*d
, *d
, e
->d
);
3670 assert(e
->x
.p
->type
== polynomial
||
3671 e
->x
.p
->type
== fractional
||
3672 e
->x
.p
->type
== flooring
);
3673 offset
= type_offset(e
->x
.p
);
3674 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3675 evalue_denom(&e
->x
.p
->arr
[i
], d
);
3678 /* Divides the evalue e by the integer n */
3679 void evalue_div(evalue
*e
, Value n
)
3683 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3686 if (value_notzero_p(e
->d
)) {
3689 value_multiply(e
->d
, e
->d
, n
);
3690 value_gcd(gc
, e
->x
.n
, e
->d
);
3691 if (value_notone_p(gc
)) {
3692 value_division(e
->d
, e
->d
, gc
);
3693 value_division(e
->x
.n
, e
->x
.n
, gc
);
3698 if (e
->x
.p
->type
== partition
) {
3699 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3700 evalue_div(&e
->x
.p
->arr
[2*i
+1], n
);
3703 offset
= type_offset(e
->x
.p
);
3704 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3705 evalue_div(&e
->x
.p
->arr
[i
], n
);
3708 /* Multiplies the evalue e by the integer n */
3709 void evalue_mul(evalue
*e
, Value n
)
3713 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3716 if (value_notzero_p(e
->d
)) {
3719 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3720 value_gcd(gc
, e
->x
.n
, e
->d
);
3721 if (value_notone_p(gc
)) {
3722 value_division(e
->d
, e
->d
, gc
);
3723 value_division(e
->x
.n
, e
->x
.n
, gc
);
3728 if (e
->x
.p
->type
== partition
) {
3729 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3730 evalue_mul(&e
->x
.p
->arr
[2*i
+1], n
);
3733 offset
= type_offset(e
->x
.p
);
3734 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3735 evalue_mul(&e
->x
.p
->arr
[i
], n
);
3738 /* Multiplies the evalue e by the n/d */
3739 void evalue_mul_div(evalue
*e
, Value n
, Value d
)
3743 if ((value_one_p(n
) && value_one_p(d
)) || EVALUE_IS_ZERO(*e
))
3746 if (value_notzero_p(e
->d
)) {
3749 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3750 value_multiply(e
->d
, e
->d
, d
);
3751 value_gcd(gc
, e
->x
.n
, e
->d
);
3752 if (value_notone_p(gc
)) {
3753 value_division(e
->d
, e
->d
, gc
);
3754 value_division(e
->x
.n
, e
->x
.n
, gc
);
3759 if (e
->x
.p
->type
== partition
) {
3760 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3761 evalue_mul_div(&e
->x
.p
->arr
[2*i
+1], n
, d
);
3764 offset
= type_offset(e
->x
.p
);
3765 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3766 evalue_mul_div(&e
->x
.p
->arr
[i
], n
, d
);
3769 void evalue_negate(evalue
*e
)
3773 if (value_notzero_p(e
->d
)) {
3774 value_oppose(e
->x
.n
, e
->x
.n
);
3777 if (e
->x
.p
->type
== partition
) {
3778 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3779 evalue_negate(&e
->x
.p
->arr
[2*i
+1]);
3782 offset
= type_offset(e
->x
.p
);
3783 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3784 evalue_negate(&e
->x
.p
->arr
[i
]);
3787 void evalue_add_constant(evalue
*e
, const Value cst
)
3791 if (value_zero_p(e
->d
)) {
3792 if (e
->x
.p
->type
== partition
) {
3793 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3794 evalue_add_constant(&e
->x
.p
->arr
[2*i
+1], cst
);
3797 if (e
->x
.p
->type
== relation
) {
3798 for (i
= 1; i
< e
->x
.p
->size
; ++i
)
3799 evalue_add_constant(&e
->x
.p
->arr
[i
], cst
);
3803 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
3804 } while (value_zero_p(e
->d
));
3806 value_addmul(e
->x
.n
, cst
, e
->d
);
3809 static void evalue_frac2polynomial_r(evalue
*e
, int *signs
, int sign
, int in_frac
)
3814 int sign_odd
= sign
;
3816 if (value_notzero_p(e
->d
)) {
3817 if (in_frac
&& sign
* value_sign(e
->x
.n
) < 0) {
3818 value_set_si(e
->x
.n
, 0);
3819 value_set_si(e
->d
, 1);
3824 if (e
->x
.p
->type
== relation
) {
3825 for (i
= e
->x
.p
->size
-1; i
>= 1; --i
)
3826 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
, sign
, in_frac
);
3830 if (e
->x
.p
->type
== polynomial
)
3831 sign_odd
*= signs
[e
->x
.p
->pos
-1];
3832 offset
= type_offset(e
->x
.p
);
3833 evalue_frac2polynomial_r(&e
->x
.p
->arr
[offset
], signs
, sign
, in_frac
);
3834 in_frac
|= e
->x
.p
->type
== fractional
;
3835 for (i
= e
->x
.p
->size
-1; i
> offset
; --i
)
3836 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
,
3837 (i
- offset
) % 2 ? sign_odd
: sign
, in_frac
);
3839 if (e
->x
.p
->type
!= fractional
)
3842 /* replace { a/m } by (m-1)/m if sign != 0
3843 * and by (m-1)/(2m) if sign == 0
3847 evalue_denom(&e
->x
.p
->arr
[0], &d
);
3848 free_evalue_refs(&e
->x
.p
->arr
[0]);
3849 value_init(e
->x
.p
->arr
[0].d
);
3850 value_init(e
->x
.p
->arr
[0].x
.n
);
3852 value_addto(e
->x
.p
->arr
[0].d
, d
, d
);
3854 value_assign(e
->x
.p
->arr
[0].d
, d
);
3855 value_decrement(e
->x
.p
->arr
[0].x
.n
, d
);
3859 reorder_terms_about(p
, &p
->arr
[0]);
3865 /* Approximate the evalue in fractional representation by a polynomial.
3866 * If sign > 0, the result is an upper bound;
3867 * if sign < 0, the result is a lower bound;
3868 * if sign = 0, the result is an intermediate approximation.
3870 void evalue_frac2polynomial(evalue
*e
, int sign
, unsigned MaxRays
)
3875 if (value_notzero_p(e
->d
))
3877 assert(e
->x
.p
->type
== partition
);
3878 /* make sure all variables in the domains have a fixed sign */
3880 evalue_split_domains_into_orthants(e
, MaxRays
);
3881 if (EVALUE_IS_ZERO(*e
))
3885 assert(e
->x
.p
->size
>= 2);
3886 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3888 signs
= alloca(sizeof(int) * dim
);
3891 for (i
= 0; i
< dim
; ++i
)
3893 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3895 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3896 evalue_frac2polynomial_r(&e
->x
.p
->arr
[2*i
+1], signs
, sign
, 0);
3900 /* Split the domains of e (which is assumed to be a partition)
3901 * such that each resulting domain lies entirely in one orthant.
3903 void evalue_split_domains_into_orthants(evalue
*e
, unsigned MaxRays
)
3906 assert(value_zero_p(e
->d
));
3907 assert(e
->x
.p
->type
== partition
);
3908 assert(e
->x
.p
->size
>= 2);
3909 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3911 for (i
= 0; i
< dim
; ++i
) {
3914 C
= Matrix_Alloc(1, 1 + dim
+ 1);
3915 value_set_si(C
->p
[0][0], 1);
3916 value_init(split
.d
);
3917 value_set_si(split
.d
, 0);
3918 split
.x
.p
= new_enode(partition
, 4, dim
);
3919 value_set_si(C
->p
[0][1+i
], 1);
3920 C2
= Matrix_Copy(C
);
3921 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[0], Constraints2Polyhedron(C2
, MaxRays
));
3923 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
3924 value_set_si(C
->p
[0][1+i
], -1);
3925 value_set_si(C
->p
[0][1+dim
], -1);
3926 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[2], Constraints2Polyhedron(C
, MaxRays
));
3927 evalue_set_si(&split
.x
.p
->arr
[3], 1, 1);
3929 free_evalue_refs(&split
);
3935 static evalue
*find_fractional_with_max_periods(evalue
*e
, Polyhedron
*D
,
3938 Value
*min
, Value
*max
)
3945 if (value_notzero_p(e
->d
))
3948 if (e
->x
.p
->type
== fractional
) {
3953 I
= polynomial_projection(e
->x
.p
, D
, &d
, &T
);
3954 bounded
= line_minmax(I
, min
, max
); /* frees I */
3958 value_set_si(mp
, max_periods
);
3959 mpz_fdiv_q(*min
, *min
, d
);
3960 mpz_fdiv_q(*max
, *max
, d
);
3961 value_assign(T
->p
[1][D
->Dimension
], d
);
3962 value_subtract(d
, *max
, *min
);
3963 if (value_ge(d
, mp
))
3966 f
= evalue_dup(&e
->x
.p
->arr
[0]);
3977 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
3978 if ((f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[i
], D
, max_periods
,
3985 static void replace_fract_by_affine(evalue
*e
, evalue
*f
, Value val
)
3989 if (value_notzero_p(e
->d
))
3992 offset
= type_offset(e
->x
.p
);
3993 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3994 replace_fract_by_affine(&e
->x
.p
->arr
[i
], f
, val
);
3996 if (e
->x
.p
->type
!= fractional
)
3999 if (!eequal(&e
->x
.p
->arr
[0], f
))
4002 replace_by_affine(e
, val
);
4005 /* Look for fractional parts that can be removed by splitting the corresponding
4006 * domain into at most max_periods parts.
4007 * We use a very simply strategy that looks for the first fractional part
4008 * that satisfies the condition, performs the split and then continues
4009 * looking for other fractional parts in the split domains until no
4010 * such fractional part can be found anymore.
4012 void evalue_split_periods(evalue
*e
, int max_periods
, unsigned int MaxRays
)
4019 if (EVALUE_IS_ZERO(*e
))
4021 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
4023 "WARNING: evalue_split_periods called on incorrect evalue type\n");
4031 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
4036 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
4038 f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[2*i
+1], D
, max_periods
,
4043 M
= Matrix_Alloc(2, 2+D
->Dimension
);
4045 value_subtract(d
, max
, min
);
4046 n
= VALUE_TO_INT(d
)+1;
4048 value_set_si(M
->p
[0][0], 1);
4049 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
4050 value_multiply(d
, max
, T
->p
[1][D
->Dimension
]);
4051 value_subtract(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
], d
);
4052 value_set_si(d
, -1);
4053 value_set_si(M
->p
[1][0], 1);
4054 Vector_Scale(T
->p
[0], M
->p
[1]+1, d
, D
->Dimension
+1);
4055 value_addmul(M
->p
[1][1+D
->Dimension
], max
, T
->p
[1][D
->Dimension
]);
4056 value_addto(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4057 T
->p
[1][D
->Dimension
]);
4058 value_decrement(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
]);
4060 p
= new_enode(partition
, e
->x
.p
->size
+ (n
-1)*2, e
->x
.p
->pos
);
4061 for (j
= 0; j
< 2*i
; ++j
) {
4062 value_clear(p
->arr
[j
].d
);
4063 p
->arr
[j
] = e
->x
.p
->arr
[j
];
4065 for (j
= 2*i
+2; j
< e
->x
.p
->size
; ++j
) {
4066 value_clear(p
->arr
[j
+2*(n
-1)].d
);
4067 p
->arr
[j
+2*(n
-1)] = e
->x
.p
->arr
[j
];
4069 for (j
= n
-1; j
>= 0; --j
) {
4071 value_clear(p
->arr
[2*i
+1].d
);
4072 p
->arr
[2*i
+1] = e
->x
.p
->arr
[2*i
+1];
4074 evalue_copy(&p
->arr
[2*(i
+j
)+1], &e
->x
.p
->arr
[2*i
+1]);
4076 value_subtract(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4077 T
->p
[1][D
->Dimension
]);
4078 value_addto(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
],
4079 T
->p
[1][D
->Dimension
]);
4081 replace_fract_by_affine(&p
->arr
[2*(i
+j
)+1], f
, max
);
4082 E
= DomainAddConstraints(D
, M
, MaxRays
);
4083 EVALUE_SET_DOMAIN(p
->arr
[2*(i
+j
)], E
);
4084 if (evalue_range_reduction_in_domain(&p
->arr
[2*(i
+j
)+1], E
))
4085 reduce_evalue(&p
->arr
[2*(i
+j
)+1]);
4086 value_decrement(max
, max
);
4088 value_clear(e
->x
.p
->arr
[2*i
].d
);
4103 void evalue_extract_affine(const evalue
*e
, Value
*coeff
, Value
*cst
, Value
*d
)
4105 value_set_si(*d
, 1);
4107 for ( ; value_zero_p(e
->d
); e
= &e
->x
.p
->arr
[0]) {
4109 assert(e
->x
.p
->type
== polynomial
);
4110 assert(e
->x
.p
->size
== 2);
4111 c
= &e
->x
.p
->arr
[1];
4112 value_multiply(coeff
[e
->x
.p
->pos
-1], *d
, c
->x
.n
);
4113 value_division(coeff
[e
->x
.p
->pos
-1], coeff
[e
->x
.p
->pos
-1], c
->d
);
4115 value_multiply(*cst
, *d
, e
->x
.n
);
4116 value_division(*cst
, *cst
, e
->d
);
4119 /* returns an evalue that corresponds to
4123 static evalue
*term(int param
, Value c
, Value den
)
4125 evalue
*EP
= ALLOC(evalue
);
4127 value_set_si(EP
->d
,0);
4128 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
4129 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
4130 value_init(EP
->x
.p
->arr
[1].x
.n
);
4131 value_assign(EP
->x
.p
->arr
[1].d
, den
);
4132 value_assign(EP
->x
.p
->arr
[1].x
.n
, c
);
4136 evalue
*affine2evalue(Value
*coeff
, Value denom
, int nvar
)
4139 evalue
*E
= ALLOC(evalue
);
4141 evalue_set(E
, coeff
[nvar
], denom
);
4142 for (i
= 0; i
< nvar
; ++i
) {
4144 if (value_zero_p(coeff
[i
]))
4146 t
= term(i
, coeff
[i
], denom
);
4153 void evalue_substitute(evalue
*e
, evalue
**subs
)
4159 if (value_notzero_p(e
->d
))
4163 assert(p
->type
!= partition
);
4165 for (i
= 0; i
< p
->size
; ++i
)
4166 evalue_substitute(&p
->arr
[i
], subs
);
4168 if (p
->type
== relation
) {
4169 /* For relation a ? b : c, compute (a' ? 1) * b' + (a' ? 0 : 1) * c' */
4173 value_set_si(v
->d
, 0);
4174 v
->x
.p
= new_enode(relation
, 3, 0);
4175 evalue_copy(&v
->x
.p
->arr
[0], &p
->arr
[0]);
4176 evalue_set_si(&v
->x
.p
->arr
[1], 0, 1);
4177 evalue_set_si(&v
->x
.p
->arr
[2], 1, 1);
4178 emul(v
, &p
->arr
[2]);
4183 value_set_si(v
->d
, 0);
4184 v
->x
.p
= new_enode(relation
, 2, 0);
4185 value_clear(v
->x
.p
->arr
[0].d
);
4186 v
->x
.p
->arr
[0] = p
->arr
[0];
4187 evalue_set_si(&v
->x
.p
->arr
[1], 1, 1);
4188 emul(v
, &p
->arr
[1]);
4191 eadd(&p
->arr
[2], &p
->arr
[1]);
4192 free_evalue_refs(&p
->arr
[2]);
4200 if (p
->type
== polynomial
)
4205 value_set_si(v
->d
, 0);
4206 v
->x
.p
= new_enode(p
->type
, 3, -1);
4207 value_clear(v
->x
.p
->arr
[0].d
);
4208 v
->x
.p
->arr
[0] = p
->arr
[0];
4209 evalue_set_si(&v
->x
.p
->arr
[1], 0, 1);
4210 evalue_set_si(&v
->x
.p
->arr
[2], 1, 1);
4213 offset
= type_offset(p
);
4215 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
4216 emul(v
, &p
->arr
[i
]);
4217 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
4218 free_evalue_refs(&(p
->arr
[i
]));
4221 if (p
->type
!= polynomial
)
4225 *e
= p
->arr
[offset
];
4229 /* evalue e is given in terms of "new" parameter; CP maps the new
4230 * parameters back to the old parameters.
4231 * Transforms e such that it refers back to the old parameters and
4232 * adds appropriate constraints to the domain.
4233 * In particular, if CP maps the new parameters onto an affine
4234 * subspace of the old parameters, then the corresponding equalities
4235 * are added to the domain.
4236 * Also, if any of the new parameters was a rational combination
4237 * of the old parameters $p' = (<a, p> + c)/m$, then modulo
4238 * constraints ${<a, p> + c)/m} = 0$ are added to ensure
4239 * the new evalue remains non-zero only for integer parameters
4240 * of the new parameters (which have been removed by the substitution).
4242 void evalue_backsubstitute(evalue
*e
, Matrix
*CP
, unsigned MaxRays
)
4249 unsigned nparam
= CP
->NbColumns
-1;
4253 if (EVALUE_IS_ZERO(*e
))
4256 assert(value_zero_p(e
->d
));
4258 assert(p
->type
== partition
);
4260 inv
= left_inverse(CP
, &eq
);
4261 subs
= ALLOCN(evalue
*, nparam
);
4262 for (i
= 0; i
< nparam
; ++i
)
4263 subs
[i
] = affine2evalue(inv
->p
[i
], inv
->p
[nparam
][inv
->NbColumns
-1],
4266 CEq
= Constraints2Polyhedron(eq
, MaxRays
);
4267 addeliminatedparams_partition(p
, inv
, CEq
, inv
->NbColumns
-1, MaxRays
);
4268 Polyhedron_Free(CEq
);
4270 for (i
= 0; i
< p
->size
/2; ++i
)
4271 evalue_substitute(&p
->arr
[2*i
+1], subs
);
4273 for (i
= 0; i
< nparam
; ++i
)
4274 evalue_free(subs
[i
]);
4278 for (i
= 0; i
< inv
->NbRows
-1; ++i
) {
4279 Vector_Gcd(inv
->p
[i
], inv
->NbColumns
, &gcd
);
4280 value_gcd(gcd
, gcd
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]);
4281 if (value_eq(gcd
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]))
4283 Vector_AntiScale(inv
->p
[i
], inv
->p
[i
], gcd
, inv
->NbColumns
);
4284 value_divexact(gcd
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1], gcd
);
4286 for (j
= 0; j
< p
->size
/2; ++j
) {
4287 evalue
*arg
= affine2evalue(inv
->p
[i
], gcd
, inv
->NbColumns
-1);
4292 value_set_si(rel
.d
, 0);
4293 rel
.x
.p
= new_enode(relation
, 2, 0);
4294 value_clear(rel
.x
.p
->arr
[1].d
);
4295 rel
.x
.p
->arr
[1] = p
->arr
[2*j
+1];
4296 ev
= &rel
.x
.p
->arr
[0];
4297 value_set_si(ev
->d
, 0);
4298 ev
->x
.p
= new_enode(fractional
, 3, -1);
4299 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
4300 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
4301 value_clear(ev
->x
.p
->arr
[0].d
);
4302 ev
->x
.p
->arr
[0] = *arg
;
4305 p
->arr
[2*j
+1] = rel
;
4316 * \sum_{i=0}^n c_i/d X^i
4318 * where d is the last element in the vector c.
4320 evalue
*evalue_polynomial(Vector
*c
, const evalue
* X
)
4322 unsigned dim
= c
->Size
-2;
4324 evalue
*EP
= ALLOC(evalue
);
4329 if (EVALUE_IS_ZERO(*X
) || dim
== 0) {
4330 evalue_set(EP
, c
->p
[0], c
->p
[dim
+1]);
4331 reduce_constant(EP
);
4335 evalue_set(EP
, c
->p
[dim
], c
->p
[dim
+1]);
4338 evalue_set(&EC
, c
->p
[dim
], c
->p
[dim
+1]);
4340 for (i
= dim
-1; i
>= 0; --i
) {
4342 value_assign(EC
.x
.n
, c
->p
[i
]);
4345 free_evalue_refs(&EC
);
4349 /* Create an evalue from an array of pairs of domains and evalues. */
4350 evalue
*evalue_from_section_array(struct evalue_section
*s
, int n
)
4355 res
= ALLOC(evalue
);
4359 evalue_set_si(res
, 0, 1);
4361 value_set_si(res
->d
, 0);
4362 res
->x
.p
= new_enode(partition
, 2*n
, s
[0].D
->Dimension
);
4363 for (i
= 0; i
< n
; ++i
) {
4364 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
], s
[i
].D
);
4365 value_clear(res
->x
.p
->arr
[2*i
+1].d
);
4366 res
->x
.p
->arr
[2*i
+1] = *s
[i
].E
;