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);
55 evalue
*EP
= ALLOC(evalue
);
57 value_set_si(EP
->d
, -2);
62 /* returns an evalue that corresponds to
66 evalue
*evalue_var(int var
)
68 evalue
*EP
= ALLOC(evalue
);
70 value_set_si(EP
->d
,0);
71 EP
->x
.p
= new_enode(polynomial
, 2, var
+ 1);
72 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
73 evalue_set_si(&EP
->x
.p
->arr
[1], 1, 1);
77 void aep_evalue(evalue
*e
, int *ref
) {
82 if (value_notzero_p(e
->d
))
83 return; /* a rational number, its already reduced */
85 return; /* hum... an overflow probably occured */
87 /* First check the components of p */
88 for (i
=0;i
<p
->size
;i
++)
89 aep_evalue(&p
->arr
[i
],ref
);
96 p
->pos
= ref
[p
->pos
-1]+1;
102 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
108 if (value_notzero_p(e
->d
))
109 return; /* a rational number, its already reduced */
111 return; /* hum... an overflow probably occured */
114 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
115 for(i
=0;i
<CT
->NbRows
-1;i
++)
116 for(j
=0;j
<CT
->NbColumns
;j
++)
117 if(value_notzero_p(CT
->p
[i
][j
])) {
122 /* Transform the references in e, using ref */
126 } /* addeliminatedparams_evalue */
128 static void addeliminatedparams_partition(enode
*p
, Matrix
*CT
, Polyhedron
*CEq
,
129 unsigned nparam
, unsigned MaxRays
)
132 assert(p
->type
== partition
);
135 for (i
= 0; i
< p
->size
/2; i
++) {
136 Polyhedron
*D
= EVALUE_DOMAIN(p
->arr
[2*i
]);
137 Polyhedron
*T
= DomainPreimage(D
, CT
, MaxRays
);
141 T
= DomainIntersection(D
, CEq
, MaxRays
);
144 EVALUE_SET_DOMAIN(p
->arr
[2*i
], T
);
148 void addeliminatedparams_enum(evalue
*e
, Matrix
*CT
, Polyhedron
*CEq
,
149 unsigned MaxRays
, unsigned nparam
)
154 if (CT
->NbRows
== CT
->NbColumns
)
157 if (EVALUE_IS_ZERO(*e
))
160 if (value_notzero_p(e
->d
)) {
163 value_set_si(res
.d
, 0);
164 res
.x
.p
= new_enode(partition
, 2, nparam
);
165 EVALUE_SET_DOMAIN(res
.x
.p
->arr
[0],
166 DomainConstraintSimplify(Polyhedron_Copy(CEq
), MaxRays
));
167 value_clear(res
.x
.p
->arr
[1].d
);
168 res
.x
.p
->arr
[1] = *e
;
176 addeliminatedparams_partition(p
, CT
, CEq
, nparam
, MaxRays
);
177 for (i
= 0; i
< p
->size
/2; i
++)
178 addeliminatedparams_evalue(&p
->arr
[2*i
+1], CT
);
181 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
189 assert(value_notzero_p(e1
->d
));
190 assert(value_notzero_p(e2
->d
));
191 value_multiply(m
, e1
->x
.n
, e2
->d
);
192 value_multiply(m2
, e2
->x
.n
, e1
->d
);
195 else if (value_gt(m
, m2
))
205 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
207 if (value_notzero_p(e1
->d
)) {
209 if (value_zero_p(e2
->d
))
211 r
= mod_rational_smaller(e1
, e2
);
212 return r
== -1 ? 0 : r
;
214 if (value_notzero_p(e2
->d
))
216 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
218 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
221 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
223 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
228 static int mod_term_smaller(const evalue
*e1
, const evalue
*e2
)
230 assert(value_zero_p(e1
->d
));
231 assert(value_zero_p(e2
->d
));
232 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
233 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
234 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
237 static void check_order(const evalue
*e
)
242 if (value_notzero_p(e
->d
))
245 switch (e
->x
.p
->type
) {
247 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
248 check_order(&e
->x
.p
->arr
[2*i
+1]);
251 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
253 if (value_notzero_p(a
->d
))
255 switch (a
->x
.p
->type
) {
257 assert(mod_term_smaller(&e
->x
.p
->arr
[0], &a
->x
.p
->arr
[0]));
266 for (i
= 0; i
< e
->x
.p
->size
; ++i
) {
268 if (value_notzero_p(a
->d
))
270 switch (a
->x
.p
->type
) {
272 assert(e
->x
.p
->pos
< a
->x
.p
->pos
);
283 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
285 if (value_notzero_p(a
->d
))
287 switch (a
->x
.p
->type
) {
298 /* Negative pos means inequality */
299 /* s is negative of substitution if m is not zero */
308 struct fixed_param
*fixed
;
313 static int relations_depth(evalue
*e
)
318 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
319 e
= &e
->x
.p
->arr
[1], ++d
);
323 static void poly_denom_not_constant(evalue
**pp
, Value
*d
)
328 while (value_zero_p(p
->d
)) {
329 assert(p
->x
.p
->type
== polynomial
);
330 assert(p
->x
.p
->size
== 2);
331 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
332 value_lcm(*d
, *d
, p
->x
.p
->arr
[1].d
);
338 static void poly_denom(evalue
*p
, Value
*d
)
340 poly_denom_not_constant(&p
, d
);
341 value_lcm(*d
, *d
, p
->d
);
344 static void realloc_substitution(struct subst
*s
, int d
)
346 struct fixed_param
*n
;
349 for (i
= 0; i
< s
->n
; ++i
)
356 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
362 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
365 /* May have been reduced already */
366 if (value_notzero_p(m
->d
))
369 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
370 assert(m
->x
.p
->size
== 3);
372 /* fractional was inverted during reduction
373 * invert it back and move constant in
375 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
376 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
377 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
378 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
379 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
380 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
381 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
382 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
383 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
384 value_set_si(m
->x
.p
->arr
[1].d
, 1);
387 /* Oops. Nested identical relations. */
388 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
391 if (s
->n
>= s
->max
) {
392 int d
= relations_depth(r
);
393 realloc_substitution(s
, d
);
397 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
398 assert(p
->x
.p
->size
== 2);
401 assert(value_pos_p(f
->x
.n
));
403 value_init(s
->fixed
[s
->n
].m
);
404 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
405 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
406 value_init(s
->fixed
[s
->n
].d
);
407 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
408 value_init(s
->fixed
[s
->n
].s
.d
);
409 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
415 static int type_offset(enode
*p
)
417 return p
->type
== fractional
? 1 :
418 p
->type
== flooring
? 1 :
419 p
->type
== relation
? 1 : 0;
422 static void reorder_terms_about(enode
*p
, evalue
*v
)
425 int offset
= type_offset(p
);
427 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
429 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
430 free_evalue_refs(&(p
->arr
[i
]));
436 static void reorder_terms(evalue
*e
)
441 assert(value_zero_p(e
->d
));
443 assert(p
->type
== fractional
); /* for now */
446 value_set_si(f
.d
, 0);
447 f
.x
.p
= new_enode(fractional
, 3, -1);
448 value_clear(f
.x
.p
->arr
[0].d
);
449 f
.x
.p
->arr
[0] = p
->arr
[0];
450 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
451 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
452 reorder_terms_about(p
, &f
);
458 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
464 if (value_notzero_p(e
->d
)) {
466 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
467 return; /* a rational number, its already reduced */
471 return; /* hum... an overflow probably occured */
473 /* First reduce the components of p */
474 add
= p
->type
== relation
;
475 for (i
=0; i
<p
->size
; i
++) {
477 add
= add_modulo_substitution(s
, e
);
479 if (i
== 0 && p
->type
==fractional
)
480 _reduce_evalue(&p
->arr
[i
], s
, 1);
482 _reduce_evalue(&p
->arr
[i
], s
, fract
);
484 if (add
&& i
== p
->size
-1) {
486 value_clear(s
->fixed
[s
->n
].m
);
487 value_clear(s
->fixed
[s
->n
].d
);
488 free_evalue_refs(&s
->fixed
[s
->n
].s
);
489 } else if (add
&& i
== 1)
490 s
->fixed
[s
->n
-1].pos
*= -1;
493 if (p
->type
==periodic
) {
495 /* Try to reduce the period */
496 for (i
=1; i
<=(p
->size
)/2; i
++) {
497 if ((p
->size
% i
)==0) {
499 /* Can we reduce the size to i ? */
501 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
502 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
505 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
509 you_lose
: /* OK, lets not do it */
514 /* Try to reduce its strength */
517 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
521 else if (p
->type
==polynomial
) {
522 for (k
= 0; s
&& k
< s
->n
; ++k
) {
523 if (s
->fixed
[k
].pos
== p
->pos
) {
524 int divide
= value_notone_p(s
->fixed
[k
].d
);
527 if (value_notzero_p(s
->fixed
[k
].m
)) {
530 assert(p
->size
== 2);
531 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
533 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
540 value_assign(d
.d
, s
->fixed
[k
].d
);
542 if (value_notzero_p(s
->fixed
[k
].m
))
543 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
545 value_set_si(d
.x
.n
, 1);
548 for (i
=p
->size
-1;i
>=1;i
--) {
549 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
551 emul(&d
, &p
->arr
[i
]);
552 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
553 free_evalue_refs(&(p
->arr
[i
]));
556 _reduce_evalue(&p
->arr
[0], s
, fract
);
559 free_evalue_refs(&d
);
565 /* Try to reduce the degree */
566 for (i
=p
->size
-1;i
>=1;i
--) {
567 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
569 /* Zero coefficient */
570 free_evalue_refs(&(p
->arr
[i
]));
575 /* Try to reduce its strength */
578 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
582 else if (p
->type
==fractional
) {
586 if (value_notzero_p(p
->arr
[0].d
)) {
588 value_assign(v
.d
, p
->arr
[0].d
);
590 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
595 evalue
*pp
= &p
->arr
[0];
596 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
597 assert(pp
->x
.p
->size
== 2);
599 /* search for exact duplicate among the modulo inequalities */
601 f
= &pp
->x
.p
->arr
[1];
602 for (k
= 0; s
&& k
< s
->n
; ++k
) {
603 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
604 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
605 value_eq(s
->fixed
[k
].m
, f
->d
) &&
606 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
613 /* replace { E/m } by { (E-1)/m } + 1/m */
618 evalue_set_si(&extra
, 1, 1);
619 value_assign(extra
.d
, g
);
620 eadd(&extra
, &v
.x
.p
->arr
[1]);
621 free_evalue_refs(&extra
);
623 /* We've been going in circles; stop now */
624 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
625 free_evalue_refs(&v
);
627 evalue_set_si(&v
, 0, 1);
632 value_set_si(v
.d
, 0);
633 v
.x
.p
= new_enode(fractional
, 3, -1);
634 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
635 value_assign(v
.x
.p
->arr
[1].d
, g
);
636 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
637 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
640 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
643 value_division(f
->d
, g
, f
->d
);
644 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
645 value_assign(f
->d
, g
);
646 value_decrement(f
->x
.n
, f
->x
.n
);
647 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
649 value_gcd(g
, f
->d
, f
->x
.n
);
650 value_division(f
->d
, f
->d
, g
);
651 value_division(f
->x
.n
, f
->x
.n
, g
);
660 /* reduction may have made this fractional arg smaller */
661 i
= reorder
? p
->size
: 1;
662 for ( ; i
< p
->size
; ++i
)
663 if (value_zero_p(p
->arr
[i
].d
) &&
664 p
->arr
[i
].x
.p
->type
== fractional
&&
665 !mod_term_smaller(e
, &p
->arr
[i
]))
669 value_set_si(v
.d
, 0);
670 v
.x
.p
= new_enode(fractional
, 3, -1);
671 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
672 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
673 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
681 evalue
*pp
= &p
->arr
[0];
684 poly_denom_not_constant(&pp
, &m
);
685 mpz_fdiv_r(r
, m
, pp
->d
);
686 if (value_notzero_p(r
)) {
688 value_set_si(v
.d
, 0);
689 v
.x
.p
= new_enode(fractional
, 3, -1);
691 value_multiply(r
, m
, pp
->x
.n
);
692 value_multiply(v
.x
.p
->arr
[1].d
, m
, pp
->d
);
693 value_init(v
.x
.p
->arr
[1].x
.n
);
694 mpz_fdiv_r(v
.x
.p
->arr
[1].x
.n
, r
, pp
->d
);
695 mpz_fdiv_q(r
, r
, pp
->d
);
697 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
698 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
700 while (value_zero_p(pp
->d
))
701 pp
= &pp
->x
.p
->arr
[0];
703 value_assign(pp
->d
, m
);
704 value_assign(pp
->x
.n
, r
);
706 value_gcd(r
, pp
->d
, pp
->x
.n
);
707 value_division(pp
->d
, pp
->d
, r
);
708 value_division(pp
->x
.n
, pp
->x
.n
, r
);
721 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
722 pp
= &pp
->x
.p
->arr
[0]) {
723 f
= &pp
->x
.p
->arr
[1];
724 assert(value_pos_p(f
->d
));
725 mpz_mul_ui(twice
, f
->x
.n
, 2);
726 if (value_lt(twice
, f
->d
))
728 if (value_eq(twice
, f
->d
))
736 value_set_si(v
.d
, 0);
737 v
.x
.p
= new_enode(fractional
, 3, -1);
738 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
739 poly_denom(&p
->arr
[0], &twice
);
740 value_assign(v
.x
.p
->arr
[1].d
, twice
);
741 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
742 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
743 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
745 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
746 pp
= &pp
->x
.p
->arr
[0]) {
747 f
= &pp
->x
.p
->arr
[1];
748 value_oppose(f
->x
.n
, f
->x
.n
);
749 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
751 value_division(pp
->d
, twice
, pp
->d
);
752 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
753 value_assign(pp
->d
, twice
);
754 value_oppose(pp
->x
.n
, pp
->x
.n
);
755 value_decrement(pp
->x
.n
, pp
->x
.n
);
756 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
758 /* Maybe we should do this during reduction of
761 value_gcd(twice
, pp
->d
, pp
->x
.n
);
762 value_division(pp
->d
, pp
->d
, twice
);
763 value_division(pp
->x
.n
, pp
->x
.n
, twice
);
773 reorder_terms_about(p
, &v
);
774 _reduce_evalue(&p
->arr
[1], s
, fract
);
777 /* Try to reduce the degree */
778 for (i
=p
->size
-1;i
>=2;i
--) {
779 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
781 /* Zero coefficient */
782 free_evalue_refs(&(p
->arr
[i
]));
787 /* Try to reduce its strength */
790 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
791 free_evalue_refs(&(p
->arr
[0]));
795 else if (p
->type
== flooring
) {
796 /* Replace floor(constant) by its value */
797 if (value_notzero_p(p
->arr
[0].d
)) {
800 value_set_si(v
.d
, 1);
802 mpz_fdiv_q(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
803 reorder_terms_about(p
, &v
);
804 _reduce_evalue(&p
->arr
[1], s
, fract
);
806 /* Try to reduce the degree */
807 for (i
=p
->size
-1;i
>=2;i
--) {
808 if (!EVALUE_IS_ZERO(p
->arr
[i
]))
810 /* Zero coefficient */
811 free_evalue_refs(&(p
->arr
[i
]));
816 /* Try to reduce its strength */
819 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
820 free_evalue_refs(&(p
->arr
[0]));
824 else if (p
->type
== relation
) {
825 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
826 free_evalue_refs(&(p
->arr
[2]));
827 free_evalue_refs(&(p
->arr
[0]));
834 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
835 free_evalue_refs(&(p
->arr
[2]));
838 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
839 free_evalue_refs(&(p
->arr
[1]));
840 free_evalue_refs(&(p
->arr
[0]));
841 evalue_set_si(e
, 0, 1);
848 /* Relation was reduced by means of an identical
849 * inequality => remove
851 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
854 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
855 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
857 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
859 free_evalue_refs(&(p
->arr
[2]));
863 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
865 evalue_set_si(e
, 0, 1);
866 free_evalue_refs(&(p
->arr
[1]));
868 free_evalue_refs(&(p
->arr
[0]));
874 } /* reduce_evalue */
876 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
881 for (k
= 0; k
< dim
; ++k
)
882 if (value_notzero_p(row
[k
+1]))
885 Vector_Normalize_Positive(row
+1, dim
+1, k
);
886 assert(s
->n
< s
->max
);
887 value_init(s
->fixed
[s
->n
].d
);
888 value_init(s
->fixed
[s
->n
].m
);
889 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
890 s
->fixed
[s
->n
].pos
= k
+1;
891 value_set_si(s
->fixed
[s
->n
].m
, 0);
892 r
= &s
->fixed
[s
->n
].s
;
894 for (l
= k
+1; l
< dim
; ++l
)
895 if (value_notzero_p(row
[l
+1])) {
896 value_set_si(r
->d
, 0);
897 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
898 value_init(r
->x
.p
->arr
[1].x
.n
);
899 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
900 value_set_si(r
->x
.p
->arr
[1].d
, 1);
904 value_oppose(r
->x
.n
, row
[dim
+1]);
905 value_set_si(r
->d
, 1);
909 static void _reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
, struct subst
*s
)
912 Polyhedron
*orig
= D
;
917 D
= DomainConvex(D
, 0);
918 /* We don't perform any substitutions if the domain is a union.
919 * We may therefore miss out on some possible simplifications,
920 * e.g., if a variable is always even in the whole union,
921 * while there is a relation in the evalue that evaluates
922 * to zero for even values of the variable.
924 if (!D
->next
&& D
->NbEq
) {
928 realloc_substitution(s
, dim
);
930 int d
= relations_depth(e
);
932 NALLOC(s
->fixed
, s
->max
);
935 for (j
= 0; j
< D
->NbEq
; ++j
)
936 add_substitution(s
, D
->Constraint
[j
], dim
);
940 _reduce_evalue(e
, s
, 0);
943 for (j
= 0; j
< s
->n
; ++j
) {
944 value_clear(s
->fixed
[j
].d
);
945 value_clear(s
->fixed
[j
].m
);
946 free_evalue_refs(&s
->fixed
[j
].s
);
951 void reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
)
953 struct subst s
= { NULL
, 0, 0 };
955 if (EVALUE_IS_ZERO(*e
))
959 evalue_set_si(e
, 0, 1);
962 _reduce_evalue_in_domain(e
, D
, &s
);
967 void reduce_evalue (evalue
*e
) {
968 struct subst s
= { NULL
, 0, 0 };
970 if (value_notzero_p(e
->d
))
971 return; /* a rational number, its already reduced */
973 if (e
->x
.p
->type
== partition
) {
976 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
977 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
979 /* This shouldn't really happen;
980 * Empty domains should not be added.
982 POL_ENSURE_VERTICES(D
);
984 _reduce_evalue_in_domain(&e
->x
.p
->arr
[2*i
+1], D
, &s
);
986 if (emptyQ(D
) || EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
987 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
988 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
989 value_clear(e
->x
.p
->arr
[2*i
].d
);
991 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
992 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
996 if (e
->x
.p
->size
== 0) {
998 evalue_set_si(e
, 0, 1);
1001 _reduce_evalue(e
, &s
, 0);
1006 static void print_evalue_r(FILE *DST
, const evalue
*e
, const char *const *pname
)
1008 if (EVALUE_IS_NAN(*e
)) {
1009 fprintf(DST
, "NaN");
1013 if(value_notzero_p(e
->d
)) {
1014 if(value_notone_p(e
->d
)) {
1015 value_print(DST
,VALUE_FMT
,e
->x
.n
);
1017 value_print(DST
,VALUE_FMT
,e
->d
);
1020 value_print(DST
,VALUE_FMT
,e
->x
.n
);
1024 print_enode(DST
,e
->x
.p
,pname
);
1026 } /* print_evalue */
1028 void print_evalue(FILE *DST
, const evalue
*e
, const char * const *pname
)
1030 print_evalue_r(DST
, e
, pname
);
1031 if (value_notzero_p(e
->d
))
1035 void print_enode(FILE *DST
, enode
*p
, const char *const *pname
)
1040 fprintf(DST
, "NULL");
1046 for (i
=0; i
<p
->size
; i
++) {
1047 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1051 fprintf(DST
, " }\n");
1055 for (i
=p
->size
-1; i
>=0; i
--) {
1056 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1057 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
1059 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
1061 fprintf(DST
, " )\n");
1065 for (i
=0; i
<p
->size
; i
++) {
1066 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1067 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
1069 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
1074 for (i
=p
->size
-1; i
>=1; i
--) {
1075 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1077 fprintf(DST
, " * ");
1078 fprintf(DST
, p
->type
== flooring
? "[" : "{");
1079 print_evalue_r(DST
, &p
->arr
[0], pname
);
1080 fprintf(DST
, p
->type
== flooring
? "]" : "}");
1082 fprintf(DST
, "^%d + ", i
-1);
1084 fprintf(DST
, " + ");
1087 fprintf(DST
, " )\n");
1091 print_evalue_r(DST
, &p
->arr
[0], pname
);
1092 fprintf(DST
, "= 0 ] * \n");
1093 print_evalue_r(DST
, &p
->arr
[1], pname
);
1095 fprintf(DST
, " +\n [ ");
1096 print_evalue_r(DST
, &p
->arr
[0], pname
);
1097 fprintf(DST
, "!= 0 ] * \n");
1098 print_evalue_r(DST
, &p
->arr
[2], pname
);
1102 char **new_names
= NULL
;
1103 const char *const *names
= pname
;
1104 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
1105 if (!pname
|| p
->pos
< maxdim
) {
1106 new_names
= ALLOCN(char *, maxdim
);
1107 for (i
= 0; i
< p
->pos
; ++i
) {
1109 new_names
[i
] = (char *)pname
[i
];
1111 new_names
[i
] = ALLOCN(char, 10);
1112 snprintf(new_names
[i
], 10, "%c", 'P'+i
);
1115 for ( ; i
< maxdim
; ++i
) {
1116 new_names
[i
] = ALLOCN(char, 10);
1117 snprintf(new_names
[i
], 10, "_p%d", i
);
1119 names
= (const char**)new_names
;
1122 for (i
=0; i
<p
->size
/2; i
++) {
1123 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
1124 print_evalue_r(DST
, &p
->arr
[2*i
+1], names
);
1125 if (value_notzero_p(p
->arr
[2*i
+1].d
))
1129 if (!pname
|| p
->pos
< maxdim
) {
1130 for (i
= pname
? p
->pos
: 0; i
< maxdim
; ++i
)
1144 * 0 if toplevels of e1 and e2 are at the same level
1145 * <0 if toplevel of e1 should be outside of toplevel of e2
1146 * >0 if toplevel of e2 should be outside of toplevel of e1
1148 static int evalue_level_cmp(const evalue
*e1
, const evalue
*e2
)
1150 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
))
1152 if (value_notzero_p(e1
->d
))
1154 if (value_notzero_p(e2
->d
))
1156 if (e1
->x
.p
->type
== partition
&& e2
->x
.p
->type
== partition
)
1158 if (e1
->x
.p
->type
== partition
)
1160 if (e2
->x
.p
->type
== partition
)
1162 if (e1
->x
.p
->type
== relation
&& e2
->x
.p
->type
== relation
) {
1163 if (eequal(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]))
1165 if (mod_term_smaller(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]))
1170 if (e1
->x
.p
->type
== relation
)
1172 if (e2
->x
.p
->type
== relation
)
1174 if (e1
->x
.p
->type
== polynomial
&& e2
->x
.p
->type
== polynomial
)
1175 return e1
->x
.p
->pos
- e2
->x
.p
->pos
;
1176 if (e1
->x
.p
->type
== polynomial
)
1178 if (e2
->x
.p
->type
== polynomial
)
1180 if (e1
->x
.p
->type
== periodic
&& e2
->x
.p
->type
== periodic
)
1181 return e1
->x
.p
->pos
- e2
->x
.p
->pos
;
1182 assert(e1
->x
.p
->type
!= periodic
);
1183 assert(e2
->x
.p
->type
!= periodic
);
1184 assert(e1
->x
.p
->type
== e2
->x
.p
->type
);
1185 if (eequal(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]))
1187 if (mod_term_smaller(e1
, e2
))
1193 static void eadd_rev(const evalue
*e1
, evalue
*res
)
1197 evalue_copy(&ev
, e1
);
1199 free_evalue_refs(res
);
1203 static void eadd_rev_cst(const evalue
*e1
, evalue
*res
)
1207 evalue_copy(&ev
, e1
);
1208 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
1209 free_evalue_refs(res
);
1213 struct section
{ Polyhedron
* D
; evalue E
; };
1215 void eadd_partitions(const evalue
*e1
, evalue
*res
)
1220 s
= (struct section
*)
1221 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
1222 sizeof(struct section
));
1224 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1225 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1226 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1229 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1230 assert(res
->x
.p
->size
>= 2);
1231 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1232 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
1234 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
1236 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1241 fd
= DomainConstraintSimplify(fd
, 0);
1246 value_init(s
[n
].E
.d
);
1247 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
1251 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1252 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
1253 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1255 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1256 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1257 d
= DomainConstraintSimplify(d
, 0);
1263 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
1264 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1266 value_init(s
[n
].E
.d
);
1267 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1268 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1273 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
1277 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1280 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1281 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1282 value_clear(res
->x
.p
->arr
[2*i
].d
);
1287 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1288 for (j
= 0; j
< n
; ++j
) {
1289 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1290 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1291 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1297 static void explicit_complement(evalue
*res
)
1299 enode
*rel
= new_enode(relation
, 3, 0);
1301 value_clear(rel
->arr
[0].d
);
1302 rel
->arr
[0] = res
->x
.p
->arr
[0];
1303 value_clear(rel
->arr
[1].d
);
1304 rel
->arr
[1] = res
->x
.p
->arr
[1];
1305 value_set_si(rel
->arr
[2].d
, 1);
1306 value_init(rel
->arr
[2].x
.n
);
1307 value_set_si(rel
->arr
[2].x
.n
, 0);
1312 static void reduce_constant(evalue
*e
)
1317 value_gcd(g
, e
->x
.n
, e
->d
);
1318 if (value_notone_p(g
)) {
1319 value_division(e
->d
, e
->d
,g
);
1320 value_division(e
->x
.n
, e
->x
.n
,g
);
1325 /* Add two rational numbers */
1326 static void eadd_rationals(const evalue
*e1
, evalue
*res
)
1328 if (value_eq(e1
->d
, res
->d
))
1329 value_addto(res
->x
.n
, res
->x
.n
, e1
->x
.n
);
1331 value_multiply(res
->x
.n
, res
->x
.n
, e1
->d
);
1332 value_addmul(res
->x
.n
, e1
->x
.n
, res
->d
);
1333 value_multiply(res
->d
,e1
->d
,res
->d
);
1335 reduce_constant(res
);
1338 static void eadd_relations(const evalue
*e1
, evalue
*res
)
1342 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1343 explicit_complement(res
);
1344 for (i
= 1; i
< e1
->x
.p
->size
; ++i
)
1345 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1348 static void eadd_arrays(const evalue
*e1
, evalue
*res
, int n
)
1352 // add any element in e1 to the corresponding element in res
1353 i
= type_offset(res
->x
.p
);
1355 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1357 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1360 static void eadd_poly(const evalue
*e1
, evalue
*res
)
1362 if (e1
->x
.p
->size
> res
->x
.p
->size
)
1365 eadd_arrays(e1
, res
, e1
->x
.p
->size
);
1369 * Product or sum of two periodics of the same parameter
1370 * and different periods
1372 static void combine_periodics(const evalue
*e1
, evalue
*res
,
1373 void (*op
)(const evalue
*, evalue
*))
1381 value_set_si(es
, e1
->x
.p
->size
);
1382 value_set_si(rs
, res
->x
.p
->size
);
1383 value_lcm(rs
, es
, rs
);
1384 size
= (int)mpz_get_si(rs
);
1387 p
= new_enode(periodic
, size
, e1
->x
.p
->pos
);
1388 for (i
= 0; i
< res
->x
.p
->size
; i
++) {
1389 value_clear(p
->arr
[i
].d
);
1390 p
->arr
[i
] = res
->x
.p
->arr
[i
];
1392 for (i
= res
->x
.p
->size
; i
< size
; i
++)
1393 evalue_copy(&p
->arr
[i
], &res
->x
.p
->arr
[i
% res
->x
.p
->size
]);
1394 for (i
= 0; i
< size
; i
++)
1395 op(&e1
->x
.p
->arr
[i
% e1
->x
.p
->size
], &p
->arr
[i
]);
1400 static void eadd_periodics(const evalue
*e1
, evalue
*res
)
1406 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1407 eadd_arrays(e1
, res
, e1
->x
.p
->size
);
1411 combine_periodics(e1
, res
, eadd
);
1414 void evalue_assign(evalue
*dst
, const evalue
*src
)
1416 if (value_pos_p(dst
->d
) && value_pos_p(src
->d
)) {
1417 value_assign(dst
->d
, src
->d
);
1418 value_assign(dst
->x
.n
, src
->x
.n
);
1421 free_evalue_refs(dst
);
1423 evalue_copy(dst
, src
);
1426 void eadd(const evalue
*e1
, evalue
*res
)
1430 if (EVALUE_IS_ZERO(*e1
))
1433 if (EVALUE_IS_NAN(*res
))
1436 if (EVALUE_IS_NAN(*e1
)) {
1437 evalue_assign(res
, e1
);
1441 if (EVALUE_IS_ZERO(*res
)) {
1442 evalue_assign(res
, e1
);
1446 cmp
= evalue_level_cmp(res
, e1
);
1448 switch (e1
->x
.p
->type
) {
1452 eadd_rev_cst(e1
, res
);
1457 } else if (cmp
== 0) {
1458 if (value_notzero_p(e1
->d
)) {
1459 eadd_rationals(e1
, res
);
1461 switch (e1
->x
.p
->type
) {
1463 eadd_partitions(e1
, res
);
1466 eadd_relations(e1
, res
);
1469 assert(e1
->x
.p
->size
== res
->x
.p
->size
);
1476 eadd_periodics(e1
, res
);
1484 switch (res
->x
.p
->type
) {
1488 /* Add to the constant term of a polynomial */
1489 eadd(e1
, &res
->x
.p
->arr
[type_offset(res
->x
.p
)]);
1492 /* Add to all elements of a periodic number */
1493 for (i
= 0; i
< res
->x
.p
->size
; i
++)
1494 eadd(e1
, &res
->x
.p
->arr
[i
]);
1497 fprintf(stderr
, "eadd: cannot add const with vector\n");
1502 /* Create (zero) complement if needed */
1503 if (res
->x
.p
->size
< 3)
1504 explicit_complement(res
);
1505 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1506 eadd(e1
, &res
->x
.p
->arr
[i
]);
1514 static void emul_rev(const evalue
*e1
, evalue
*res
)
1518 evalue_copy(&ev
, e1
);
1520 free_evalue_refs(res
);
1524 static void emul_poly(const evalue
*e1
, evalue
*res
)
1526 int i
, j
, offset
= type_offset(res
->x
.p
);
1529 int size
= (e1
->x
.p
->size
+ res
->x
.p
->size
- offset
- 1);
1531 p
= new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1533 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1534 if (!EVALUE_IS_ZERO(e1
->x
.p
->arr
[i
]))
1537 /* special case pure power */
1538 if (i
== e1
->x
.p
->size
-1) {
1540 value_clear(p
->arr
[0].d
);
1541 p
->arr
[0] = res
->x
.p
->arr
[0];
1543 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1544 evalue_set_si(&p
->arr
[i
], 0, 1);
1545 for (i
= offset
; i
< res
->x
.p
->size
; ++i
) {
1546 value_clear(p
->arr
[i
+e1
->x
.p
->size
-offset
-1].d
);
1547 p
->arr
[i
+e1
->x
.p
->size
-offset
-1] = res
->x
.p
->arr
[i
];
1548 emul(&e1
->x
.p
->arr
[e1
->x
.p
->size
-1],
1549 &p
->arr
[i
+e1
->x
.p
->size
-offset
-1]);
1557 value_set_si(tmp
.d
,0);
1560 evalue_copy(&p
->arr
[0], &e1
->x
.p
->arr
[0]);
1561 for (i
= offset
; i
< e1
->x
.p
->size
; i
++) {
1562 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1563 emul(&res
->x
.p
->arr
[offset
], &tmp
.x
.p
->arr
[i
]);
1566 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1567 for (i
= offset
+1; i
<res
->x
.p
->size
; i
++)
1568 for (j
= offset
; j
<e1
->x
.p
->size
; j
++) {
1571 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1572 emul(&res
->x
.p
->arr
[i
], &ev
);
1573 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-offset
]);
1574 free_evalue_refs(&ev
);
1576 free_evalue_refs(res
);
1580 void emul_partitions(const evalue
*e1
, evalue
*res
)
1585 s
= (struct section
*)
1586 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1587 sizeof(struct section
));
1589 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1590 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1591 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1594 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1595 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1596 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1597 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1598 d
= DomainConstraintSimplify(d
, 0);
1604 /* This code is only needed because the partitions
1605 are not true partitions.
1607 for (k
= 0; k
< n
; ++k
) {
1608 if (DomainIncludes(s
[k
].D
, d
))
1610 if (DomainIncludes(d
, s
[k
].D
)) {
1611 Domain_Free(s
[k
].D
);
1612 free_evalue_refs(&s
[k
].E
);
1623 value_init(s
[n
].E
.d
);
1624 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1625 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1629 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1630 value_clear(res
->x
.p
->arr
[2*i
].d
);
1631 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1636 evalue_set_si(res
, 0, 1);
1638 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1639 for (j
= 0; j
< n
; ++j
) {
1640 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1641 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1642 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1649 /* Product of two rational numbers */
1650 static void emul_rationals(const evalue
*e1
, evalue
*res
)
1652 value_multiply(res
->d
, e1
->d
, res
->d
);
1653 value_multiply(res
->x
.n
, e1
->x
.n
, res
->x
.n
);
1654 reduce_constant(res
);
1657 static void emul_relations(const evalue
*e1
, evalue
*res
)
1661 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3) {
1662 free_evalue_refs(&res
->x
.p
->arr
[2]);
1665 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1666 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1669 static void emul_periodics(const evalue
*e1
, evalue
*res
)
1676 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1677 /* Product of two periodics of the same parameter and period */
1678 for (i
= 0; i
< res
->x
.p
->size
; i
++)
1679 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1683 combine_periodics(e1
, res
, emul
);
1686 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1688 static void emul_fractionals(const evalue
*e1
, evalue
*res
)
1692 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1693 if (!value_two_p(d
.d
))
1698 value_set_si(d
.x
.n
, 1);
1699 /* { x }^2 == { x }/2 */
1700 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1701 assert(e1
->x
.p
->size
== 3);
1702 assert(res
->x
.p
->size
== 3);
1704 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1706 eadd(&res
->x
.p
->arr
[1], &tmp
);
1707 emul(&e1
->x
.p
->arr
[2], &tmp
);
1708 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[1]);
1709 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[2]);
1710 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1711 free_evalue_refs(&tmp
);
1717 /* Computes the product of two evalues "e1" and "res" and puts
1718 * the result in "res". You need to make a copy of "res"
1719 * before calling this function if you still need it afterward.
1720 * The vector type of evalues is not treated here
1722 void emul(const evalue
*e1
, evalue
*res
)
1726 assert(!(value_zero_p(e1
->d
) && e1
->x
.p
->type
== evector
));
1727 assert(!(value_zero_p(res
->d
) && res
->x
.p
->type
== evector
));
1729 if (EVALUE_IS_ZERO(*res
))
1732 if (EVALUE_IS_ONE(*e1
))
1735 if (EVALUE_IS_ZERO(*e1
)) {
1736 evalue_assign(res
, e1
);
1740 if (EVALUE_IS_NAN(*res
))
1743 if (EVALUE_IS_NAN(*e1
)) {
1744 evalue_assign(res
, e1
);
1748 cmp
= evalue_level_cmp(res
, e1
);
1751 } else if (cmp
== 0) {
1752 if (value_notzero_p(e1
->d
)) {
1753 emul_rationals(e1
, res
);
1755 switch (e1
->x
.p
->type
) {
1757 emul_partitions(e1
, res
);
1760 emul_relations(e1
, res
);
1767 emul_periodics(e1
, res
);
1770 emul_fractionals(e1
, res
);
1776 switch (res
->x
.p
->type
) {
1778 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1779 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1786 for (i
= type_offset(res
->x
.p
); i
< res
->x
.p
->size
; ++i
)
1787 emul(e1
, &res
->x
.p
->arr
[i
]);
1793 /* Frees mask content ! */
1794 void emask(evalue
*mask
, evalue
*res
) {
1801 if (EVALUE_IS_ZERO(*res
)) {
1802 free_evalue_refs(mask
);
1806 assert(value_zero_p(mask
->d
));
1807 assert(mask
->x
.p
->type
== partition
);
1808 assert(value_zero_p(res
->d
));
1809 assert(res
->x
.p
->type
== partition
);
1810 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1811 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1812 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1813 pos
= res
->x
.p
->pos
;
1815 s
= (struct section
*)
1816 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1817 sizeof(struct section
));
1821 evalue_set_si(&mone
, -1, 1);
1824 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1825 assert(mask
->x
.p
->size
>= 2);
1826 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1827 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1829 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1831 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1840 value_init(s
[n
].E
.d
);
1841 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1845 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1846 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1849 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1850 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1851 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1852 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1854 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1855 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1861 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1862 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1864 value_init(s
[n
].E
.d
);
1865 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1866 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1872 /* Just ignore; this may have been previously masked off */
1874 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1878 free_evalue_refs(&mone
);
1879 free_evalue_refs(mask
);
1880 free_evalue_refs(res
);
1883 evalue_set_si(res
, 0, 1);
1885 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1886 for (j
= 0; j
< n
; ++j
) {
1887 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1888 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1889 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1896 void evalue_copy(evalue
*dst
, const evalue
*src
)
1898 value_assign(dst
->d
, src
->d
);
1899 if (EVALUE_IS_NAN(*dst
)) {
1903 if (value_pos_p(src
->d
)) {
1904 value_init(dst
->x
.n
);
1905 value_assign(dst
->x
.n
, src
->x
.n
);
1907 dst
->x
.p
= ecopy(src
->x
.p
);
1910 evalue
*evalue_dup(const evalue
*e
)
1912 evalue
*res
= ALLOC(evalue
);
1914 evalue_copy(res
, e
);
1918 enode
*new_enode(enode_type type
,int size
,int pos
) {
1924 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1927 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1931 for(i
=0; i
<size
; i
++) {
1932 value_init(res
->arr
[i
].d
);
1933 value_set_si(res
->arr
[i
].d
,0);
1934 res
->arr
[i
].x
.p
= 0;
1939 enode
*ecopy(enode
*e
) {
1944 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1945 for(i
=0;i
<e
->size
;++i
) {
1946 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1947 if(value_zero_p(res
->arr
[i
].d
))
1948 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1949 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1950 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1952 value_init(res
->arr
[i
].x
.n
);
1953 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1959 int ecmp(const evalue
*e1
, const evalue
*e2
)
1965 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
1969 value_multiply(m
, e1
->x
.n
, e2
->d
);
1970 value_multiply(m2
, e2
->x
.n
, e1
->d
);
1972 if (value_lt(m
, m2
))
1974 else if (value_gt(m
, m2
))
1984 if (value_notzero_p(e1
->d
))
1986 if (value_notzero_p(e2
->d
))
1992 if (p1
->type
!= p2
->type
)
1993 return p1
->type
- p2
->type
;
1994 if (p1
->pos
!= p2
->pos
)
1995 return p1
->pos
- p2
->pos
;
1996 if (p1
->size
!= p2
->size
)
1997 return p1
->size
- p2
->size
;
1999 for (i
= p1
->size
-1; i
>= 0; --i
)
2000 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
2006 int eequal(const evalue
*e1
, const evalue
*e2
)
2011 if (value_ne(e1
->d
,e2
->d
))
2014 if (EVALUE_IS_DOMAIN(*e1
))
2015 return PolyhedronIncludes(EVALUE_DOMAIN(*e2
), EVALUE_DOMAIN(*e1
)) &&
2016 PolyhedronIncludes(EVALUE_DOMAIN(*e1
), EVALUE_DOMAIN(*e2
));
2018 if (EVALUE_IS_NAN(*e1
))
2021 assert(value_posz_p(e1
->d
));
2023 /* e1->d == e2->d */
2024 if (value_notzero_p(e1
->d
)) {
2025 if (value_ne(e1
->x
.n
,e2
->x
.n
))
2028 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2032 /* e1->d == e2->d == 0 */
2035 if (p1
->type
!= p2
->type
) return 0;
2036 if (p1
->size
!= p2
->size
) return 0;
2037 if (p1
->pos
!= p2
->pos
) return 0;
2038 for (i
=0; i
<p1
->size
; i
++)
2039 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
2044 void free_evalue_refs(evalue
*e
) {
2049 if (EVALUE_IS_NAN(*e
)) {
2054 if (EVALUE_IS_DOMAIN(*e
)) {
2055 Domain_Free(EVALUE_DOMAIN(*e
));
2058 } else if (value_pos_p(e
->d
)) {
2060 /* 'e' stores a constant */
2062 value_clear(e
->x
.n
);
2065 assert(value_zero_p(e
->d
));
2068 if (!p
) return; /* null pointer */
2069 for (i
=0; i
<p
->size
; i
++) {
2070 free_evalue_refs(&(p
->arr
[i
]));
2074 } /* free_evalue_refs */
2076 void evalue_free(evalue
*e
)
2078 free_evalue_refs(e
);
2082 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
2083 Vector
* val
, evalue
*res
)
2085 unsigned nparam
= periods
->Size
;
2088 double d
= compute_evalue(e
, val
->p
);
2089 d
*= VALUE_TO_DOUBLE(m
);
2094 value_assign(res
->d
, m
);
2095 value_init(res
->x
.n
);
2096 value_set_double(res
->x
.n
, d
);
2097 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
2100 if (value_one_p(periods
->p
[p
]))
2101 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
2106 value_assign(tmp
, periods
->p
[p
]);
2107 value_set_si(res
->d
, 0);
2108 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
2110 value_decrement(tmp
, tmp
);
2111 value_assign(val
->p
[p
], tmp
);
2112 mod2table_r(e
, periods
, m
, p
+1, val
,
2113 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
2114 } while (value_pos_p(tmp
));
2120 static void rel2table(evalue
*e
, int zero
)
2122 if (value_pos_p(e
->d
)) {
2123 if (value_zero_p(e
->x
.n
) == zero
)
2124 value_set_si(e
->x
.n
, 1);
2126 value_set_si(e
->x
.n
, 0);
2127 value_set_si(e
->d
, 1);
2130 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
2131 rel2table(&e
->x
.p
->arr
[i
], zero
);
2135 void evalue_mod2table(evalue
*e
, int nparam
)
2140 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2143 for (i
=0; i
<p
->size
; i
++) {
2144 evalue_mod2table(&(p
->arr
[i
]), nparam
);
2146 if (p
->type
== relation
) {
2151 evalue_copy(©
, &p
->arr
[0]);
2153 rel2table(&p
->arr
[0], 1);
2154 emul(&p
->arr
[0], &p
->arr
[1]);
2156 rel2table(©
, 0);
2157 emul(©
, &p
->arr
[2]);
2158 eadd(&p
->arr
[2], &p
->arr
[1]);
2159 free_evalue_refs(&p
->arr
[2]);
2160 free_evalue_refs(©
);
2162 free_evalue_refs(&p
->arr
[0]);
2166 } else if (p
->type
== fractional
) {
2167 Vector
*periods
= Vector_Alloc(nparam
);
2168 Vector
*val
= Vector_Alloc(nparam
);
2174 value_set_si(tmp
, 1);
2175 Vector_Set(periods
->p
, 1, nparam
);
2176 Vector_Set(val
->p
, 0, nparam
);
2177 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
2180 assert(p
->type
== polynomial
);
2181 assert(p
->size
== 2);
2182 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
2183 value_lcm(tmp
, tmp
, p
->arr
[1].d
);
2185 value_lcm(tmp
, tmp
, ev
->d
);
2187 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
2190 evalue_set_si(&res
, 0, 1);
2191 /* Compute the polynomial using Horner's rule */
2192 for (i
=p
->size
-1;i
>1;i
--) {
2193 eadd(&p
->arr
[i
], &res
);
2196 eadd(&p
->arr
[1], &res
);
2198 free_evalue_refs(e
);
2199 free_evalue_refs(&EP
);
2204 Vector_Free(periods
);
2206 } /* evalue_mod2table */
2208 /********************************************************/
2209 /* function in domain */
2210 /* check if the parameters in list_args */
2211 /* verifies the constraints of Domain P */
2212 /********************************************************/
2213 int in_domain(Polyhedron
*P
, Value
*list_args
)
2216 Value v
; /* value of the constraint of a row when
2217 parameters are instantiated*/
2221 for (row
= 0; row
< P
->NbConstraints
; row
++) {
2222 Inner_Product(P
->Constraint
[row
]+1, list_args
, P
->Dimension
, &v
);
2223 value_addto(v
, v
, P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
2224 if (value_neg_p(v
) ||
2225 value_zero_p(P
->Constraint
[row
][0]) && value_notzero_p(v
)) {
2232 return in
|| (P
->next
&& in_domain(P
->next
, list_args
));
2235 /****************************************************/
2236 /* function compute enode */
2237 /* compute the value of enode p with parameters */
2238 /* list "list_args */
2239 /* compute the polynomial or the periodic */
2240 /****************************************************/
2242 static double compute_enode(enode
*p
, Value
*list_args
) {
2254 if (p
->type
== polynomial
) {
2256 value_assign(param
,list_args
[p
->pos
-1]);
2258 /* Compute the polynomial using Horner's rule */
2259 for (i
=p
->size
-1;i
>0;i
--) {
2260 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2261 res
*=VALUE_TO_DOUBLE(param
);
2263 res
+=compute_evalue(&p
->arr
[0],list_args
);
2265 else if (p
->type
== fractional
) {
2266 double d
= compute_evalue(&p
->arr
[0], list_args
);
2267 d
-= floor(d
+1e-10);
2269 /* Compute the polynomial using Horner's rule */
2270 for (i
=p
->size
-1;i
>1;i
--) {
2271 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2274 res
+=compute_evalue(&p
->arr
[1],list_args
);
2276 else if (p
->type
== flooring
) {
2277 double d
= compute_evalue(&p
->arr
[0], list_args
);
2280 /* Compute the polynomial using Horner's rule */
2281 for (i
=p
->size
-1;i
>1;i
--) {
2282 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2285 res
+=compute_evalue(&p
->arr
[1],list_args
);
2287 else if (p
->type
== periodic
) {
2288 value_assign(m
,list_args
[p
->pos
-1]);
2290 /* Choose the right element of the periodic */
2291 value_set_si(param
,p
->size
);
2292 value_pmodulus(m
,m
,param
);
2293 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2295 else if (p
->type
== relation
) {
2296 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2297 res
= compute_evalue(&p
->arr
[1], list_args
);
2298 else if (p
->size
> 2)
2299 res
= compute_evalue(&p
->arr
[2], list_args
);
2301 else if (p
->type
== partition
) {
2302 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2303 Value
*vals
= list_args
;
2306 for (i
= 0; i
< dim
; ++i
) {
2307 value_init(vals
[i
]);
2309 value_assign(vals
[i
], list_args
[i
]);
2312 for (i
= 0; i
< p
->size
/2; ++i
)
2313 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2314 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2318 for (i
= 0; i
< dim
; ++i
)
2319 value_clear(vals
[i
]);
2328 } /* compute_enode */
2330 /*************************************************/
2331 /* return the value of Ehrhart Polynomial */
2332 /* It returns a double, because since it is */
2333 /* a recursive function, some intermediate value */
2334 /* might not be integral */
2335 /*************************************************/
2337 double compute_evalue(const evalue
*e
, Value
*list_args
)
2341 if (value_notzero_p(e
->d
)) {
2342 if (value_notone_p(e
->d
))
2343 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2345 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2348 res
= compute_enode(e
->x
.p
,list_args
);
2350 } /* compute_evalue */
2353 /****************************************************/
2354 /* function compute_poly : */
2355 /* Check for the good validity domain */
2356 /* return the number of point in the Polyhedron */
2357 /* in allocated memory */
2358 /* Using the Ehrhart pseudo-polynomial */
2359 /****************************************************/
2360 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2363 /* double d; int i; */
2365 tmp
= (Value
*) malloc (sizeof(Value
));
2366 assert(tmp
!= NULL
);
2368 value_set_si(*tmp
,0);
2371 return(tmp
); /* no ehrhart polynomial */
2372 if(en
->ValidityDomain
) {
2373 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2374 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2379 return(tmp
); /* no Validity Domain */
2381 if(in_domain(en
->ValidityDomain
,list_args
)) {
2383 #ifdef EVAL_EHRHART_DEBUG
2384 Print_Domain(stdout
,en
->ValidityDomain
);
2385 print_evalue(stdout
,&en
->EP
);
2388 /* d = compute_evalue(&en->EP,list_args);
2390 printf("(double)%lf = %d\n", d, i ); */
2391 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2397 value_set_si(*tmp
,0);
2398 return(tmp
); /* no compatible domain with the arguments */
2399 } /* compute_poly */
2401 static evalue
*eval_polynomial(const enode
*p
, int offset
,
2402 evalue
*base
, Value
*values
)
2407 res
= evalue_zero();
2408 for (i
= p
->size
-1; i
> offset
; --i
) {
2409 c
= evalue_eval(&p
->arr
[i
], values
);
2414 c
= evalue_eval(&p
->arr
[offset
], values
);
2421 evalue
*evalue_eval(const evalue
*e
, Value
*values
)
2428 if (value_notzero_p(e
->d
)) {
2429 res
= ALLOC(evalue
);
2431 evalue_copy(res
, e
);
2434 switch (e
->x
.p
->type
) {
2436 value_init(param
.x
.n
);
2437 value_assign(param
.x
.n
, values
[e
->x
.p
->pos
-1]);
2438 value_init(param
.d
);
2439 value_set_si(param
.d
, 1);
2441 res
= eval_polynomial(e
->x
.p
, 0, ¶m
, values
);
2442 free_evalue_refs(¶m
);
2445 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2446 mpz_fdiv_r(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2448 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2449 evalue_free(param2
);
2452 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2453 mpz_fdiv_q(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2454 value_set_si(param2
->d
, 1);
2456 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2457 evalue_free(param2
);
2460 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2461 if (value_zero_p(param2
->x
.n
))
2462 res
= evalue_eval(&e
->x
.p
->arr
[1], values
);
2463 else if (e
->x
.p
->size
> 2)
2464 res
= evalue_eval(&e
->x
.p
->arr
[2], values
);
2466 res
= evalue_zero();
2467 evalue_free(param2
);
2470 assert(e
->x
.p
->pos
== EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
);
2471 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2472 if (in_domain(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), values
)) {
2473 res
= evalue_eval(&e
->x
.p
->arr
[2*i
+1], values
);
2477 res
= evalue_zero();
2485 size_t value_size(Value v
) {
2486 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2487 * sizeof(v
[0]._mp_d
[0]);
2490 size_t domain_size(Polyhedron
*D
)
2493 size_t s
= sizeof(*D
);
2495 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2496 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2497 s
+= value_size(D
->Constraint
[i
][j
]);
2500 for (i = 0; i < D->NbRays; ++i)
2501 for (j = 0; j < D->Dimension+2; ++j)
2502 s += value_size(D->Ray[i][j]);
2505 return D
->next
? s
+domain_size(D
->next
) : s
;
2508 size_t enode_size(enode
*p
) {
2509 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2512 if (p
->type
== partition
)
2513 for (i
= 0; i
< p
->size
/2; ++i
) {
2514 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2515 s
+= evalue_size(&p
->arr
[2*i
+1]);
2518 for (i
= 0; i
< p
->size
; ++i
) {
2519 s
+= evalue_size(&p
->arr
[i
]);
2524 size_t evalue_size(evalue
*e
)
2526 size_t s
= sizeof(*e
);
2527 s
+= value_size(e
->d
);
2528 if (value_notzero_p(e
->d
))
2529 s
+= value_size(e
->x
.n
);
2531 s
+= enode_size(e
->x
.p
);
2535 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2537 evalue
*found
= NULL
;
2542 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2545 value_init(offset
.d
);
2546 value_init(offset
.x
.n
);
2547 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2548 value_lcm(offset
.d
, m
, offset
.d
);
2549 value_set_si(offset
.x
.n
, 1);
2552 evalue_copy(©
, cst
);
2555 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2557 if (eequal(base
, &e
->x
.p
->arr
[0]))
2558 found
= &e
->x
.p
->arr
[0];
2560 value_set_si(offset
.x
.n
, -2);
2563 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2565 if (eequal(base
, &e
->x
.p
->arr
[0]))
2568 free_evalue_refs(cst
);
2569 free_evalue_refs(&offset
);
2572 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2573 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2578 static evalue
*find_relation_pair(evalue
*e
)
2581 evalue
*found
= NULL
;
2583 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2586 if (e
->x
.p
->type
== fractional
) {
2591 poly_denom(&e
->x
.p
->arr
[0], &m
);
2593 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2594 cst
= &cst
->x
.p
->arr
[0])
2597 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2598 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2603 i
= e
->x
.p
->type
== relation
;
2604 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2605 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2610 void evalue_mod2relation(evalue
*e
) {
2613 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2616 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2617 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2618 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2619 value_clear(e
->x
.p
->arr
[2*i
].d
);
2620 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2622 if (2*i
< e
->x
.p
->size
) {
2623 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2624 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2629 if (e
->x
.p
->size
== 0) {
2631 evalue_set_si(e
, 0, 1);
2637 while ((d
= find_relation_pair(e
)) != NULL
) {
2641 value_init(split
.d
);
2642 value_set_si(split
.d
, 0);
2643 split
.x
.p
= new_enode(relation
, 3, 0);
2644 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2645 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2647 ev
= &split
.x
.p
->arr
[0];
2648 value_set_si(ev
->d
, 0);
2649 ev
->x
.p
= new_enode(fractional
, 3, -1);
2650 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2651 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2652 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2658 free_evalue_refs(&split
);
2662 static int evalue_comp(const void * a
, const void * b
)
2664 const evalue
*e1
= *(const evalue
**)a
;
2665 const evalue
*e2
= *(const evalue
**)b
;
2666 return ecmp(e1
, e2
);
2669 void evalue_combine(evalue
*e
)
2676 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2679 NALLOC(evs
, e
->x
.p
->size
/2);
2680 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2681 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2682 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2683 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2684 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2685 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2686 value_clear(p
->arr
[2*k
].d
);
2687 value_clear(p
->arr
[2*k
+1].d
);
2688 p
->arr
[2*k
] = *(evs
[i
]-1);
2689 p
->arr
[2*k
+1] = *(evs
[i
]);
2692 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2695 value_clear((evs
[i
]-1)->d
);
2699 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2700 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2701 free_evalue_refs(evs
[i
]);
2705 for (i
= 2*k
; i
< p
->size
; ++i
)
2706 value_clear(p
->arr
[i
].d
);
2713 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2715 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2717 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2720 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2721 Polyhedron
*D
, *N
, **P
;
2724 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2731 if (D
->NbEq
<= H
->NbEq
) {
2737 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2738 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2739 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2740 reduce_evalue(&tmp
);
2741 if (value_notzero_p(tmp
.d
) ||
2742 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2745 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2746 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2749 free_evalue_refs(&tmp
);
2755 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2757 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2759 value_clear(e
->x
.p
->arr
[2*i
].d
);
2760 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2762 if (2*i
< e
->x
.p
->size
) {
2763 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2764 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2771 H
= DomainConvex(D
, 0);
2772 E
= DomainDifference(H
, D
, 0);
2774 D
= DomainDifference(H
, E
, 0);
2777 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2781 /* Use smallest representative for coefficients in affine form in
2782 * argument of fractional.
2783 * Since any change will make the argument non-standard,
2784 * the containing evalue will have to be reduced again afterward.
2786 static void fractional_minimal_coefficients(enode
*p
)
2792 assert(p
->type
== fractional
);
2794 while (value_zero_p(pp
->d
)) {
2795 assert(pp
->x
.p
->type
== polynomial
);
2796 assert(pp
->x
.p
->size
== 2);
2797 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2798 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2799 if (value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2800 value_subtract(pp
->x
.p
->arr
[1].x
.n
,
2801 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2802 pp
= &pp
->x
.p
->arr
[0];
2808 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2813 unsigned dim
= D
->Dimension
;
2814 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2817 assert(p
->type
== fractional
|| p
->type
== flooring
);
2818 value_set_si(T
->p
[1][dim
], 1);
2819 evalue_extract_affine(&p
->arr
[0], T
->p
[0], &T
->p
[0][dim
], d
);
2820 I
= DomainImage(D
, T
, 0);
2821 H
= DomainConvex(I
, 0);
2831 static void replace_by_affine(evalue
*e
, Value offset
)
2838 value_init(inc
.x
.n
);
2839 value_set_si(inc
.d
, 1);
2840 value_oppose(inc
.x
.n
, offset
);
2841 eadd(&inc
, &p
->arr
[0]);
2842 reorder_terms_about(p
, &p
->arr
[0]); /* frees arr[0] */
2846 free_evalue_refs(&inc
);
2849 int evalue_range_reduction_in_domain(evalue
*e
, Polyhedron
*D
)
2858 if (value_notzero_p(e
->d
))
2863 if (p
->type
== relation
) {
2870 fractional_minimal_coefficients(p
->arr
[0].x
.p
);
2871 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
);
2872 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2873 equal
= value_eq(min
, max
);
2874 mpz_cdiv_q(min
, min
, d
);
2875 mpz_fdiv_q(max
, max
, d
);
2877 if (bounded
&& value_gt(min
, max
)) {
2883 evalue_set_si(e
, 0, 1);
2886 free_evalue_refs(&(p
->arr
[1]));
2887 free_evalue_refs(&(p
->arr
[0]));
2893 return r
? r
: evalue_range_reduction_in_domain(e
, D
);
2894 } else if (bounded
&& equal
) {
2897 free_evalue_refs(&(p
->arr
[2]));
2900 free_evalue_refs(&(p
->arr
[0]));
2906 return evalue_range_reduction_in_domain(e
, D
);
2907 } else if (bounded
&& value_eq(min
, max
)) {
2908 /* zero for a single value */
2910 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2911 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2912 value_multiply(min
, min
, d
);
2913 value_subtract(M
->p
[0][D
->Dimension
+1],
2914 M
->p
[0][D
->Dimension
+1], min
);
2915 E
= DomainAddConstraints(D
, M
, 0);
2921 r
= evalue_range_reduction_in_domain(&p
->arr
[1], E
);
2923 r
|= evalue_range_reduction_in_domain(&p
->arr
[2], D
);
2925 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2933 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2936 i
= p
->type
== relation
? 1 :
2937 p
->type
== fractional
? 1 : 0;
2938 for (; i
<p
->size
; i
++)
2939 r
|= evalue_range_reduction_in_domain(&p
->arr
[i
], D
);
2941 if (p
->type
!= fractional
) {
2942 if (r
&& p
->type
== polynomial
) {
2945 value_set_si(f
.d
, 0);
2946 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2947 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2948 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2949 reorder_terms_about(p
, &f
);
2960 fractional_minimal_coefficients(p
);
2961 I
= polynomial_projection(p
, D
, &d
, NULL
);
2962 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2963 mpz_fdiv_q(min
, min
, d
);
2964 mpz_fdiv_q(max
, max
, d
);
2965 value_subtract(d
, max
, min
);
2967 if (bounded
&& value_eq(min
, max
)) {
2968 replace_by_affine(e
, min
);
2970 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
2971 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
2972 * See pages 199-200 of PhD thesis.
2980 value_set_si(rem
.d
, 0);
2981 rem
.x
.p
= new_enode(fractional
, 3, -1);
2982 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
2983 value_clear(rem
.x
.p
->arr
[1].d
);
2984 value_clear(rem
.x
.p
->arr
[2].d
);
2985 rem
.x
.p
->arr
[1] = p
->arr
[1];
2986 rem
.x
.p
->arr
[2] = p
->arr
[2];
2987 for (i
= 3; i
< p
->size
; ++i
)
2988 p
->arr
[i
-2] = p
->arr
[i
];
2992 value_init(inc
.x
.n
);
2993 value_set_si(inc
.d
, 1);
2994 value_oppose(inc
.x
.n
, min
);
2997 evalue_copy(&t
, &p
->arr
[0]);
3001 value_set_si(f
.d
, 0);
3002 f
.x
.p
= new_enode(fractional
, 3, -1);
3003 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3004 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3005 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
3007 value_init(factor
.d
);
3008 evalue_set_si(&factor
, -1, 1);
3014 value_clear(f
.x
.p
->arr
[1].x
.n
);
3015 value_clear(f
.x
.p
->arr
[2].x
.n
);
3016 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3017 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3021 reorder_terms(&rem
);
3028 free_evalue_refs(&inc
);
3029 free_evalue_refs(&t
);
3030 free_evalue_refs(&f
);
3031 free_evalue_refs(&factor
);
3032 free_evalue_refs(&rem
);
3034 evalue_range_reduction_in_domain(e
, D
);
3038 _reduce_evalue(&p
->arr
[0], 0, 1);
3050 void evalue_range_reduction(evalue
*e
)
3053 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3056 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3057 if (evalue_range_reduction_in_domain(&e
->x
.p
->arr
[2*i
+1],
3058 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
3059 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3061 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
3062 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
3063 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3064 value_clear(e
->x
.p
->arr
[2*i
].d
);
3066 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
3067 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
3075 Enumeration
* partition2enumeration(evalue
*EP
)
3078 Enumeration
*en
, *res
= NULL
;
3080 if (EVALUE_IS_ZERO(*EP
)) {
3085 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
3086 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
3087 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
3090 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
3091 value_clear(EP
->x
.p
->arr
[2*i
].d
);
3092 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
3100 int evalue_frac2floor_in_domain3(evalue
*e
, Polyhedron
*D
, int shift
)
3109 if (value_notzero_p(e
->d
))
3114 i
= p
->type
== relation
? 1 :
3115 p
->type
== fractional
? 1 : 0;
3116 for (; i
<p
->size
; i
++)
3117 r
|= evalue_frac2floor_in_domain3(&p
->arr
[i
], D
, shift
);
3119 if (p
->type
!= fractional
) {
3120 if (r
&& p
->type
== polynomial
) {
3123 value_set_si(f
.d
, 0);
3124 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3125 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3126 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3127 reorder_terms_about(p
, &f
);
3137 I
= polynomial_projection(p
, D
, &d
, NULL
);
3140 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3143 assert(I
->NbEq
== 0); /* Should have been reduced */
3146 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3147 if (value_pos_p(I
->Constraint
[i
][1]))
3150 if (i
< I
->NbConstraints
) {
3152 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3153 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3154 if (value_neg_p(min
)) {
3156 mpz_fdiv_q(min
, min
, d
);
3157 value_init(offset
.d
);
3158 value_set_si(offset
.d
, 1);
3159 value_init(offset
.x
.n
);
3160 value_oppose(offset
.x
.n
, min
);
3161 eadd(&offset
, &p
->arr
[0]);
3162 free_evalue_refs(&offset
);
3172 value_set_si(fl
.d
, 0);
3173 fl
.x
.p
= new_enode(flooring
, 3, -1);
3174 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
3175 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
3176 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
3178 eadd(&fl
, &p
->arr
[0]);
3179 reorder_terms_about(p
, &p
->arr
[0]);
3183 free_evalue_refs(&fl
);
3188 int evalue_frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
3190 return evalue_frac2floor_in_domain3(e
, D
, 1);
3193 void evalue_frac2floor2(evalue
*e
, int shift
)
3196 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
3198 if (evalue_frac2floor_in_domain3(e
, NULL
, 0))
3204 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3205 if (evalue_frac2floor_in_domain3(&e
->x
.p
->arr
[2*i
+1],
3206 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), shift
))
3207 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3210 void evalue_frac2floor(evalue
*e
)
3212 evalue_frac2floor2(e
, 1);
3215 /* Add a new variable with lower bound 1 and upper bound specified
3216 * by row. If negative is true, then the new variable has upper
3217 * bound -1 and lower bound specified by row.
3219 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
3220 Vector
*row
, int negative
)
3224 int nparam
= D
->Dimension
- nvar
;
3227 nr
= D
->NbConstraints
+ 2;
3228 nc
= D
->Dimension
+ 2 + 1;
3229 C
= Matrix_Alloc(nr
, nc
);
3230 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
3231 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
3232 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3233 D
->Dimension
+ 1 - nvar
);
3238 nc
= C
->NbColumns
+ 1;
3239 C
= Matrix_Alloc(nr
, nc
);
3240 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
3241 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
3242 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3243 oldC
->NbColumns
- 1 - nvar
);
3246 value_set_si(C
->p
[nr
-2][0], 1);
3248 value_set_si(C
->p
[nr
-2][1 + nvar
], -1);
3250 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
3251 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
3253 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
3254 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
3260 static void floor2frac_r(evalue
*e
, int nvar
)
3267 if (value_notzero_p(e
->d
))
3272 assert(p
->type
== flooring
);
3273 for (i
= 1; i
< p
->size
; i
++)
3274 floor2frac_r(&p
->arr
[i
], nvar
);
3276 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
3277 assert(pp
->x
.p
->type
== polynomial
);
3278 pp
->x
.p
->pos
-= nvar
;
3282 value_set_si(f
.d
, 0);
3283 f
.x
.p
= new_enode(fractional
, 3, -1);
3284 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3285 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3286 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3288 eadd(&f
, &p
->arr
[0]);
3289 reorder_terms_about(p
, &p
->arr
[0]);
3293 free_evalue_refs(&f
);
3296 /* Convert flooring back to fractional and shift position
3297 * of the parameters by nvar
3299 static void floor2frac(evalue
*e
, int nvar
)
3301 floor2frac_r(e
, nvar
);
3305 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
3308 int nparam
= D
->Dimension
- nvar
;
3312 D
= Constraints2Polyhedron(C
, 0);
3316 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
3318 /* Double check that D was not unbounded. */
3319 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
3327 static evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
3328 int *signs
, Matrix
*C
, unsigned MaxRays
)
3334 evalue
*factor
= NULL
;
3338 if (EVALUE_IS_ZERO(*e
))
3342 Polyhedron
*DD
= Disjoint_Domain(D
, 0, MaxRays
);
3349 res
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3352 for (Q
= DD
; Q
; Q
= DD
) {
3358 t
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3371 if (value_notzero_p(e
->d
)) {
3374 t
= esum_over_domain_cst(nvar
, D
, C
);
3376 if (!EVALUE_IS_ONE(*e
))
3382 switch (e
->x
.p
->type
) {
3384 evalue
*pp
= &e
->x
.p
->arr
[0];
3386 if (pp
->x
.p
->pos
> nvar
) {
3387 /* remainder is independent of the summated vars */
3393 floor2frac(&f
, nvar
);
3395 t
= esum_over_domain_cst(nvar
, D
, C
);
3399 free_evalue_refs(&f
);
3404 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3405 poly_denom(pp
, &row
->p
[1 + nvar
]);
3406 value_set_si(row
->p
[0], 1);
3407 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3408 pp
= &pp
->x
.p
->arr
[0]) {
3410 assert(pp
->x
.p
->type
== polynomial
);
3412 if (pos
>= 1 + nvar
)
3414 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3415 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3416 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3418 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3419 value_division(row
->p
[1 + D
->Dimension
+ 1],
3420 row
->p
[1 + D
->Dimension
+ 1],
3422 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3423 row
->p
[1 + D
->Dimension
+ 1],
3425 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3429 int pos
= e
->x
.p
->pos
;
3432 factor
= ALLOC(evalue
);
3433 value_init(factor
->d
);
3434 value_set_si(factor
->d
, 0);
3435 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3436 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3437 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3441 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3442 negative
= signs
[pos
-1] < 0;
3443 value_set_si(row
->p
[0], 1);
3445 value_set_si(row
->p
[pos
], -1);
3446 value_set_si(row
->p
[1 + nvar
], 1);
3448 value_set_si(row
->p
[pos
], 1);
3449 value_set_si(row
->p
[1 + nvar
], -1);
3457 offset
= type_offset(e
->x
.p
);
3459 res
= esum_over_domain(&e
->x
.p
->arr
[offset
], nvar
, D
, signs
, C
, MaxRays
);
3463 evalue_copy(&cum
, factor
);
3467 for (i
= 1; offset
+i
< e
->x
.p
->size
; ++i
) {
3471 C
= esum_add_constraint(nvar
, D
, C
, row
, negative
);
3477 Vector_Print(stderr, P_VALUE_FMT, row);
3479 Matrix_Print(stderr, P_VALUE_FMT, C);
3481 t
= esum_over_domain(&e
->x
.p
->arr
[offset
+i
], nvar
, D
, signs
, C
, MaxRays
);
3486 if (negative
&& (i
% 2))
3496 if (factor
&& offset
+i
+1 < e
->x
.p
->size
)
3503 free_evalue_refs(&cum
);
3504 evalue_free(factor
);
3515 static void domain_signs(Polyhedron
*D
, int *signs
)
3519 POL_ENSURE_VERTICES(D
);
3520 for (j
= 0; j
< D
->Dimension
; ++j
) {
3522 for (k
= 0; k
< D
->NbRays
; ++k
) {
3523 signs
[j
] = value_sign(D
->Ray
[k
][1+j
]);
3530 static void shift_floor_in_domain(evalue
*e
, Polyhedron
*D
)
3537 if (value_notzero_p(e
->d
))
3542 for (i
= type_offset(p
); i
< p
->size
; ++i
)
3543 shift_floor_in_domain(&p
->arr
[i
], D
);
3545 if (p
->type
!= flooring
)
3551 I
= polynomial_projection(p
, D
, &d
, NULL
);
3552 assert(I
->NbEq
== 0); /* Should have been reduced */
3554 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3555 if (value_pos_p(I
->Constraint
[i
][1]))
3557 assert(i
< I
->NbConstraints
);
3558 if (i
< I
->NbConstraints
) {
3559 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3560 mpz_fdiv_q(m
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3561 if (value_neg_p(m
)) {
3562 /* replace [e] by [e-m]+m such that e-m >= 0 */
3567 value_set_si(f
.d
, 1);
3568 value_oppose(f
.x
.n
, m
);
3569 eadd(&f
, &p
->arr
[0]);
3572 value_set_si(f
.d
, 0);
3573 f
.x
.p
= new_enode(flooring
, 3, -1);
3574 value_clear(f
.x
.p
->arr
[0].d
);
3575 f
.x
.p
->arr
[0] = p
->arr
[0];
3576 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
3577 value_set_si(f
.x
.p
->arr
[1].d
, 1);
3578 value_init(f
.x
.p
->arr
[1].x
.n
);
3579 value_assign(f
.x
.p
->arr
[1].x
.n
, m
);
3580 reorder_terms_about(p
, &f
);
3591 /* Make arguments of all floors non-negative */
3592 static void shift_floor_arguments(evalue
*e
)
3596 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3599 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3600 shift_floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
3601 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3604 evalue
*evalue_sum(evalue
*e
, int nvar
, unsigned MaxRays
)
3608 evalue
*res
= ALLOC(evalue
);
3612 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3613 evalue_copy(res
, e
);
3617 evalue_split_domains_into_orthants(e
, MaxRays
);
3619 evalue_frac2floor2(e
, 0);
3620 evalue_set_si(res
, 0, 1);
3622 assert(value_zero_p(e
->d
));
3623 assert(e
->x
.p
->type
== partition
);
3624 shift_floor_arguments(e
);
3626 assert(e
->x
.p
->size
>= 2);
3627 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3629 signs
= alloca(sizeof(int) * dim
);
3631 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3633 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3634 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3635 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
, 0,
3646 evalue
*esum(evalue
*e
, int nvar
)
3648 return evalue_sum(e
, nvar
, 0);
3651 /* Initial silly implementation */
3652 void eor(evalue
*e1
, evalue
*res
)
3658 evalue_set_si(&mone
, -1, 1);
3660 evalue_copy(&E
, res
);
3666 free_evalue_refs(&E
);
3667 free_evalue_refs(&mone
);
3670 /* computes denominator of polynomial evalue
3671 * d should point to a value initialized to 1
3673 void evalue_denom(const evalue
*e
, Value
*d
)
3677 if (value_notzero_p(e
->d
)) {
3678 value_lcm(*d
, *d
, e
->d
);
3681 assert(e
->x
.p
->type
== polynomial
||
3682 e
->x
.p
->type
== fractional
||
3683 e
->x
.p
->type
== flooring
);
3684 offset
= type_offset(e
->x
.p
);
3685 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3686 evalue_denom(&e
->x
.p
->arr
[i
], d
);
3689 /* Divides the evalue e by the integer n */
3690 void evalue_div(evalue
*e
, Value n
)
3694 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3697 if (value_notzero_p(e
->d
)) {
3700 value_multiply(e
->d
, e
->d
, n
);
3701 value_gcd(gc
, e
->x
.n
, e
->d
);
3702 if (value_notone_p(gc
)) {
3703 value_division(e
->d
, e
->d
, gc
);
3704 value_division(e
->x
.n
, e
->x
.n
, gc
);
3709 if (e
->x
.p
->type
== partition
) {
3710 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3711 evalue_div(&e
->x
.p
->arr
[2*i
+1], n
);
3714 offset
= type_offset(e
->x
.p
);
3715 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3716 evalue_div(&e
->x
.p
->arr
[i
], n
);
3719 /* Multiplies the evalue e by the integer n */
3720 void evalue_mul(evalue
*e
, Value n
)
3724 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3727 if (value_notzero_p(e
->d
)) {
3730 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3731 value_gcd(gc
, e
->x
.n
, e
->d
);
3732 if (value_notone_p(gc
)) {
3733 value_division(e
->d
, e
->d
, gc
);
3734 value_division(e
->x
.n
, e
->x
.n
, gc
);
3739 if (e
->x
.p
->type
== partition
) {
3740 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3741 evalue_mul(&e
->x
.p
->arr
[2*i
+1], n
);
3744 offset
= type_offset(e
->x
.p
);
3745 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3746 evalue_mul(&e
->x
.p
->arr
[i
], n
);
3749 /* Multiplies the evalue e by the n/d */
3750 void evalue_mul_div(evalue
*e
, Value n
, Value d
)
3754 if ((value_one_p(n
) && value_one_p(d
)) || EVALUE_IS_ZERO(*e
))
3757 if (value_notzero_p(e
->d
)) {
3760 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3761 value_multiply(e
->d
, e
->d
, d
);
3762 value_gcd(gc
, e
->x
.n
, e
->d
);
3763 if (value_notone_p(gc
)) {
3764 value_division(e
->d
, e
->d
, gc
);
3765 value_division(e
->x
.n
, e
->x
.n
, gc
);
3770 if (e
->x
.p
->type
== partition
) {
3771 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3772 evalue_mul_div(&e
->x
.p
->arr
[2*i
+1], n
, d
);
3775 offset
= type_offset(e
->x
.p
);
3776 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3777 evalue_mul_div(&e
->x
.p
->arr
[i
], n
, d
);
3780 void evalue_negate(evalue
*e
)
3784 if (value_notzero_p(e
->d
)) {
3785 value_oppose(e
->x
.n
, e
->x
.n
);
3788 if (e
->x
.p
->type
== partition
) {
3789 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3790 evalue_negate(&e
->x
.p
->arr
[2*i
+1]);
3793 offset
= type_offset(e
->x
.p
);
3794 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3795 evalue_negate(&e
->x
.p
->arr
[i
]);
3798 void evalue_add_constant(evalue
*e
, const Value cst
)
3802 if (value_zero_p(e
->d
)) {
3803 if (e
->x
.p
->type
== partition
) {
3804 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3805 evalue_add_constant(&e
->x
.p
->arr
[2*i
+1], cst
);
3808 if (e
->x
.p
->type
== relation
) {
3809 for (i
= 1; i
< e
->x
.p
->size
; ++i
)
3810 evalue_add_constant(&e
->x
.p
->arr
[i
], cst
);
3814 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
3815 } while (value_zero_p(e
->d
));
3817 value_addmul(e
->x
.n
, cst
, e
->d
);
3820 static void evalue_frac2polynomial_r(evalue
*e
, int *signs
, int sign
, int in_frac
)
3825 int sign_odd
= sign
;
3827 if (value_notzero_p(e
->d
)) {
3828 if (in_frac
&& sign
* value_sign(e
->x
.n
) < 0) {
3829 value_set_si(e
->x
.n
, 0);
3830 value_set_si(e
->d
, 1);
3835 if (e
->x
.p
->type
== relation
) {
3836 for (i
= e
->x
.p
->size
-1; i
>= 1; --i
)
3837 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
, sign
, in_frac
);
3841 if (e
->x
.p
->type
== polynomial
)
3842 sign_odd
*= signs
[e
->x
.p
->pos
-1];
3843 offset
= type_offset(e
->x
.p
);
3844 evalue_frac2polynomial_r(&e
->x
.p
->arr
[offset
], signs
, sign
, in_frac
);
3845 in_frac
|= e
->x
.p
->type
== fractional
;
3846 for (i
= e
->x
.p
->size
-1; i
> offset
; --i
)
3847 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
,
3848 (i
- offset
) % 2 ? sign_odd
: sign
, in_frac
);
3850 if (e
->x
.p
->type
!= fractional
)
3853 /* replace { a/m } by (m-1)/m if sign != 0
3854 * and by (m-1)/(2m) if sign == 0
3858 evalue_denom(&e
->x
.p
->arr
[0], &d
);
3859 free_evalue_refs(&e
->x
.p
->arr
[0]);
3860 value_init(e
->x
.p
->arr
[0].d
);
3861 value_init(e
->x
.p
->arr
[0].x
.n
);
3863 value_addto(e
->x
.p
->arr
[0].d
, d
, d
);
3865 value_assign(e
->x
.p
->arr
[0].d
, d
);
3866 value_decrement(e
->x
.p
->arr
[0].x
.n
, d
);
3870 reorder_terms_about(p
, &p
->arr
[0]);
3876 /* Approximate the evalue in fractional representation by a polynomial.
3877 * If sign > 0, the result is an upper bound;
3878 * if sign < 0, the result is a lower bound;
3879 * if sign = 0, the result is an intermediate approximation.
3881 void evalue_frac2polynomial(evalue
*e
, int sign
, unsigned MaxRays
)
3886 if (value_notzero_p(e
->d
))
3888 assert(e
->x
.p
->type
== partition
);
3889 /* make sure all variables in the domains have a fixed sign */
3891 evalue_split_domains_into_orthants(e
, MaxRays
);
3892 if (EVALUE_IS_ZERO(*e
))
3896 assert(e
->x
.p
->size
>= 2);
3897 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3899 signs
= alloca(sizeof(int) * dim
);
3902 for (i
= 0; i
< dim
; ++i
)
3904 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3906 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3907 evalue_frac2polynomial_r(&e
->x
.p
->arr
[2*i
+1], signs
, sign
, 0);
3911 /* Split the domains of e (which is assumed to be a partition)
3912 * such that each resulting domain lies entirely in one orthant.
3914 void evalue_split_domains_into_orthants(evalue
*e
, unsigned MaxRays
)
3917 assert(value_zero_p(e
->d
));
3918 assert(e
->x
.p
->type
== partition
);
3919 assert(e
->x
.p
->size
>= 2);
3920 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3922 for (i
= 0; i
< dim
; ++i
) {
3925 C
= Matrix_Alloc(1, 1 + dim
+ 1);
3926 value_set_si(C
->p
[0][0], 1);
3927 value_init(split
.d
);
3928 value_set_si(split
.d
, 0);
3929 split
.x
.p
= new_enode(partition
, 4, dim
);
3930 value_set_si(C
->p
[0][1+i
], 1);
3931 C2
= Matrix_Copy(C
);
3932 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[0], Constraints2Polyhedron(C2
, MaxRays
));
3934 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
3935 value_set_si(C
->p
[0][1+i
], -1);
3936 value_set_si(C
->p
[0][1+dim
], -1);
3937 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[2], Constraints2Polyhedron(C
, MaxRays
));
3938 evalue_set_si(&split
.x
.p
->arr
[3], 1, 1);
3940 free_evalue_refs(&split
);
3945 static evalue
*find_fractional_with_max_periods(evalue
*e
, Polyhedron
*D
,
3948 Value
*min
, Value
*max
)
3955 if (value_notzero_p(e
->d
))
3958 if (e
->x
.p
->type
== fractional
) {
3963 I
= polynomial_projection(e
->x
.p
, D
, &d
, &T
);
3964 bounded
= line_minmax(I
, min
, max
); /* frees I */
3968 value_set_si(mp
, max_periods
);
3969 mpz_fdiv_q(*min
, *min
, d
);
3970 mpz_fdiv_q(*max
, *max
, d
);
3971 value_assign(T
->p
[1][D
->Dimension
], d
);
3972 value_subtract(d
, *max
, *min
);
3973 if (value_ge(d
, mp
))
3976 f
= evalue_dup(&e
->x
.p
->arr
[0]);
3987 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
3988 if ((f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[i
], D
, max_periods
,
3995 static void replace_fract_by_affine(evalue
*e
, evalue
*f
, Value val
)
3999 if (value_notzero_p(e
->d
))
4002 offset
= type_offset(e
->x
.p
);
4003 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
4004 replace_fract_by_affine(&e
->x
.p
->arr
[i
], f
, val
);
4006 if (e
->x
.p
->type
!= fractional
)
4009 if (!eequal(&e
->x
.p
->arr
[0], f
))
4012 replace_by_affine(e
, val
);
4015 /* Look for fractional parts that can be removed by splitting the corresponding
4016 * domain into at most max_periods parts.
4017 * We use a very simply strategy that looks for the first fractional part
4018 * that satisfies the condition, performs the split and then continues
4019 * looking for other fractional parts in the split domains until no
4020 * such fractional part can be found anymore.
4022 void evalue_split_periods(evalue
*e
, int max_periods
, unsigned int MaxRays
)
4029 if (EVALUE_IS_ZERO(*e
))
4031 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
4033 "WARNING: evalue_split_periods called on incorrect evalue type\n");
4041 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
4046 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
4048 f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[2*i
+1], D
, max_periods
,
4053 M
= Matrix_Alloc(2, 2+D
->Dimension
);
4055 value_subtract(d
, max
, min
);
4056 n
= VALUE_TO_INT(d
)+1;
4058 value_set_si(M
->p
[0][0], 1);
4059 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
4060 value_multiply(d
, max
, T
->p
[1][D
->Dimension
]);
4061 value_subtract(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
], d
);
4062 value_set_si(d
, -1);
4063 value_set_si(M
->p
[1][0], 1);
4064 Vector_Scale(T
->p
[0], M
->p
[1]+1, d
, D
->Dimension
+1);
4065 value_addmul(M
->p
[1][1+D
->Dimension
], max
, T
->p
[1][D
->Dimension
]);
4066 value_addto(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4067 T
->p
[1][D
->Dimension
]);
4068 value_decrement(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
]);
4070 p
= new_enode(partition
, e
->x
.p
->size
+ (n
-1)*2, e
->x
.p
->pos
);
4071 for (j
= 0; j
< 2*i
; ++j
) {
4072 value_clear(p
->arr
[j
].d
);
4073 p
->arr
[j
] = e
->x
.p
->arr
[j
];
4075 for (j
= 2*i
+2; j
< e
->x
.p
->size
; ++j
) {
4076 value_clear(p
->arr
[j
+2*(n
-1)].d
);
4077 p
->arr
[j
+2*(n
-1)] = e
->x
.p
->arr
[j
];
4079 for (j
= n
-1; j
>= 0; --j
) {
4081 value_clear(p
->arr
[2*i
+1].d
);
4082 p
->arr
[2*i
+1] = e
->x
.p
->arr
[2*i
+1];
4084 evalue_copy(&p
->arr
[2*(i
+j
)+1], &e
->x
.p
->arr
[2*i
+1]);
4086 value_subtract(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4087 T
->p
[1][D
->Dimension
]);
4088 value_addto(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
],
4089 T
->p
[1][D
->Dimension
]);
4091 replace_fract_by_affine(&p
->arr
[2*(i
+j
)+1], f
, max
);
4092 E
= DomainAddConstraints(D
, M
, MaxRays
);
4093 EVALUE_SET_DOMAIN(p
->arr
[2*(i
+j
)], E
);
4094 if (evalue_range_reduction_in_domain(&p
->arr
[2*(i
+j
)+1], E
))
4095 reduce_evalue(&p
->arr
[2*(i
+j
)+1]);
4096 value_decrement(max
, max
);
4098 value_clear(e
->x
.p
->arr
[2*i
].d
);
4113 void evalue_extract_affine(const evalue
*e
, Value
*coeff
, Value
*cst
, Value
*d
)
4115 value_set_si(*d
, 1);
4117 for ( ; value_zero_p(e
->d
); e
= &e
->x
.p
->arr
[0]) {
4119 assert(e
->x
.p
->type
== polynomial
);
4120 assert(e
->x
.p
->size
== 2);
4121 c
= &e
->x
.p
->arr
[1];
4122 value_multiply(coeff
[e
->x
.p
->pos
-1], *d
, c
->x
.n
);
4123 value_division(coeff
[e
->x
.p
->pos
-1], coeff
[e
->x
.p
->pos
-1], c
->d
);
4125 value_multiply(*cst
, *d
, e
->x
.n
);
4126 value_division(*cst
, *cst
, e
->d
);
4129 /* returns an evalue that corresponds to
4133 static evalue
*term(int param
, Value c
, Value den
)
4135 evalue
*EP
= ALLOC(evalue
);
4137 value_set_si(EP
->d
,0);
4138 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
4139 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
4140 value_init(EP
->x
.p
->arr
[1].x
.n
);
4141 value_assign(EP
->x
.p
->arr
[1].d
, den
);
4142 value_assign(EP
->x
.p
->arr
[1].x
.n
, c
);
4146 evalue
*affine2evalue(Value
*coeff
, Value denom
, int nvar
)
4149 evalue
*E
= ALLOC(evalue
);
4151 evalue_set(E
, coeff
[nvar
], denom
);
4152 for (i
= 0; i
< nvar
; ++i
) {
4154 if (value_zero_p(coeff
[i
]))
4156 t
= term(i
, coeff
[i
], denom
);
4163 void evalue_substitute(evalue
*e
, evalue
**subs
)
4169 if (value_notzero_p(e
->d
))
4173 assert(p
->type
!= partition
);
4175 for (i
= 0; i
< p
->size
; ++i
)
4176 evalue_substitute(&p
->arr
[i
], subs
);
4178 if (p
->type
== relation
) {
4179 /* For relation a ? b : c, compute (a' ? 1) * b' + (a' ? 0 : 1) * c' */
4183 value_set_si(v
->d
, 0);
4184 v
->x
.p
= new_enode(relation
, 3, 0);
4185 evalue_copy(&v
->x
.p
->arr
[0], &p
->arr
[0]);
4186 evalue_set_si(&v
->x
.p
->arr
[1], 0, 1);
4187 evalue_set_si(&v
->x
.p
->arr
[2], 1, 1);
4188 emul(v
, &p
->arr
[2]);
4193 value_set_si(v
->d
, 0);
4194 v
->x
.p
= new_enode(relation
, 2, 0);
4195 value_clear(v
->x
.p
->arr
[0].d
);
4196 v
->x
.p
->arr
[0] = p
->arr
[0];
4197 evalue_set_si(&v
->x
.p
->arr
[1], 1, 1);
4198 emul(v
, &p
->arr
[1]);
4201 eadd(&p
->arr
[2], &p
->arr
[1]);
4202 free_evalue_refs(&p
->arr
[2]);
4210 if (p
->type
== polynomial
)
4215 value_set_si(v
->d
, 0);
4216 v
->x
.p
= new_enode(p
->type
, 3, -1);
4217 value_clear(v
->x
.p
->arr
[0].d
);
4218 v
->x
.p
->arr
[0] = p
->arr
[0];
4219 evalue_set_si(&v
->x
.p
->arr
[1], 0, 1);
4220 evalue_set_si(&v
->x
.p
->arr
[2], 1, 1);
4223 offset
= type_offset(p
);
4225 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
4226 emul(v
, &p
->arr
[i
]);
4227 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
4228 free_evalue_refs(&(p
->arr
[i
]));
4231 if (p
->type
!= polynomial
)
4235 *e
= p
->arr
[offset
];
4239 /* evalue e is given in terms of "new" parameter; CP maps the new
4240 * parameters back to the old parameters.
4241 * Transforms e such that it refers back to the old parameters and
4242 * adds appropriate constraints to the domain.
4243 * In particular, if CP maps the new parameters onto an affine
4244 * subspace of the old parameters, then the corresponding equalities
4245 * are added to the domain.
4246 * Also, if any of the new parameters was a rational combination
4247 * of the old parameters $p' = (<a, p> + c)/m$, then modulo
4248 * constraints ${<a, p> + c)/m} = 0$ are added to ensure
4249 * the new evalue remains non-zero only for integer parameters
4250 * of the new parameters (which have been removed by the substitution).
4252 void evalue_backsubstitute(evalue
*e
, Matrix
*CP
, unsigned MaxRays
)
4259 unsigned nparam
= CP
->NbColumns
-1;
4263 if (EVALUE_IS_ZERO(*e
))
4266 assert(value_zero_p(e
->d
));
4268 assert(p
->type
== partition
);
4270 inv
= left_inverse(CP
, &eq
);
4271 subs
= ALLOCN(evalue
*, nparam
);
4272 for (i
= 0; i
< nparam
; ++i
)
4273 subs
[i
] = affine2evalue(inv
->p
[i
], inv
->p
[nparam
][inv
->NbColumns
-1],
4276 CEq
= Constraints2Polyhedron(eq
, MaxRays
);
4277 addeliminatedparams_partition(p
, inv
, CEq
, inv
->NbColumns
-1, MaxRays
);
4278 Polyhedron_Free(CEq
);
4280 for (i
= 0; i
< p
->size
/2; ++i
)
4281 evalue_substitute(&p
->arr
[2*i
+1], subs
);
4283 for (i
= 0; i
< nparam
; ++i
)
4284 evalue_free(subs
[i
]);
4288 for (i
= 0; i
< inv
->NbRows
-1; ++i
) {
4289 Vector_Gcd(inv
->p
[i
], inv
->NbColumns
, &gcd
);
4290 value_gcd(gcd
, gcd
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]);
4291 if (value_eq(gcd
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]))
4293 Vector_AntiScale(inv
->p
[i
], inv
->p
[i
], gcd
, inv
->NbColumns
);
4294 value_divexact(gcd
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1], gcd
);
4296 for (j
= 0; j
< p
->size
/2; ++j
) {
4297 evalue
*arg
= affine2evalue(inv
->p
[i
], gcd
, inv
->NbColumns
-1);
4302 value_set_si(rel
.d
, 0);
4303 rel
.x
.p
= new_enode(relation
, 2, 0);
4304 value_clear(rel
.x
.p
->arr
[1].d
);
4305 rel
.x
.p
->arr
[1] = p
->arr
[2*j
+1];
4306 ev
= &rel
.x
.p
->arr
[0];
4307 value_set_si(ev
->d
, 0);
4308 ev
->x
.p
= new_enode(fractional
, 3, -1);
4309 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
4310 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
4311 value_clear(ev
->x
.p
->arr
[0].d
);
4312 ev
->x
.p
->arr
[0] = *arg
;
4315 p
->arr
[2*j
+1] = rel
;
4326 * \sum_{i=0}^n c_i/d X^i
4328 * where d is the last element in the vector c.
4330 evalue
*evalue_polynomial(Vector
*c
, const evalue
* X
)
4332 unsigned dim
= c
->Size
-2;
4334 evalue
*EP
= ALLOC(evalue
);
4339 if (EVALUE_IS_ZERO(*X
) || dim
== 0) {
4340 evalue_set(EP
, c
->p
[0], c
->p
[dim
+1]);
4341 reduce_constant(EP
);
4345 evalue_set(EP
, c
->p
[dim
], c
->p
[dim
+1]);
4348 evalue_set(&EC
, c
->p
[dim
], c
->p
[dim
+1]);
4350 for (i
= dim
-1; i
>= 0; --i
) {
4352 value_assign(EC
.x
.n
, c
->p
[i
]);
4355 free_evalue_refs(&EC
);
4359 /* Create an evalue from an array of pairs of domains and evalues. */
4360 evalue
*evalue_from_section_array(struct evalue_section
*s
, int n
)
4365 res
= ALLOC(evalue
);
4369 evalue_set_si(res
, 0, 1);
4371 value_set_si(res
->d
, 0);
4372 res
->x
.p
= new_enode(partition
, 2*n
, s
[0].D
->Dimension
);
4373 for (i
= 0; i
< n
; ++i
) {
4374 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
], s
[i
].D
);
4375 value_clear(res
->x
.p
->arr
[2*i
+1].d
);
4376 res
->x
.p
->arr
[2*i
+1] = *s
[i
].E
;
4383 /* shift variables in polynomial n up */
4384 void evalue_shift_variables(evalue
*e
, int n
)
4387 if (value_notzero_p(e
->d
))
4389 assert(e
->x
.p
->type
== polynomial
|| e
->x
.p
->type
== fractional
);
4390 if (e
->x
.p
->type
== polynomial
) {
4391 assert(e
->x
.p
->pos
+ n
>= 1);
4394 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
4395 evalue_shift_variables(&e
->x
.p
->arr
[i
], n
);