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 /***********************************************************************/
15 #include <barvinok/evalue.h>
16 #include <barvinok/barvinok.h>
17 #include <barvinok/util.h>
19 #ifndef value_pmodulus
20 #define value_pmodulus(ref,val1,val2) (mpz_fdiv_r((ref),(val1),(val2)))
23 #define ALLOC(type) (type*)malloc(sizeof(type))
24 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
27 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
29 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
32 void evalue_set_si(evalue
*ev
, int n
, int d
) {
33 value_set_si(ev
->d
, d
);
35 value_set_si(ev
->x
.n
, n
);
38 void evalue_set(evalue
*ev
, Value n
, Value d
) {
39 value_assign(ev
->d
, d
);
41 value_assign(ev
->x
.n
, n
);
46 evalue
*EP
= ALLOC(evalue
);
48 evalue_set_si(EP
, 0, 1);
52 void aep_evalue(evalue
*e
, int *ref
) {
57 if (value_notzero_p(e
->d
))
58 return; /* a rational number, its already reduced */
60 return; /* hum... an overflow probably occured */
62 /* First check the components of p */
63 for (i
=0;i
<p
->size
;i
++)
64 aep_evalue(&p
->arr
[i
],ref
);
71 p
->pos
= ref
[p
->pos
-1]+1;
77 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
83 if (value_notzero_p(e
->d
))
84 return; /* a rational number, its already reduced */
86 return; /* hum... an overflow probably occured */
89 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
90 for(i
=0;i
<CT
->NbRows
-1;i
++)
91 for(j
=0;j
<CT
->NbColumns
;j
++)
92 if(value_notzero_p(CT
->p
[i
][j
])) {
97 /* Transform the references in e, using ref */
101 } /* addeliminatedparams_evalue */
103 static void addeliminatedparams_partition(enode
*p
, Matrix
*CT
, Polyhedron
*CEq
,
104 unsigned nparam
, unsigned MaxRays
)
107 assert(p
->type
== partition
);
110 for (i
= 0; i
< p
->size
/2; i
++) {
111 Polyhedron
*D
= EVALUE_DOMAIN(p
->arr
[2*i
]);
112 Polyhedron
*T
= DomainPreimage(D
, CT
, MaxRays
);
116 T
= DomainIntersection(D
, CEq
, MaxRays
);
119 EVALUE_SET_DOMAIN(p
->arr
[2*i
], T
);
123 void addeliminatedparams_enum(evalue
*e
, Matrix
*CT
, Polyhedron
*CEq
,
124 unsigned MaxRays
, unsigned nparam
)
129 if (CT
->NbRows
== CT
->NbColumns
)
132 if (EVALUE_IS_ZERO(*e
))
135 if (value_notzero_p(e
->d
)) {
138 value_set_si(res
.d
, 0);
139 res
.x
.p
= new_enode(partition
, 2, nparam
);
140 EVALUE_SET_DOMAIN(res
.x
.p
->arr
[0],
141 DomainConstraintSimplify(Polyhedron_Copy(CEq
), MaxRays
));
142 value_clear(res
.x
.p
->arr
[1].d
);
143 res
.x
.p
->arr
[1] = *e
;
151 addeliminatedparams_partition(p
, CT
, CEq
, nparam
, MaxRays
);
152 for (i
= 0; i
< p
->size
/2; i
++)
153 addeliminatedparams_evalue(&p
->arr
[2*i
+1], CT
);
156 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
164 assert(value_notzero_p(e1
->d
));
165 assert(value_notzero_p(e2
->d
));
166 value_multiply(m
, e1
->x
.n
, e2
->d
);
167 value_multiply(m2
, e2
->x
.n
, e1
->d
);
170 else if (value_gt(m
, m2
))
180 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
182 if (value_notzero_p(e1
->d
)) {
184 if (value_zero_p(e2
->d
))
186 r
= mod_rational_smaller(e1
, e2
);
187 return r
== -1 ? 0 : r
;
189 if (value_notzero_p(e2
->d
))
191 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
193 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
196 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
198 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
203 static int mod_term_smaller(const evalue
*e1
, const evalue
*e2
)
205 assert(value_zero_p(e1
->d
));
206 assert(value_zero_p(e2
->d
));
207 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
208 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
209 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
212 static void check_order(const evalue
*e
)
217 if (value_notzero_p(e
->d
))
220 switch (e
->x
.p
->type
) {
222 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
223 check_order(&e
->x
.p
->arr
[2*i
+1]);
226 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
228 if (value_notzero_p(a
->d
))
230 switch (a
->x
.p
->type
) {
232 assert(mod_term_smaller(&e
->x
.p
->arr
[0], &a
->x
.p
->arr
[0]));
241 for (i
= 0; i
< e
->x
.p
->size
; ++i
) {
243 if (value_notzero_p(a
->d
))
245 switch (a
->x
.p
->type
) {
247 assert(e
->x
.p
->pos
< a
->x
.p
->pos
);
258 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
260 if (value_notzero_p(a
->d
))
262 switch (a
->x
.p
->type
) {
273 /* Negative pos means inequality */
274 /* s is negative of substitution if m is not zero */
283 struct fixed_param
*fixed
;
288 static int relations_depth(evalue
*e
)
293 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
294 e
= &e
->x
.p
->arr
[1], ++d
);
298 static void poly_denom_not_constant(evalue
**pp
, Value
*d
)
303 while (value_zero_p(p
->d
)) {
304 assert(p
->x
.p
->type
== polynomial
);
305 assert(p
->x
.p
->size
== 2);
306 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
307 value_lcm(*d
, p
->x
.p
->arr
[1].d
, d
);
313 static void poly_denom(evalue
*p
, Value
*d
)
315 poly_denom_not_constant(&p
, d
);
316 value_lcm(*d
, p
->d
, d
);
319 static void realloc_substitution(struct subst
*s
, int d
)
321 struct fixed_param
*n
;
324 for (i
= 0; i
< s
->n
; ++i
)
331 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
337 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
340 /* May have been reduced already */
341 if (value_notzero_p(m
->d
))
344 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
345 assert(m
->x
.p
->size
== 3);
347 /* fractional was inverted during reduction
348 * invert it back and move constant in
350 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
351 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
352 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
353 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
354 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
355 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
356 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
357 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
358 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
359 value_set_si(m
->x
.p
->arr
[1].d
, 1);
362 /* Oops. Nested identical relations. */
363 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
366 if (s
->n
>= s
->max
) {
367 int d
= relations_depth(r
);
368 realloc_substitution(s
, d
);
372 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
373 assert(p
->x
.p
->size
== 2);
376 assert(value_pos_p(f
->x
.n
));
378 value_init(s
->fixed
[s
->n
].m
);
379 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
380 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
381 value_init(s
->fixed
[s
->n
].d
);
382 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
383 value_init(s
->fixed
[s
->n
].s
.d
);
384 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
390 static int type_offset(enode
*p
)
392 return p
->type
== fractional
? 1 :
393 p
->type
== flooring
? 1 : 0;
396 static void reorder_terms_about(enode
*p
, evalue
*v
)
399 int offset
= type_offset(p
);
401 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
403 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
404 free_evalue_refs(&(p
->arr
[i
]));
410 static void reorder_terms(evalue
*e
)
415 assert(value_zero_p(e
->d
));
417 assert(p
->type
== fractional
); /* for now */
420 value_set_si(f
.d
, 0);
421 f
.x
.p
= new_enode(fractional
, 3, -1);
422 value_clear(f
.x
.p
->arr
[0].d
);
423 f
.x
.p
->arr
[0] = p
->arr
[0];
424 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
425 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
426 reorder_terms_about(p
, &f
);
432 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
438 if (value_notzero_p(e
->d
)) {
440 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
441 return; /* a rational number, its already reduced */
445 return; /* hum... an overflow probably occured */
447 /* First reduce the components of p */
448 add
= p
->type
== relation
;
449 for (i
=0; i
<p
->size
; i
++) {
451 add
= add_modulo_substitution(s
, e
);
453 if (i
== 0 && p
->type
==fractional
)
454 _reduce_evalue(&p
->arr
[i
], s
, 1);
456 _reduce_evalue(&p
->arr
[i
], s
, fract
);
458 if (add
&& i
== p
->size
-1) {
460 value_clear(s
->fixed
[s
->n
].m
);
461 value_clear(s
->fixed
[s
->n
].d
);
462 free_evalue_refs(&s
->fixed
[s
->n
].s
);
463 } else if (add
&& i
== 1)
464 s
->fixed
[s
->n
-1].pos
*= -1;
467 if (p
->type
==periodic
) {
469 /* Try to reduce the period */
470 for (i
=1; i
<=(p
->size
)/2; i
++) {
471 if ((p
->size
% i
)==0) {
473 /* Can we reduce the size to i ? */
475 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
476 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
479 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
483 you_lose
: /* OK, lets not do it */
488 /* Try to reduce its strength */
491 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
495 else if (p
->type
==polynomial
) {
496 for (k
= 0; s
&& k
< s
->n
; ++k
) {
497 if (s
->fixed
[k
].pos
== p
->pos
) {
498 int divide
= value_notone_p(s
->fixed
[k
].d
);
501 if (value_notzero_p(s
->fixed
[k
].m
)) {
504 assert(p
->size
== 2);
505 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
507 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
514 value_assign(d
.d
, s
->fixed
[k
].d
);
516 if (value_notzero_p(s
->fixed
[k
].m
))
517 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
519 value_set_si(d
.x
.n
, 1);
522 for (i
=p
->size
-1;i
>=1;i
--) {
523 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
525 emul(&d
, &p
->arr
[i
]);
526 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
527 free_evalue_refs(&(p
->arr
[i
]));
530 _reduce_evalue(&p
->arr
[0], s
, fract
);
533 free_evalue_refs(&d
);
539 /* Try to reduce the degree */
540 for (i
=p
->size
-1;i
>=1;i
--) {
541 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
543 /* Zero coefficient */
544 free_evalue_refs(&(p
->arr
[i
]));
549 /* Try to reduce its strength */
552 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
556 else if (p
->type
==fractional
) {
560 if (value_notzero_p(p
->arr
[0].d
)) {
562 value_assign(v
.d
, p
->arr
[0].d
);
564 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
569 evalue
*pp
= &p
->arr
[0];
570 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
571 assert(pp
->x
.p
->size
== 2);
573 /* search for exact duplicate among the modulo inequalities */
575 f
= &pp
->x
.p
->arr
[1];
576 for (k
= 0; s
&& k
< s
->n
; ++k
) {
577 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
578 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
579 value_eq(s
->fixed
[k
].m
, f
->d
) &&
580 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
587 /* replace { E/m } by { (E-1)/m } + 1/m */
592 evalue_set_si(&extra
, 1, 1);
593 value_assign(extra
.d
, g
);
594 eadd(&extra
, &v
.x
.p
->arr
[1]);
595 free_evalue_refs(&extra
);
597 /* We've been going in circles; stop now */
598 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
599 free_evalue_refs(&v
);
601 evalue_set_si(&v
, 0, 1);
606 value_set_si(v
.d
, 0);
607 v
.x
.p
= new_enode(fractional
, 3, -1);
608 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
609 value_assign(v
.x
.p
->arr
[1].d
, g
);
610 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
611 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
614 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
617 value_division(f
->d
, g
, f
->d
);
618 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
619 value_assign(f
->d
, g
);
620 value_decrement(f
->x
.n
, f
->x
.n
);
621 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
623 Gcd(f
->d
, f
->x
.n
, &g
);
624 value_division(f
->d
, f
->d
, g
);
625 value_division(f
->x
.n
, f
->x
.n
, g
);
634 /* reduction may have made this fractional arg smaller */
635 i
= reorder
? p
->size
: 1;
636 for ( ; i
< p
->size
; ++i
)
637 if (value_zero_p(p
->arr
[i
].d
) &&
638 p
->arr
[i
].x
.p
->type
== fractional
&&
639 !mod_term_smaller(e
, &p
->arr
[i
]))
643 value_set_si(v
.d
, 0);
644 v
.x
.p
= new_enode(fractional
, 3, -1);
645 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
646 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
647 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
655 evalue
*pp
= &p
->arr
[0];
658 poly_denom_not_constant(&pp
, &m
);
659 mpz_fdiv_r(r
, m
, pp
->d
);
660 if (value_notzero_p(r
)) {
662 value_set_si(v
.d
, 0);
663 v
.x
.p
= new_enode(fractional
, 3, -1);
665 value_multiply(r
, m
, pp
->x
.n
);
666 value_multiply(v
.x
.p
->arr
[1].d
, m
, pp
->d
);
667 value_init(v
.x
.p
->arr
[1].x
.n
);
668 mpz_fdiv_r(v
.x
.p
->arr
[1].x
.n
, r
, pp
->d
);
669 mpz_fdiv_q(r
, r
, pp
->d
);
671 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
672 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
674 while (value_zero_p(pp
->d
))
675 pp
= &pp
->x
.p
->arr
[0];
677 value_assign(pp
->d
, m
);
678 value_assign(pp
->x
.n
, r
);
680 Gcd(pp
->d
, pp
->x
.n
, &r
);
681 value_division(pp
->d
, pp
->d
, r
);
682 value_division(pp
->x
.n
, pp
->x
.n
, r
);
695 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
696 pp
= &pp
->x
.p
->arr
[0]) {
697 f
= &pp
->x
.p
->arr
[1];
698 assert(value_pos_p(f
->d
));
699 mpz_mul_ui(twice
, f
->x
.n
, 2);
700 if (value_lt(twice
, f
->d
))
702 if (value_eq(twice
, f
->d
))
710 value_set_si(v
.d
, 0);
711 v
.x
.p
= new_enode(fractional
, 3, -1);
712 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
713 poly_denom(&p
->arr
[0], &twice
);
714 value_assign(v
.x
.p
->arr
[1].d
, twice
);
715 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
716 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
717 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
719 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
720 pp
= &pp
->x
.p
->arr
[0]) {
721 f
= &pp
->x
.p
->arr
[1];
722 value_oppose(f
->x
.n
, f
->x
.n
);
723 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
725 value_division(pp
->d
, twice
, pp
->d
);
726 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
727 value_assign(pp
->d
, twice
);
728 value_oppose(pp
->x
.n
, pp
->x
.n
);
729 value_decrement(pp
->x
.n
, pp
->x
.n
);
730 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
732 /* Maybe we should do this during reduction of
735 Gcd(pp
->d
, pp
->x
.n
, &twice
);
736 value_division(pp
->d
, pp
->d
, twice
);
737 value_division(pp
->x
.n
, pp
->x
.n
, twice
);
747 reorder_terms_about(p
, &v
);
748 _reduce_evalue(&p
->arr
[1], s
, fract
);
751 /* Try to reduce the degree */
752 for (i
=p
->size
-1;i
>=2;i
--) {
753 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
755 /* Zero coefficient */
756 free_evalue_refs(&(p
->arr
[i
]));
761 /* Try to reduce its strength */
764 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
765 free_evalue_refs(&(p
->arr
[0]));
769 else if (p
->type
== flooring
) {
770 /* Try to reduce the degree */
771 for (i
=p
->size
-1;i
>=2;i
--) {
772 if (!EVALUE_IS_ZERO(p
->arr
[i
]))
774 /* Zero coefficient */
775 free_evalue_refs(&(p
->arr
[i
]));
780 /* Try to reduce its strength */
783 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
784 free_evalue_refs(&(p
->arr
[0]));
788 else if (p
->type
== relation
) {
789 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
790 free_evalue_refs(&(p
->arr
[2]));
791 free_evalue_refs(&(p
->arr
[0]));
798 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
799 free_evalue_refs(&(p
->arr
[2]));
802 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
803 free_evalue_refs(&(p
->arr
[1]));
804 free_evalue_refs(&(p
->arr
[0]));
805 evalue_set_si(e
, 0, 1);
812 /* Relation was reduced by means of an identical
813 * inequality => remove
815 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
818 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
819 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
821 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
823 free_evalue_refs(&(p
->arr
[2]));
827 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
829 evalue_set_si(e
, 0, 1);
830 free_evalue_refs(&(p
->arr
[1]));
832 free_evalue_refs(&(p
->arr
[0]));
838 } /* reduce_evalue */
840 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
845 for (k
= 0; k
< dim
; ++k
)
846 if (value_notzero_p(row
[k
+1]))
849 Vector_Normalize_Positive(row
+1, dim
+1, k
);
850 assert(s
->n
< s
->max
);
851 value_init(s
->fixed
[s
->n
].d
);
852 value_init(s
->fixed
[s
->n
].m
);
853 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
854 s
->fixed
[s
->n
].pos
= k
+1;
855 value_set_si(s
->fixed
[s
->n
].m
, 0);
856 r
= &s
->fixed
[s
->n
].s
;
858 for (l
= k
+1; l
< dim
; ++l
)
859 if (value_notzero_p(row
[l
+1])) {
860 value_set_si(r
->d
, 0);
861 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
862 value_init(r
->x
.p
->arr
[1].x
.n
);
863 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
864 value_set_si(r
->x
.p
->arr
[1].d
, 1);
868 value_oppose(r
->x
.n
, row
[dim
+1]);
869 value_set_si(r
->d
, 1);
873 static void _reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
, struct subst
*s
)
876 Polyhedron
*orig
= D
;
881 D
= DomainConvex(D
, 0);
882 if (!D
->next
&& D
->NbEq
) {
886 realloc_substitution(s
, dim
);
888 int d
= relations_depth(e
);
890 NALLOC(s
->fixed
, s
->max
);
893 for (j
= 0; j
< D
->NbEq
; ++j
)
894 add_substitution(s
, D
->Constraint
[j
], dim
);
898 _reduce_evalue(e
, s
, 0);
901 for (j
= 0; j
< s
->n
; ++j
) {
902 value_clear(s
->fixed
[j
].d
);
903 value_clear(s
->fixed
[j
].m
);
904 free_evalue_refs(&s
->fixed
[j
].s
);
909 void reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
)
911 struct subst s
= { NULL
, 0, 0 };
913 if (EVALUE_IS_ZERO(*e
))
917 evalue_set_si(e
, 0, 1);
920 _reduce_evalue_in_domain(e
, D
, &s
);
925 void reduce_evalue (evalue
*e
) {
926 struct subst s
= { NULL
, 0, 0 };
928 if (value_notzero_p(e
->d
))
929 return; /* a rational number, its already reduced */
931 if (e
->x
.p
->type
== partition
) {
934 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
935 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
937 /* This shouldn't really happen;
938 * Empty domains should not be added.
940 POL_ENSURE_VERTICES(D
);
942 _reduce_evalue_in_domain(&e
->x
.p
->arr
[2*i
+1], D
, &s
);
944 if (emptyQ(D
) || EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
945 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
946 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
947 value_clear(e
->x
.p
->arr
[2*i
].d
);
949 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
950 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
954 if (e
->x
.p
->size
== 0) {
956 evalue_set_si(e
, 0, 1);
959 _reduce_evalue(e
, &s
, 0);
964 static void print_evalue_r(FILE *DST
, const evalue
*e
, const char *const *pname
)
966 if(value_notzero_p(e
->d
)) {
967 if(value_notone_p(e
->d
)) {
968 value_print(DST
,VALUE_FMT
,e
->x
.n
);
970 value_print(DST
,VALUE_FMT
,e
->d
);
973 value_print(DST
,VALUE_FMT
,e
->x
.n
);
977 print_enode(DST
,e
->x
.p
,pname
);
981 void print_evalue(FILE *DST
, const evalue
*e
, const char * const *pname
)
983 print_evalue_r(DST
, e
, pname
);
984 if (value_notzero_p(e
->d
))
988 void print_enode(FILE *DST
, enode
*p
, const char *const *pname
)
993 fprintf(DST
, "NULL");
999 for (i
=0; i
<p
->size
; i
++) {
1000 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1004 fprintf(DST
, " }\n");
1008 for (i
=p
->size
-1; i
>=0; i
--) {
1009 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1010 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
1012 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
1014 fprintf(DST
, " )\n");
1018 for (i
=0; i
<p
->size
; i
++) {
1019 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1020 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
1022 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
1027 for (i
=p
->size
-1; i
>=1; i
--) {
1028 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1030 fprintf(DST
, " * ");
1031 fprintf(DST
, p
->type
== flooring
? "[" : "{");
1032 print_evalue_r(DST
, &p
->arr
[0], pname
);
1033 fprintf(DST
, p
->type
== flooring
? "]" : "}");
1035 fprintf(DST
, "^%d + ", i
-1);
1037 fprintf(DST
, " + ");
1040 fprintf(DST
, " )\n");
1044 print_evalue_r(DST
, &p
->arr
[0], pname
);
1045 fprintf(DST
, "= 0 ] * \n");
1046 print_evalue_r(DST
, &p
->arr
[1], pname
);
1048 fprintf(DST
, " +\n [ ");
1049 print_evalue_r(DST
, &p
->arr
[0], pname
);
1050 fprintf(DST
, "!= 0 ] * \n");
1051 print_evalue_r(DST
, &p
->arr
[2], pname
);
1055 char **new_names
= NULL
;
1056 const char *const *names
= pname
;
1057 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
1058 if (!pname
|| p
->pos
< maxdim
) {
1059 new_names
= ALLOCN(char *, maxdim
);
1060 for (i
= 0; i
< p
->pos
; ++i
) {
1062 new_names
[i
] = (char *)pname
[i
];
1064 new_names
[i
] = ALLOCN(char, 10);
1065 snprintf(new_names
[i
], 10, "%c", 'P'+i
);
1068 for ( ; i
< maxdim
; ++i
) {
1069 new_names
[i
] = ALLOCN(char, 10);
1070 snprintf(new_names
[i
], 10, "_p%d", i
);
1072 names
= (const char**)new_names
;
1075 for (i
=0; i
<p
->size
/2; i
++) {
1076 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
1077 print_evalue_r(DST
, &p
->arr
[2*i
+1], names
);
1078 if (value_notzero_p(p
->arr
[2*i
+1].d
))
1082 if (!pname
|| p
->pos
< maxdim
) {
1083 for (i
= pname
? p
->pos
: 0; i
< maxdim
; ++i
)
1096 static void eadd_rev(const evalue
*e1
, evalue
*res
)
1100 evalue_copy(&ev
, e1
);
1102 free_evalue_refs(res
);
1106 static void eadd_rev_cst(const evalue
*e1
, evalue
*res
)
1110 evalue_copy(&ev
, e1
);
1111 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
1112 free_evalue_refs(res
);
1116 static int is_zero_on(evalue
*e
, Polyhedron
*D
)
1121 tmp
.x
.p
= new_enode(partition
, 2, D
->Dimension
);
1122 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Domain_Copy(D
));
1123 evalue_copy(&tmp
.x
.p
->arr
[1], e
);
1124 reduce_evalue(&tmp
);
1125 is_zero
= EVALUE_IS_ZERO(tmp
);
1126 free_evalue_refs(&tmp
);
1130 struct section
{ Polyhedron
* D
; evalue E
; };
1132 void eadd_partitions(const evalue
*e1
, evalue
*res
)
1137 s
= (struct section
*)
1138 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
1139 sizeof(struct section
));
1141 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1142 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1143 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1146 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1147 assert(res
->x
.p
->size
>= 2);
1148 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1149 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
1151 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
1153 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1162 /* See if we can extend one of the domains in res to cover fd */
1163 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1164 if (is_zero_on(&res
->x
.p
->arr
[2*i
+1], fd
))
1166 if (i
< res
->x
.p
->size
/2) {
1167 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
1168 DomainConcat(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
])));
1171 value_init(s
[n
].E
.d
);
1172 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
1176 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1177 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
1178 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1180 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1181 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1187 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
1188 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1190 value_init(s
[n
].E
.d
);
1191 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1192 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1193 if (!emptyQ(fd
) && is_zero_on(&e1
->x
.p
->arr
[2*j
+1], fd
)) {
1194 d
= DomainConcat(fd
, d
);
1195 fd
= Empty_Polyhedron(fd
->Dimension
);
1201 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
1205 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1208 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1209 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1210 value_clear(res
->x
.p
->arr
[2*i
].d
);
1215 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1216 for (j
= 0; j
< n
; ++j
) {
1217 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1218 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1219 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1220 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1226 static void explicit_complement(evalue
*res
)
1228 enode
*rel
= new_enode(relation
, 3, 0);
1230 value_clear(rel
->arr
[0].d
);
1231 rel
->arr
[0] = res
->x
.p
->arr
[0];
1232 value_clear(rel
->arr
[1].d
);
1233 rel
->arr
[1] = res
->x
.p
->arr
[1];
1234 value_set_si(rel
->arr
[2].d
, 1);
1235 value_init(rel
->arr
[2].x
.n
);
1236 value_set_si(rel
->arr
[2].x
.n
, 0);
1241 void eadd(const evalue
*e1
, evalue
*res
)
1244 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1245 /* Add two rational numbers */
1251 value_multiply(m1
,e1
->x
.n
,res
->d
);
1252 value_multiply(m2
,res
->x
.n
,e1
->d
);
1253 value_addto(res
->x
.n
,m1
,m2
);
1254 value_multiply(res
->d
,e1
->d
,res
->d
);
1255 Gcd(res
->x
.n
,res
->d
,&g
);
1256 if (value_notone_p(g
)) {
1257 value_division(res
->d
,res
->d
,g
);
1258 value_division(res
->x
.n
,res
->x
.n
,g
);
1260 value_clear(g
); value_clear(m1
); value_clear(m2
);
1263 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
1264 switch (res
->x
.p
->type
) {
1266 /* Add the constant to the constant term of a polynomial*/
1267 eadd(e1
, &res
->x
.p
->arr
[0]);
1270 /* Add the constant to all elements of a periodic number */
1271 for (i
=0; i
<res
->x
.p
->size
; i
++) {
1272 eadd(e1
, &res
->x
.p
->arr
[i
]);
1276 fprintf(stderr
, "eadd: cannot add const with vector\n");
1280 eadd(e1
, &res
->x
.p
->arr
[1]);
1283 assert(EVALUE_IS_ZERO(*e1
));
1284 break; /* Do nothing */
1286 /* Create (zero) complement if needed */
1287 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
1288 explicit_complement(res
);
1289 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1290 eadd(e1
, &res
->x
.p
->arr
[i
]);
1296 /* add polynomial or periodic to constant
1297 * you have to exchange e1 and res, before doing addition */
1299 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1303 else { // ((e1->d==0) && (res->d==0))
1304 assert(!((e1
->x
.p
->type
== partition
) ^
1305 (res
->x
.p
->type
== partition
)));
1306 if (e1
->x
.p
->type
== partition
) {
1307 eadd_partitions(e1
, res
);
1310 if (e1
->x
.p
->type
== relation
&&
1311 (res
->x
.p
->type
!= relation
||
1312 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1316 if (res
->x
.p
->type
== relation
) {
1317 if (e1
->x
.p
->type
== relation
&&
1318 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1319 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1320 explicit_complement(res
);
1321 for (i
= 1; i
< e1
->x
.p
->size
; ++i
)
1322 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1325 if (res
->x
.p
->size
< 3)
1326 explicit_complement(res
);
1327 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1328 eadd(e1
, &res
->x
.p
->arr
[i
]);
1331 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
1332 /* adding to evalues of different type. two cases are possible
1333 * res is periodic and e1 is polynomial, you have to exchange
1334 * e1 and res then to add e1 to the constant term of res */
1335 if (e1
->x
.p
->type
== polynomial
) {
1336 eadd_rev_cst(e1
, res
);
1338 else if (res
->x
.p
->type
== polynomial
) {
1339 /* res is polynomial and e1 is periodic,
1340 add e1 to the constant term of res */
1342 eadd(e1
,&res
->x
.p
->arr
[0]);
1348 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
1349 ((res
->x
.p
->type
== fractional
||
1350 res
->x
.p
->type
== flooring
) &&
1351 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1352 /* adding evalues of different position (i.e function of different unknowns
1353 * to case are possible */
1355 switch (res
->x
.p
->type
) {
1358 if (mod_term_smaller(res
, e1
))
1359 eadd(e1
,&res
->x
.p
->arr
[1]);
1361 eadd_rev_cst(e1
, res
);
1363 case polynomial
: // res and e1 are polynomials
1364 // add e1 to the constant term of res
1366 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1367 eadd(e1
,&res
->x
.p
->arr
[0]);
1369 eadd_rev_cst(e1
, res
);
1370 // value_clear(g); value_clear(m1); value_clear(m2);
1372 case periodic
: // res and e1 are pointers to periodic numbers
1373 //add e1 to all elements of res
1375 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1376 for (i
=0;i
<res
->x
.p
->size
;i
++) {
1377 eadd(e1
,&res
->x
.p
->arr
[i
]);
1388 //same type , same pos and same size
1389 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1390 // add any element in e1 to the corresponding element in res
1391 i
= type_offset(res
->x
.p
);
1393 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1394 for (; i
<res
->x
.p
->size
; i
++) {
1395 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1400 /* Sizes are different */
1401 switch(res
->x
.p
->type
) {
1405 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1406 /* new enode and add res to that new node. If you do not do */
1407 /* that, you lose the the upper weight part of e1 ! */
1409 if(e1
->x
.p
->size
> res
->x
.p
->size
)
1412 i
= type_offset(res
->x
.p
);
1414 assert(eequal(&e1
->x
.p
->arr
[0],
1415 &res
->x
.p
->arr
[0]));
1416 for (; i
<e1
->x
.p
->size
; i
++) {
1417 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1424 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1427 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1428 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1429 to add periodicaly elements of e1 to elements of ne, and finaly to
1434 value_init(ex
); value_init(ey
);value_init(ep
);
1437 value_set_si(ex
,e1
->x
.p
->size
);
1438 value_set_si(ey
,res
->x
.p
->size
);
1439 value_assign (ep
,*Lcm(ex
,ey
));
1440 p
=(int)mpz_get_si(ep
);
1441 ne
= (evalue
*) malloc (sizeof(evalue
));
1443 value_set_si( ne
->d
,0);
1445 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1447 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1450 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1453 value_assign(res
->d
, ne
->d
);
1459 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
1468 static void emul_rev(const evalue
*e1
, evalue
*res
)
1472 evalue_copy(&ev
, e1
);
1474 free_evalue_refs(res
);
1478 static void emul_poly(const evalue
*e1
, evalue
*res
)
1480 int i
, j
, offset
= type_offset(res
->x
.p
);
1483 int size
= (e1
->x
.p
->size
+ res
->x
.p
->size
- offset
- 1);
1485 p
= new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1487 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1488 if (!EVALUE_IS_ZERO(e1
->x
.p
->arr
[i
]))
1491 /* special case pure power */
1492 if (i
== e1
->x
.p
->size
-1) {
1494 value_clear(p
->arr
[0].d
);
1495 p
->arr
[0] = res
->x
.p
->arr
[0];
1497 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1498 evalue_set_si(&p
->arr
[i
], 0, 1);
1499 for (i
= offset
; i
< res
->x
.p
->size
; ++i
) {
1500 value_clear(p
->arr
[i
+e1
->x
.p
->size
-offset
-1].d
);
1501 p
->arr
[i
+e1
->x
.p
->size
-offset
-1] = res
->x
.p
->arr
[i
];
1502 emul(&e1
->x
.p
->arr
[e1
->x
.p
->size
-1],
1503 &p
->arr
[i
+e1
->x
.p
->size
-offset
-1]);
1511 value_set_si(tmp
.d
,0);
1514 evalue_copy(&p
->arr
[0], &e1
->x
.p
->arr
[0]);
1515 for (i
= offset
; i
< e1
->x
.p
->size
; i
++) {
1516 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1517 emul(&res
->x
.p
->arr
[offset
], &tmp
.x
.p
->arr
[i
]);
1520 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1521 for (i
= offset
+1; i
<res
->x
.p
->size
; i
++)
1522 for (j
= offset
; j
<e1
->x
.p
->size
; j
++) {
1525 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1526 emul(&res
->x
.p
->arr
[i
], &ev
);
1527 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-offset
]);
1528 free_evalue_refs(&ev
);
1530 free_evalue_refs(res
);
1534 void emul_partitions(const evalue
*e1
, evalue
*res
)
1539 s
= (struct section
*)
1540 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1541 sizeof(struct section
));
1543 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1544 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1545 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1548 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1549 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1550 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1551 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1557 /* This code is only needed because the partitions
1558 are not true partitions.
1560 for (k
= 0; k
< n
; ++k
) {
1561 if (DomainIncludes(s
[k
].D
, d
))
1563 if (DomainIncludes(d
, s
[k
].D
)) {
1564 Domain_Free(s
[k
].D
);
1565 free_evalue_refs(&s
[k
].E
);
1576 value_init(s
[n
].E
.d
);
1577 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1578 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1582 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1583 value_clear(res
->x
.p
->arr
[2*i
].d
);
1584 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1589 evalue_set_si(res
, 0, 1);
1591 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1592 for (j
= 0; j
< n
; ++j
) {
1593 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1594 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1595 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1596 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1603 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1605 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1606 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1607 * evalues is not treated here */
1609 void emul(const evalue
*e1
, evalue
*res
)
1613 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1614 fprintf(stderr
, "emul: do not proced on evector type !\n");
1618 if (EVALUE_IS_ZERO(*res
))
1621 if (EVALUE_IS_ONE(*e1
))
1624 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1625 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1626 emul_partitions(e1
, res
);
1629 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1630 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1631 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1633 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1634 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1635 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1636 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3) {
1637 free_evalue_refs(&res
->x
.p
->arr
[2]);
1640 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1641 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1644 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1645 emul(e1
, &res
->x
.p
->arr
[i
]);
1647 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1648 switch(e1
->x
.p
->type
) {
1650 switch(res
->x
.p
->type
) {
1652 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1653 /* Product of two polynomials of the same variable */
1658 /* Product of two polynomials of different variables */
1660 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1661 for( i
=0; i
<res
->x
.p
->size
; i
++)
1662 emul(e1
, &res
->x
.p
->arr
[i
]);
1671 /* Product of a polynomial and a periodic or fractional */
1678 switch(res
->x
.p
->type
) {
1680 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1681 /* Product of two periodics of the same parameter and period */
1683 for(i
=0; i
<res
->x
.p
->size
;i
++)
1684 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1689 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1690 /* Product of two periodics of the same parameter and different periods */
1694 value_init(x
); value_init(y
);value_init(z
);
1697 value_set_si(x
,e1
->x
.p
->size
);
1698 value_set_si(y
,res
->x
.p
->size
);
1699 value_assign (z
,*Lcm(x
,y
));
1700 lcm
=(int)mpz_get_si(z
);
1701 newp
= (evalue
*) malloc (sizeof(evalue
));
1702 value_init(newp
->d
);
1703 value_set_si( newp
->d
,0);
1704 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1705 for(i
=0;i
<lcm
;i
++) {
1706 evalue_copy(&newp
->x
.p
->arr
[i
],
1707 &res
->x
.p
->arr
[i
%iy
]);
1710 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1712 value_assign(res
->d
,newp
->d
);
1715 value_clear(x
); value_clear(y
);value_clear(z
);
1719 /* Product of two periodics of different parameters */
1721 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1722 for(i
=0; i
<res
->x
.p
->size
; i
++)
1723 emul(e1
, &(res
->x
.p
->arr
[i
]));
1731 /* Product of a periodic and a polynomial */
1733 for(i
=0; i
<res
->x
.p
->size
; i
++)
1734 emul(e1
, &(res
->x
.p
->arr
[i
]));
1741 switch(res
->x
.p
->type
) {
1743 for(i
=0; i
<res
->x
.p
->size
; i
++)
1744 emul(e1
, &(res
->x
.p
->arr
[i
]));
1751 assert(e1
->x
.p
->type
== res
->x
.p
->type
);
1752 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1753 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1756 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1757 if (e1
->x
.p
->type
!= fractional
|| !value_two_p(d
.d
))
1762 value_set_si(d
.x
.n
, 1);
1763 /* { x }^2 == { x }/2 */
1764 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1765 assert(e1
->x
.p
->size
== 3);
1766 assert(res
->x
.p
->size
== 3);
1768 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1770 eadd(&res
->x
.p
->arr
[1], &tmp
);
1771 emul(&e1
->x
.p
->arr
[2], &tmp
);
1772 emul(&e1
->x
.p
->arr
[1], res
);
1773 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1774 free_evalue_refs(&tmp
);
1779 if(mod_term_smaller(res
, e1
))
1780 for(i
=1; i
<res
->x
.p
->size
; i
++)
1781 emul(e1
, &(res
->x
.p
->arr
[i
]));
1796 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1797 /* Product of two rational numbers */
1801 value_multiply(res
->d
,e1
->d
,res
->d
);
1802 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1803 Gcd(res
->x
.n
, res
->d
,&g
);
1804 if (value_notone_p(g
)) {
1805 value_division(res
->d
,res
->d
,g
);
1806 value_division(res
->x
.n
,res
->x
.n
,g
);
1812 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1813 /* Product of an expression (polynomial or peririodic) and a rational number */
1819 /* Product of a rationel number and an expression (polynomial or peririodic) */
1821 i
= type_offset(res
->x
.p
);
1822 for (; i
<res
->x
.p
->size
; i
++)
1823 emul(e1
, &res
->x
.p
->arr
[i
]);
1833 /* Frees mask content ! */
1834 void emask(evalue
*mask
, evalue
*res
) {
1841 if (EVALUE_IS_ZERO(*res
)) {
1842 free_evalue_refs(mask
);
1846 assert(value_zero_p(mask
->d
));
1847 assert(mask
->x
.p
->type
== partition
);
1848 assert(value_zero_p(res
->d
));
1849 assert(res
->x
.p
->type
== partition
);
1850 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1851 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1852 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1853 pos
= res
->x
.p
->pos
;
1855 s
= (struct section
*)
1856 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1857 sizeof(struct section
));
1861 evalue_set_si(&mone
, -1, 1);
1864 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1865 assert(mask
->x
.p
->size
>= 2);
1866 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1867 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1869 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1871 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1880 value_init(s
[n
].E
.d
);
1881 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1885 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1886 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1889 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1890 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1891 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1892 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1894 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1895 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1901 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1902 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1904 value_init(s
[n
].E
.d
);
1905 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1906 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1912 /* Just ignore; this may have been previously masked off */
1914 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1918 free_evalue_refs(&mone
);
1919 free_evalue_refs(mask
);
1920 free_evalue_refs(res
);
1923 evalue_set_si(res
, 0, 1);
1925 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1926 for (j
= 0; j
< n
; ++j
) {
1927 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1928 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1929 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1936 void evalue_copy(evalue
*dst
, const evalue
*src
)
1938 value_assign(dst
->d
, src
->d
);
1939 if(value_notzero_p(src
->d
)) {
1940 value_init(dst
->x
.n
);
1941 value_assign(dst
->x
.n
, src
->x
.n
);
1943 dst
->x
.p
= ecopy(src
->x
.p
);
1946 evalue
*evalue_dup(const evalue
*e
)
1948 evalue
*res
= ALLOC(evalue
);
1950 evalue_copy(res
, e
);
1954 enode
*new_enode(enode_type type
,int size
,int pos
) {
1960 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1963 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1967 for(i
=0; i
<size
; i
++) {
1968 value_init(res
->arr
[i
].d
);
1969 value_set_si(res
->arr
[i
].d
,0);
1970 res
->arr
[i
].x
.p
= 0;
1975 enode
*ecopy(enode
*e
) {
1980 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1981 for(i
=0;i
<e
->size
;++i
) {
1982 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1983 if(value_zero_p(res
->arr
[i
].d
))
1984 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1985 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1986 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1988 value_init(res
->arr
[i
].x
.n
);
1989 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1995 int ecmp(const evalue
*e1
, const evalue
*e2
)
2001 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
2005 value_multiply(m
, e1
->x
.n
, e2
->d
);
2006 value_multiply(m2
, e2
->x
.n
, e1
->d
);
2008 if (value_lt(m
, m2
))
2010 else if (value_gt(m
, m2
))
2020 if (value_notzero_p(e1
->d
))
2022 if (value_notzero_p(e2
->d
))
2028 if (p1
->type
!= p2
->type
)
2029 return p1
->type
- p2
->type
;
2030 if (p1
->pos
!= p2
->pos
)
2031 return p1
->pos
- p2
->pos
;
2032 if (p1
->size
!= p2
->size
)
2033 return p1
->size
- p2
->size
;
2035 for (i
= p1
->size
-1; i
>= 0; --i
)
2036 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
2042 int eequal(const evalue
*e1
, const evalue
*e2
)
2047 if (value_ne(e1
->d
,e2
->d
))
2050 /* e1->d == e2->d */
2051 if (value_notzero_p(e1
->d
)) {
2052 if (value_ne(e1
->x
.n
,e2
->x
.n
))
2055 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2059 /* e1->d == e2->d == 0 */
2062 if (p1
->type
!= p2
->type
) return 0;
2063 if (p1
->size
!= p2
->size
) return 0;
2064 if (p1
->pos
!= p2
->pos
) return 0;
2065 for (i
=0; i
<p1
->size
; i
++)
2066 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
2071 void free_evalue_refs(evalue
*e
) {
2076 if (EVALUE_IS_DOMAIN(*e
)) {
2077 Domain_Free(EVALUE_DOMAIN(*e
));
2080 } else if (value_pos_p(e
->d
)) {
2082 /* 'e' stores a constant */
2084 value_clear(e
->x
.n
);
2087 assert(value_zero_p(e
->d
));
2090 if (!p
) return; /* null pointer */
2091 for (i
=0; i
<p
->size
; i
++) {
2092 free_evalue_refs(&(p
->arr
[i
]));
2096 } /* free_evalue_refs */
2098 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
2099 Vector
* val
, evalue
*res
)
2101 unsigned nparam
= periods
->Size
;
2104 double d
= compute_evalue(e
, val
->p
);
2105 d
*= VALUE_TO_DOUBLE(m
);
2110 value_assign(res
->d
, m
);
2111 value_init(res
->x
.n
);
2112 value_set_double(res
->x
.n
, d
);
2113 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
2116 if (value_one_p(periods
->p
[p
]))
2117 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
2122 value_assign(tmp
, periods
->p
[p
]);
2123 value_set_si(res
->d
, 0);
2124 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
2126 value_decrement(tmp
, tmp
);
2127 value_assign(val
->p
[p
], tmp
);
2128 mod2table_r(e
, periods
, m
, p
+1, val
,
2129 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
2130 } while (value_pos_p(tmp
));
2136 static void rel2table(evalue
*e
, int zero
)
2138 if (value_pos_p(e
->d
)) {
2139 if (value_zero_p(e
->x
.n
) == zero
)
2140 value_set_si(e
->x
.n
, 1);
2142 value_set_si(e
->x
.n
, 0);
2143 value_set_si(e
->d
, 1);
2146 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
2147 rel2table(&e
->x
.p
->arr
[i
], zero
);
2151 void evalue_mod2table(evalue
*e
, int nparam
)
2156 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2159 for (i
=0; i
<p
->size
; i
++) {
2160 evalue_mod2table(&(p
->arr
[i
]), nparam
);
2162 if (p
->type
== relation
) {
2167 evalue_copy(©
, &p
->arr
[0]);
2169 rel2table(&p
->arr
[0], 1);
2170 emul(&p
->arr
[0], &p
->arr
[1]);
2172 rel2table(©
, 0);
2173 emul(©
, &p
->arr
[2]);
2174 eadd(&p
->arr
[2], &p
->arr
[1]);
2175 free_evalue_refs(&p
->arr
[2]);
2176 free_evalue_refs(©
);
2178 free_evalue_refs(&p
->arr
[0]);
2182 } else if (p
->type
== fractional
) {
2183 Vector
*periods
= Vector_Alloc(nparam
);
2184 Vector
*val
= Vector_Alloc(nparam
);
2190 value_set_si(tmp
, 1);
2191 Vector_Set(periods
->p
, 1, nparam
);
2192 Vector_Set(val
->p
, 0, nparam
);
2193 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
2196 assert(p
->type
== polynomial
);
2197 assert(p
->size
== 2);
2198 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
2199 value_lcm(tmp
, p
->arr
[1].d
, &tmp
);
2201 value_lcm(tmp
, ev
->d
, &tmp
);
2203 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
2206 evalue_set_si(&res
, 0, 1);
2207 /* Compute the polynomial using Horner's rule */
2208 for (i
=p
->size
-1;i
>1;i
--) {
2209 eadd(&p
->arr
[i
], &res
);
2212 eadd(&p
->arr
[1], &res
);
2214 free_evalue_refs(e
);
2215 free_evalue_refs(&EP
);
2220 Vector_Free(periods
);
2222 } /* evalue_mod2table */
2224 /********************************************************/
2225 /* function in domain */
2226 /* check if the parameters in list_args */
2227 /* verifies the constraints of Domain P */
2228 /********************************************************/
2229 int in_domain(Polyhedron
*P
, Value
*list_args
)
2232 Value v
; /* value of the constraint of a row when
2233 parameters are instantiated*/
2237 for (row
= 0; row
< P
->NbConstraints
; row
++) {
2238 Inner_Product(P
->Constraint
[row
]+1, list_args
, P
->Dimension
, &v
);
2239 value_addto(v
, v
, P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
2240 if (value_neg_p(v
) ||
2241 value_zero_p(P
->Constraint
[row
][0]) && value_notzero_p(v
)) {
2248 return in
|| (P
->next
&& in_domain(P
->next
, list_args
));
2251 /****************************************************/
2252 /* function compute enode */
2253 /* compute the value of enode p with parameters */
2254 /* list "list_args */
2255 /* compute the polynomial or the periodic */
2256 /****************************************************/
2258 static double compute_enode(enode
*p
, Value
*list_args
) {
2270 if (p
->type
== polynomial
) {
2272 value_assign(param
,list_args
[p
->pos
-1]);
2274 /* Compute the polynomial using Horner's rule */
2275 for (i
=p
->size
-1;i
>0;i
--) {
2276 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2277 res
*=VALUE_TO_DOUBLE(param
);
2279 res
+=compute_evalue(&p
->arr
[0],list_args
);
2281 else if (p
->type
== fractional
) {
2282 double d
= compute_evalue(&p
->arr
[0], list_args
);
2283 d
-= floor(d
+1e-10);
2285 /* Compute the polynomial using Horner's rule */
2286 for (i
=p
->size
-1;i
>1;i
--) {
2287 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2290 res
+=compute_evalue(&p
->arr
[1],list_args
);
2292 else if (p
->type
== flooring
) {
2293 double d
= compute_evalue(&p
->arr
[0], list_args
);
2296 /* Compute the polynomial using Horner's rule */
2297 for (i
=p
->size
-1;i
>1;i
--) {
2298 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2301 res
+=compute_evalue(&p
->arr
[1],list_args
);
2303 else if (p
->type
== periodic
) {
2304 value_assign(m
,list_args
[p
->pos
-1]);
2306 /* Choose the right element of the periodic */
2307 value_set_si(param
,p
->size
);
2308 value_pmodulus(m
,m
,param
);
2309 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2311 else if (p
->type
== relation
) {
2312 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2313 res
= compute_evalue(&p
->arr
[1], list_args
);
2314 else if (p
->size
> 2)
2315 res
= compute_evalue(&p
->arr
[2], list_args
);
2317 else if (p
->type
== partition
) {
2318 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2319 Value
*vals
= list_args
;
2322 for (i
= 0; i
< dim
; ++i
) {
2323 value_init(vals
[i
]);
2325 value_assign(vals
[i
], list_args
[i
]);
2328 for (i
= 0; i
< p
->size
/2; ++i
)
2329 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2330 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2334 for (i
= 0; i
< dim
; ++i
)
2335 value_clear(vals
[i
]);
2344 } /* compute_enode */
2346 /*************************************************/
2347 /* return the value of Ehrhart Polynomial */
2348 /* It returns a double, because since it is */
2349 /* a recursive function, some intermediate value */
2350 /* might not be integral */
2351 /*************************************************/
2353 double compute_evalue(const evalue
*e
, Value
*list_args
)
2357 if (value_notzero_p(e
->d
)) {
2358 if (value_notone_p(e
->d
))
2359 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2361 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2364 res
= compute_enode(e
->x
.p
,list_args
);
2366 } /* compute_evalue */
2369 /****************************************************/
2370 /* function compute_poly : */
2371 /* Check for the good validity domain */
2372 /* return the number of point in the Polyhedron */
2373 /* in allocated memory */
2374 /* Using the Ehrhart pseudo-polynomial */
2375 /****************************************************/
2376 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2379 /* double d; int i; */
2381 tmp
= (Value
*) malloc (sizeof(Value
));
2382 assert(tmp
!= NULL
);
2384 value_set_si(*tmp
,0);
2387 return(tmp
); /* no ehrhart polynomial */
2388 if(en
->ValidityDomain
) {
2389 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2390 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2395 return(tmp
); /* no Validity Domain */
2397 if(in_domain(en
->ValidityDomain
,list_args
)) {
2399 #ifdef EVAL_EHRHART_DEBUG
2400 Print_Domain(stdout
,en
->ValidityDomain
);
2401 print_evalue(stdout
,&en
->EP
);
2404 /* d = compute_evalue(&en->EP,list_args);
2406 printf("(double)%lf = %d\n", d, i ); */
2407 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2413 value_set_si(*tmp
,0);
2414 return(tmp
); /* no compatible domain with the arguments */
2415 } /* compute_poly */
2417 static evalue
*eval_polynomial(const enode
*p
, int offset
,
2418 evalue
*base
, Value
*values
)
2423 res
= evalue_zero();
2424 for (i
= p
->size
-1; i
> offset
; --i
) {
2425 c
= evalue_eval(&p
->arr
[i
], values
);
2427 free_evalue_refs(c
);
2431 c
= evalue_eval(&p
->arr
[offset
], values
);
2433 free_evalue_refs(c
);
2439 evalue
*evalue_eval(const evalue
*e
, Value
*values
)
2446 if (value_notzero_p(e
->d
)) {
2447 res
= ALLOC(evalue
);
2449 evalue_copy(res
, e
);
2452 switch (e
->x
.p
->type
) {
2454 value_init(param
.x
.n
);
2455 value_assign(param
.x
.n
, values
[e
->x
.p
->pos
-1]);
2456 value_init(param
.d
);
2457 value_set_si(param
.d
, 1);
2459 res
= eval_polynomial(e
->x
.p
, 0, ¶m
, values
);
2460 free_evalue_refs(¶m
);
2463 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2464 mpz_fdiv_r(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2466 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2467 free_evalue_refs(param2
);
2471 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2472 mpz_fdiv_q(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2473 value_set_si(param2
->d
, 1);
2475 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2476 free_evalue_refs(param2
);
2480 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2481 if (value_zero_p(param2
->x
.n
))
2482 res
= evalue_eval(&e
->x
.p
->arr
[1], values
);
2483 else if (e
->x
.p
->size
> 2)
2484 res
= evalue_eval(&e
->x
.p
->arr
[2], values
);
2486 res
= evalue_zero();
2487 free_evalue_refs(param2
);
2491 assert(e
->x
.p
->pos
== EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
);
2492 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2493 if (in_domain(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), values
)) {
2494 res
= evalue_eval(&e
->x
.p
->arr
[2*i
+1], values
);
2498 res
= evalue_zero();
2506 size_t value_size(Value v
) {
2507 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2508 * sizeof(v
[0]._mp_d
[0]);
2511 size_t domain_size(Polyhedron
*D
)
2514 size_t s
= sizeof(*D
);
2516 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2517 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2518 s
+= value_size(D
->Constraint
[i
][j
]);
2521 for (i = 0; i < D->NbRays; ++i)
2522 for (j = 0; j < D->Dimension+2; ++j)
2523 s += value_size(D->Ray[i][j]);
2526 return D
->next
? s
+domain_size(D
->next
) : s
;
2529 size_t enode_size(enode
*p
) {
2530 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2533 if (p
->type
== partition
)
2534 for (i
= 0; i
< p
->size
/2; ++i
) {
2535 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2536 s
+= evalue_size(&p
->arr
[2*i
+1]);
2539 for (i
= 0; i
< p
->size
; ++i
) {
2540 s
+= evalue_size(&p
->arr
[i
]);
2545 size_t evalue_size(evalue
*e
)
2547 size_t s
= sizeof(*e
);
2548 s
+= value_size(e
->d
);
2549 if (value_notzero_p(e
->d
))
2550 s
+= value_size(e
->x
.n
);
2552 s
+= enode_size(e
->x
.p
);
2556 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2558 evalue
*found
= NULL
;
2563 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2566 value_init(offset
.d
);
2567 value_init(offset
.x
.n
);
2568 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2569 value_lcm(m
, offset
.d
, &offset
.d
);
2570 value_set_si(offset
.x
.n
, 1);
2573 evalue_copy(©
, cst
);
2576 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2578 if (eequal(base
, &e
->x
.p
->arr
[0]))
2579 found
= &e
->x
.p
->arr
[0];
2581 value_set_si(offset
.x
.n
, -2);
2584 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2586 if (eequal(base
, &e
->x
.p
->arr
[0]))
2589 free_evalue_refs(cst
);
2590 free_evalue_refs(&offset
);
2593 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2594 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2599 static evalue
*find_relation_pair(evalue
*e
)
2602 evalue
*found
= NULL
;
2604 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2607 if (e
->x
.p
->type
== fractional
) {
2612 poly_denom(&e
->x
.p
->arr
[0], &m
);
2614 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2615 cst
= &cst
->x
.p
->arr
[0])
2618 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2619 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2624 i
= e
->x
.p
->type
== relation
;
2625 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2626 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2631 void evalue_mod2relation(evalue
*e
) {
2634 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2637 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2638 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2639 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2640 value_clear(e
->x
.p
->arr
[2*i
].d
);
2641 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2643 if (2*i
< e
->x
.p
->size
) {
2644 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2645 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2650 if (e
->x
.p
->size
== 0) {
2652 evalue_set_si(e
, 0, 1);
2658 while ((d
= find_relation_pair(e
)) != NULL
) {
2662 value_init(split
.d
);
2663 value_set_si(split
.d
, 0);
2664 split
.x
.p
= new_enode(relation
, 3, 0);
2665 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2666 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2668 ev
= &split
.x
.p
->arr
[0];
2669 value_set_si(ev
->d
, 0);
2670 ev
->x
.p
= new_enode(fractional
, 3, -1);
2671 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2672 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2673 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2679 free_evalue_refs(&split
);
2683 static int evalue_comp(const void * a
, const void * b
)
2685 const evalue
*e1
= *(const evalue
**)a
;
2686 const evalue
*e2
= *(const evalue
**)b
;
2687 return ecmp(e1
, e2
);
2690 void evalue_combine(evalue
*e
)
2697 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2700 NALLOC(evs
, e
->x
.p
->size
/2);
2701 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2702 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2703 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2704 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2705 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2706 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2707 value_clear(p
->arr
[2*k
].d
);
2708 value_clear(p
->arr
[2*k
+1].d
);
2709 p
->arr
[2*k
] = *(evs
[i
]-1);
2710 p
->arr
[2*k
+1] = *(evs
[i
]);
2713 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2716 value_clear((evs
[i
]-1)->d
);
2720 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2721 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2722 free_evalue_refs(evs
[i
]);
2726 for (i
= 2*k
; i
< p
->size
; ++i
)
2727 value_clear(p
->arr
[i
].d
);
2734 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2736 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2738 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2741 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2742 Polyhedron
*D
, *N
, **P
;
2745 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2752 if (D
->NbEq
<= H
->NbEq
) {
2758 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2759 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2760 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2761 reduce_evalue(&tmp
);
2762 if (value_notzero_p(tmp
.d
) ||
2763 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2766 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2767 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2770 free_evalue_refs(&tmp
);
2776 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2778 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2780 value_clear(e
->x
.p
->arr
[2*i
].d
);
2781 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2783 if (2*i
< e
->x
.p
->size
) {
2784 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2785 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2792 H
= DomainConvex(D
, 0);
2793 E
= DomainDifference(H
, D
, 0);
2795 D
= DomainDifference(H
, E
, 0);
2798 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2802 /* Use smallest representative for coefficients in affine form in
2803 * argument of fractional.
2804 * Since any change will make the argument non-standard,
2805 * the containing evalue will have to be reduced again afterward.
2807 static void fractional_minimal_coefficients(enode
*p
)
2813 assert(p
->type
== fractional
);
2815 while (value_zero_p(pp
->d
)) {
2816 assert(pp
->x
.p
->type
== polynomial
);
2817 assert(pp
->x
.p
->size
== 2);
2818 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2819 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2820 if (value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2821 value_subtract(pp
->x
.p
->arr
[1].x
.n
,
2822 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2823 pp
= &pp
->x
.p
->arr
[0];
2829 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2834 unsigned dim
= D
->Dimension
;
2835 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2838 assert(p
->type
== fractional
|| p
->type
== flooring
);
2839 value_set_si(T
->p
[1][dim
], 1);
2840 evalue_extract_affine(&p
->arr
[0], T
->p
[0], &T
->p
[0][dim
], d
);
2841 I
= DomainImage(D
, T
, 0);
2842 H
= DomainConvex(I
, 0);
2852 static void replace_by_affine(evalue
*e
, Value offset
)
2859 value_init(inc
.x
.n
);
2860 value_set_si(inc
.d
, 1);
2861 value_oppose(inc
.x
.n
, offset
);
2862 eadd(&inc
, &p
->arr
[0]);
2863 reorder_terms_about(p
, &p
->arr
[0]); /* frees arr[0] */
2867 free_evalue_refs(&inc
);
2870 int evalue_range_reduction_in_domain(evalue
*e
, Polyhedron
*D
)
2879 if (value_notzero_p(e
->d
))
2884 if (p
->type
== relation
) {
2891 fractional_minimal_coefficients(p
->arr
[0].x
.p
);
2892 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
);
2893 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2894 equal
= value_eq(min
, max
);
2895 mpz_cdiv_q(min
, min
, d
);
2896 mpz_fdiv_q(max
, max
, d
);
2898 if (bounded
&& value_gt(min
, max
)) {
2904 evalue_set_si(e
, 0, 1);
2907 free_evalue_refs(&(p
->arr
[1]));
2908 free_evalue_refs(&(p
->arr
[0]));
2914 return r
? r
: evalue_range_reduction_in_domain(e
, D
);
2915 } else if (bounded
&& equal
) {
2918 free_evalue_refs(&(p
->arr
[2]));
2921 free_evalue_refs(&(p
->arr
[0]));
2927 return evalue_range_reduction_in_domain(e
, D
);
2928 } else if (bounded
&& value_eq(min
, max
)) {
2929 /* zero for a single value */
2931 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2932 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2933 value_multiply(min
, min
, d
);
2934 value_subtract(M
->p
[0][D
->Dimension
+1],
2935 M
->p
[0][D
->Dimension
+1], min
);
2936 E
= DomainAddConstraints(D
, M
, 0);
2942 r
= evalue_range_reduction_in_domain(&p
->arr
[1], E
);
2944 r
|= evalue_range_reduction_in_domain(&p
->arr
[2], D
);
2946 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2954 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2957 i
= p
->type
== relation
? 1 :
2958 p
->type
== fractional
? 1 : 0;
2959 for (; i
<p
->size
; i
++)
2960 r
|= evalue_range_reduction_in_domain(&p
->arr
[i
], D
);
2962 if (p
->type
!= fractional
) {
2963 if (r
&& p
->type
== polynomial
) {
2966 value_set_si(f
.d
, 0);
2967 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2968 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2969 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2970 reorder_terms_about(p
, &f
);
2981 fractional_minimal_coefficients(p
);
2982 I
= polynomial_projection(p
, D
, &d
, NULL
);
2983 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2984 mpz_fdiv_q(min
, min
, d
);
2985 mpz_fdiv_q(max
, max
, d
);
2986 value_subtract(d
, max
, min
);
2988 if (bounded
&& value_eq(min
, max
)) {
2989 replace_by_affine(e
, min
);
2991 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
2992 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
2993 * See pages 199-200 of PhD thesis.
3001 value_set_si(rem
.d
, 0);
3002 rem
.x
.p
= new_enode(fractional
, 3, -1);
3003 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
3004 value_clear(rem
.x
.p
->arr
[1].d
);
3005 value_clear(rem
.x
.p
->arr
[2].d
);
3006 rem
.x
.p
->arr
[1] = p
->arr
[1];
3007 rem
.x
.p
->arr
[2] = p
->arr
[2];
3008 for (i
= 3; i
< p
->size
; ++i
)
3009 p
->arr
[i
-2] = p
->arr
[i
];
3013 value_init(inc
.x
.n
);
3014 value_set_si(inc
.d
, 1);
3015 value_oppose(inc
.x
.n
, min
);
3018 evalue_copy(&t
, &p
->arr
[0]);
3022 value_set_si(f
.d
, 0);
3023 f
.x
.p
= new_enode(fractional
, 3, -1);
3024 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3025 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3026 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
3028 value_init(factor
.d
);
3029 evalue_set_si(&factor
, -1, 1);
3035 value_clear(f
.x
.p
->arr
[1].x
.n
);
3036 value_clear(f
.x
.p
->arr
[2].x
.n
);
3037 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3038 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3042 reorder_terms(&rem
);
3049 free_evalue_refs(&inc
);
3050 free_evalue_refs(&t
);
3051 free_evalue_refs(&f
);
3052 free_evalue_refs(&factor
);
3053 free_evalue_refs(&rem
);
3055 evalue_range_reduction_in_domain(e
, D
);
3059 _reduce_evalue(&p
->arr
[0], 0, 1);
3071 void evalue_range_reduction(evalue
*e
)
3074 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3077 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3078 if (evalue_range_reduction_in_domain(&e
->x
.p
->arr
[2*i
+1],
3079 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
3080 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3082 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
3083 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
3084 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3085 value_clear(e
->x
.p
->arr
[2*i
].d
);
3087 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
3088 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
3096 Enumeration
* partition2enumeration(evalue
*EP
)
3099 Enumeration
*en
, *res
= NULL
;
3101 if (EVALUE_IS_ZERO(*EP
)) {
3106 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
3107 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
3108 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
3111 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
3112 value_clear(EP
->x
.p
->arr
[2*i
].d
);
3113 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
3121 int evalue_frac2floor_in_domain3(evalue
*e
, Polyhedron
*D
, int shift
)
3130 if (value_notzero_p(e
->d
))
3135 i
= p
->type
== relation
? 1 :
3136 p
->type
== fractional
? 1 : 0;
3137 for (; i
<p
->size
; i
++)
3138 r
|= evalue_frac2floor_in_domain3(&p
->arr
[i
], D
, shift
);
3140 if (p
->type
!= fractional
) {
3141 if (r
&& p
->type
== polynomial
) {
3144 value_set_si(f
.d
, 0);
3145 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3146 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3147 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3148 reorder_terms_about(p
, &f
);
3158 I
= polynomial_projection(p
, D
, &d
, NULL
);
3161 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3164 assert(I
->NbEq
== 0); /* Should have been reduced */
3167 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3168 if (value_pos_p(I
->Constraint
[i
][1]))
3171 if (i
< I
->NbConstraints
) {
3173 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3174 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3175 if (value_neg_p(min
)) {
3177 mpz_fdiv_q(min
, min
, d
);
3178 value_init(offset
.d
);
3179 value_set_si(offset
.d
, 1);
3180 value_init(offset
.x
.n
);
3181 value_oppose(offset
.x
.n
, min
);
3182 eadd(&offset
, &p
->arr
[0]);
3183 free_evalue_refs(&offset
);
3193 value_set_si(fl
.d
, 0);
3194 fl
.x
.p
= new_enode(flooring
, 3, -1);
3195 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
3196 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
3197 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
3199 eadd(&fl
, &p
->arr
[0]);
3200 reorder_terms_about(p
, &p
->arr
[0]);
3204 free_evalue_refs(&fl
);
3209 int evalue_frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
3211 return evalue_frac2floor_in_domain3(e
, D
, 1);
3214 void evalue_frac2floor2(evalue
*e
, int shift
)
3217 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
3219 if (evalue_frac2floor_in_domain3(e
, NULL
, 0))
3225 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3226 if (evalue_frac2floor_in_domain3(&e
->x
.p
->arr
[2*i
+1],
3227 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), shift
))
3228 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3231 void evalue_frac2floor(evalue
*e
)
3233 evalue_frac2floor2(e
, 1);
3236 /* Add a new variable with lower bound 1 and upper bound specified
3237 * by row. If negative is true, then the new variable has upper
3238 * bound -1 and lower bound specified by row.
3240 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
3241 Vector
*row
, int negative
)
3245 int nparam
= D
->Dimension
- nvar
;
3248 nr
= D
->NbConstraints
+ 2;
3249 nc
= D
->Dimension
+ 2 + 1;
3250 C
= Matrix_Alloc(nr
, nc
);
3251 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
3252 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
3253 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3254 D
->Dimension
+ 1 - nvar
);
3259 nc
= C
->NbColumns
+ 1;
3260 C
= Matrix_Alloc(nr
, nc
);
3261 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
3262 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
3263 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3264 oldC
->NbColumns
- 1 - nvar
);
3267 value_set_si(C
->p
[nr
-2][0], 1);
3269 value_set_si(C
->p
[nr
-2][1 + nvar
], -1);
3271 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
3272 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
3274 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
3275 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
3281 static void floor2frac_r(evalue
*e
, int nvar
)
3288 if (value_notzero_p(e
->d
))
3293 assert(p
->type
== flooring
);
3294 for (i
= 1; i
< p
->size
; i
++)
3295 floor2frac_r(&p
->arr
[i
], nvar
);
3297 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
3298 assert(pp
->x
.p
->type
== polynomial
);
3299 pp
->x
.p
->pos
-= nvar
;
3303 value_set_si(f
.d
, 0);
3304 f
.x
.p
= new_enode(fractional
, 3, -1);
3305 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3306 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3307 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3309 eadd(&f
, &p
->arr
[0]);
3310 reorder_terms_about(p
, &p
->arr
[0]);
3314 free_evalue_refs(&f
);
3317 /* Convert flooring back to fractional and shift position
3318 * of the parameters by nvar
3320 static void floor2frac(evalue
*e
, int nvar
)
3322 floor2frac_r(e
, nvar
);
3326 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
3329 int nparam
= D
->Dimension
- nvar
;
3333 D
= Constraints2Polyhedron(C
, 0);
3337 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
3339 /* Double check that D was not unbounded. */
3340 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
3348 static evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
3349 int *signs
, Matrix
*C
, unsigned MaxRays
)
3355 evalue
*factor
= NULL
;
3359 if (EVALUE_IS_ZERO(*e
))
3363 Polyhedron
*DD
= Disjoint_Domain(D
, 0, MaxRays
);
3370 res
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3373 for (Q
= DD
; Q
; Q
= DD
) {
3379 t
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3386 free_evalue_refs(t
);
3393 if (value_notzero_p(e
->d
)) {
3396 t
= esum_over_domain_cst(nvar
, D
, C
);
3398 if (!EVALUE_IS_ONE(*e
))
3404 switch (e
->x
.p
->type
) {
3406 evalue
*pp
= &e
->x
.p
->arr
[0];
3408 if (pp
->x
.p
->pos
> nvar
) {
3409 /* remainder is independent of the summated vars */
3415 floor2frac(&f
, nvar
);
3417 t
= esum_over_domain_cst(nvar
, D
, C
);
3421 free_evalue_refs(&f
);
3426 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3427 poly_denom(pp
, &row
->p
[1 + nvar
]);
3428 value_set_si(row
->p
[0], 1);
3429 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3430 pp
= &pp
->x
.p
->arr
[0]) {
3432 assert(pp
->x
.p
->type
== polynomial
);
3434 if (pos
>= 1 + nvar
)
3436 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3437 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3438 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3440 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3441 value_division(row
->p
[1 + D
->Dimension
+ 1],
3442 row
->p
[1 + D
->Dimension
+ 1],
3444 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3445 row
->p
[1 + D
->Dimension
+ 1],
3447 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3451 int pos
= e
->x
.p
->pos
;
3454 factor
= ALLOC(evalue
);
3455 value_init(factor
->d
);
3456 value_set_si(factor
->d
, 0);
3457 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3458 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3459 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3463 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3464 negative
= signs
[pos
-1] < 0;
3465 value_set_si(row
->p
[0], 1);
3467 value_set_si(row
->p
[pos
], -1);
3468 value_set_si(row
->p
[1 + nvar
], 1);
3470 value_set_si(row
->p
[pos
], 1);
3471 value_set_si(row
->p
[1 + nvar
], -1);
3479 offset
= type_offset(e
->x
.p
);
3481 res
= esum_over_domain(&e
->x
.p
->arr
[offset
], nvar
, D
, signs
, C
, MaxRays
);
3485 evalue_copy(&cum
, factor
);
3489 for (i
= 1; offset
+i
< e
->x
.p
->size
; ++i
) {
3493 C
= esum_add_constraint(nvar
, D
, C
, row
, negative
);
3499 Vector_Print(stderr, P_VALUE_FMT, row);
3501 Matrix_Print(stderr, P_VALUE_FMT, C);
3503 t
= esum_over_domain(&e
->x
.p
->arr
[offset
+i
], nvar
, D
, signs
, C
, MaxRays
);
3508 if (negative
&& (i
% 2))
3516 free_evalue_refs(t
);
3519 if (factor
&& offset
+i
+1 < e
->x
.p
->size
)
3526 free_evalue_refs(factor
);
3527 free_evalue_refs(&cum
);
3539 static void domain_signs(Polyhedron
*D
, int *signs
)
3543 POL_ENSURE_VERTICES(D
);
3544 for (j
= 0; j
< D
->Dimension
; ++j
) {
3546 for (k
= 0; k
< D
->NbRays
; ++k
) {
3547 signs
[j
] = value_sign(D
->Ray
[k
][1+j
]);
3554 static void shift_floor_in_domain(evalue
*e
, Polyhedron
*D
)
3561 if (value_notzero_p(e
->d
))
3566 for (i
= type_offset(p
); i
< p
->size
; ++i
)
3567 shift_floor_in_domain(&p
->arr
[i
], D
);
3569 if (p
->type
!= flooring
)
3575 I
= polynomial_projection(p
, D
, &d
, NULL
);
3576 assert(I
->NbEq
== 0); /* Should have been reduced */
3578 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3579 if (value_pos_p(I
->Constraint
[i
][1]))
3581 assert(i
< I
->NbConstraints
);
3582 if (i
< I
->NbConstraints
) {
3583 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3584 mpz_fdiv_q(m
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3585 if (value_neg_p(m
)) {
3586 /* replace [e] by [e-m]+m such that e-m >= 0 */
3591 value_set_si(f
.d
, 1);
3592 value_oppose(f
.x
.n
, m
);
3593 eadd(&f
, &p
->arr
[0]);
3596 value_set_si(f
.d
, 0);
3597 f
.x
.p
= new_enode(flooring
, 3, -1);
3598 value_clear(f
.x
.p
->arr
[0].d
);
3599 f
.x
.p
->arr
[0] = p
->arr
[0];
3600 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
3601 value_set_si(f
.x
.p
->arr
[1].d
, 1);
3602 value_init(f
.x
.p
->arr
[1].x
.n
);
3603 value_assign(f
.x
.p
->arr
[1].x
.n
, m
);
3604 reorder_terms_about(p
, &f
);
3615 /* Make arguments of all floors non-negative */
3616 static void shift_floor_arguments(evalue
*e
)
3620 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3623 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3624 shift_floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
3625 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3628 evalue
*evalue_sum(evalue
*e
, int nvar
, unsigned MaxRays
)
3632 evalue
*res
= ALLOC(evalue
);
3636 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3637 evalue_copy(res
, e
);
3641 evalue_split_domains_into_orthants(e
, MaxRays
);
3642 evalue_frac2floor2(e
, 0);
3643 evalue_set_si(res
, 0, 1);
3645 assert(value_zero_p(e
->d
));
3646 assert(e
->x
.p
->type
== partition
);
3647 shift_floor_arguments(e
);
3649 assert(e
->x
.p
->size
>= 2);
3650 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3652 signs
= alloca(sizeof(int) * dim
);
3654 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3656 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3657 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3658 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
, 0,
3661 free_evalue_refs(t
);
3670 evalue
*esum(evalue
*e
, int nvar
)
3672 return evalue_sum(e
, nvar
, 0);
3675 /* Initial silly implementation */
3676 void eor(evalue
*e1
, evalue
*res
)
3682 evalue_set_si(&mone
, -1, 1);
3684 evalue_copy(&E
, res
);
3690 free_evalue_refs(&E
);
3691 free_evalue_refs(&mone
);
3694 /* computes denominator of polynomial evalue
3695 * d should point to a value initialized to 1
3697 void evalue_denom(const evalue
*e
, Value
*d
)
3701 if (value_notzero_p(e
->d
)) {
3702 value_lcm(*d
, e
->d
, d
);
3705 assert(e
->x
.p
->type
== polynomial
||
3706 e
->x
.p
->type
== fractional
||
3707 e
->x
.p
->type
== flooring
);
3708 offset
= type_offset(e
->x
.p
);
3709 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3710 evalue_denom(&e
->x
.p
->arr
[i
], d
);
3713 /* Divides the evalue e by the integer n */
3714 void evalue_div(evalue
*e
, Value n
)
3718 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3721 if (value_notzero_p(e
->d
)) {
3724 value_multiply(e
->d
, e
->d
, n
);
3725 Gcd(e
->x
.n
, e
->d
, &gc
);
3726 if (value_notone_p(gc
)) {
3727 value_division(e
->d
, e
->d
, gc
);
3728 value_division(e
->x
.n
, e
->x
.n
, gc
);
3733 if (e
->x
.p
->type
== partition
) {
3734 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3735 evalue_div(&e
->x
.p
->arr
[2*i
+1], n
);
3738 offset
= type_offset(e
->x
.p
);
3739 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3740 evalue_div(&e
->x
.p
->arr
[i
], n
);
3743 /* Multiplies the evalue e by the integer n */
3744 void evalue_mul(evalue
*e
, Value n
)
3748 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3751 if (value_notzero_p(e
->d
)) {
3754 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3755 Gcd(e
->x
.n
, e
->d
, &gc
);
3756 if (value_notone_p(gc
)) {
3757 value_division(e
->d
, e
->d
, gc
);
3758 value_division(e
->x
.n
, e
->x
.n
, gc
);
3763 if (e
->x
.p
->type
== partition
) {
3764 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3765 evalue_mul(&e
->x
.p
->arr
[2*i
+1], n
);
3768 offset
= type_offset(e
->x
.p
);
3769 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3770 evalue_mul(&e
->x
.p
->arr
[i
], n
);
3773 /* Multiplies the evalue e by the n/d */
3774 void evalue_mul_div(evalue
*e
, Value n
, Value d
)
3778 if ((value_one_p(n
) && value_one_p(d
)) || EVALUE_IS_ZERO(*e
))
3781 if (value_notzero_p(e
->d
)) {
3784 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3785 value_multiply(e
->d
, e
->d
, d
);
3786 Gcd(e
->x
.n
, e
->d
, &gc
);
3787 if (value_notone_p(gc
)) {
3788 value_division(e
->d
, e
->d
, gc
);
3789 value_division(e
->x
.n
, e
->x
.n
, gc
);
3794 if (e
->x
.p
->type
== partition
) {
3795 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3796 evalue_mul_div(&e
->x
.p
->arr
[2*i
+1], n
, d
);
3799 offset
= type_offset(e
->x
.p
);
3800 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3801 evalue_mul_div(&e
->x
.p
->arr
[i
], n
, d
);
3804 void evalue_negate(evalue
*e
)
3808 if (value_notzero_p(e
->d
)) {
3809 value_oppose(e
->x
.n
, e
->x
.n
);
3812 if (e
->x
.p
->type
== partition
) {
3813 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3814 evalue_negate(&e
->x
.p
->arr
[2*i
+1]);
3817 offset
= type_offset(e
->x
.p
);
3818 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3819 evalue_negate(&e
->x
.p
->arr
[i
]);
3822 void evalue_add_constant(evalue
*e
, const Value cst
)
3826 if (value_zero_p(e
->d
)) {
3827 if (e
->x
.p
->type
== partition
) {
3828 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3829 evalue_add_constant(&e
->x
.p
->arr
[2*i
+1], cst
);
3832 if (e
->x
.p
->type
== relation
) {
3833 for (i
= 1; i
< e
->x
.p
->size
; ++i
)
3834 evalue_add_constant(&e
->x
.p
->arr
[i
], cst
);
3838 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
3839 } while (value_zero_p(e
->d
));
3841 value_addmul(e
->x
.n
, cst
, e
->d
);
3844 static void evalue_frac2polynomial_r(evalue
*e
, int *signs
, int sign
, int in_frac
)
3849 int sign_odd
= sign
;
3851 if (value_notzero_p(e
->d
)) {
3852 if (in_frac
&& sign
* value_sign(e
->x
.n
) < 0) {
3853 value_set_si(e
->x
.n
, 0);
3854 value_set_si(e
->d
, 1);
3859 if (e
->x
.p
->type
== relation
) {
3860 for (i
= e
->x
.p
->size
-1; i
>= 1; --i
)
3861 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
, sign
, in_frac
);
3865 if (e
->x
.p
->type
== polynomial
)
3866 sign_odd
*= signs
[e
->x
.p
->pos
-1];
3867 offset
= type_offset(e
->x
.p
);
3868 evalue_frac2polynomial_r(&e
->x
.p
->arr
[offset
], signs
, sign
, in_frac
);
3869 in_frac
|= e
->x
.p
->type
== fractional
;
3870 for (i
= e
->x
.p
->size
-1; i
> offset
; --i
)
3871 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
,
3872 (i
- offset
) % 2 ? sign_odd
: sign
, in_frac
);
3874 if (e
->x
.p
->type
!= fractional
)
3877 /* replace { a/m } by (m-1)/m if sign != 0
3878 * and by (m-1)/(2m) if sign == 0
3882 evalue_denom(&e
->x
.p
->arr
[0], &d
);
3883 free_evalue_refs(&e
->x
.p
->arr
[0]);
3884 value_init(e
->x
.p
->arr
[0].d
);
3885 value_init(e
->x
.p
->arr
[0].x
.n
);
3887 value_addto(e
->x
.p
->arr
[0].d
, d
, d
);
3889 value_assign(e
->x
.p
->arr
[0].d
, d
);
3890 value_decrement(e
->x
.p
->arr
[0].x
.n
, d
);
3894 reorder_terms_about(p
, &p
->arr
[0]);
3900 /* Approximate the evalue in fractional representation by a polynomial.
3901 * If sign > 0, the result is an upper bound;
3902 * if sign < 0, the result is a lower bound;
3903 * if sign = 0, the result is an intermediate approximation.
3905 void evalue_frac2polynomial(evalue
*e
, int sign
, unsigned MaxRays
)
3910 if (value_notzero_p(e
->d
))
3912 assert(e
->x
.p
->type
== partition
);
3913 /* make sure all variables in the domains have a fixed sign */
3915 evalue_split_domains_into_orthants(e
, MaxRays
);
3916 if (EVALUE_IS_ZERO(*e
))
3920 assert(e
->x
.p
->size
>= 2);
3921 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3923 signs
= alloca(sizeof(int) * dim
);
3926 for (i
= 0; i
< dim
; ++i
)
3928 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3930 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3931 evalue_frac2polynomial_r(&e
->x
.p
->arr
[2*i
+1], signs
, sign
, 0);
3935 /* Split the domains of e (which is assumed to be a partition)
3936 * such that each resulting domain lies entirely in one orthant.
3938 void evalue_split_domains_into_orthants(evalue
*e
, unsigned MaxRays
)
3941 assert(value_zero_p(e
->d
));
3942 assert(e
->x
.p
->type
== partition
);
3943 assert(e
->x
.p
->size
>= 2);
3944 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3946 for (i
= 0; i
< dim
; ++i
) {
3949 C
= Matrix_Alloc(1, 1 + dim
+ 1);
3950 value_set_si(C
->p
[0][0], 1);
3951 value_init(split
.d
);
3952 value_set_si(split
.d
, 0);
3953 split
.x
.p
= new_enode(partition
, 4, dim
);
3954 value_set_si(C
->p
[0][1+i
], 1);
3955 C2
= Matrix_Copy(C
);
3956 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[0], Constraints2Polyhedron(C2
, MaxRays
));
3958 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
3959 value_set_si(C
->p
[0][1+i
], -1);
3960 value_set_si(C
->p
[0][1+dim
], -1);
3961 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[2], Constraints2Polyhedron(C
, MaxRays
));
3962 evalue_set_si(&split
.x
.p
->arr
[3], 1, 1);
3964 free_evalue_refs(&split
);
3970 static evalue
*find_fractional_with_max_periods(evalue
*e
, Polyhedron
*D
,
3973 Value
*min
, Value
*max
)
3980 if (value_notzero_p(e
->d
))
3983 if (e
->x
.p
->type
== fractional
) {
3988 I
= polynomial_projection(e
->x
.p
, D
, &d
, &T
);
3989 bounded
= line_minmax(I
, min
, max
); /* frees I */
3993 value_set_si(mp
, max_periods
);
3994 mpz_fdiv_q(*min
, *min
, d
);
3995 mpz_fdiv_q(*max
, *max
, d
);
3996 value_assign(T
->p
[1][D
->Dimension
], d
);
3997 value_subtract(d
, *max
, *min
);
3998 if (value_ge(d
, mp
))
4001 f
= evalue_dup(&e
->x
.p
->arr
[0]);
4012 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
4013 if ((f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[i
], D
, max_periods
,
4020 static void replace_fract_by_affine(evalue
*e
, evalue
*f
, Value val
)
4024 if (value_notzero_p(e
->d
))
4027 offset
= type_offset(e
->x
.p
);
4028 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
4029 replace_fract_by_affine(&e
->x
.p
->arr
[i
], f
, val
);
4031 if (e
->x
.p
->type
!= fractional
)
4034 if (!eequal(&e
->x
.p
->arr
[0], f
))
4037 replace_by_affine(e
, val
);
4040 /* Look for fractional parts that can be removed by splitting the corresponding
4041 * domain into at most max_periods parts.
4042 * We use a very simply strategy that looks for the first fractional part
4043 * that satisfies the condition, performs the split and then continues
4044 * looking for other fractional parts in the split domains until no
4045 * such fractional part can be found anymore.
4047 void evalue_split_periods(evalue
*e
, int max_periods
, unsigned int MaxRays
)
4054 if (EVALUE_IS_ZERO(*e
))
4056 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
4058 "WARNING: evalue_split_periods called on incorrect evalue type\n");
4066 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
4071 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
4073 f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[2*i
+1], D
, max_periods
,
4078 M
= Matrix_Alloc(2, 2+D
->Dimension
);
4080 value_subtract(d
, max
, min
);
4081 n
= VALUE_TO_INT(d
)+1;
4083 value_set_si(M
->p
[0][0], 1);
4084 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
4085 value_multiply(d
, max
, T
->p
[1][D
->Dimension
]);
4086 value_subtract(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
], d
);
4087 value_set_si(d
, -1);
4088 value_set_si(M
->p
[1][0], 1);
4089 Vector_Scale(T
->p
[0], M
->p
[1]+1, d
, D
->Dimension
+1);
4090 value_addmul(M
->p
[1][1+D
->Dimension
], max
, T
->p
[1][D
->Dimension
]);
4091 value_addto(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4092 T
->p
[1][D
->Dimension
]);
4093 value_decrement(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
]);
4095 p
= new_enode(partition
, e
->x
.p
->size
+ (n
-1)*2, e
->x
.p
->pos
);
4096 for (j
= 0; j
< 2*i
; ++j
) {
4097 value_clear(p
->arr
[j
].d
);
4098 p
->arr
[j
] = e
->x
.p
->arr
[j
];
4100 for (j
= 2*i
+2; j
< e
->x
.p
->size
; ++j
) {
4101 value_clear(p
->arr
[j
+2*(n
-1)].d
);
4102 p
->arr
[j
+2*(n
-1)] = e
->x
.p
->arr
[j
];
4104 for (j
= n
-1; j
>= 0; --j
) {
4106 value_clear(p
->arr
[2*i
+1].d
);
4107 p
->arr
[2*i
+1] = e
->x
.p
->arr
[2*i
+1];
4109 evalue_copy(&p
->arr
[2*(i
+j
)+1], &e
->x
.p
->arr
[2*i
+1]);
4111 value_subtract(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4112 T
->p
[1][D
->Dimension
]);
4113 value_addto(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
],
4114 T
->p
[1][D
->Dimension
]);
4116 replace_fract_by_affine(&p
->arr
[2*(i
+j
)+1], f
, max
);
4117 E
= DomainAddConstraints(D
, M
, MaxRays
);
4118 EVALUE_SET_DOMAIN(p
->arr
[2*(i
+j
)], E
);
4119 if (evalue_range_reduction_in_domain(&p
->arr
[2*(i
+j
)+1], E
))
4120 reduce_evalue(&p
->arr
[2*(i
+j
)+1]);
4121 value_decrement(max
, max
);
4123 value_clear(e
->x
.p
->arr
[2*i
].d
);
4127 free_evalue_refs(f
);
4139 void evalue_extract_affine(const evalue
*e
, Value
*coeff
, Value
*cst
, Value
*d
)
4141 value_set_si(*d
, 1);
4143 for ( ; value_zero_p(e
->d
); e
= &e
->x
.p
->arr
[0]) {
4145 assert(e
->x
.p
->type
== polynomial
);
4146 assert(e
->x
.p
->size
== 2);
4147 c
= &e
->x
.p
->arr
[1];
4148 value_multiply(coeff
[e
->x
.p
->pos
-1], *d
, c
->x
.n
);
4149 value_division(coeff
[e
->x
.p
->pos
-1], coeff
[e
->x
.p
->pos
-1], c
->d
);
4151 value_multiply(*cst
, *d
, e
->x
.n
);
4152 value_division(*cst
, *cst
, e
->d
);
4155 /* returns an evalue that corresponds to
4159 static evalue
*term(int param
, Value c
, Value den
)
4161 evalue
*EP
= ALLOC(evalue
);
4163 value_set_si(EP
->d
,0);
4164 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
4165 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
4166 value_init(EP
->x
.p
->arr
[1].x
.n
);
4167 value_assign(EP
->x
.p
->arr
[1].d
, den
);
4168 value_assign(EP
->x
.p
->arr
[1].x
.n
, c
);
4172 evalue
*affine2evalue(Value
*coeff
, Value denom
, int nvar
)
4175 evalue
*E
= ALLOC(evalue
);
4177 evalue_set(E
, coeff
[nvar
], denom
);
4178 for (i
= 0; i
< nvar
; ++i
) {
4180 if (value_zero_p(coeff
[i
]))
4182 t
= term(i
, coeff
[i
], denom
);
4184 free_evalue_refs(t
);
4190 void evalue_substitute(evalue
*e
, evalue
**subs
)
4196 if (value_notzero_p(e
->d
))
4200 assert(p
->type
!= partition
);
4202 for (i
= 0; i
< p
->size
; ++i
)
4203 evalue_substitute(&p
->arr
[i
], subs
);
4205 if (p
->type
== polynomial
)
4210 value_set_si(v
->d
, 0);
4211 v
->x
.p
= new_enode(p
->type
, 3, -1);
4212 value_clear(v
->x
.p
->arr
[0].d
);
4213 v
->x
.p
->arr
[0] = p
->arr
[0];
4214 evalue_set_si(&v
->x
.p
->arr
[1], 0, 1);
4215 evalue_set_si(&v
->x
.p
->arr
[2], 1, 1);
4218 offset
= type_offset(p
);
4220 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
4221 emul(v
, &p
->arr
[i
]);
4222 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
4223 free_evalue_refs(&(p
->arr
[i
]));
4226 if (p
->type
!= polynomial
) {
4227 free_evalue_refs(v
);
4232 *e
= p
->arr
[offset
];
4236 /* evalue e is given in terms of "new" parameter; CP maps the new
4237 * parameters back to the old parameters.
4238 * Transforms e such that it refers back to the old parameters.
4240 void evalue_backsubstitute(evalue
*e
, Matrix
*CP
, unsigned MaxRays
)
4247 unsigned nparam
= CP
->NbColumns
-1;
4250 if (EVALUE_IS_ZERO(*e
))
4253 assert(value_zero_p(e
->d
));
4255 assert(p
->type
== partition
);
4257 inv
= left_inverse(CP
, &eq
);
4258 subs
= ALLOCN(evalue
*, nparam
);
4259 for (i
= 0; i
< nparam
; ++i
)
4260 subs
[i
] = affine2evalue(inv
->p
[i
], inv
->p
[nparam
][inv
->NbColumns
-1],
4263 CEq
= Constraints2Polyhedron(eq
, MaxRays
);
4264 addeliminatedparams_partition(p
, inv
, CEq
, inv
->NbColumns
-1, MaxRays
);
4265 Polyhedron_Free(CEq
);
4267 for (i
= 0; i
< p
->size
/2; ++i
)
4268 evalue_substitute(&p
->arr
[2*i
+1], subs
);
4270 for (i
= 0; i
< nparam
; ++i
) {
4271 free_evalue_refs(subs
[i
]);
4281 * \sum_{i=0}^n c_i/d X^i
4283 * where d is the last element in the vector c.
4285 evalue
*evalue_polynomial(Vector
*c
, const evalue
* X
)
4287 unsigned dim
= c
->Size
-2;
4289 evalue
*EP
= ALLOC(evalue
);
4293 evalue_set(&EC
, c
->p
[dim
], c
->p
[dim
+1]);
4296 evalue_set(EP
, c
->p
[dim
], c
->p
[dim
+1]);
4298 for (i
= dim
-1; i
>= 0; --i
) {
4300 value_assign(EC
.x
.n
, c
->p
[i
]);
4303 free_evalue_refs(&EC
);
4307 /* Create an evalue from an array of pairs of domains and evalues. */
4308 evalue
*evalue_from_section_array(struct evalue_section
*s
, int n
)
4313 res
= ALLOC(evalue
);
4317 evalue_set_si(res
, 0, 1);
4319 value_set_si(res
->d
, 0);
4320 res
->x
.p
= new_enode(partition
, 2*n
, s
[0].D
->Dimension
);
4321 for (i
= 0; i
< n
; ++i
) {
4322 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
], s
[i
].D
);
4323 value_clear(res
->x
.p
->arr
[2*i
+1].d
);
4324 res
->x
.p
->arr
[2*i
+1] = *s
[i
].E
;