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
, char **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
, char **pname
)
983 print_evalue_r(DST
, e
, pname
);
984 if (value_notzero_p(e
->d
))
988 void print_enode(FILE *DST
,enode
*p
,char **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 **names
= pname
;
1056 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
1057 if (!pname
|| p
->pos
< maxdim
) {
1058 NALLOC(names
, maxdim
);
1059 for (i
= 0; i
< p
->pos
; ++i
) {
1061 names
[i
] = pname
[i
];
1063 NALLOC(names
[i
], 10);
1064 snprintf(names
[i
], 10, "%c", 'P'+i
);
1067 for ( ; i
< maxdim
; ++i
) {
1068 NALLOC(names
[i
], 10);
1069 snprintf(names
[i
], 10, "_p%d", i
);
1073 for (i
=0; i
<p
->size
/2; i
++) {
1074 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
1075 print_evalue_r(DST
, &p
->arr
[2*i
+1], names
);
1078 if (!pname
|| p
->pos
< maxdim
) {
1079 for (i
= pname
? p
->pos
: 0; i
< maxdim
; ++i
)
1092 static void eadd_rev(const evalue
*e1
, evalue
*res
)
1096 evalue_copy(&ev
, e1
);
1098 free_evalue_refs(res
);
1102 static void eadd_rev_cst(const evalue
*e1
, evalue
*res
)
1106 evalue_copy(&ev
, e1
);
1107 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
1108 free_evalue_refs(res
);
1112 static int is_zero_on(evalue
*e
, Polyhedron
*D
)
1117 tmp
.x
.p
= new_enode(partition
, 2, D
->Dimension
);
1118 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Domain_Copy(D
));
1119 evalue_copy(&tmp
.x
.p
->arr
[1], e
);
1120 reduce_evalue(&tmp
);
1121 is_zero
= EVALUE_IS_ZERO(tmp
);
1122 free_evalue_refs(&tmp
);
1126 struct section
{ Polyhedron
* D
; evalue E
; };
1128 void eadd_partitions(const evalue
*e1
, evalue
*res
)
1133 s
= (struct section
*)
1134 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
1135 sizeof(struct section
));
1137 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1138 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1139 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1142 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1143 assert(res
->x
.p
->size
>= 2);
1144 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1145 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
1147 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
1149 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1158 /* See if we can extend one of the domains in res to cover fd */
1159 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1160 if (is_zero_on(&res
->x
.p
->arr
[2*i
+1], fd
))
1162 if (i
< res
->x
.p
->size
/2) {
1163 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
1164 DomainConcat(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
])));
1167 value_init(s
[n
].E
.d
);
1168 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
1172 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1173 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
1174 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1176 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1177 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1183 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
1184 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1186 value_init(s
[n
].E
.d
);
1187 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1188 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1189 if (!emptyQ(fd
) && is_zero_on(&e1
->x
.p
->arr
[2*j
+1], fd
)) {
1190 d
= DomainConcat(fd
, d
);
1191 fd
= Empty_Polyhedron(fd
->Dimension
);
1197 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
1201 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1204 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1205 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1206 value_clear(res
->x
.p
->arr
[2*i
].d
);
1211 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1212 for (j
= 0; j
< n
; ++j
) {
1213 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1214 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1215 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1216 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1222 static void explicit_complement(evalue
*res
)
1224 enode
*rel
= new_enode(relation
, 3, 0);
1226 value_clear(rel
->arr
[0].d
);
1227 rel
->arr
[0] = res
->x
.p
->arr
[0];
1228 value_clear(rel
->arr
[1].d
);
1229 rel
->arr
[1] = res
->x
.p
->arr
[1];
1230 value_set_si(rel
->arr
[2].d
, 1);
1231 value_init(rel
->arr
[2].x
.n
);
1232 value_set_si(rel
->arr
[2].x
.n
, 0);
1237 void eadd(const evalue
*e1
, evalue
*res
)
1240 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1241 /* Add two rational numbers */
1247 value_multiply(m1
,e1
->x
.n
,res
->d
);
1248 value_multiply(m2
,res
->x
.n
,e1
->d
);
1249 value_addto(res
->x
.n
,m1
,m2
);
1250 value_multiply(res
->d
,e1
->d
,res
->d
);
1251 Gcd(res
->x
.n
,res
->d
,&g
);
1252 if (value_notone_p(g
)) {
1253 value_division(res
->d
,res
->d
,g
);
1254 value_division(res
->x
.n
,res
->x
.n
,g
);
1256 value_clear(g
); value_clear(m1
); value_clear(m2
);
1259 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
1260 switch (res
->x
.p
->type
) {
1262 /* Add the constant to the constant term of a polynomial*/
1263 eadd(e1
, &res
->x
.p
->arr
[0]);
1266 /* Add the constant to all elements of a periodic number */
1267 for (i
=0; i
<res
->x
.p
->size
; i
++) {
1268 eadd(e1
, &res
->x
.p
->arr
[i
]);
1272 fprintf(stderr
, "eadd: cannot add const with vector\n");
1276 eadd(e1
, &res
->x
.p
->arr
[1]);
1279 assert(EVALUE_IS_ZERO(*e1
));
1280 break; /* Do nothing */
1282 /* Create (zero) complement if needed */
1283 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
1284 explicit_complement(res
);
1285 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1286 eadd(e1
, &res
->x
.p
->arr
[i
]);
1292 /* add polynomial or periodic to constant
1293 * you have to exchange e1 and res, before doing addition */
1295 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1299 else { // ((e1->d==0) && (res->d==0))
1300 assert(!((e1
->x
.p
->type
== partition
) ^
1301 (res
->x
.p
->type
== partition
)));
1302 if (e1
->x
.p
->type
== partition
) {
1303 eadd_partitions(e1
, res
);
1306 if (e1
->x
.p
->type
== relation
&&
1307 (res
->x
.p
->type
!= relation
||
1308 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1312 if (res
->x
.p
->type
== relation
) {
1313 if (e1
->x
.p
->type
== relation
&&
1314 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1315 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1316 explicit_complement(res
);
1317 for (i
= 1; i
< e1
->x
.p
->size
; ++i
)
1318 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1321 if (res
->x
.p
->size
< 3)
1322 explicit_complement(res
);
1323 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1324 eadd(e1
, &res
->x
.p
->arr
[i
]);
1327 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
1328 /* adding to evalues of different type. two cases are possible
1329 * res is periodic and e1 is polynomial, you have to exchange
1330 * e1 and res then to add e1 to the constant term of res */
1331 if (e1
->x
.p
->type
== polynomial
) {
1332 eadd_rev_cst(e1
, res
);
1334 else if (res
->x
.p
->type
== polynomial
) {
1335 /* res is polynomial and e1 is periodic,
1336 add e1 to the constant term of res */
1338 eadd(e1
,&res
->x
.p
->arr
[0]);
1344 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
1345 ((res
->x
.p
->type
== fractional
||
1346 res
->x
.p
->type
== flooring
) &&
1347 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1348 /* adding evalues of different position (i.e function of different unknowns
1349 * to case are possible */
1351 switch (res
->x
.p
->type
) {
1354 if (mod_term_smaller(res
, e1
))
1355 eadd(e1
,&res
->x
.p
->arr
[1]);
1357 eadd_rev_cst(e1
, res
);
1359 case polynomial
: // res and e1 are polynomials
1360 // add e1 to the constant term of res
1362 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1363 eadd(e1
,&res
->x
.p
->arr
[0]);
1365 eadd_rev_cst(e1
, res
);
1366 // value_clear(g); value_clear(m1); value_clear(m2);
1368 case periodic
: // res and e1 are pointers to periodic numbers
1369 //add e1 to all elements of res
1371 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1372 for (i
=0;i
<res
->x
.p
->size
;i
++) {
1373 eadd(e1
,&res
->x
.p
->arr
[i
]);
1384 //same type , same pos and same size
1385 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1386 // add any element in e1 to the corresponding element in res
1387 i
= type_offset(res
->x
.p
);
1389 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1390 for (; i
<res
->x
.p
->size
; i
++) {
1391 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1396 /* Sizes are different */
1397 switch(res
->x
.p
->type
) {
1401 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1402 /* new enode and add res to that new node. If you do not do */
1403 /* that, you lose the the upper weight part of e1 ! */
1405 if(e1
->x
.p
->size
> res
->x
.p
->size
)
1408 i
= type_offset(res
->x
.p
);
1410 assert(eequal(&e1
->x
.p
->arr
[0],
1411 &res
->x
.p
->arr
[0]));
1412 for (; i
<e1
->x
.p
->size
; i
++) {
1413 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1420 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1423 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1424 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1425 to add periodicaly elements of e1 to elements of ne, and finaly to
1430 value_init(ex
); value_init(ey
);value_init(ep
);
1433 value_set_si(ex
,e1
->x
.p
->size
);
1434 value_set_si(ey
,res
->x
.p
->size
);
1435 value_assign (ep
,*Lcm(ex
,ey
));
1436 p
=(int)mpz_get_si(ep
);
1437 ne
= (evalue
*) malloc (sizeof(evalue
));
1439 value_set_si( ne
->d
,0);
1441 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1443 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1446 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1449 value_assign(res
->d
, ne
->d
);
1455 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
1464 static void emul_rev (evalue
*e1
, evalue
*res
)
1468 evalue_copy(&ev
, e1
);
1470 free_evalue_refs(res
);
1474 static void emul_poly(evalue
*e1
, evalue
*res
)
1476 int i
, j
, offset
= type_offset(res
->x
.p
);
1479 int size
= (e1
->x
.p
->size
+ res
->x
.p
->size
- offset
- 1);
1481 p
= new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1483 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1484 if (!EVALUE_IS_ZERO(e1
->x
.p
->arr
[i
]))
1487 /* special case pure power */
1488 if (i
== e1
->x
.p
->size
-1) {
1490 value_clear(p
->arr
[0].d
);
1491 p
->arr
[0] = res
->x
.p
->arr
[0];
1493 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1494 evalue_set_si(&p
->arr
[i
], 0, 1);
1495 for (i
= offset
; i
< res
->x
.p
->size
; ++i
) {
1496 value_clear(p
->arr
[i
+e1
->x
.p
->size
-offset
-1].d
);
1497 p
->arr
[i
+e1
->x
.p
->size
-offset
-1] = res
->x
.p
->arr
[i
];
1498 emul(&e1
->x
.p
->arr
[e1
->x
.p
->size
-1],
1499 &p
->arr
[i
+e1
->x
.p
->size
-offset
-1]);
1507 value_set_si(tmp
.d
,0);
1510 evalue_copy(&p
->arr
[0], &e1
->x
.p
->arr
[0]);
1511 for (i
= offset
; i
< e1
->x
.p
->size
; i
++) {
1512 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1513 emul(&res
->x
.p
->arr
[offset
], &tmp
.x
.p
->arr
[i
]);
1516 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1517 for (i
= offset
+1; i
<res
->x
.p
->size
; i
++)
1518 for (j
= offset
; j
<e1
->x
.p
->size
; j
++) {
1521 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1522 emul(&res
->x
.p
->arr
[i
], &ev
);
1523 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-offset
]);
1524 free_evalue_refs(&ev
);
1526 free_evalue_refs(res
);
1530 void emul_partitions (evalue
*e1
,evalue
*res
)
1535 s
= (struct section
*)
1536 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1537 sizeof(struct section
));
1539 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1540 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1541 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1544 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1545 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1546 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1547 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1553 /* This code is only needed because the partitions
1554 are not true partitions.
1556 for (k
= 0; k
< n
; ++k
) {
1557 if (DomainIncludes(s
[k
].D
, d
))
1559 if (DomainIncludes(d
, s
[k
].D
)) {
1560 Domain_Free(s
[k
].D
);
1561 free_evalue_refs(&s
[k
].E
);
1572 value_init(s
[n
].E
.d
);
1573 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1574 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1578 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1579 value_clear(res
->x
.p
->arr
[2*i
].d
);
1580 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1585 evalue_set_si(res
, 0, 1);
1587 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1588 for (j
= 0; j
< n
; ++j
) {
1589 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1590 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1591 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1592 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1599 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1601 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1602 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1603 * evalues is not treated here */
1605 void emul (evalue
*e1
, evalue
*res
){
1608 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1609 fprintf(stderr
, "emul: do not proced on evector type !\n");
1613 if (EVALUE_IS_ZERO(*res
))
1616 if (EVALUE_IS_ONE(*e1
))
1619 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1620 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1621 emul_partitions(e1
, res
);
1624 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1625 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1626 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1628 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1629 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1630 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1631 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1632 explicit_complement(res
);
1633 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1634 explicit_complement(e1
);
1635 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1636 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1639 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1640 emul(e1
, &res
->x
.p
->arr
[i
]);
1642 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1643 switch(e1
->x
.p
->type
) {
1645 switch(res
->x
.p
->type
) {
1647 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1648 /* Product of two polynomials of the same variable */
1653 /* Product of two polynomials of different variables */
1655 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1656 for( i
=0; i
<res
->x
.p
->size
; i
++)
1657 emul(e1
, &res
->x
.p
->arr
[i
]);
1666 /* Product of a polynomial and a periodic or fractional */
1673 switch(res
->x
.p
->type
) {
1675 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1676 /* 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
]));
1684 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1685 /* Product of two periodics of the same parameter and different periods */
1689 value_init(x
); value_init(y
);value_init(z
);
1692 value_set_si(x
,e1
->x
.p
->size
);
1693 value_set_si(y
,res
->x
.p
->size
);
1694 value_assign (z
,*Lcm(x
,y
));
1695 lcm
=(int)mpz_get_si(z
);
1696 newp
= (evalue
*) malloc (sizeof(evalue
));
1697 value_init(newp
->d
);
1698 value_set_si( newp
->d
,0);
1699 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1700 for(i
=0;i
<lcm
;i
++) {
1701 evalue_copy(&newp
->x
.p
->arr
[i
],
1702 &res
->x
.p
->arr
[i
%iy
]);
1705 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1707 value_assign(res
->d
,newp
->d
);
1710 value_clear(x
); value_clear(y
);value_clear(z
);
1714 /* Product of two periodics of different parameters */
1716 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1717 for(i
=0; i
<res
->x
.p
->size
; i
++)
1718 emul(e1
, &(res
->x
.p
->arr
[i
]));
1726 /* Product of a periodic and a polynomial */
1728 for(i
=0; i
<res
->x
.p
->size
; i
++)
1729 emul(e1
, &(res
->x
.p
->arr
[i
]));
1736 switch(res
->x
.p
->type
) {
1738 for(i
=0; i
<res
->x
.p
->size
; i
++)
1739 emul(e1
, &(res
->x
.p
->arr
[i
]));
1746 assert(e1
->x
.p
->type
== res
->x
.p
->type
);
1747 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1748 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1751 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1752 if (e1
->x
.p
->type
!= fractional
|| !value_two_p(d
.d
))
1757 value_set_si(d
.x
.n
, 1);
1758 /* { x }^2 == { x }/2 */
1759 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1760 assert(e1
->x
.p
->size
== 3);
1761 assert(res
->x
.p
->size
== 3);
1763 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1765 eadd(&res
->x
.p
->arr
[1], &tmp
);
1766 emul(&e1
->x
.p
->arr
[2], &tmp
);
1767 emul(&e1
->x
.p
->arr
[1], res
);
1768 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1769 free_evalue_refs(&tmp
);
1774 if(mod_term_smaller(res
, e1
))
1775 for(i
=1; i
<res
->x
.p
->size
; i
++)
1776 emul(e1
, &(res
->x
.p
->arr
[i
]));
1791 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1792 /* Product of two rational numbers */
1796 value_multiply(res
->d
,e1
->d
,res
->d
);
1797 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1798 Gcd(res
->x
.n
, res
->d
,&g
);
1799 if (value_notone_p(g
)) {
1800 value_division(res
->d
,res
->d
,g
);
1801 value_division(res
->x
.n
,res
->x
.n
,g
);
1807 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1808 /* Product of an expression (polynomial or peririodic) and a rational number */
1814 /* Product of a rationel number and an expression (polynomial or peririodic) */
1816 i
= type_offset(res
->x
.p
);
1817 for (; i
<res
->x
.p
->size
; i
++)
1818 emul(e1
, &res
->x
.p
->arr
[i
]);
1828 /* Frees mask content ! */
1829 void emask(evalue
*mask
, evalue
*res
) {
1836 if (EVALUE_IS_ZERO(*res
)) {
1837 free_evalue_refs(mask
);
1841 assert(value_zero_p(mask
->d
));
1842 assert(mask
->x
.p
->type
== partition
);
1843 assert(value_zero_p(res
->d
));
1844 assert(res
->x
.p
->type
== partition
);
1845 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1846 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1847 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1848 pos
= res
->x
.p
->pos
;
1850 s
= (struct section
*)
1851 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1852 sizeof(struct section
));
1856 evalue_set_si(&mone
, -1, 1);
1859 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1860 assert(mask
->x
.p
->size
>= 2);
1861 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1862 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1864 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1866 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1875 value_init(s
[n
].E
.d
);
1876 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1880 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1881 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1884 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1885 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1886 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1887 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1889 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1890 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1896 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1897 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1899 value_init(s
[n
].E
.d
);
1900 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1901 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1907 /* Just ignore; this may have been previously masked off */
1909 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1913 free_evalue_refs(&mone
);
1914 free_evalue_refs(mask
);
1915 free_evalue_refs(res
);
1918 evalue_set_si(res
, 0, 1);
1920 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1921 for (j
= 0; j
< n
; ++j
) {
1922 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1923 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1924 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1931 void evalue_copy(evalue
*dst
, const evalue
*src
)
1933 value_assign(dst
->d
, src
->d
);
1934 if(value_notzero_p(src
->d
)) {
1935 value_init(dst
->x
.n
);
1936 value_assign(dst
->x
.n
, src
->x
.n
);
1938 dst
->x
.p
= ecopy(src
->x
.p
);
1941 enode
*new_enode(enode_type type
,int size
,int pos
) {
1947 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1950 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1954 for(i
=0; i
<size
; i
++) {
1955 value_init(res
->arr
[i
].d
);
1956 value_set_si(res
->arr
[i
].d
,0);
1957 res
->arr
[i
].x
.p
= 0;
1962 enode
*ecopy(enode
*e
) {
1967 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1968 for(i
=0;i
<e
->size
;++i
) {
1969 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1970 if(value_zero_p(res
->arr
[i
].d
))
1971 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1972 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1973 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1975 value_init(res
->arr
[i
].x
.n
);
1976 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1982 int ecmp(const evalue
*e1
, const evalue
*e2
)
1988 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
1992 value_multiply(m
, e1
->x
.n
, e2
->d
);
1993 value_multiply(m2
, e2
->x
.n
, e1
->d
);
1995 if (value_lt(m
, m2
))
1997 else if (value_gt(m
, m2
))
2007 if (value_notzero_p(e1
->d
))
2009 if (value_notzero_p(e2
->d
))
2015 if (p1
->type
!= p2
->type
)
2016 return p1
->type
- p2
->type
;
2017 if (p1
->pos
!= p2
->pos
)
2018 return p1
->pos
- p2
->pos
;
2019 if (p1
->size
!= p2
->size
)
2020 return p1
->size
- p2
->size
;
2022 for (i
= p1
->size
-1; i
>= 0; --i
)
2023 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
2029 int eequal(const evalue
*e1
, const evalue
*e2
)
2034 if (value_ne(e1
->d
,e2
->d
))
2037 /* e1->d == e2->d */
2038 if (value_notzero_p(e1
->d
)) {
2039 if (value_ne(e1
->x
.n
,e2
->x
.n
))
2042 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2046 /* e1->d == e2->d == 0 */
2049 if (p1
->type
!= p2
->type
) return 0;
2050 if (p1
->size
!= p2
->size
) return 0;
2051 if (p1
->pos
!= p2
->pos
) return 0;
2052 for (i
=0; i
<p1
->size
; i
++)
2053 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
2058 void free_evalue_refs(evalue
*e
) {
2063 if (EVALUE_IS_DOMAIN(*e
)) {
2064 Domain_Free(EVALUE_DOMAIN(*e
));
2067 } else if (value_pos_p(e
->d
)) {
2069 /* 'e' stores a constant */
2071 value_clear(e
->x
.n
);
2074 assert(value_zero_p(e
->d
));
2077 if (!p
) return; /* null pointer */
2078 for (i
=0; i
<p
->size
; i
++) {
2079 free_evalue_refs(&(p
->arr
[i
]));
2083 } /* free_evalue_refs */
2085 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
2086 Vector
* val
, evalue
*res
)
2088 unsigned nparam
= periods
->Size
;
2091 double d
= compute_evalue(e
, val
->p
);
2092 d
*= VALUE_TO_DOUBLE(m
);
2097 value_assign(res
->d
, m
);
2098 value_init(res
->x
.n
);
2099 value_set_double(res
->x
.n
, d
);
2100 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
2103 if (value_one_p(periods
->p
[p
]))
2104 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
2109 value_assign(tmp
, periods
->p
[p
]);
2110 value_set_si(res
->d
, 0);
2111 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
2113 value_decrement(tmp
, tmp
);
2114 value_assign(val
->p
[p
], tmp
);
2115 mod2table_r(e
, periods
, m
, p
+1, val
,
2116 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
2117 } while (value_pos_p(tmp
));
2123 static void rel2table(evalue
*e
, int zero
)
2125 if (value_pos_p(e
->d
)) {
2126 if (value_zero_p(e
->x
.n
) == zero
)
2127 value_set_si(e
->x
.n
, 1);
2129 value_set_si(e
->x
.n
, 0);
2130 value_set_si(e
->d
, 1);
2133 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
2134 rel2table(&e
->x
.p
->arr
[i
], zero
);
2138 void evalue_mod2table(evalue
*e
, int nparam
)
2143 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2146 for (i
=0; i
<p
->size
; i
++) {
2147 evalue_mod2table(&(p
->arr
[i
]), nparam
);
2149 if (p
->type
== relation
) {
2154 evalue_copy(©
, &p
->arr
[0]);
2156 rel2table(&p
->arr
[0], 1);
2157 emul(&p
->arr
[0], &p
->arr
[1]);
2159 rel2table(©
, 0);
2160 emul(©
, &p
->arr
[2]);
2161 eadd(&p
->arr
[2], &p
->arr
[1]);
2162 free_evalue_refs(&p
->arr
[2]);
2163 free_evalue_refs(©
);
2165 free_evalue_refs(&p
->arr
[0]);
2169 } else if (p
->type
== fractional
) {
2170 Vector
*periods
= Vector_Alloc(nparam
);
2171 Vector
*val
= Vector_Alloc(nparam
);
2177 value_set_si(tmp
, 1);
2178 Vector_Set(periods
->p
, 1, nparam
);
2179 Vector_Set(val
->p
, 0, nparam
);
2180 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
2183 assert(p
->type
== polynomial
);
2184 assert(p
->size
== 2);
2185 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
2186 value_lcm(tmp
, p
->arr
[1].d
, &tmp
);
2188 value_lcm(tmp
, ev
->d
, &tmp
);
2190 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
2193 evalue_set_si(&res
, 0, 1);
2194 /* Compute the polynomial using Horner's rule */
2195 for (i
=p
->size
-1;i
>1;i
--) {
2196 eadd(&p
->arr
[i
], &res
);
2199 eadd(&p
->arr
[1], &res
);
2201 free_evalue_refs(e
);
2202 free_evalue_refs(&EP
);
2207 Vector_Free(periods
);
2209 } /* evalue_mod2table */
2211 /********************************************************/
2212 /* function in domain */
2213 /* check if the parameters in list_args */
2214 /* verifies the constraints of Domain P */
2215 /********************************************************/
2216 int in_domain(Polyhedron
*P
, Value
*list_args
)
2219 Value v
; /* value of the constraint of a row when
2220 parameters are instantiated*/
2224 for (row
= 0; row
< P
->NbConstraints
; row
++) {
2225 Inner_Product(P
->Constraint
[row
]+1, list_args
, P
->Dimension
, &v
);
2226 value_addto(v
, v
, P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
2227 if (value_neg_p(v
) ||
2228 value_zero_p(P
->Constraint
[row
][0]) && value_notzero_p(v
)) {
2235 return in
|| (P
->next
&& in_domain(P
->next
, list_args
));
2238 /****************************************************/
2239 /* function compute enode */
2240 /* compute the value of enode p with parameters */
2241 /* list "list_args */
2242 /* compute the polynomial or the periodic */
2243 /****************************************************/
2245 static double compute_enode(enode
*p
, Value
*list_args
) {
2257 if (p
->type
== polynomial
) {
2259 value_assign(param
,list_args
[p
->pos
-1]);
2261 /* Compute the polynomial using Horner's rule */
2262 for (i
=p
->size
-1;i
>0;i
--) {
2263 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2264 res
*=VALUE_TO_DOUBLE(param
);
2266 res
+=compute_evalue(&p
->arr
[0],list_args
);
2268 else if (p
->type
== fractional
) {
2269 double d
= compute_evalue(&p
->arr
[0], list_args
);
2270 d
-= floor(d
+1e-10);
2272 /* Compute the polynomial using Horner's rule */
2273 for (i
=p
->size
-1;i
>1;i
--) {
2274 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2277 res
+=compute_evalue(&p
->arr
[1],list_args
);
2279 else if (p
->type
== flooring
) {
2280 double d
= compute_evalue(&p
->arr
[0], list_args
);
2283 /* Compute the polynomial using Horner's rule */
2284 for (i
=p
->size
-1;i
>1;i
--) {
2285 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2288 res
+=compute_evalue(&p
->arr
[1],list_args
);
2290 else if (p
->type
== periodic
) {
2291 value_assign(m
,list_args
[p
->pos
-1]);
2293 /* Choose the right element of the periodic */
2294 value_set_si(param
,p
->size
);
2295 value_pmodulus(m
,m
,param
);
2296 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2298 else if (p
->type
== relation
) {
2299 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2300 res
= compute_evalue(&p
->arr
[1], list_args
);
2301 else if (p
->size
> 2)
2302 res
= compute_evalue(&p
->arr
[2], list_args
);
2304 else if (p
->type
== partition
) {
2305 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2306 Value
*vals
= list_args
;
2309 for (i
= 0; i
< dim
; ++i
) {
2310 value_init(vals
[i
]);
2312 value_assign(vals
[i
], list_args
[i
]);
2315 for (i
= 0; i
< p
->size
/2; ++i
)
2316 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2317 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2321 for (i
= 0; i
< dim
; ++i
)
2322 value_clear(vals
[i
]);
2331 } /* compute_enode */
2333 /*************************************************/
2334 /* return the value of Ehrhart Polynomial */
2335 /* It returns a double, because since it is */
2336 /* a recursive function, some intermediate value */
2337 /* might not be integral */
2338 /*************************************************/
2340 double compute_evalue(const evalue
*e
, Value
*list_args
)
2344 if (value_notzero_p(e
->d
)) {
2345 if (value_notone_p(e
->d
))
2346 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2348 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2351 res
= compute_enode(e
->x
.p
,list_args
);
2353 } /* compute_evalue */
2356 /****************************************************/
2357 /* function compute_poly : */
2358 /* Check for the good validity domain */
2359 /* return the number of point in the Polyhedron */
2360 /* in allocated memory */
2361 /* Using the Ehrhart pseudo-polynomial */
2362 /****************************************************/
2363 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2366 /* double d; int i; */
2368 tmp
= (Value
*) malloc (sizeof(Value
));
2369 assert(tmp
!= NULL
);
2371 value_set_si(*tmp
,0);
2374 return(tmp
); /* no ehrhart polynomial */
2375 if(en
->ValidityDomain
) {
2376 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2377 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2382 return(tmp
); /* no Validity Domain */
2384 if(in_domain(en
->ValidityDomain
,list_args
)) {
2386 #ifdef EVAL_EHRHART_DEBUG
2387 Print_Domain(stdout
,en
->ValidityDomain
);
2388 print_evalue(stdout
,&en
->EP
);
2391 /* d = compute_evalue(&en->EP,list_args);
2393 printf("(double)%lf = %d\n", d, i ); */
2394 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2400 value_set_si(*tmp
,0);
2401 return(tmp
); /* no compatible domain with the arguments */
2402 } /* compute_poly */
2404 static evalue
*eval_polynomial(const enode
*p
, int offset
,
2405 evalue
*base
, Value
*values
)
2410 res
= evalue_zero();
2411 for (i
= p
->size
-1; i
> offset
; --i
) {
2412 c
= evalue_eval(&p
->arr
[i
], values
);
2414 free_evalue_refs(c
);
2418 c
= evalue_eval(&p
->arr
[offset
], values
);
2420 free_evalue_refs(c
);
2426 evalue
*evalue_eval(const evalue
*e
, Value
*values
)
2433 if (value_notzero_p(e
->d
)) {
2434 res
= ALLOC(evalue
);
2436 evalue_copy(res
, e
);
2439 switch (e
->x
.p
->type
) {
2441 value_init(param
.x
.n
);
2442 value_assign(param
.x
.n
, values
[e
->x
.p
->pos
-1]);
2443 value_init(param
.d
);
2444 value_set_si(param
.d
, 1);
2446 res
= eval_polynomial(e
->x
.p
, 0, ¶m
, values
);
2447 free_evalue_refs(¶m
);
2450 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2451 mpz_fdiv_r(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2453 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2454 free_evalue_refs(param2
);
2458 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2459 mpz_fdiv_q(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2460 value_set_si(param2
->d
, 1);
2462 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2463 free_evalue_refs(param2
);
2467 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2468 if (value_zero_p(param2
->x
.n
))
2469 res
= evalue_eval(&e
->x
.p
->arr
[1], values
);
2470 else if (e
->x
.p
->size
> 2)
2471 res
= evalue_eval(&e
->x
.p
->arr
[2], values
);
2473 res
= evalue_zero();
2474 free_evalue_refs(param2
);
2478 assert(e
->x
.p
->pos
== EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
);
2479 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2480 if (in_domain(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), values
)) {
2481 res
= evalue_eval(&e
->x
.p
->arr
[2*i
+1], values
);
2485 res
= evalue_zero();
2493 size_t value_size(Value v
) {
2494 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2495 * sizeof(v
[0]._mp_d
[0]);
2498 size_t domain_size(Polyhedron
*D
)
2501 size_t s
= sizeof(*D
);
2503 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2504 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2505 s
+= value_size(D
->Constraint
[i
][j
]);
2508 for (i = 0; i < D->NbRays; ++i)
2509 for (j = 0; j < D->Dimension+2; ++j)
2510 s += value_size(D->Ray[i][j]);
2513 return D
->next
? s
+domain_size(D
->next
) : s
;
2516 size_t enode_size(enode
*p
) {
2517 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2520 if (p
->type
== partition
)
2521 for (i
= 0; i
< p
->size
/2; ++i
) {
2522 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2523 s
+= evalue_size(&p
->arr
[2*i
+1]);
2526 for (i
= 0; i
< p
->size
; ++i
) {
2527 s
+= evalue_size(&p
->arr
[i
]);
2532 size_t evalue_size(evalue
*e
)
2534 size_t s
= sizeof(*e
);
2535 s
+= value_size(e
->d
);
2536 if (value_notzero_p(e
->d
))
2537 s
+= value_size(e
->x
.n
);
2539 s
+= enode_size(e
->x
.p
);
2543 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2545 evalue
*found
= NULL
;
2550 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2553 value_init(offset
.d
);
2554 value_init(offset
.x
.n
);
2555 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2556 value_lcm(m
, offset
.d
, &offset
.d
);
2557 value_set_si(offset
.x
.n
, 1);
2560 evalue_copy(©
, cst
);
2563 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2565 if (eequal(base
, &e
->x
.p
->arr
[0]))
2566 found
= &e
->x
.p
->arr
[0];
2568 value_set_si(offset
.x
.n
, -2);
2571 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2573 if (eequal(base
, &e
->x
.p
->arr
[0]))
2576 free_evalue_refs(cst
);
2577 free_evalue_refs(&offset
);
2580 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2581 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2586 static evalue
*find_relation_pair(evalue
*e
)
2589 evalue
*found
= NULL
;
2591 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2594 if (e
->x
.p
->type
== fractional
) {
2599 poly_denom(&e
->x
.p
->arr
[0], &m
);
2601 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2602 cst
= &cst
->x
.p
->arr
[0])
2605 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2606 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2611 i
= e
->x
.p
->type
== relation
;
2612 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2613 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2618 void evalue_mod2relation(evalue
*e
) {
2621 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2624 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2625 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2626 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2627 value_clear(e
->x
.p
->arr
[2*i
].d
);
2628 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2630 if (2*i
< e
->x
.p
->size
) {
2631 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2632 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2637 if (e
->x
.p
->size
== 0) {
2639 evalue_set_si(e
, 0, 1);
2645 while ((d
= find_relation_pair(e
)) != NULL
) {
2649 value_init(split
.d
);
2650 value_set_si(split
.d
, 0);
2651 split
.x
.p
= new_enode(relation
, 3, 0);
2652 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2653 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2655 ev
= &split
.x
.p
->arr
[0];
2656 value_set_si(ev
->d
, 0);
2657 ev
->x
.p
= new_enode(fractional
, 3, -1);
2658 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2659 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2660 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2666 free_evalue_refs(&split
);
2670 static int evalue_comp(const void * a
, const void * b
)
2672 const evalue
*e1
= *(const evalue
**)a
;
2673 const evalue
*e2
= *(const evalue
**)b
;
2674 return ecmp(e1
, e2
);
2677 void evalue_combine(evalue
*e
)
2684 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2687 NALLOC(evs
, e
->x
.p
->size
/2);
2688 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2689 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2690 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2691 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2692 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2693 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2694 value_clear(p
->arr
[2*k
].d
);
2695 value_clear(p
->arr
[2*k
+1].d
);
2696 p
->arr
[2*k
] = *(evs
[i
]-1);
2697 p
->arr
[2*k
+1] = *(evs
[i
]);
2700 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2703 value_clear((evs
[i
]-1)->d
);
2707 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2708 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2709 free_evalue_refs(evs
[i
]);
2713 for (i
= 2*k
; i
< p
->size
; ++i
)
2714 value_clear(p
->arr
[i
].d
);
2721 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2723 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2725 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2728 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2729 Polyhedron
*D
, *N
, **P
;
2732 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2739 if (D
->NbEq
<= H
->NbEq
) {
2745 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2746 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2747 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2748 reduce_evalue(&tmp
);
2749 if (value_notzero_p(tmp
.d
) ||
2750 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2753 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2754 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2757 free_evalue_refs(&tmp
);
2763 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2765 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2767 value_clear(e
->x
.p
->arr
[2*i
].d
);
2768 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2770 if (2*i
< e
->x
.p
->size
) {
2771 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2772 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2779 H
= DomainConvex(D
, 0);
2780 E
= DomainDifference(H
, D
, 0);
2782 D
= DomainDifference(H
, E
, 0);
2785 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2789 /* Use smallest representative for coefficients in affine form in
2790 * argument of fractional.
2791 * Since any change will make the argument non-standard,
2792 * the containing evalue will have to be reduced again afterward.
2794 static void fractional_minimal_coefficients(enode
*p
)
2800 assert(p
->type
== fractional
);
2802 while (value_zero_p(pp
->d
)) {
2803 assert(pp
->x
.p
->type
== polynomial
);
2804 assert(pp
->x
.p
->size
== 2);
2805 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2806 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2807 if (value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2808 value_subtract(pp
->x
.p
->arr
[1].x
.n
,
2809 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2810 pp
= &pp
->x
.p
->arr
[0];
2816 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2821 unsigned dim
= D
->Dimension
;
2822 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2825 assert(p
->type
== fractional
);
2827 value_set_si(T
->p
[1][dim
], 1);
2829 while (value_zero_p(pp
->d
)) {
2830 assert(pp
->x
.p
->type
== polynomial
);
2831 assert(pp
->x
.p
->size
== 2);
2832 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2833 value_division(T
->p
[0][pp
->x
.p
->pos
-1], *d
, pp
->x
.p
->arr
[1].d
);
2834 value_multiply(T
->p
[0][pp
->x
.p
->pos
-1],
2835 T
->p
[0][pp
->x
.p
->pos
-1], pp
->x
.p
->arr
[1].x
.n
);
2836 pp
= &pp
->x
.p
->arr
[0];
2838 value_division(T
->p
[0][dim
], *d
, pp
->d
);
2839 value_multiply(T
->p
[0][dim
], T
->p
[0][dim
], pp
->x
.n
);
2840 I
= DomainImage(D
, T
, 0);
2841 H
= DomainConvex(I
, 0);
2851 int evalue_range_reduction_in_domain(evalue
*e
, Polyhedron
*D
)
2860 if (value_notzero_p(e
->d
))
2865 if (p
->type
== relation
) {
2872 fractional_minimal_coefficients(p
->arr
[0].x
.p
);
2873 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
);
2874 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2875 equal
= value_eq(min
, max
);
2876 mpz_cdiv_q(min
, min
, d
);
2877 mpz_fdiv_q(max
, max
, d
);
2879 if (bounded
&& value_gt(min
, max
)) {
2885 evalue_set_si(e
, 0, 1);
2888 free_evalue_refs(&(p
->arr
[1]));
2889 free_evalue_refs(&(p
->arr
[0]));
2895 return r
? r
: evalue_range_reduction_in_domain(e
, D
);
2896 } else if (bounded
&& equal
) {
2899 free_evalue_refs(&(p
->arr
[2]));
2902 free_evalue_refs(&(p
->arr
[0]));
2908 return evalue_range_reduction_in_domain(e
, D
);
2909 } else if (bounded
&& value_eq(min
, max
)) {
2910 /* zero for a single value */
2912 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2913 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2914 value_multiply(min
, min
, d
);
2915 value_subtract(M
->p
[0][D
->Dimension
+1],
2916 M
->p
[0][D
->Dimension
+1], min
);
2917 E
= DomainAddConstraints(D
, M
, 0);
2923 r
= evalue_range_reduction_in_domain(&p
->arr
[1], E
);
2925 r
|= evalue_range_reduction_in_domain(&p
->arr
[2], D
);
2927 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2935 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2938 i
= p
->type
== relation
? 1 :
2939 p
->type
== fractional
? 1 : 0;
2940 for (; i
<p
->size
; i
++)
2941 r
|= evalue_range_reduction_in_domain(&p
->arr
[i
], D
);
2943 if (p
->type
!= fractional
) {
2944 if (r
&& p
->type
== polynomial
) {
2947 value_set_si(f
.d
, 0);
2948 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2949 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2950 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2951 reorder_terms_about(p
, &f
);
2962 fractional_minimal_coefficients(p
);
2963 I
= polynomial_projection(p
, D
, &d
, NULL
);
2964 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2965 mpz_fdiv_q(min
, min
, d
);
2966 mpz_fdiv_q(max
, max
, d
);
2967 value_subtract(d
, max
, min
);
2969 if (bounded
&& value_eq(min
, max
)) {
2972 value_init(inc
.x
.n
);
2973 value_set_si(inc
.d
, 1);
2974 value_oppose(inc
.x
.n
, min
);
2975 eadd(&inc
, &p
->arr
[0]);
2976 reorder_terms_about(p
, &p
->arr
[0]); /* frees arr[0] */
2980 free_evalue_refs(&inc
);
2982 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
2983 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
2984 * See pages 199-200 of PhD thesis.
2992 value_set_si(rem
.d
, 0);
2993 rem
.x
.p
= new_enode(fractional
, 3, -1);
2994 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
2995 value_clear(rem
.x
.p
->arr
[1].d
);
2996 value_clear(rem
.x
.p
->arr
[2].d
);
2997 rem
.x
.p
->arr
[1] = p
->arr
[1];
2998 rem
.x
.p
->arr
[2] = p
->arr
[2];
2999 for (i
= 3; i
< p
->size
; ++i
)
3000 p
->arr
[i
-2] = p
->arr
[i
];
3004 value_init(inc
.x
.n
);
3005 value_set_si(inc
.d
, 1);
3006 value_oppose(inc
.x
.n
, min
);
3009 evalue_copy(&t
, &p
->arr
[0]);
3013 value_set_si(f
.d
, 0);
3014 f
.x
.p
= new_enode(fractional
, 3, -1);
3015 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3016 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3017 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
3019 value_init(factor
.d
);
3020 evalue_set_si(&factor
, -1, 1);
3026 value_clear(f
.x
.p
->arr
[1].x
.n
);
3027 value_clear(f
.x
.p
->arr
[2].x
.n
);
3028 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3029 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3033 reorder_terms(&rem
);
3040 free_evalue_refs(&inc
);
3041 free_evalue_refs(&t
);
3042 free_evalue_refs(&f
);
3043 free_evalue_refs(&factor
);
3044 free_evalue_refs(&rem
);
3046 evalue_range_reduction_in_domain(e
, D
);
3050 _reduce_evalue(&p
->arr
[0], 0, 1);
3062 void evalue_range_reduction(evalue
*e
)
3065 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3068 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3069 if (evalue_range_reduction_in_domain(&e
->x
.p
->arr
[2*i
+1],
3070 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
3071 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3073 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
3074 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
3075 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3076 value_clear(e
->x
.p
->arr
[2*i
].d
);
3078 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
3079 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
3087 Enumeration
* partition2enumeration(evalue
*EP
)
3090 Enumeration
*en
, *res
= NULL
;
3092 if (EVALUE_IS_ZERO(*EP
)) {
3097 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
3098 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
3099 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
3102 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
3103 value_clear(EP
->x
.p
->arr
[2*i
].d
);
3104 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
3112 int evalue_frac2floor_in_domain3(evalue
*e
, Polyhedron
*D
, int shift
)
3121 if (value_notzero_p(e
->d
))
3126 i
= p
->type
== relation
? 1 :
3127 p
->type
== fractional
? 1 : 0;
3128 for (; i
<p
->size
; i
++)
3129 r
|= evalue_frac2floor_in_domain3(&p
->arr
[i
], D
, shift
);
3131 if (p
->type
!= fractional
) {
3132 if (r
&& p
->type
== polynomial
) {
3135 value_set_si(f
.d
, 0);
3136 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3137 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3138 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3139 reorder_terms_about(p
, &f
);
3149 I
= polynomial_projection(p
, D
, &d
, NULL
);
3152 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3155 assert(I
->NbEq
== 0); /* Should have been reduced */
3158 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3159 if (value_pos_p(I
->Constraint
[i
][1]))
3162 if (i
< I
->NbConstraints
) {
3164 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3165 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3166 if (value_neg_p(min
)) {
3168 mpz_fdiv_q(min
, min
, d
);
3169 value_init(offset
.d
);
3170 value_set_si(offset
.d
, 1);
3171 value_init(offset
.x
.n
);
3172 value_oppose(offset
.x
.n
, min
);
3173 eadd(&offset
, &p
->arr
[0]);
3174 free_evalue_refs(&offset
);
3184 value_set_si(fl
.d
, 0);
3185 fl
.x
.p
= new_enode(flooring
, 3, -1);
3186 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
3187 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
3188 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
3190 eadd(&fl
, &p
->arr
[0]);
3191 reorder_terms_about(p
, &p
->arr
[0]);
3195 free_evalue_refs(&fl
);
3200 int evalue_frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
3202 return evalue_frac2floor_in_domain3(e
, D
, 1);
3205 void evalue_frac2floor2(evalue
*e
, int shift
)
3208 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
3210 if (evalue_frac2floor_in_domain3(e
, NULL
, 0))
3216 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3217 if (evalue_frac2floor_in_domain3(&e
->x
.p
->arr
[2*i
+1],
3218 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), shift
))
3219 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3222 void evalue_frac2floor(evalue
*e
)
3224 evalue_frac2floor2(e
, 1);
3227 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
3232 int nparam
= D
->Dimension
- nvar
;
3235 nr
= D
->NbConstraints
+ 2;
3236 nc
= D
->Dimension
+ 2 + 1;
3237 C
= Matrix_Alloc(nr
, nc
);
3238 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
3239 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
3240 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3241 D
->Dimension
+ 1 - nvar
);
3246 nc
= C
->NbColumns
+ 1;
3247 C
= Matrix_Alloc(nr
, nc
);
3248 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
3249 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
3250 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3251 oldC
->NbColumns
- 1 - nvar
);
3254 value_set_si(C
->p
[nr
-2][0], 1);
3255 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
3256 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
3258 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
3259 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
3265 static void floor2frac_r(evalue
*e
, int nvar
)
3272 if (value_notzero_p(e
->d
))
3277 assert(p
->type
== flooring
);
3278 for (i
= 1; i
< p
->size
; i
++)
3279 floor2frac_r(&p
->arr
[i
], nvar
);
3281 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
3282 assert(pp
->x
.p
->type
== polynomial
);
3283 pp
->x
.p
->pos
-= nvar
;
3287 value_set_si(f
.d
, 0);
3288 f
.x
.p
= new_enode(fractional
, 3, -1);
3289 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3290 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3291 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3293 eadd(&f
, &p
->arr
[0]);
3294 reorder_terms_about(p
, &p
->arr
[0]);
3298 free_evalue_refs(&f
);
3301 /* Convert flooring back to fractional and shift position
3302 * of the parameters by nvar
3304 static void floor2frac(evalue
*e
, int nvar
)
3306 floor2frac_r(e
, nvar
);
3310 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
3313 int nparam
= D
->Dimension
- nvar
;
3317 D
= Constraints2Polyhedron(C
, 0);
3321 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
3323 /* Double check that D was not unbounded. */
3324 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
3332 evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
3339 evalue
*factor
= NULL
;
3342 if (EVALUE_IS_ZERO(*e
))
3346 Polyhedron
*DD
= Disjoint_Domain(D
, 0, 0);
3353 res
= esum_over_domain(e
, nvar
, Q
, C
);
3356 for (Q
= DD
; Q
; Q
= DD
) {
3362 t
= esum_over_domain(e
, nvar
, Q
, C
);
3369 free_evalue_refs(t
);
3376 if (value_notzero_p(e
->d
)) {
3379 t
= esum_over_domain_cst(nvar
, D
, C
);
3381 if (!EVALUE_IS_ONE(*e
))
3387 switch (e
->x
.p
->type
) {
3389 evalue
*pp
= &e
->x
.p
->arr
[0];
3391 if (pp
->x
.p
->pos
> nvar
) {
3392 /* remainder is independent of the summated vars */
3398 floor2frac(&f
, nvar
);
3400 t
= esum_over_domain_cst(nvar
, D
, C
);
3404 free_evalue_refs(&f
);
3409 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3410 poly_denom(pp
, &row
->p
[1 + nvar
]);
3411 value_set_si(row
->p
[0], 1);
3412 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3413 pp
= &pp
->x
.p
->arr
[0]) {
3415 assert(pp
->x
.p
->type
== polynomial
);
3417 if (pos
>= 1 + nvar
)
3419 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3420 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3421 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3423 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3424 value_division(row
->p
[1 + D
->Dimension
+ 1],
3425 row
->p
[1 + D
->Dimension
+ 1],
3427 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3428 row
->p
[1 + D
->Dimension
+ 1],
3430 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3434 int pos
= e
->x
.p
->pos
;
3437 factor
= ALLOC(evalue
);
3438 value_init(factor
->d
);
3439 value_set_si(factor
->d
, 0);
3440 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3441 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3442 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3446 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3447 for (i
= 0; i
< D
->NbRays
; ++i
)
3448 if (value_notzero_p(D
->Ray
[i
][pos
]))
3450 assert(i
< D
->NbRays
);
3451 if (value_neg_p(D
->Ray
[i
][pos
])) {
3452 factor
= ALLOC(evalue
);
3453 value_init(factor
->d
);
3454 evalue_set_si(factor
, -1, 1);
3456 value_set_si(row
->p
[0], 1);
3457 value_set_si(row
->p
[pos
], 1);
3458 value_set_si(row
->p
[1 + nvar
], -1);
3465 i
= type_offset(e
->x
.p
);
3467 res
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
3472 evalue_copy(&cum
, factor
);
3476 for (; i
< e
->x
.p
->size
; ++i
) {
3480 C
= esum_add_constraint(nvar
, D
, C
, row
);
3486 Vector_Print(stderr, P_VALUE_FMT, row);
3488 Matrix_Print(stderr, P_VALUE_FMT, C);
3490 t
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
3499 free_evalue_refs(t
);
3502 if (factor
&& i
+1 < e
->x
.p
->size
)
3509 free_evalue_refs(factor
);
3510 free_evalue_refs(&cum
);
3522 evalue
*esum(evalue
*e
, int nvar
)
3525 evalue
*res
= ALLOC(evalue
);
3529 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3530 evalue_copy(res
, e
);
3534 evalue_set_si(res
, 0, 1);
3536 assert(value_zero_p(e
->d
));
3537 assert(e
->x
.p
->type
== partition
);
3539 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3541 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3542 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
3544 free_evalue_refs(t
);
3553 /* Initial silly implementation */
3554 void eor(evalue
*e1
, evalue
*res
)
3560 evalue_set_si(&mone
, -1, 1);
3562 evalue_copy(&E
, res
);
3568 free_evalue_refs(&E
);
3569 free_evalue_refs(&mone
);
3572 /* computes denominator of polynomial evalue
3573 * d should point to a value initialized to 1
3575 void evalue_denom(const evalue
*e
, Value
*d
)
3579 if (value_notzero_p(e
->d
)) {
3580 value_lcm(*d
, e
->d
, d
);
3583 assert(e
->x
.p
->type
== polynomial
||
3584 e
->x
.p
->type
== fractional
||
3585 e
->x
.p
->type
== flooring
);
3586 offset
= type_offset(e
->x
.p
);
3587 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3588 evalue_denom(&e
->x
.p
->arr
[i
], d
);
3591 /* Divides the evalue e by the integer n */
3592 void evalue_div(evalue
* e
, Value n
)
3596 if (value_notzero_p(e
->d
)) {
3599 value_multiply(e
->d
, e
->d
, n
);
3600 Gcd(e
->x
.n
, e
->d
, &gc
);
3601 if (value_notone_p(gc
)) {
3602 value_division(e
->d
, e
->d
, gc
);
3603 value_division(e
->x
.n
, e
->x
.n
, gc
);
3608 if (e
->x
.p
->type
== partition
) {
3609 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3610 evalue_div(&e
->x
.p
->arr
[2*i
+1], n
);
3613 offset
= type_offset(e
->x
.p
);
3614 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3615 evalue_div(&e
->x
.p
->arr
[i
], n
);
3618 void evalue_add_constant(evalue
*e
, const Value cst
)
3622 if (value_zero_p(e
->d
)) {
3623 if (e
->x
.p
->type
== partition
) {
3624 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3625 evalue_add_constant(&e
->x
.p
->arr
[2*i
+1], cst
);
3628 if (e
->x
.p
->type
== relation
) {
3629 for (i
= 1; i
< e
->x
.p
->size
; ++i
)
3630 evalue_add_constant(&e
->x
.p
->arr
[i
], cst
);
3634 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
3635 } while (value_zero_p(e
->d
));
3637 value_addmul(e
->x
.n
, cst
, e
->d
);
3640 static void evalue_frac2polynomial_r(evalue
*e
, int *signs
, int sign
, int in_frac
)
3645 int sign_odd
= sign
;
3647 if (value_notzero_p(e
->d
)) {
3648 if (in_frac
&& sign
* value_sign(e
->x
.n
) < 0) {
3649 value_set_si(e
->x
.n
, 0);
3650 value_set_si(e
->d
, 1);
3655 if (e
->x
.p
->type
== relation
) {
3656 for (i
= e
->x
.p
->size
-1; i
>= 1; --i
)
3657 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
, sign
, in_frac
);
3661 if (e
->x
.p
->type
== polynomial
)
3662 sign_odd
*= signs
[e
->x
.p
->pos
-1];
3663 offset
= type_offset(e
->x
.p
);
3664 evalue_frac2polynomial_r(&e
->x
.p
->arr
[offset
], signs
, sign
, in_frac
);
3665 in_frac
|= e
->x
.p
->type
== fractional
;
3666 for (i
= e
->x
.p
->size
-1; i
> offset
; --i
)
3667 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
,
3668 (i
- offset
) % 2 ? sign_odd
: sign
, in_frac
);
3670 if (e
->x
.p
->type
!= fractional
)
3673 /* replace { a/m } by (m-1)/m if sign != 0
3674 * and by (m-1)/(2m) if sign == 0
3678 evalue_denom(&e
->x
.p
->arr
[0], &d
);
3679 free_evalue_refs(&e
->x
.p
->arr
[0]);
3680 value_init(e
->x
.p
->arr
[0].d
);
3681 value_init(e
->x
.p
->arr
[0].x
.n
);
3683 value_addto(e
->x
.p
->arr
[0].d
, d
, d
);
3685 value_assign(e
->x
.p
->arr
[0].d
, d
);
3686 value_decrement(e
->x
.p
->arr
[0].x
.n
, d
);
3690 reorder_terms_about(p
, &p
->arr
[0]);
3696 /* Approximate the evalue in fractional representation by a polynomial.
3697 * If sign > 0, the result is an upper bound;
3698 * if sign < 0, the result is a lower bound;
3699 * if sign = 0, the result is an intermediate approximation.
3701 void evalue_frac2polynomial(evalue
*e
, int sign
, unsigned MaxRays
)
3706 if (value_notzero_p(e
->d
))
3708 assert(e
->x
.p
->type
== partition
);
3709 /* make sure all variables in the domains have a fixed sign */
3711 evalue_split_domains_into_orthants(e
, MaxRays
);
3713 assert(e
->x
.p
->size
>= 2);
3714 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3716 signs
= alloca(sizeof(int) * dim
);
3718 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3719 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
3720 POL_ENSURE_VERTICES(D
);
3721 for (j
= 0; j
< dim
; ++j
) {
3725 for (k
= 0; k
< D
->NbRays
; ++k
) {
3726 signs
[j
] = value_sign(D
->Ray
[k
][1+j
]);
3731 evalue_frac2polynomial_r(&e
->x
.p
->arr
[2*i
+1], signs
, sign
, 0);
3735 /* Split the domains of e (which is assumed to be a partition)
3736 * such that each resulting domain lies entirely in one orthant.
3738 void evalue_split_domains_into_orthants(evalue
*e
, unsigned MaxRays
)
3741 assert(value_zero_p(e
->d
));
3742 assert(e
->x
.p
->type
== partition
);
3743 assert(e
->x
.p
->size
>= 2);
3744 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3746 for (i
= 0; i
< dim
; ++i
) {
3749 C
= Matrix_Alloc(1, 1 + dim
+ 1);
3750 value_set_si(C
->p
[0][0], 1);
3751 value_init(split
.d
);
3752 value_set_si(split
.d
, 0);
3753 split
.x
.p
= new_enode(partition
, 4, dim
);
3754 value_set_si(C
->p
[0][1+i
], 1);
3755 C2
= Matrix_Copy(C
);
3756 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[0], Constraints2Polyhedron(C2
, MaxRays
));
3758 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
3759 value_set_si(C
->p
[0][1+i
], -1);
3760 value_set_si(C
->p
[0][1+dim
], -1);
3761 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[2], Constraints2Polyhedron(C
, MaxRays
));
3762 evalue_set_si(&split
.x
.p
->arr
[3], 1, 1);
3764 free_evalue_refs(&split
);
3769 static Matrix
*find_fractional_with_max_periods(evalue
*e
, Polyhedron
*D
,
3771 Value
*min
, Value
*max
)
3777 if (value_notzero_p(e
->d
))
3780 if (e
->x
.p
->type
== fractional
) {
3785 I
= polynomial_projection(e
->x
.p
, D
, &d
, &T
);
3786 bounded
= line_minmax(I
, min
, max
); /* frees I */
3790 value_set_si(mp
, max_periods
);
3791 mpz_fdiv_q(*min
, *min
, d
);
3792 mpz_fdiv_q(*max
, *max
, d
);
3793 value_assign(T
->p
[1][D
->Dimension
], d
);
3794 value_subtract(d
, *max
, *min
);
3795 if (value_ge(d
, mp
)) {
3809 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
3810 if ((T
= find_fractional_with_max_periods(&e
->x
.p
->arr
[i
], D
, max_periods
,
3817 /* Look for fractional parts that can be removed by splitting the corresponding
3818 * domain into at most max_periods parts.
3819 * We use a very simply strategy that looks for the first fractional part
3820 * that satisfies the condition, performs the split and then continues
3821 * looking for other fractional parts in the split domains until no
3822 * such fractional part can be found anymore.
3824 void evalue_split_periods(evalue
*e
, int max_periods
, unsigned int MaxRays
)
3831 if (EVALUE_IS_ZERO(*e
))
3833 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
3835 "WARNING: evalue_split_periods called on incorrect evalue type\n");
3843 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3847 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
3849 T
= find_fractional_with_max_periods(&e
->x
.p
->arr
[2*i
+1], D
, max_periods
,
3854 M
= Matrix_Alloc(2, 2+D
->Dimension
);
3856 value_subtract(d
, max
, min
);
3857 n
= VALUE_TO_INT(d
)+1;
3859 value_set_si(M
->p
[0][0], 1);
3860 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
3861 value_multiply(d
, max
, T
->p
[1][D
->Dimension
]);
3862 value_subtract(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
], d
);
3863 value_set_si(d
, -1);
3864 value_set_si(M
->p
[1][0], 1);
3865 Vector_Scale(T
->p
[0], M
->p
[1]+1, d
, D
->Dimension
+1);
3866 value_addmul(M
->p
[1][1+D
->Dimension
], max
, T
->p
[1][D
->Dimension
]);
3867 value_addto(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
3868 T
->p
[1][D
->Dimension
]);
3869 value_decrement(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
]);
3871 p
= new_enode(partition
, e
->x
.p
->size
+ (n
-1)*2, e
->x
.p
->pos
);
3872 for (j
= 0; j
< 2*i
; ++j
) {
3873 value_clear(p
->arr
[j
].d
);
3874 p
->arr
[j
] = e
->x
.p
->arr
[j
];
3876 for (j
= 2*i
+2; j
< e
->x
.p
->size
; ++j
) {
3877 value_clear(p
->arr
[j
+2*(n
-1)].d
);
3878 p
->arr
[j
+2*(n
-1)] = e
->x
.p
->arr
[j
];
3880 for (j
= n
-1; j
>= 0; --j
) {
3882 value_clear(p
->arr
[2*i
+1].d
);
3883 p
->arr
[2*i
+1] = e
->x
.p
->arr
[2*i
+1];
3885 evalue_copy(&p
->arr
[2*(i
+j
)+1], &e
->x
.p
->arr
[2*i
+1]);
3887 value_subtract(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
3888 T
->p
[1][D
->Dimension
]);
3889 value_addto(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
],
3890 T
->p
[1][D
->Dimension
]);
3892 E
= DomainAddConstraints(D
, M
, MaxRays
);
3893 EVALUE_SET_DOMAIN(p
->arr
[2*(i
+j
)], E
);
3894 if (evalue_range_reduction_in_domain(&p
->arr
[2*(i
+j
)+1], E
))
3895 reduce_evalue(&p
->arr
[2*(i
+j
)+1]);
3897 value_clear(e
->x
.p
->arr
[2*i
].d
);
3911 void evalue_extract_affine(const evalue
*e
, Value
*coeff
, Value
*cst
, Value
*d
)
3913 value_set_si(*d
, 1);
3915 for ( ; value_zero_p(e
->d
); e
= &e
->x
.p
->arr
[0]) {
3917 assert(e
->x
.p
->type
== polynomial
);
3918 assert(e
->x
.p
->size
== 2);
3919 c
= &e
->x
.p
->arr
[1];
3920 value_multiply(coeff
[e
->x
.p
->pos
-1], *d
, c
->x
.n
);
3921 value_division(coeff
[e
->x
.p
->pos
-1], coeff
[e
->x
.p
->pos
-1], c
->d
);
3923 value_multiply(*cst
, *d
, e
->x
.n
);
3924 value_division(*cst
, *cst
, e
->d
);
3927 /* returns an evalue that corresponds to
3931 static evalue
*term(int param
, Value c
, Value den
)
3933 evalue
*EP
= ALLOC(evalue
);
3935 value_set_si(EP
->d
,0);
3936 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
3937 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
3938 value_init(EP
->x
.p
->arr
[1].x
.n
);
3939 value_assign(EP
->x
.p
->arr
[1].d
, den
);
3940 value_assign(EP
->x
.p
->arr
[1].x
.n
, c
);
3944 evalue
*affine2evalue(Value
*coeff
, Value denom
, int nvar
)
3947 evalue
*E
= ALLOC(evalue
);
3949 evalue_set(E
, coeff
[nvar
], denom
);
3950 for (i
= 0; i
< nvar
; ++i
) {
3951 if (value_zero_p(coeff
[i
]))
3953 evalue
*t
= term(i
, coeff
[i
], denom
);
3955 free_evalue_refs(t
);
3961 void evalue_substitute(evalue
*e
, evalue
**subs
)
3967 if (value_notzero_p(e
->d
))
3971 assert(p
->type
!= partition
);
3973 for (i
= 0; i
< p
->size
; ++i
)
3974 evalue_substitute(&p
->arr
[i
], subs
);
3976 if (p
->type
== polynomial
)
3981 value_set_si(v
->d
, 0);
3982 v
->x
.p
= new_enode(p
->type
, 3, -1);
3983 value_clear(v
->x
.p
->arr
[0].d
);
3984 v
->x
.p
->arr
[0] = p
->arr
[0];
3985 evalue_set_si(&v
->x
.p
->arr
[1], 0, 1);
3986 evalue_set_si(&v
->x
.p
->arr
[2], 1, 1);
3989 offset
= type_offset(p
);
3991 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
3992 emul(v
, &p
->arr
[i
]);
3993 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
3994 free_evalue_refs(&(p
->arr
[i
]));
3997 if (p
->type
!= polynomial
) {
3998 free_evalue_refs(v
);
4003 *e
= p
->arr
[offset
];
4007 /* evalue e is given in terms of "new" parameter; CP maps the new
4008 * parameters back to the old parameters.
4009 * Transforms e such that it refers back to the old parameters.
4011 void evalue_backsubstitute(evalue
*e
, Matrix
*CP
, unsigned MaxRays
)
4018 unsigned nparam
= CP
->NbColumns
-1;
4021 if (EVALUE_IS_ZERO(*e
))
4024 assert(value_zero_p(e
->d
));
4026 assert(p
->type
== partition
);
4028 inv
= left_inverse(CP
, &eq
);
4029 subs
= ALLOCN(evalue
*, nparam
);
4030 for (i
= 0; i
< nparam
; ++i
)
4031 subs
[i
] = affine2evalue(inv
->p
[i
], inv
->p
[nparam
][inv
->NbColumns
-1],
4034 CEq
= Constraints2Polyhedron(eq
, MaxRays
);
4035 addeliminatedparams_partition(p
, inv
, CEq
, inv
->NbColumns
-1, MaxRays
);
4036 Polyhedron_Free(CEq
);
4038 for (i
= 0; i
< p
->size
/2; ++i
)
4039 evalue_substitute(&p
->arr
[2*i
+1], subs
);
4041 for (i
= 0; i
< nparam
; ++i
) {
4042 free_evalue_refs(subs
[i
]);