1 /***********************************************************************/
2 /* copyright 1997, Doran Wilde */
3 /* copyright 1997-2000, Vincent Loechner */
4 /* copyright 2003-2006, Sven Verdoolaege */
5 /* Permission is granted to copy, use, and distribute */
6 /* for any commercial or noncommercial purpose under the terms */
7 /* of the GNU General Public license, version 2, June 1991 */
8 /* (see file : LICENSE). */
9 /***********************************************************************/
16 #include <barvinok/evalue.h>
17 #include <barvinok/barvinok.h>
18 #include <barvinok/util.h>
21 #ifndef value_pmodulus
22 #define value_pmodulus(ref,val1,val2) (mpz_fdiv_r((ref),(val1),(val2)))
25 #define ALLOC(type) (type*)malloc(sizeof(type))
26 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
29 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
31 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
34 void evalue_set_si(evalue
*ev
, int n
, int d
) {
35 value_set_si(ev
->d
, d
);
37 value_set_si(ev
->x
.n
, n
);
40 void evalue_set(evalue
*ev
, Value n
, Value d
) {
41 value_assign(ev
->d
, d
);
43 value_assign(ev
->x
.n
, n
);
46 void evalue_set_reduce(evalue
*ev
, Value n
, Value d
) {
48 value_gcd(ev
->x
.n
, n
, d
);
49 value_divexact(ev
->d
, d
, ev
->x
.n
);
50 value_divexact(ev
->x
.n
, n
, ev
->x
.n
);
55 evalue
*EP
= ALLOC(evalue
);
57 evalue_set_si(EP
, 0, 1);
63 evalue
*EP
= ALLOC(evalue
);
65 value_set_si(EP
->d
, -2);
70 /* returns an evalue that corresponds to
74 evalue
*evalue_var(int var
)
76 evalue
*EP
= ALLOC(evalue
);
78 value_set_si(EP
->d
,0);
79 EP
->x
.p
= new_enode(polynomial
, 2, var
+ 1);
80 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
81 evalue_set_si(&EP
->x
.p
->arr
[1], 1, 1);
85 void aep_evalue(evalue
*e
, int *ref
) {
90 if (value_notzero_p(e
->d
))
91 return; /* a rational number, its already reduced */
93 return; /* hum... an overflow probably occured */
95 /* First check the components of p */
96 for (i
=0;i
<p
->size
;i
++)
97 aep_evalue(&p
->arr
[i
],ref
);
104 p
->pos
= ref
[p
->pos
-1]+1;
110 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
116 if (value_notzero_p(e
->d
))
117 return; /* a rational number, its already reduced */
119 return; /* hum... an overflow probably occured */
122 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
123 for(i
=0;i
<CT
->NbRows
-1;i
++)
124 for(j
=0;j
<CT
->NbColumns
;j
++)
125 if(value_notzero_p(CT
->p
[i
][j
])) {
130 /* Transform the references in e, using ref */
134 } /* addeliminatedparams_evalue */
136 static void addeliminatedparams_partition(enode
*p
, Matrix
*CT
, Polyhedron
*CEq
,
137 unsigned nparam
, unsigned MaxRays
)
140 assert(p
->type
== partition
);
143 for (i
= 0; i
< p
->size
/2; i
++) {
144 Polyhedron
*D
= EVALUE_DOMAIN(p
->arr
[2*i
]);
145 Polyhedron
*T
= DomainPreimage(D
, CT
, MaxRays
);
149 T
= DomainIntersection(D
, CEq
, MaxRays
);
152 EVALUE_SET_DOMAIN(p
->arr
[2*i
], T
);
156 void addeliminatedparams_enum(evalue
*e
, Matrix
*CT
, Polyhedron
*CEq
,
157 unsigned MaxRays
, unsigned nparam
)
162 if (CT
->NbRows
== CT
->NbColumns
)
165 if (EVALUE_IS_ZERO(*e
))
168 if (value_notzero_p(e
->d
)) {
171 value_set_si(res
.d
, 0);
172 res
.x
.p
= new_enode(partition
, 2, nparam
);
173 EVALUE_SET_DOMAIN(res
.x
.p
->arr
[0],
174 DomainConstraintSimplify(Polyhedron_Copy(CEq
), MaxRays
));
175 value_clear(res
.x
.p
->arr
[1].d
);
176 res
.x
.p
->arr
[1] = *e
;
184 addeliminatedparams_partition(p
, CT
, CEq
, nparam
, MaxRays
);
185 for (i
= 0; i
< p
->size
/2; i
++)
186 addeliminatedparams_evalue(&p
->arr
[2*i
+1], CT
);
189 static int mod_rational_cmp(evalue
*e1
, evalue
*e2
)
197 assert(value_notzero_p(e1
->d
));
198 assert(value_notzero_p(e2
->d
));
199 value_multiply(m
, e1
->x
.n
, e2
->d
);
200 value_multiply(m2
, e2
->x
.n
, e1
->d
);
203 else if (value_gt(m
, m2
))
213 static int mod_term_cmp_r(evalue
*e1
, evalue
*e2
)
215 if (value_notzero_p(e1
->d
)) {
217 if (value_zero_p(e2
->d
))
219 return mod_rational_cmp(e1
, e2
);
221 if (value_notzero_p(e2
->d
))
223 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
225 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
228 int r
= mod_rational_cmp(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
230 ? mod_term_cmp_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
235 static int mod_term_cmp(const evalue
*e1
, const evalue
*e2
)
237 assert(value_zero_p(e1
->d
));
238 assert(value_zero_p(e2
->d
));
239 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
240 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
241 return mod_term_cmp_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
244 static void check_order(const evalue
*e
)
249 if (value_notzero_p(e
->d
))
252 switch (e
->x
.p
->type
) {
254 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
255 check_order(&e
->x
.p
->arr
[2*i
+1]);
258 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
260 if (value_notzero_p(a
->d
))
262 switch (a
->x
.p
->type
) {
264 assert(mod_term_cmp(&e
->x
.p
->arr
[0], &a
->x
.p
->arr
[0]) < 0);
273 for (i
= 0; i
< e
->x
.p
->size
; ++i
) {
275 if (value_notzero_p(a
->d
))
277 switch (a
->x
.p
->type
) {
279 assert(e
->x
.p
->pos
< a
->x
.p
->pos
);
290 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
292 if (value_notzero_p(a
->d
))
294 switch (a
->x
.p
->type
) {
305 /* Negative pos means inequality */
306 /* s is negative of substitution if m is not zero */
315 struct fixed_param
*fixed
;
320 static int relations_depth(evalue
*e
)
325 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
326 e
= &e
->x
.p
->arr
[1], ++d
);
330 static void poly_denom_not_constant(evalue
**pp
, Value
*d
)
335 while (value_zero_p(p
->d
)) {
336 assert(p
->x
.p
->type
== polynomial
);
337 assert(p
->x
.p
->size
== 2);
338 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
339 value_lcm(*d
, *d
, p
->x
.p
->arr
[1].d
);
345 static void poly_denom(evalue
*p
, Value
*d
)
347 poly_denom_not_constant(&p
, d
);
348 value_lcm(*d
, *d
, p
->d
);
351 static void realloc_substitution(struct subst
*s
, int d
)
353 struct fixed_param
*n
;
356 for (i
= 0; i
< s
->n
; ++i
)
363 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
369 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
372 /* May have been reduced already */
373 if (value_notzero_p(m
->d
))
376 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
377 assert(m
->x
.p
->size
== 3);
379 /* fractional was inverted during reduction
380 * invert it back and move constant in
382 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
383 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
384 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
385 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
386 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
387 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
388 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
389 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
390 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
391 value_set_si(m
->x
.p
->arr
[1].d
, 1);
394 /* Oops. Nested identical relations. */
395 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
398 if (s
->n
>= s
->max
) {
399 int d
= relations_depth(r
);
400 realloc_substitution(s
, d
);
404 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
405 assert(p
->x
.p
->size
== 2);
408 assert(value_pos_p(f
->x
.n
));
410 value_init(s
->fixed
[s
->n
].m
);
411 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
412 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
413 value_init(s
->fixed
[s
->n
].d
);
414 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
415 value_init(s
->fixed
[s
->n
].s
.d
);
416 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
422 static int type_offset(enode
*p
)
424 return p
->type
== fractional
? 1 :
425 p
->type
== flooring
? 1 :
426 p
->type
== relation
? 1 : 0;
429 static void reorder_terms_about(enode
*p
, evalue
*v
)
432 int offset
= type_offset(p
);
434 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
436 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
437 free_evalue_refs(&(p
->arr
[i
]));
443 void evalue_reorder_terms(evalue
*e
)
449 assert(value_zero_p(e
->d
));
451 assert(p
->type
== fractional
||
452 p
->type
== flooring
||
453 p
->type
== polynomial
); /* for now */
455 offset
= type_offset(p
);
457 value_set_si(f
.d
, 0);
458 f
.x
.p
= new_enode(p
->type
, offset
+2, p
->pos
);
460 value_clear(f
.x
.p
->arr
[0].d
);
461 f
.x
.p
->arr
[0] = p
->arr
[0];
463 evalue_set_si(&f
.x
.p
->arr
[offset
], 0, 1);
464 evalue_set_si(&f
.x
.p
->arr
[offset
+1], 1, 1);
465 reorder_terms_about(p
, &f
);
471 static void evalue_reduce_size(evalue
*e
)
475 assert(value_zero_p(e
->d
));
478 offset
= type_offset(p
);
480 /* Try to reduce the degree */
481 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
482 if (!EVALUE_IS_ZERO(p
->arr
[i
]))
484 free_evalue_refs(&p
->arr
[i
]);
489 /* Try to reduce its strength */
490 if (p
->type
== relation
) {
492 free_evalue_refs(&p
->arr
[0]);
493 evalue_set_si(e
, 0, 1);
496 } else if (p
->size
== offset
+1) {
498 memcpy(e
, &p
->arr
[offset
], sizeof(evalue
));
500 free_evalue_refs(&p
->arr
[0]);
505 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
511 if (value_notzero_p(e
->d
)) {
513 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
514 return; /* a rational number, its already reduced */
518 return; /* hum... an overflow probably occured */
520 /* First reduce the components of p */
521 add
= p
->type
== relation
;
522 for (i
=0; i
<p
->size
; i
++) {
524 add
= add_modulo_substitution(s
, e
);
526 if (i
== 0 && p
->type
==fractional
)
527 _reduce_evalue(&p
->arr
[i
], s
, 1);
529 _reduce_evalue(&p
->arr
[i
], s
, fract
);
531 if (add
&& i
== p
->size
-1) {
533 value_clear(s
->fixed
[s
->n
].m
);
534 value_clear(s
->fixed
[s
->n
].d
);
535 free_evalue_refs(&s
->fixed
[s
->n
].s
);
536 } else if (add
&& i
== 1)
537 s
->fixed
[s
->n
-1].pos
*= -1;
540 if (p
->type
==periodic
) {
542 /* Try to reduce the period */
543 for (i
=1; i
<=(p
->size
)/2; i
++) {
544 if ((p
->size
% i
)==0) {
546 /* Can we reduce the size to i ? */
548 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
549 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
552 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
556 you_lose
: /* OK, lets not do it */
561 /* Try to reduce its strength */
564 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
568 else if (p
->type
==polynomial
) {
569 for (k
= 0; s
&& k
< s
->n
; ++k
) {
570 if (s
->fixed
[k
].pos
== p
->pos
) {
571 int divide
= value_notone_p(s
->fixed
[k
].d
);
574 if (value_notzero_p(s
->fixed
[k
].m
)) {
577 assert(p
->size
== 2);
578 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
580 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
587 value_assign(d
.d
, s
->fixed
[k
].d
);
589 if (value_notzero_p(s
->fixed
[k
].m
))
590 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
592 value_set_si(d
.x
.n
, 1);
595 for (i
=p
->size
-1;i
>=1;i
--) {
596 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
598 emul(&d
, &p
->arr
[i
]);
599 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
600 free_evalue_refs(&(p
->arr
[i
]));
603 _reduce_evalue(&p
->arr
[0], s
, fract
);
606 free_evalue_refs(&d
);
612 evalue_reduce_size(e
);
614 else if (p
->type
==fractional
) {
618 if (value_notzero_p(p
->arr
[0].d
)) {
620 value_assign(v
.d
, p
->arr
[0].d
);
622 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
627 evalue
*pp
= &p
->arr
[0];
628 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
629 assert(pp
->x
.p
->size
== 2);
631 /* search for exact duplicate among the modulo inequalities */
633 f
= &pp
->x
.p
->arr
[1];
634 for (k
= 0; s
&& k
< s
->n
; ++k
) {
635 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
636 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
637 value_eq(s
->fixed
[k
].m
, f
->d
) &&
638 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
645 /* replace { E/m } by { (E-1)/m } + 1/m */
650 evalue_set_si(&extra
, 1, 1);
651 value_assign(extra
.d
, g
);
652 eadd(&extra
, &v
.x
.p
->arr
[1]);
653 free_evalue_refs(&extra
);
655 /* We've been going in circles; stop now */
656 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
657 free_evalue_refs(&v
);
659 evalue_set_si(&v
, 0, 1);
664 value_set_si(v
.d
, 0);
665 v
.x
.p
= new_enode(fractional
, 3, -1);
666 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
667 value_assign(v
.x
.p
->arr
[1].d
, g
);
668 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
669 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
672 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
675 value_division(f
->d
, g
, f
->d
);
676 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
677 value_assign(f
->d
, g
);
678 value_decrement(f
->x
.n
, f
->x
.n
);
679 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
681 value_gcd(g
, f
->d
, f
->x
.n
);
682 value_division(f
->d
, f
->d
, g
);
683 value_division(f
->x
.n
, f
->x
.n
, g
);
692 /* reduction may have made this fractional arg smaller */
693 i
= reorder
? p
->size
: 1;
694 for ( ; i
< p
->size
; ++i
)
695 if (value_zero_p(p
->arr
[i
].d
) &&
696 p
->arr
[i
].x
.p
->type
== fractional
&&
697 mod_term_cmp(e
, &p
->arr
[i
]) >= 0)
701 value_set_si(v
.d
, 0);
702 v
.x
.p
= new_enode(fractional
, 3, -1);
703 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
704 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
705 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
713 evalue
*pp
= &p
->arr
[0];
716 poly_denom_not_constant(&pp
, &m
);
717 mpz_fdiv_r(r
, m
, pp
->d
);
718 if (value_notzero_p(r
)) {
720 value_set_si(v
.d
, 0);
721 v
.x
.p
= new_enode(fractional
, 3, -1);
723 value_multiply(r
, m
, pp
->x
.n
);
724 value_multiply(v
.x
.p
->arr
[1].d
, m
, pp
->d
);
725 value_init(v
.x
.p
->arr
[1].x
.n
);
726 mpz_fdiv_r(v
.x
.p
->arr
[1].x
.n
, r
, pp
->d
);
727 mpz_fdiv_q(r
, r
, pp
->d
);
729 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
730 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
732 while (value_zero_p(pp
->d
))
733 pp
= &pp
->x
.p
->arr
[0];
735 value_assign(pp
->d
, m
);
736 value_assign(pp
->x
.n
, r
);
738 value_gcd(r
, pp
->d
, pp
->x
.n
);
739 value_division(pp
->d
, pp
->d
, r
);
740 value_division(pp
->x
.n
, pp
->x
.n
, r
);
753 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
754 pp
= &pp
->x
.p
->arr
[0]) {
755 f
= &pp
->x
.p
->arr
[1];
756 assert(value_pos_p(f
->d
));
757 mpz_mul_ui(twice
, f
->x
.n
, 2);
758 if (value_lt(twice
, f
->d
))
760 if (value_eq(twice
, f
->d
))
768 value_set_si(v
.d
, 0);
769 v
.x
.p
= new_enode(fractional
, 3, -1);
770 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
771 poly_denom(&p
->arr
[0], &twice
);
772 value_assign(v
.x
.p
->arr
[1].d
, twice
);
773 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
774 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
775 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
777 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
778 pp
= &pp
->x
.p
->arr
[0]) {
779 f
= &pp
->x
.p
->arr
[1];
780 value_oppose(f
->x
.n
, f
->x
.n
);
781 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
783 value_division(pp
->d
, twice
, pp
->d
);
784 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
785 value_assign(pp
->d
, twice
);
786 value_oppose(pp
->x
.n
, pp
->x
.n
);
787 value_decrement(pp
->x
.n
, pp
->x
.n
);
788 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
790 /* Maybe we should do this during reduction of
793 value_gcd(twice
, pp
->d
, pp
->x
.n
);
794 value_division(pp
->d
, pp
->d
, twice
);
795 value_division(pp
->x
.n
, pp
->x
.n
, twice
);
805 reorder_terms_about(p
, &v
);
806 _reduce_evalue(&p
->arr
[1], s
, fract
);
809 evalue_reduce_size(e
);
811 else if (p
->type
== flooring
) {
812 /* Replace floor(constant) by its value */
813 if (value_notzero_p(p
->arr
[0].d
)) {
816 value_set_si(v
.d
, 1);
818 mpz_fdiv_q(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
819 reorder_terms_about(p
, &v
);
820 _reduce_evalue(&p
->arr
[1], s
, fract
);
822 evalue_reduce_size(e
);
824 else if (p
->type
== relation
) {
825 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
826 free_evalue_refs(&(p
->arr
[2]));
827 free_evalue_refs(&(p
->arr
[0]));
834 evalue_reduce_size(e
);
835 if (value_notzero_p(e
->d
) || p
!= e
->x
.p
)
842 /* Relation was reduced by means of an identical
843 * inequality => remove
845 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
848 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
849 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
851 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
853 free_evalue_refs(&(p
->arr
[2]));
857 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
859 evalue_set_si(e
, 0, 1);
860 free_evalue_refs(&(p
->arr
[1]));
862 free_evalue_refs(&(p
->arr
[0]));
868 } /* reduce_evalue */
870 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
875 for (k
= 0; k
< dim
; ++k
)
876 if (value_notzero_p(row
[k
+1]))
879 Vector_Normalize_Positive(row
+1, dim
+1, k
);
880 assert(s
->n
< s
->max
);
881 value_init(s
->fixed
[s
->n
].d
);
882 value_init(s
->fixed
[s
->n
].m
);
883 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
884 s
->fixed
[s
->n
].pos
= k
+1;
885 value_set_si(s
->fixed
[s
->n
].m
, 0);
886 r
= &s
->fixed
[s
->n
].s
;
888 for (l
= k
+1; l
< dim
; ++l
)
889 if (value_notzero_p(row
[l
+1])) {
890 value_set_si(r
->d
, 0);
891 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
892 value_init(r
->x
.p
->arr
[1].x
.n
);
893 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
894 value_set_si(r
->x
.p
->arr
[1].d
, 1);
898 value_oppose(r
->x
.n
, row
[dim
+1]);
899 value_set_si(r
->d
, 1);
903 static void _reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
, struct subst
*s
)
906 Polyhedron
*orig
= D
;
911 D
= DomainConvex(D
, 0);
912 /* We don't perform any substitutions if the domain is a union.
913 * We may therefore miss out on some possible simplifications,
914 * e.g., if a variable is always even in the whole union,
915 * while there is a relation in the evalue that evaluates
916 * to zero for even values of the variable.
918 if (!D
->next
&& D
->NbEq
) {
922 realloc_substitution(s
, dim
);
924 int d
= relations_depth(e
);
926 NALLOC(s
->fixed
, s
->max
);
929 for (j
= 0; j
< D
->NbEq
; ++j
)
930 add_substitution(s
, D
->Constraint
[j
], dim
);
934 _reduce_evalue(e
, s
, 0);
937 for (j
= 0; j
< s
->n
; ++j
) {
938 value_clear(s
->fixed
[j
].d
);
939 value_clear(s
->fixed
[j
].m
);
940 free_evalue_refs(&s
->fixed
[j
].s
);
945 void reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
)
947 struct subst s
= { NULL
, 0, 0 };
948 POL_ENSURE_VERTICES(D
);
950 if (EVALUE_IS_ZERO(*e
))
954 evalue_set_si(e
, 0, 1);
957 _reduce_evalue_in_domain(e
, D
, &s
);
962 void reduce_evalue (evalue
*e
) {
963 struct subst s
= { NULL
, 0, 0 };
965 if (value_notzero_p(e
->d
))
966 return; /* a rational number, its already reduced */
968 if (e
->x
.p
->type
== partition
) {
971 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
972 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
974 /* This shouldn't really happen;
975 * Empty domains should not be added.
977 POL_ENSURE_VERTICES(D
);
979 _reduce_evalue_in_domain(&e
->x
.p
->arr
[2*i
+1], D
, &s
);
981 if (emptyQ(D
) || EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
982 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
983 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
984 value_clear(e
->x
.p
->arr
[2*i
].d
);
986 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
987 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
991 if (e
->x
.p
->size
== 0) {
993 evalue_set_si(e
, 0, 1);
996 _reduce_evalue(e
, &s
, 0);
1001 static void print_evalue_r(FILE *DST
, const evalue
*e
, const char **pname
)
1003 if (EVALUE_IS_NAN(*e
)) {
1004 fprintf(DST
, "NaN");
1008 if(value_notzero_p(e
->d
)) {
1009 if(value_notone_p(e
->d
)) {
1010 value_print(DST
,VALUE_FMT
,e
->x
.n
);
1012 value_print(DST
,VALUE_FMT
,e
->d
);
1015 value_print(DST
,VALUE_FMT
,e
->x
.n
);
1019 print_enode(DST
,e
->x
.p
,pname
);
1021 } /* print_evalue */
1023 void print_evalue(FILE *DST
, const evalue
*e
, const char **pname
)
1025 print_evalue_r(DST
, e
, pname
);
1026 if (value_notzero_p(e
->d
))
1030 void print_enode(FILE *DST
, enode
*p
, const char **pname
)
1035 fprintf(DST
, "NULL");
1041 for (i
=0; i
<p
->size
; i
++) {
1042 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1046 fprintf(DST
, " }\n");
1050 for (i
=p
->size
-1; i
>=0; i
--) {
1051 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1052 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
1054 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
1056 fprintf(DST
, " )\n");
1060 for (i
=0; i
<p
->size
; i
++) {
1061 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1062 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
1064 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
1069 for (i
=p
->size
-1; i
>=1; i
--) {
1070 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1072 fprintf(DST
, " * ");
1073 fprintf(DST
, p
->type
== flooring
? "[" : "{");
1074 print_evalue_r(DST
, &p
->arr
[0], pname
);
1075 fprintf(DST
, p
->type
== flooring
? "]" : "}");
1077 fprintf(DST
, "^%d + ", i
-1);
1079 fprintf(DST
, " + ");
1082 fprintf(DST
, " )\n");
1086 print_evalue_r(DST
, &p
->arr
[0], pname
);
1087 fprintf(DST
, "= 0 ] * \n");
1088 print_evalue_r(DST
, &p
->arr
[1], pname
);
1090 fprintf(DST
, " +\n [ ");
1091 print_evalue_r(DST
, &p
->arr
[0], pname
);
1092 fprintf(DST
, "!= 0 ] * \n");
1093 print_evalue_r(DST
, &p
->arr
[2], pname
);
1097 char **new_names
= NULL
;
1098 const char **names
= pname
;
1099 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
1100 if (!pname
|| p
->pos
< maxdim
) {
1101 new_names
= ALLOCN(char *, maxdim
);
1102 for (i
= 0; i
< p
->pos
; ++i
) {
1104 new_names
[i
] = (char *)pname
[i
];
1106 new_names
[i
] = ALLOCN(char, 10);
1107 snprintf(new_names
[i
], 10, "%c", 'P'+i
);
1110 for ( ; i
< maxdim
; ++i
) {
1111 new_names
[i
] = ALLOCN(char, 10);
1112 snprintf(new_names
[i
], 10, "_p%d", i
);
1114 names
= (const char**)new_names
;
1117 for (i
=0; i
<p
->size
/2; i
++) {
1118 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
1119 print_evalue_r(DST
, &p
->arr
[2*i
+1], names
);
1120 if (value_notzero_p(p
->arr
[2*i
+1].d
))
1124 if (!pname
|| p
->pos
< maxdim
) {
1125 for (i
= pname
? p
->pos
: 0; i
< maxdim
; ++i
)
1139 * 0 if toplevels of e1 and e2 are at the same level
1140 * <0 if toplevel of e1 should be outside of toplevel of e2
1141 * >0 if toplevel of e2 should be outside of toplevel of e1
1143 static int evalue_level_cmp(const evalue
*e1
, const evalue
*e2
)
1145 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
))
1147 if (value_notzero_p(e1
->d
))
1149 if (value_notzero_p(e2
->d
))
1151 if (e1
->x
.p
->type
== partition
&& e2
->x
.p
->type
== partition
)
1153 if (e1
->x
.p
->type
== partition
)
1155 if (e2
->x
.p
->type
== partition
)
1157 if (e1
->x
.p
->type
== relation
&& e2
->x
.p
->type
== relation
) {
1158 if (eequal(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]))
1160 return mod_term_cmp(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
1162 if (e1
->x
.p
->type
== relation
)
1164 if (e2
->x
.p
->type
== relation
)
1166 if (e1
->x
.p
->type
== polynomial
&& e2
->x
.p
->type
== polynomial
)
1167 return e1
->x
.p
->pos
- e2
->x
.p
->pos
;
1168 if (e1
->x
.p
->type
== polynomial
)
1170 if (e2
->x
.p
->type
== polynomial
)
1172 if (e1
->x
.p
->type
== periodic
&& e2
->x
.p
->type
== periodic
)
1173 return e1
->x
.p
->pos
- e2
->x
.p
->pos
;
1174 assert(e1
->x
.p
->type
!= periodic
);
1175 assert(e2
->x
.p
->type
!= periodic
);
1176 assert(e1
->x
.p
->type
== e2
->x
.p
->type
);
1177 if (eequal(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]))
1179 return mod_term_cmp(e1
, e2
);
1182 static void eadd_rev(const evalue
*e1
, evalue
*res
)
1186 evalue_copy(&ev
, e1
);
1188 free_evalue_refs(res
);
1192 static void eadd_rev_cst(const evalue
*e1
, evalue
*res
)
1196 evalue_copy(&ev
, e1
);
1197 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
1198 free_evalue_refs(res
);
1202 struct section
{ Polyhedron
* D
; evalue E
; };
1204 void eadd_partitions(const evalue
*e1
, evalue
*res
)
1209 s
= (struct section
*)
1210 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
1211 sizeof(struct section
));
1213 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1214 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1215 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1218 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1219 assert(res
->x
.p
->size
>= 2);
1220 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1221 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
1223 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
1225 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1230 fd
= DomainConstraintSimplify(fd
, 0);
1235 value_init(s
[n
].E
.d
);
1236 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
1240 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1241 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
1242 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1244 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1245 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1246 d
= DomainConstraintSimplify(d
, 0);
1252 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
1253 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1255 value_init(s
[n
].E
.d
);
1256 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1257 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1262 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
1266 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1269 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1270 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1271 value_clear(res
->x
.p
->arr
[2*i
].d
);
1276 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1277 for (j
= 0; j
< n
; ++j
) {
1278 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1279 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1280 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1286 static void explicit_complement(evalue
*res
)
1288 enode
*rel
= new_enode(relation
, 3, 0);
1290 value_clear(rel
->arr
[0].d
);
1291 rel
->arr
[0] = res
->x
.p
->arr
[0];
1292 value_clear(rel
->arr
[1].d
);
1293 rel
->arr
[1] = res
->x
.p
->arr
[1];
1294 value_set_si(rel
->arr
[2].d
, 1);
1295 value_init(rel
->arr
[2].x
.n
);
1296 value_set_si(rel
->arr
[2].x
.n
, 0);
1301 static void reduce_constant(evalue
*e
)
1306 value_gcd(g
, e
->x
.n
, e
->d
);
1307 if (value_notone_p(g
)) {
1308 value_division(e
->d
, e
->d
,g
);
1309 value_division(e
->x
.n
, e
->x
.n
,g
);
1314 /* Add two rational numbers */
1315 static void eadd_rationals(const evalue
*e1
, evalue
*res
)
1317 if (value_eq(e1
->d
, res
->d
))
1318 value_addto(res
->x
.n
, res
->x
.n
, e1
->x
.n
);
1320 value_multiply(res
->x
.n
, res
->x
.n
, e1
->d
);
1321 value_addmul(res
->x
.n
, e1
->x
.n
, res
->d
);
1322 value_multiply(res
->d
,e1
->d
,res
->d
);
1324 reduce_constant(res
);
1327 static void eadd_relations(const evalue
*e1
, evalue
*res
)
1331 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1332 explicit_complement(res
);
1333 for (i
= 1; i
< e1
->x
.p
->size
; ++i
)
1334 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1337 static void eadd_arrays(const evalue
*e1
, evalue
*res
, int n
)
1341 // add any element in e1 to the corresponding element in res
1342 i
= type_offset(res
->x
.p
);
1344 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1346 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1349 static void eadd_poly(const evalue
*e1
, evalue
*res
)
1351 if (e1
->x
.p
->size
> res
->x
.p
->size
)
1354 eadd_arrays(e1
, res
, e1
->x
.p
->size
);
1358 * Product or sum of two periodics of the same parameter
1359 * and different periods
1361 static void combine_periodics(const evalue
*e1
, evalue
*res
,
1362 void (*op
)(const evalue
*, evalue
*))
1370 value_set_si(es
, e1
->x
.p
->size
);
1371 value_set_si(rs
, res
->x
.p
->size
);
1372 value_lcm(rs
, es
, rs
);
1373 size
= (int)mpz_get_si(rs
);
1376 p
= new_enode(periodic
, size
, e1
->x
.p
->pos
);
1377 for (i
= 0; i
< res
->x
.p
->size
; i
++) {
1378 value_clear(p
->arr
[i
].d
);
1379 p
->arr
[i
] = res
->x
.p
->arr
[i
];
1381 for (i
= res
->x
.p
->size
; i
< size
; i
++)
1382 evalue_copy(&p
->arr
[i
], &res
->x
.p
->arr
[i
% res
->x
.p
->size
]);
1383 for (i
= 0; i
< size
; i
++)
1384 op(&e1
->x
.p
->arr
[i
% e1
->x
.p
->size
], &p
->arr
[i
]);
1389 static void eadd_periodics(const evalue
*e1
, evalue
*res
)
1395 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1396 eadd_arrays(e1
, res
, e1
->x
.p
->size
);
1400 combine_periodics(e1
, res
, eadd
);
1403 void evalue_assign(evalue
*dst
, const evalue
*src
)
1405 if (value_pos_p(dst
->d
) && value_pos_p(src
->d
)) {
1406 value_assign(dst
->d
, src
->d
);
1407 value_assign(dst
->x
.n
, src
->x
.n
);
1410 free_evalue_refs(dst
);
1412 evalue_copy(dst
, src
);
1415 void eadd(const evalue
*e1
, evalue
*res
)
1419 if (EVALUE_IS_ZERO(*e1
))
1422 if (EVALUE_IS_NAN(*res
))
1425 if (EVALUE_IS_NAN(*e1
)) {
1426 evalue_assign(res
, e1
);
1430 if (EVALUE_IS_ZERO(*res
)) {
1431 evalue_assign(res
, e1
);
1435 cmp
= evalue_level_cmp(res
, e1
);
1437 switch (e1
->x
.p
->type
) {
1441 eadd_rev_cst(e1
, res
);
1446 } else if (cmp
== 0) {
1447 if (value_notzero_p(e1
->d
)) {
1448 eadd_rationals(e1
, res
);
1450 switch (e1
->x
.p
->type
) {
1452 eadd_partitions(e1
, res
);
1455 eadd_relations(e1
, res
);
1458 assert(e1
->x
.p
->size
== res
->x
.p
->size
);
1465 eadd_periodics(e1
, res
);
1473 switch (res
->x
.p
->type
) {
1477 /* Add to the constant term of a polynomial */
1478 eadd(e1
, &res
->x
.p
->arr
[type_offset(res
->x
.p
)]);
1481 /* Add to all elements of a periodic number */
1482 for (i
= 0; i
< res
->x
.p
->size
; i
++)
1483 eadd(e1
, &res
->x
.p
->arr
[i
]);
1486 fprintf(stderr
, "eadd: cannot add const with vector\n");
1491 /* Create (zero) complement if needed */
1492 if (res
->x
.p
->size
< 3)
1493 explicit_complement(res
);
1494 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1495 eadd(e1
, &res
->x
.p
->arr
[i
]);
1503 static void emul_rev(const evalue
*e1
, evalue
*res
)
1507 evalue_copy(&ev
, e1
);
1509 free_evalue_refs(res
);
1513 static void emul_poly(const evalue
*e1
, evalue
*res
)
1515 int i
, j
, offset
= type_offset(res
->x
.p
);
1518 int size
= (e1
->x
.p
->size
+ res
->x
.p
->size
- offset
- 1);
1520 p
= new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1522 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1523 if (!EVALUE_IS_ZERO(e1
->x
.p
->arr
[i
]))
1526 /* special case pure power */
1527 if (i
== e1
->x
.p
->size
-1) {
1529 value_clear(p
->arr
[0].d
);
1530 p
->arr
[0] = res
->x
.p
->arr
[0];
1532 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1533 evalue_set_si(&p
->arr
[i
], 0, 1);
1534 for (i
= offset
; i
< res
->x
.p
->size
; ++i
) {
1535 value_clear(p
->arr
[i
+e1
->x
.p
->size
-offset
-1].d
);
1536 p
->arr
[i
+e1
->x
.p
->size
-offset
-1] = res
->x
.p
->arr
[i
];
1537 emul(&e1
->x
.p
->arr
[e1
->x
.p
->size
-1],
1538 &p
->arr
[i
+e1
->x
.p
->size
-offset
-1]);
1546 value_set_si(tmp
.d
,0);
1549 evalue_copy(&p
->arr
[0], &e1
->x
.p
->arr
[0]);
1550 for (i
= offset
; i
< e1
->x
.p
->size
; i
++) {
1551 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1552 emul(&res
->x
.p
->arr
[offset
], &tmp
.x
.p
->arr
[i
]);
1555 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1556 for (i
= offset
+1; i
<res
->x
.p
->size
; i
++)
1557 for (j
= offset
; j
<e1
->x
.p
->size
; j
++) {
1560 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1561 emul(&res
->x
.p
->arr
[i
], &ev
);
1562 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-offset
]);
1563 free_evalue_refs(&ev
);
1565 free_evalue_refs(res
);
1569 void emul_partitions(const evalue
*e1
, evalue
*res
)
1574 s
= (struct section
*)
1575 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1576 sizeof(struct section
));
1578 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1579 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1580 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1583 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1584 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1585 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1586 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1587 d
= DomainConstraintSimplify(d
, 0);
1593 /* This code is only needed because the partitions
1594 are not true partitions.
1596 for (k
= 0; k
< n
; ++k
) {
1597 if (DomainIncludes(s
[k
].D
, d
))
1599 if (DomainIncludes(d
, s
[k
].D
)) {
1600 Domain_Free(s
[k
].D
);
1601 free_evalue_refs(&s
[k
].E
);
1612 value_init(s
[n
].E
.d
);
1613 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1614 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1618 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1619 value_clear(res
->x
.p
->arr
[2*i
].d
);
1620 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1625 evalue_set_si(res
, 0, 1);
1627 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1628 for (j
= 0; j
< n
; ++j
) {
1629 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1630 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1631 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1638 /* Product of two rational numbers */
1639 static void emul_rationals(const evalue
*e1
, evalue
*res
)
1641 value_multiply(res
->d
, e1
->d
, res
->d
);
1642 value_multiply(res
->x
.n
, e1
->x
.n
, res
->x
.n
);
1643 reduce_constant(res
);
1646 static void emul_relations(const evalue
*e1
, evalue
*res
)
1650 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3) {
1651 free_evalue_refs(&res
->x
.p
->arr
[2]);
1654 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1655 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1658 static void emul_periodics(const evalue
*e1
, evalue
*res
)
1665 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1666 /* Product of two periodics of the same parameter and period */
1667 for (i
= 0; i
< res
->x
.p
->size
; i
++)
1668 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1672 combine_periodics(e1
, res
, emul
);
1675 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1677 static void emul_fractionals(const evalue
*e1
, evalue
*res
)
1681 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1682 if (!value_two_p(d
.d
))
1687 value_set_si(d
.x
.n
, 1);
1688 /* { x }^2 == { x }/2 */
1689 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1690 assert(e1
->x
.p
->size
== 3);
1691 assert(res
->x
.p
->size
== 3);
1693 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1695 eadd(&res
->x
.p
->arr
[1], &tmp
);
1696 emul(&e1
->x
.p
->arr
[2], &tmp
);
1697 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[1]);
1698 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[2]);
1699 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1700 free_evalue_refs(&tmp
);
1706 /* Computes the product of two evalues "e1" and "res" and puts
1707 * the result in "res". You need to make a copy of "res"
1708 * before calling this function if you still need it afterward.
1709 * The vector type of evalues is not treated here
1711 void emul(const evalue
*e1
, evalue
*res
)
1715 assert(!(value_zero_p(e1
->d
) && e1
->x
.p
->type
== evector
));
1716 assert(!(value_zero_p(res
->d
) && res
->x
.p
->type
== evector
));
1718 if (EVALUE_IS_ZERO(*res
))
1721 if (EVALUE_IS_ONE(*e1
))
1724 if (EVALUE_IS_ZERO(*e1
)) {
1725 evalue_assign(res
, e1
);
1729 if (EVALUE_IS_NAN(*res
))
1732 if (EVALUE_IS_NAN(*e1
)) {
1733 evalue_assign(res
, e1
);
1737 cmp
= evalue_level_cmp(res
, e1
);
1740 } else if (cmp
== 0) {
1741 if (value_notzero_p(e1
->d
)) {
1742 emul_rationals(e1
, res
);
1744 switch (e1
->x
.p
->type
) {
1746 emul_partitions(e1
, res
);
1749 emul_relations(e1
, res
);
1756 emul_periodics(e1
, res
);
1759 emul_fractionals(e1
, res
);
1765 switch (res
->x
.p
->type
) {
1767 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1768 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1775 for (i
= type_offset(res
->x
.p
); i
< res
->x
.p
->size
; ++i
)
1776 emul(e1
, &res
->x
.p
->arr
[i
]);
1782 /* Frees mask content ! */
1783 void emask(evalue
*mask
, evalue
*res
) {
1790 if (EVALUE_IS_ZERO(*res
)) {
1791 free_evalue_refs(mask
);
1795 assert(value_zero_p(mask
->d
));
1796 assert(mask
->x
.p
->type
== partition
);
1797 assert(value_zero_p(res
->d
));
1798 assert(res
->x
.p
->type
== partition
);
1799 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1800 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1801 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1802 pos
= res
->x
.p
->pos
;
1804 s
= (struct section
*)
1805 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1806 sizeof(struct section
));
1810 evalue_set_si(&mone
, -1, 1);
1813 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1814 assert(mask
->x
.p
->size
>= 2);
1815 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1816 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1818 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1820 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1829 value_init(s
[n
].E
.d
);
1830 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1834 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1835 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1838 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1839 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1840 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1841 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1843 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1844 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1850 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1851 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1853 value_init(s
[n
].E
.d
);
1854 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1855 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1861 /* Just ignore; this may have been previously masked off */
1863 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1867 free_evalue_refs(&mone
);
1868 free_evalue_refs(mask
);
1869 free_evalue_refs(res
);
1872 evalue_set_si(res
, 0, 1);
1874 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1875 for (j
= 0; j
< n
; ++j
) {
1876 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1877 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1878 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1885 void evalue_copy(evalue
*dst
, const evalue
*src
)
1887 value_assign(dst
->d
, src
->d
);
1888 if (EVALUE_IS_NAN(*dst
)) {
1892 if (value_pos_p(src
->d
)) {
1893 value_init(dst
->x
.n
);
1894 value_assign(dst
->x
.n
, src
->x
.n
);
1896 dst
->x
.p
= ecopy(src
->x
.p
);
1899 evalue
*evalue_dup(const evalue
*e
)
1901 evalue
*res
= ALLOC(evalue
);
1903 evalue_copy(res
, e
);
1907 enode
*new_enode(enode_type type
,int size
,int pos
) {
1913 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1916 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1920 for(i
=0; i
<size
; i
++) {
1921 value_init(res
->arr
[i
].d
);
1922 value_set_si(res
->arr
[i
].d
,0);
1923 res
->arr
[i
].x
.p
= 0;
1928 enode
*ecopy(enode
*e
) {
1933 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1934 for(i
=0;i
<e
->size
;++i
) {
1935 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1936 if(value_zero_p(res
->arr
[i
].d
))
1937 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1938 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1939 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1941 value_init(res
->arr
[i
].x
.n
);
1942 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1948 int ecmp(const evalue
*e1
, const evalue
*e2
)
1954 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
1958 value_multiply(m
, e1
->x
.n
, e2
->d
);
1959 value_multiply(m2
, e2
->x
.n
, e1
->d
);
1961 if (value_lt(m
, m2
))
1963 else if (value_gt(m
, m2
))
1973 if (value_notzero_p(e1
->d
))
1975 if (value_notzero_p(e2
->d
))
1981 if (p1
->type
!= p2
->type
)
1982 return p1
->type
- p2
->type
;
1983 if (p1
->pos
!= p2
->pos
)
1984 return p1
->pos
- p2
->pos
;
1985 if (p1
->size
!= p2
->size
)
1986 return p1
->size
- p2
->size
;
1988 for (i
= p1
->size
-1; i
>= 0; --i
)
1989 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
1995 int eequal(const evalue
*e1
, const evalue
*e2
)
2000 if (value_ne(e1
->d
,e2
->d
))
2003 if (EVALUE_IS_DOMAIN(*e1
))
2004 return PolyhedronIncludes(EVALUE_DOMAIN(*e2
), EVALUE_DOMAIN(*e1
)) &&
2005 PolyhedronIncludes(EVALUE_DOMAIN(*e1
), EVALUE_DOMAIN(*e2
));
2007 if (EVALUE_IS_NAN(*e1
))
2010 assert(value_posz_p(e1
->d
));
2012 /* e1->d == e2->d */
2013 if (value_notzero_p(e1
->d
)) {
2014 if (value_ne(e1
->x
.n
,e2
->x
.n
))
2017 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2021 /* e1->d == e2->d == 0 */
2024 if (p1
->type
!= p2
->type
) return 0;
2025 if (p1
->size
!= p2
->size
) return 0;
2026 if (p1
->pos
!= p2
->pos
) return 0;
2027 for (i
=0; i
<p1
->size
; i
++)
2028 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
2033 void free_evalue_refs(evalue
*e
) {
2038 if (EVALUE_IS_NAN(*e
)) {
2043 if (EVALUE_IS_DOMAIN(*e
)) {
2044 Domain_Free(EVALUE_DOMAIN(*e
));
2047 } else if (value_pos_p(e
->d
)) {
2049 /* 'e' stores a constant */
2051 value_clear(e
->x
.n
);
2054 assert(value_zero_p(e
->d
));
2057 if (!p
) return; /* null pointer */
2058 for (i
=0; i
<p
->size
; i
++) {
2059 free_evalue_refs(&(p
->arr
[i
]));
2063 } /* free_evalue_refs */
2065 void evalue_free(evalue
*e
)
2067 free_evalue_refs(e
);
2071 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
2072 Vector
* val
, evalue
*res
)
2074 unsigned nparam
= periods
->Size
;
2077 double d
= compute_evalue(e
, val
->p
);
2078 d
*= VALUE_TO_DOUBLE(m
);
2083 value_assign(res
->d
, m
);
2084 value_init(res
->x
.n
);
2085 value_set_double(res
->x
.n
, d
);
2086 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
2089 if (value_one_p(periods
->p
[p
]))
2090 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
2095 value_assign(tmp
, periods
->p
[p
]);
2096 value_set_si(res
->d
, 0);
2097 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
2099 value_decrement(tmp
, tmp
);
2100 value_assign(val
->p
[p
], tmp
);
2101 mod2table_r(e
, periods
, m
, p
+1, val
,
2102 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
2103 } while (value_pos_p(tmp
));
2109 static void rel2table(evalue
*e
, int zero
)
2111 if (value_pos_p(e
->d
)) {
2112 if (value_zero_p(e
->x
.n
) == zero
)
2113 value_set_si(e
->x
.n
, 1);
2115 value_set_si(e
->x
.n
, 0);
2116 value_set_si(e
->d
, 1);
2119 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
2120 rel2table(&e
->x
.p
->arr
[i
], zero
);
2124 void evalue_mod2table(evalue
*e
, int nparam
)
2129 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2132 for (i
=0; i
<p
->size
; i
++) {
2133 evalue_mod2table(&(p
->arr
[i
]), nparam
);
2135 if (p
->type
== relation
) {
2140 evalue_copy(©
, &p
->arr
[0]);
2142 rel2table(&p
->arr
[0], 1);
2143 emul(&p
->arr
[0], &p
->arr
[1]);
2145 rel2table(©
, 0);
2146 emul(©
, &p
->arr
[2]);
2147 eadd(&p
->arr
[2], &p
->arr
[1]);
2148 free_evalue_refs(&p
->arr
[2]);
2149 free_evalue_refs(©
);
2151 free_evalue_refs(&p
->arr
[0]);
2155 } else if (p
->type
== fractional
) {
2156 Vector
*periods
= Vector_Alloc(nparam
);
2157 Vector
*val
= Vector_Alloc(nparam
);
2163 value_set_si(tmp
, 1);
2164 Vector_Set(periods
->p
, 1, nparam
);
2165 Vector_Set(val
->p
, 0, nparam
);
2166 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
2169 assert(p
->type
== polynomial
);
2170 assert(p
->size
== 2);
2171 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
2172 value_lcm(tmp
, tmp
, p
->arr
[1].d
);
2174 value_lcm(tmp
, tmp
, ev
->d
);
2176 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
2179 evalue_set_si(&res
, 0, 1);
2180 /* Compute the polynomial using Horner's rule */
2181 for (i
=p
->size
-1;i
>1;i
--) {
2182 eadd(&p
->arr
[i
], &res
);
2185 eadd(&p
->arr
[1], &res
);
2187 free_evalue_refs(e
);
2188 free_evalue_refs(&EP
);
2193 Vector_Free(periods
);
2195 } /* evalue_mod2table */
2197 /********************************************************/
2198 /* function in domain */
2199 /* check if the parameters in list_args */
2200 /* verifies the constraints of Domain P */
2201 /********************************************************/
2202 int in_domain(Polyhedron
*P
, Value
*list_args
)
2205 Value v
; /* value of the constraint of a row when
2206 parameters are instantiated*/
2208 if (P
->Dimension
== 0)
2213 for (row
= 0; row
< P
->NbConstraints
; row
++) {
2214 Inner_Product(P
->Constraint
[row
]+1, list_args
, P
->Dimension
, &v
);
2215 value_addto(v
, v
, P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
2216 if (value_neg_p(v
) ||
2217 value_zero_p(P
->Constraint
[row
][0]) && value_notzero_p(v
)) {
2224 return in
|| (P
->next
&& in_domain(P
->next
, list_args
));
2227 /****************************************************/
2228 /* function compute enode */
2229 /* compute the value of enode p with parameters */
2230 /* list "list_args */
2231 /* compute the polynomial or the periodic */
2232 /****************************************************/
2234 static double compute_enode(enode
*p
, Value
*list_args
) {
2246 if (p
->type
== polynomial
) {
2248 value_assign(param
,list_args
[p
->pos
-1]);
2250 /* Compute the polynomial using Horner's rule */
2251 for (i
=p
->size
-1;i
>0;i
--) {
2252 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2253 res
*=VALUE_TO_DOUBLE(param
);
2255 res
+=compute_evalue(&p
->arr
[0],list_args
);
2257 else if (p
->type
== fractional
) {
2258 double d
= compute_evalue(&p
->arr
[0], list_args
);
2259 d
-= floor(d
+1e-10);
2261 /* Compute the polynomial using Horner's rule */
2262 for (i
=p
->size
-1;i
>1;i
--) {
2263 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2266 res
+=compute_evalue(&p
->arr
[1],list_args
);
2268 else if (p
->type
== flooring
) {
2269 double d
= compute_evalue(&p
->arr
[0], list_args
);
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
== periodic
) {
2280 value_assign(m
,list_args
[p
->pos
-1]);
2282 /* Choose the right element of the periodic */
2283 value_set_si(param
,p
->size
);
2284 value_pmodulus(m
,m
,param
);
2285 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2287 else if (p
->type
== relation
) {
2288 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2289 res
= compute_evalue(&p
->arr
[1], list_args
);
2290 else if (p
->size
> 2)
2291 res
= compute_evalue(&p
->arr
[2], list_args
);
2293 else if (p
->type
== partition
) {
2294 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2295 Value
*vals
= list_args
;
2298 for (i
= 0; i
< dim
; ++i
) {
2299 value_init(vals
[i
]);
2301 value_assign(vals
[i
], list_args
[i
]);
2304 for (i
= 0; i
< p
->size
/2; ++i
)
2305 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2306 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2310 for (i
= 0; i
< dim
; ++i
)
2311 value_clear(vals
[i
]);
2320 } /* compute_enode */
2322 /*************************************************/
2323 /* return the value of Ehrhart Polynomial */
2324 /* It returns a double, because since it is */
2325 /* a recursive function, some intermediate value */
2326 /* might not be integral */
2327 /*************************************************/
2329 double compute_evalue(const evalue
*e
, Value
*list_args
)
2333 if (value_notzero_p(e
->d
)) {
2334 if (value_notone_p(e
->d
))
2335 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2337 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2340 res
= compute_enode(e
->x
.p
,list_args
);
2342 } /* compute_evalue */
2345 /****************************************************/
2346 /* function compute_poly : */
2347 /* Check for the good validity domain */
2348 /* return the number of point in the Polyhedron */
2349 /* in allocated memory */
2350 /* Using the Ehrhart pseudo-polynomial */
2351 /****************************************************/
2352 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2355 /* double d; int i; */
2357 tmp
= (Value
*) malloc (sizeof(Value
));
2358 assert(tmp
!= NULL
);
2360 value_set_si(*tmp
,0);
2363 return(tmp
); /* no ehrhart polynomial */
2364 if(en
->ValidityDomain
) {
2365 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2366 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2371 return(tmp
); /* no Validity Domain */
2373 if(in_domain(en
->ValidityDomain
,list_args
)) {
2375 #ifdef EVAL_EHRHART_DEBUG
2376 Print_Domain(stdout
,en
->ValidityDomain
);
2377 print_evalue(stdout
,&en
->EP
);
2380 /* d = compute_evalue(&en->EP,list_args);
2382 printf("(double)%lf = %d\n", d, i ); */
2383 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2389 value_set_si(*tmp
,0);
2390 return(tmp
); /* no compatible domain with the arguments */
2391 } /* compute_poly */
2393 static evalue
*eval_polynomial(const enode
*p
, int offset
,
2394 evalue
*base
, Value
*values
)
2399 res
= evalue_zero();
2400 for (i
= p
->size
-1; i
> offset
; --i
) {
2401 c
= evalue_eval(&p
->arr
[i
], values
);
2406 c
= evalue_eval(&p
->arr
[offset
], values
);
2413 evalue
*evalue_eval(const evalue
*e
, Value
*values
)
2420 if (value_notzero_p(e
->d
)) {
2421 res
= ALLOC(evalue
);
2423 evalue_copy(res
, e
);
2426 switch (e
->x
.p
->type
) {
2428 value_init(param
.x
.n
);
2429 value_assign(param
.x
.n
, values
[e
->x
.p
->pos
-1]);
2430 value_init(param
.d
);
2431 value_set_si(param
.d
, 1);
2433 res
= eval_polynomial(e
->x
.p
, 0, ¶m
, values
);
2434 free_evalue_refs(¶m
);
2437 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2438 mpz_fdiv_r(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2440 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2441 evalue_free(param2
);
2444 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2445 mpz_fdiv_q(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2446 value_set_si(param2
->d
, 1);
2448 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2449 evalue_free(param2
);
2452 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2453 if (value_zero_p(param2
->x
.n
))
2454 res
= evalue_eval(&e
->x
.p
->arr
[1], values
);
2455 else if (e
->x
.p
->size
> 2)
2456 res
= evalue_eval(&e
->x
.p
->arr
[2], values
);
2458 res
= evalue_zero();
2459 evalue_free(param2
);
2462 assert(e
->x
.p
->pos
== EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
);
2463 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2464 if (in_domain(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), values
)) {
2465 res
= evalue_eval(&e
->x
.p
->arr
[2*i
+1], values
);
2469 res
= evalue_zero();
2477 size_t value_size(Value v
) {
2478 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2479 * sizeof(v
[0]._mp_d
[0]);
2482 size_t domain_size(Polyhedron
*D
)
2485 size_t s
= sizeof(*D
);
2487 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2488 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2489 s
+= value_size(D
->Constraint
[i
][j
]);
2492 for (i = 0; i < D->NbRays; ++i)
2493 for (j = 0; j < D->Dimension+2; ++j)
2494 s += value_size(D->Ray[i][j]);
2497 return D
->next
? s
+domain_size(D
->next
) : s
;
2500 size_t enode_size(enode
*p
) {
2501 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2504 if (p
->type
== partition
)
2505 for (i
= 0; i
< p
->size
/2; ++i
) {
2506 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2507 s
+= evalue_size(&p
->arr
[2*i
+1]);
2510 for (i
= 0; i
< p
->size
; ++i
) {
2511 s
+= evalue_size(&p
->arr
[i
]);
2516 size_t evalue_size(evalue
*e
)
2518 size_t s
= sizeof(*e
);
2519 s
+= value_size(e
->d
);
2520 if (value_notzero_p(e
->d
))
2521 s
+= value_size(e
->x
.n
);
2523 s
+= enode_size(e
->x
.p
);
2527 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2529 evalue
*found
= NULL
;
2534 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2537 value_init(offset
.d
);
2538 value_init(offset
.x
.n
);
2539 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2540 value_lcm(offset
.d
, m
, offset
.d
);
2541 value_set_si(offset
.x
.n
, 1);
2544 evalue_copy(©
, cst
);
2547 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2549 if (eequal(base
, &e
->x
.p
->arr
[0]))
2550 found
= &e
->x
.p
->arr
[0];
2552 value_set_si(offset
.x
.n
, -2);
2555 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2557 if (eequal(base
, &e
->x
.p
->arr
[0]))
2560 free_evalue_refs(cst
);
2561 free_evalue_refs(&offset
);
2564 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2565 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2570 static evalue
*find_relation_pair(evalue
*e
)
2573 evalue
*found
= NULL
;
2575 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2578 if (e
->x
.p
->type
== fractional
) {
2583 poly_denom(&e
->x
.p
->arr
[0], &m
);
2585 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2586 cst
= &cst
->x
.p
->arr
[0])
2589 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2590 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2595 i
= e
->x
.p
->type
== relation
;
2596 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2597 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2602 void evalue_mod2relation(evalue
*e
) {
2605 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2608 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2609 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2610 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2611 value_clear(e
->x
.p
->arr
[2*i
].d
);
2612 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2614 if (2*i
< e
->x
.p
->size
) {
2615 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2616 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2621 if (e
->x
.p
->size
== 0) {
2623 evalue_set_si(e
, 0, 1);
2629 while ((d
= find_relation_pair(e
)) != NULL
) {
2633 value_init(split
.d
);
2634 value_set_si(split
.d
, 0);
2635 split
.x
.p
= new_enode(relation
, 3, 0);
2636 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2637 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2639 ev
= &split
.x
.p
->arr
[0];
2640 value_set_si(ev
->d
, 0);
2641 ev
->x
.p
= new_enode(fractional
, 3, -1);
2642 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2643 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2644 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2650 free_evalue_refs(&split
);
2654 static int evalue_comp(const void * a
, const void * b
)
2656 const evalue
*e1
= *(const evalue
**)a
;
2657 const evalue
*e2
= *(const evalue
**)b
;
2658 return ecmp(e1
, e2
);
2661 void evalue_combine(evalue
*e
)
2668 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2671 NALLOC(evs
, e
->x
.p
->size
/2);
2672 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2673 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2674 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2675 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2676 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2677 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2678 value_clear(p
->arr
[2*k
].d
);
2679 value_clear(p
->arr
[2*k
+1].d
);
2680 p
->arr
[2*k
] = *(evs
[i
]-1);
2681 p
->arr
[2*k
+1] = *(evs
[i
]);
2684 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2687 value_clear((evs
[i
]-1)->d
);
2691 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2692 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2693 free_evalue_refs(evs
[i
]);
2697 for (i
= 2*k
; i
< p
->size
; ++i
)
2698 value_clear(p
->arr
[i
].d
);
2705 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2707 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2709 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2712 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2713 Polyhedron
*D
, *N
, **P
;
2716 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2723 if (D
->NbEq
<= H
->NbEq
) {
2729 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2730 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2731 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2732 reduce_evalue(&tmp
);
2733 if (value_notzero_p(tmp
.d
) ||
2734 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2737 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2738 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2741 free_evalue_refs(&tmp
);
2747 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2749 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2751 value_clear(e
->x
.p
->arr
[2*i
].d
);
2752 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2754 if (2*i
< e
->x
.p
->size
) {
2755 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2756 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2763 H
= DomainConvex(D
, 0);
2764 E
= DomainDifference(H
, D
, 0);
2766 D
= DomainDifference(H
, E
, 0);
2769 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2773 /* Use smallest representative for coefficients in affine form in
2774 * argument of fractional.
2775 * Since any change will make the argument non-standard,
2776 * the containing evalue will have to be reduced again afterward.
2778 static void fractional_minimal_coefficients(enode
*p
)
2784 assert(p
->type
== fractional
);
2786 while (value_zero_p(pp
->d
)) {
2787 assert(pp
->x
.p
->type
== polynomial
);
2788 assert(pp
->x
.p
->size
== 2);
2789 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2790 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2791 if (value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2792 value_subtract(pp
->x
.p
->arr
[1].x
.n
,
2793 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2794 pp
= &pp
->x
.p
->arr
[0];
2800 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2805 unsigned dim
= D
->Dimension
;
2806 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2809 assert(p
->type
== fractional
|| p
->type
== flooring
);
2810 value_set_si(T
->p
[1][dim
], 1);
2811 evalue_extract_affine(&p
->arr
[0], T
->p
[0], &T
->p
[0][dim
], d
);
2812 I
= DomainImage(D
, T
, 0);
2813 H
= DomainConvex(I
, 0);
2823 static void replace_by_affine(evalue
*e
, Value offset
)
2830 value_init(inc
.x
.n
);
2831 value_set_si(inc
.d
, 1);
2832 value_oppose(inc
.x
.n
, offset
);
2833 eadd(&inc
, &p
->arr
[0]);
2834 reorder_terms_about(p
, &p
->arr
[0]); /* frees arr[0] */
2838 free_evalue_refs(&inc
);
2841 int evalue_range_reduction_in_domain(evalue
*e
, Polyhedron
*D
)
2850 if (value_notzero_p(e
->d
))
2855 if (p
->type
== relation
) {
2862 fractional_minimal_coefficients(p
->arr
[0].x
.p
);
2863 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
);
2864 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2865 equal
= value_eq(min
, max
);
2866 mpz_cdiv_q(min
, min
, d
);
2867 mpz_fdiv_q(max
, max
, d
);
2869 if (bounded
&& value_gt(min
, max
)) {
2875 evalue_set_si(e
, 0, 1);
2878 free_evalue_refs(&(p
->arr
[1]));
2879 free_evalue_refs(&(p
->arr
[0]));
2885 return r
? r
: evalue_range_reduction_in_domain(e
, D
);
2886 } else if (bounded
&& equal
) {
2889 free_evalue_refs(&(p
->arr
[2]));
2892 free_evalue_refs(&(p
->arr
[0]));
2898 return evalue_range_reduction_in_domain(e
, D
);
2899 } else if (bounded
&& value_eq(min
, max
)) {
2900 /* zero for a single value */
2902 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2903 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2904 value_multiply(min
, min
, d
);
2905 value_subtract(M
->p
[0][D
->Dimension
+1],
2906 M
->p
[0][D
->Dimension
+1], min
);
2907 E
= DomainAddConstraints(D
, M
, 0);
2913 r
= evalue_range_reduction_in_domain(&p
->arr
[1], E
);
2915 r
|= evalue_range_reduction_in_domain(&p
->arr
[2], D
);
2917 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2925 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2928 i
= p
->type
== relation
? 1 :
2929 p
->type
== fractional
? 1 : 0;
2930 for (; i
<p
->size
; i
++)
2931 r
|= evalue_range_reduction_in_domain(&p
->arr
[i
], D
);
2933 if (p
->type
!= fractional
) {
2934 if (r
&& p
->type
== polynomial
) {
2937 value_set_si(f
.d
, 0);
2938 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2939 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2940 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2941 reorder_terms_about(p
, &f
);
2952 fractional_minimal_coefficients(p
);
2953 I
= polynomial_projection(p
, D
, &d
, NULL
);
2954 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2955 mpz_fdiv_q(min
, min
, d
);
2956 mpz_fdiv_q(max
, max
, d
);
2957 value_subtract(d
, max
, min
);
2959 if (bounded
&& value_eq(min
, max
)) {
2960 replace_by_affine(e
, min
);
2962 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
2963 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
2964 * See pages 199-200 of PhD thesis.
2972 value_set_si(rem
.d
, 0);
2973 rem
.x
.p
= new_enode(fractional
, 3, -1);
2974 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
2975 value_clear(rem
.x
.p
->arr
[1].d
);
2976 value_clear(rem
.x
.p
->arr
[2].d
);
2977 rem
.x
.p
->arr
[1] = p
->arr
[1];
2978 rem
.x
.p
->arr
[2] = p
->arr
[2];
2979 for (i
= 3; i
< p
->size
; ++i
)
2980 p
->arr
[i
-2] = p
->arr
[i
];
2984 value_init(inc
.x
.n
);
2985 value_set_si(inc
.d
, 1);
2986 value_oppose(inc
.x
.n
, min
);
2989 evalue_copy(&t
, &p
->arr
[0]);
2993 value_set_si(f
.d
, 0);
2994 f
.x
.p
= new_enode(fractional
, 3, -1);
2995 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
2996 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2997 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
2999 value_init(factor
.d
);
3000 evalue_set_si(&factor
, -1, 1);
3006 value_clear(f
.x
.p
->arr
[1].x
.n
);
3007 value_clear(f
.x
.p
->arr
[2].x
.n
);
3008 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3009 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3013 evalue_reorder_terms(&rem
);
3014 evalue_reorder_terms(e
);
3020 free_evalue_refs(&inc
);
3021 free_evalue_refs(&t
);
3022 free_evalue_refs(&f
);
3023 free_evalue_refs(&factor
);
3024 free_evalue_refs(&rem
);
3026 evalue_range_reduction_in_domain(e
, D
);
3030 _reduce_evalue(&p
->arr
[0], 0, 1);
3032 evalue_reorder_terms(e
);
3042 void evalue_range_reduction(evalue
*e
)
3045 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3048 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3049 if (evalue_range_reduction_in_domain(&e
->x
.p
->arr
[2*i
+1],
3050 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
3051 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3053 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
3054 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
3055 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3056 value_clear(e
->x
.p
->arr
[2*i
].d
);
3058 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
3059 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
3067 Enumeration
* partition2enumeration(evalue
*EP
)
3070 Enumeration
*en
, *res
= NULL
;
3072 if (EVALUE_IS_ZERO(*EP
)) {
3077 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
3078 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
3079 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
3082 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
3083 value_clear(EP
->x
.p
->arr
[2*i
].d
);
3084 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
3092 int evalue_frac2floor_in_domain3(evalue
*e
, Polyhedron
*D
, int shift
)
3101 if (value_notzero_p(e
->d
))
3106 i
= p
->type
== relation
? 1 :
3107 p
->type
== fractional
? 1 : 0;
3108 for (; i
<p
->size
; i
++)
3109 r
|= evalue_frac2floor_in_domain3(&p
->arr
[i
], D
, shift
);
3111 if (p
->type
!= fractional
) {
3112 if (r
&& p
->type
== polynomial
) {
3115 value_set_si(f
.d
, 0);
3116 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3117 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3118 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3119 reorder_terms_about(p
, &f
);
3129 I
= polynomial_projection(p
, D
, &d
, NULL
);
3132 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3135 assert(I
->NbEq
== 0); /* Should have been reduced */
3138 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3139 if (value_pos_p(I
->Constraint
[i
][1]))
3142 if (i
< I
->NbConstraints
) {
3144 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3145 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3146 if (value_neg_p(min
)) {
3148 mpz_fdiv_q(min
, min
, d
);
3149 value_init(offset
.d
);
3150 value_set_si(offset
.d
, 1);
3151 value_init(offset
.x
.n
);
3152 value_oppose(offset
.x
.n
, min
);
3153 eadd(&offset
, &p
->arr
[0]);
3154 free_evalue_refs(&offset
);
3164 value_set_si(fl
.d
, 0);
3165 fl
.x
.p
= new_enode(flooring
, 3, -1);
3166 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
3167 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
3168 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
3170 eadd(&fl
, &p
->arr
[0]);
3171 reorder_terms_about(p
, &p
->arr
[0]);
3175 free_evalue_refs(&fl
);
3180 int evalue_frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
3182 return evalue_frac2floor_in_domain3(e
, D
, 1);
3185 void evalue_frac2floor2(evalue
*e
, int shift
)
3188 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
3190 if (evalue_frac2floor_in_domain3(e
, NULL
, 0))
3196 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3197 if (evalue_frac2floor_in_domain3(&e
->x
.p
->arr
[2*i
+1],
3198 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), shift
))
3199 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3202 void evalue_frac2floor(evalue
*e
)
3204 evalue_frac2floor2(e
, 1);
3207 /* Add a new variable with lower bound 1 and upper bound specified
3208 * by row. If negative is true, then the new variable has upper
3209 * bound -1 and lower bound specified by row.
3211 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
3212 Vector
*row
, int negative
)
3216 int nparam
= D
->Dimension
- nvar
;
3219 nr
= D
->NbConstraints
+ 2;
3220 nc
= D
->Dimension
+ 2 + 1;
3221 C
= Matrix_Alloc(nr
, nc
);
3222 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
3223 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
3224 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3225 D
->Dimension
+ 1 - nvar
);
3230 nc
= C
->NbColumns
+ 1;
3231 C
= Matrix_Alloc(nr
, nc
);
3232 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
3233 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
3234 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3235 oldC
->NbColumns
- 1 - nvar
);
3238 value_set_si(C
->p
[nr
-2][0], 1);
3240 value_set_si(C
->p
[nr
-2][1 + nvar
], -1);
3242 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
3243 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
3245 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
3246 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
3252 static void floor2frac_r(evalue
*e
, int nvar
)
3259 if (value_notzero_p(e
->d
))
3264 assert(p
->type
== flooring
);
3265 for (i
= 1; i
< p
->size
; i
++)
3266 floor2frac_r(&p
->arr
[i
], nvar
);
3268 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
3269 assert(pp
->x
.p
->type
== polynomial
);
3270 pp
->x
.p
->pos
-= nvar
;
3274 value_set_si(f
.d
, 0);
3275 f
.x
.p
= new_enode(fractional
, 3, -1);
3276 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3277 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3278 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3280 eadd(&f
, &p
->arr
[0]);
3281 reorder_terms_about(p
, &p
->arr
[0]);
3285 free_evalue_refs(&f
);
3288 /* Convert flooring back to fractional and shift position
3289 * of the parameters by nvar
3291 static void floor2frac(evalue
*e
, int nvar
)
3293 floor2frac_r(e
, nvar
);
3297 int evalue_floor2frac(evalue
*e
)
3302 if (value_notzero_p(e
->d
))
3305 if (e
->x
.p
->type
== partition
) {
3306 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3307 if (evalue_floor2frac(&e
->x
.p
->arr
[2*i
+1]))
3308 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3312 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
3313 r
|= evalue_floor2frac(&e
->x
.p
->arr
[i
]);
3315 if (e
->x
.p
->type
== flooring
) {
3321 evalue_reorder_terms(e
);
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 void domain_signs(Polyhedron
*D
, int *signs
)
3352 POL_ENSURE_VERTICES(D
);
3353 for (j
= 0; j
< D
->Dimension
; ++j
) {
3355 for (k
= 0; k
< D
->NbRays
; ++k
) {
3356 signs
[j
] = value_sign(D
->Ray
[k
][1+j
]);
3363 static evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
3364 int *signs
, Matrix
*C
, unsigned MaxRays
)
3370 evalue
*factor
= NULL
;
3374 if (EVALUE_IS_ZERO(*e
))
3378 Polyhedron
*DD
= Disjoint_Domain(D
, 0, MaxRays
);
3385 res
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3388 for (Q
= DD
; Q
; Q
= DD
) {
3394 t
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3407 if (value_notzero_p(e
->d
)) {
3410 t
= esum_over_domain_cst(nvar
, D
, C
);
3412 if (!EVALUE_IS_ONE(*e
))
3419 signs
= alloca(sizeof(int) * D
->Dimension
);
3420 domain_signs(D
, signs
);
3423 switch (e
->x
.p
->type
) {
3425 evalue
*pp
= &e
->x
.p
->arr
[0];
3427 if (pp
->x
.p
->pos
> nvar
) {
3428 /* remainder is independent of the summated vars */
3434 floor2frac(&f
, nvar
);
3436 t
= esum_over_domain_cst(nvar
, D
, C
);
3440 free_evalue_refs(&f
);
3445 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3446 poly_denom(pp
, &row
->p
[1 + nvar
]);
3447 value_set_si(row
->p
[0], 1);
3448 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3449 pp
= &pp
->x
.p
->arr
[0]) {
3451 assert(pp
->x
.p
->type
== polynomial
);
3453 if (pos
>= 1 + nvar
)
3455 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3456 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3457 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3459 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3460 value_division(row
->p
[1 + D
->Dimension
+ 1],
3461 row
->p
[1 + D
->Dimension
+ 1],
3463 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3464 row
->p
[1 + D
->Dimension
+ 1],
3466 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3470 int pos
= e
->x
.p
->pos
;
3473 factor
= ALLOC(evalue
);
3474 value_init(factor
->d
);
3475 value_set_si(factor
->d
, 0);
3476 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3477 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3478 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3482 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3483 negative
= signs
[pos
-1] < 0;
3484 value_set_si(row
->p
[0], 1);
3486 value_set_si(row
->p
[pos
], -1);
3487 value_set_si(row
->p
[1 + nvar
], 1);
3489 value_set_si(row
->p
[pos
], 1);
3490 value_set_si(row
->p
[1 + nvar
], -1);
3498 offset
= type_offset(e
->x
.p
);
3500 res
= esum_over_domain(&e
->x
.p
->arr
[offset
], nvar
, D
, signs
, C
, MaxRays
);
3504 evalue_copy(&cum
, factor
);
3508 for (i
= 1; offset
+i
< e
->x
.p
->size
; ++i
) {
3512 C
= esum_add_constraint(nvar
, D
, C
, row
, negative
);
3518 Vector_Print(stderr, P_VALUE_FMT, row);
3520 Matrix_Print(stderr, P_VALUE_FMT, C);
3522 t
= esum_over_domain(&e
->x
.p
->arr
[offset
+i
], nvar
, D
, signs
, C
, MaxRays
);
3527 if (negative
&& (i
% 2))
3537 if (factor
&& offset
+i
+1 < e
->x
.p
->size
)
3544 free_evalue_refs(&cum
);
3545 evalue_free(factor
);
3556 static Polyhedron_Insert(Polyhedron
***next
, Polyhedron
*Q
)
3566 static Polyhedron
*Polyhedron_Split_Into_Orthants(Polyhedron
*P
,
3571 Vector
*c
= Vector_Alloc(1 + P
->Dimension
+ 1);
3572 value_set_si(c
->p
[0], 1);
3574 if (P
->Dimension
== 0)
3575 return Polyhedron_Copy(P
);
3577 for (i
= 0; i
< P
->Dimension
; ++i
) {
3578 Polyhedron
*L
= NULL
;
3579 Polyhedron
**next
= &L
;
3582 for (I
= D
; I
; I
= I
->next
) {
3584 value_set_si(c
->p
[1+i
], 1);
3585 value_set_si(c
->p
[1+P
->Dimension
], 0);
3586 Q
= AddConstraints(c
->p
, 1, I
, MaxRays
);
3587 Polyhedron_Insert(&next
, Q
);
3588 value_set_si(c
->p
[1+i
], -1);
3589 value_set_si(c
->p
[1+P
->Dimension
], -1);
3590 Q
= AddConstraints(c
->p
, 1, I
, MaxRays
);
3591 Polyhedron_Insert(&next
, Q
);
3592 value_set_si(c
->p
[1+i
], 0);
3602 /* Make arguments of all floors non-negative */
3603 static void shift_floor_in_domain(evalue
*e
, Polyhedron
*D
)
3610 if (value_notzero_p(e
->d
))
3615 for (i
= type_offset(p
); i
< p
->size
; ++i
)
3616 shift_floor_in_domain(&p
->arr
[i
], D
);
3618 if (p
->type
!= flooring
)
3624 I
= polynomial_projection(p
, D
, &d
, NULL
);
3625 assert(I
->NbEq
== 0); /* Should have been reduced */
3627 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3628 if (value_pos_p(I
->Constraint
[i
][1]))
3630 assert(i
< I
->NbConstraints
);
3631 if (i
< I
->NbConstraints
) {
3632 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3633 mpz_fdiv_q(m
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3634 if (value_neg_p(m
)) {
3635 /* replace [e] by [e-m]+m such that e-m >= 0 */
3640 value_set_si(f
.d
, 1);
3641 value_oppose(f
.x
.n
, m
);
3642 eadd(&f
, &p
->arr
[0]);
3645 value_set_si(f
.d
, 0);
3646 f
.x
.p
= new_enode(flooring
, 3, -1);
3647 value_clear(f
.x
.p
->arr
[0].d
);
3648 f
.x
.p
->arr
[0] = p
->arr
[0];
3649 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
3650 value_set_si(f
.x
.p
->arr
[1].d
, 1);
3651 value_init(f
.x
.p
->arr
[1].x
.n
);
3652 value_assign(f
.x
.p
->arr
[1].x
.n
, m
);
3653 reorder_terms_about(p
, &f
);
3664 evalue
*box_summate(Polyhedron
*P
, evalue
*E
, unsigned nvar
, unsigned MaxRays
)
3666 evalue
*sum
= evalue_zero();
3667 Polyhedron
*D
= Polyhedron_Split_Into_Orthants(P
, MaxRays
);
3668 for (P
= D
; P
; P
= P
->next
) {
3670 evalue
*fe
= evalue_dup(E
);
3671 Polyhedron
*next
= P
->next
;
3673 reduce_evalue_in_domain(fe
, P
);
3674 evalue_frac2floor2(fe
, 0);
3675 shift_floor_in_domain(fe
, P
);
3676 t
= esum_over_domain(fe
, nvar
, P
, NULL
, NULL
, MaxRays
);
3688 /* Initial silly implementation */
3689 void eor(evalue
*e1
, evalue
*res
)
3695 evalue_set_si(&mone
, -1, 1);
3697 evalue_copy(&E
, res
);
3703 free_evalue_refs(&E
);
3704 free_evalue_refs(&mone
);
3707 /* computes denominator of polynomial evalue
3708 * d should point to a value initialized to 1
3710 void evalue_denom(const evalue
*e
, Value
*d
)
3714 if (value_notzero_p(e
->d
)) {
3715 value_lcm(*d
, *d
, e
->d
);
3718 assert(e
->x
.p
->type
== polynomial
||
3719 e
->x
.p
->type
== fractional
||
3720 e
->x
.p
->type
== flooring
);
3721 offset
= type_offset(e
->x
.p
);
3722 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3723 evalue_denom(&e
->x
.p
->arr
[i
], d
);
3726 /* Divides the evalue e by the integer n */
3727 void evalue_div(evalue
*e
, Value n
)
3731 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3734 if (value_notzero_p(e
->d
)) {
3737 value_multiply(e
->d
, e
->d
, n
);
3738 value_gcd(gc
, e
->x
.n
, e
->d
);
3739 if (value_notone_p(gc
)) {
3740 value_division(e
->d
, e
->d
, gc
);
3741 value_division(e
->x
.n
, e
->x
.n
, gc
);
3746 if (e
->x
.p
->type
== partition
) {
3747 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3748 evalue_div(&e
->x
.p
->arr
[2*i
+1], n
);
3751 offset
= type_offset(e
->x
.p
);
3752 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3753 evalue_div(&e
->x
.p
->arr
[i
], n
);
3756 /* Multiplies the evalue e by the integer n */
3757 void evalue_mul(evalue
*e
, Value n
)
3761 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3764 if (value_notzero_p(e
->d
)) {
3767 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3768 value_gcd(gc
, e
->x
.n
, e
->d
);
3769 if (value_notone_p(gc
)) {
3770 value_division(e
->d
, e
->d
, gc
);
3771 value_division(e
->x
.n
, e
->x
.n
, gc
);
3776 if (e
->x
.p
->type
== partition
) {
3777 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3778 evalue_mul(&e
->x
.p
->arr
[2*i
+1], n
);
3781 offset
= type_offset(e
->x
.p
);
3782 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3783 evalue_mul(&e
->x
.p
->arr
[i
], n
);
3786 /* Multiplies the evalue e by the n/d */
3787 void evalue_mul_div(evalue
*e
, Value n
, Value d
)
3791 if ((value_one_p(n
) && value_one_p(d
)) || EVALUE_IS_ZERO(*e
))
3794 if (value_notzero_p(e
->d
)) {
3797 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3798 value_multiply(e
->d
, e
->d
, d
);
3799 value_gcd(gc
, e
->x
.n
, e
->d
);
3800 if (value_notone_p(gc
)) {
3801 value_division(e
->d
, e
->d
, gc
);
3802 value_division(e
->x
.n
, e
->x
.n
, gc
);
3807 if (e
->x
.p
->type
== partition
) {
3808 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3809 evalue_mul_div(&e
->x
.p
->arr
[2*i
+1], n
, d
);
3812 offset
= type_offset(e
->x
.p
);
3813 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3814 evalue_mul_div(&e
->x
.p
->arr
[i
], n
, d
);
3817 void evalue_negate(evalue
*e
)
3821 if (value_notzero_p(e
->d
)) {
3822 value_oppose(e
->x
.n
, e
->x
.n
);
3825 if (e
->x
.p
->type
== partition
) {
3826 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3827 evalue_negate(&e
->x
.p
->arr
[2*i
+1]);
3830 offset
= type_offset(e
->x
.p
);
3831 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3832 evalue_negate(&e
->x
.p
->arr
[i
]);
3835 void evalue_add_constant(evalue
*e
, const Value cst
)
3839 if (value_zero_p(e
->d
)) {
3840 if (e
->x
.p
->type
== partition
) {
3841 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3842 evalue_add_constant(&e
->x
.p
->arr
[2*i
+1], cst
);
3845 if (e
->x
.p
->type
== relation
) {
3846 for (i
= 1; i
< e
->x
.p
->size
; ++i
)
3847 evalue_add_constant(&e
->x
.p
->arr
[i
], cst
);
3851 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
3852 } while (value_zero_p(e
->d
));
3854 value_addmul(e
->x
.n
, cst
, e
->d
);
3857 static void evalue_frac2polynomial_r(evalue
*e
, int *signs
, int sign
, int in_frac
)
3862 int sign_odd
= sign
;
3864 if (value_notzero_p(e
->d
)) {
3865 if (in_frac
&& sign
* value_sign(e
->x
.n
) < 0) {
3866 value_set_si(e
->x
.n
, 0);
3867 value_set_si(e
->d
, 1);
3872 if (e
->x
.p
->type
== relation
) {
3873 for (i
= e
->x
.p
->size
-1; i
>= 1; --i
)
3874 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
, sign
, in_frac
);
3878 if (e
->x
.p
->type
== polynomial
)
3879 sign_odd
*= signs
[e
->x
.p
->pos
-1];
3880 offset
= type_offset(e
->x
.p
);
3881 evalue_frac2polynomial_r(&e
->x
.p
->arr
[offset
], signs
, sign
, in_frac
);
3882 in_frac
|= e
->x
.p
->type
== fractional
;
3883 for (i
= e
->x
.p
->size
-1; i
> offset
; --i
)
3884 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
,
3885 (i
- offset
) % 2 ? sign_odd
: sign
, in_frac
);
3887 if (e
->x
.p
->type
!= fractional
)
3890 /* replace { a/m } by (m-1)/m if sign != 0
3891 * and by (m-1)/(2m) if sign == 0
3895 evalue_denom(&e
->x
.p
->arr
[0], &d
);
3896 free_evalue_refs(&e
->x
.p
->arr
[0]);
3897 value_init(e
->x
.p
->arr
[0].d
);
3898 value_init(e
->x
.p
->arr
[0].x
.n
);
3900 value_addto(e
->x
.p
->arr
[0].d
, d
, d
);
3902 value_assign(e
->x
.p
->arr
[0].d
, d
);
3903 value_decrement(e
->x
.p
->arr
[0].x
.n
, d
);
3907 reorder_terms_about(p
, &p
->arr
[0]);
3913 /* Approximate the evalue in fractional representation by a polynomial.
3914 * If sign > 0, the result is an upper bound;
3915 * if sign < 0, the result is a lower bound;
3916 * if sign = 0, the result is an intermediate approximation.
3918 void evalue_frac2polynomial(evalue
*e
, int sign
, unsigned MaxRays
)
3923 if (value_notzero_p(e
->d
))
3925 assert(e
->x
.p
->type
== partition
);
3926 /* make sure all variables in the domains have a fixed sign */
3928 evalue_split_domains_into_orthants(e
, MaxRays
);
3929 if (EVALUE_IS_ZERO(*e
))
3933 assert(e
->x
.p
->size
>= 2);
3934 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3936 signs
= alloca(sizeof(int) * dim
);
3939 for (i
= 0; i
< dim
; ++i
)
3941 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3943 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3944 evalue_frac2polynomial_r(&e
->x
.p
->arr
[2*i
+1], signs
, sign
, 0);
3948 /* Split the domains of e (which is assumed to be a partition)
3949 * such that each resulting domain lies entirely in one orthant.
3951 void evalue_split_domains_into_orthants(evalue
*e
, unsigned MaxRays
)
3954 assert(value_zero_p(e
->d
));
3955 assert(e
->x
.p
->type
== partition
);
3956 assert(e
->x
.p
->size
>= 2);
3957 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3959 for (i
= 0; i
< dim
; ++i
) {
3962 C
= Matrix_Alloc(1, 1 + dim
+ 1);
3963 value_set_si(C
->p
[0][0], 1);
3964 value_init(split
.d
);
3965 value_set_si(split
.d
, 0);
3966 split
.x
.p
= new_enode(partition
, 4, dim
);
3967 value_set_si(C
->p
[0][1+i
], 1);
3968 C2
= Matrix_Copy(C
);
3969 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[0], Constraints2Polyhedron(C2
, MaxRays
));
3971 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
3972 value_set_si(C
->p
[0][1+i
], -1);
3973 value_set_si(C
->p
[0][1+dim
], -1);
3974 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[2], Constraints2Polyhedron(C
, MaxRays
));
3975 evalue_set_si(&split
.x
.p
->arr
[3], 1, 1);
3977 free_evalue_refs(&split
);
3982 static evalue
*find_fractional_with_max_periods(evalue
*e
, Polyhedron
*D
,
3985 Value
*min
, Value
*max
)
3992 if (value_notzero_p(e
->d
))
3995 if (e
->x
.p
->type
== fractional
) {
4000 I
= polynomial_projection(e
->x
.p
, D
, &d
, &T
);
4001 bounded
= line_minmax(I
, min
, max
); /* frees I */
4005 value_set_si(mp
, max_periods
);
4006 mpz_fdiv_q(*min
, *min
, d
);
4007 mpz_fdiv_q(*max
, *max
, d
);
4008 value_assign(T
->p
[1][D
->Dimension
], d
);
4009 value_subtract(d
, *max
, *min
);
4010 if (value_ge(d
, mp
))
4013 f
= evalue_dup(&e
->x
.p
->arr
[0]);
4024 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
4025 if ((f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[i
], D
, max_periods
,
4032 static void replace_fract_by_affine(evalue
*e
, evalue
*f
, Value val
)
4036 if (value_notzero_p(e
->d
))
4039 offset
= type_offset(e
->x
.p
);
4040 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
4041 replace_fract_by_affine(&e
->x
.p
->arr
[i
], f
, val
);
4043 if (e
->x
.p
->type
!= fractional
)
4046 if (!eequal(&e
->x
.p
->arr
[0], f
))
4049 replace_by_affine(e
, val
);
4052 /* Look for fractional parts that can be removed by splitting the corresponding
4053 * domain into at most max_periods parts.
4054 * We use a very simply strategy that looks for the first fractional part
4055 * that satisfies the condition, performs the split and then continues
4056 * looking for other fractional parts in the split domains until no
4057 * such fractional part can be found anymore.
4059 void evalue_split_periods(evalue
*e
, int max_periods
, unsigned int MaxRays
)
4066 if (EVALUE_IS_ZERO(*e
))
4068 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
4070 "WARNING: evalue_split_periods called on incorrect evalue type\n");
4078 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
4083 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
4085 f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[2*i
+1], D
, max_periods
,
4090 M
= Matrix_Alloc(2, 2+D
->Dimension
);
4092 value_subtract(d
, max
, min
);
4093 n
= VALUE_TO_INT(d
)+1;
4095 value_set_si(M
->p
[0][0], 1);
4096 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
4097 value_multiply(d
, max
, T
->p
[1][D
->Dimension
]);
4098 value_subtract(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
], d
);
4099 value_set_si(d
, -1);
4100 value_set_si(M
->p
[1][0], 1);
4101 Vector_Scale(T
->p
[0], M
->p
[1]+1, d
, D
->Dimension
+1);
4102 value_addmul(M
->p
[1][1+D
->Dimension
], max
, T
->p
[1][D
->Dimension
]);
4103 value_addto(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4104 T
->p
[1][D
->Dimension
]);
4105 value_decrement(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
]);
4107 p
= new_enode(partition
, e
->x
.p
->size
+ (n
-1)*2, e
->x
.p
->pos
);
4108 for (j
= 0; j
< 2*i
; ++j
) {
4109 value_clear(p
->arr
[j
].d
);
4110 p
->arr
[j
] = e
->x
.p
->arr
[j
];
4112 for (j
= 2*i
+2; j
< e
->x
.p
->size
; ++j
) {
4113 value_clear(p
->arr
[j
+2*(n
-1)].d
);
4114 p
->arr
[j
+2*(n
-1)] = e
->x
.p
->arr
[j
];
4116 for (j
= n
-1; j
>= 0; --j
) {
4118 value_clear(p
->arr
[2*i
+1].d
);
4119 p
->arr
[2*i
+1] = e
->x
.p
->arr
[2*i
+1];
4121 evalue_copy(&p
->arr
[2*(i
+j
)+1], &e
->x
.p
->arr
[2*i
+1]);
4123 value_subtract(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4124 T
->p
[1][D
->Dimension
]);
4125 value_addto(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
],
4126 T
->p
[1][D
->Dimension
]);
4128 replace_fract_by_affine(&p
->arr
[2*(i
+j
)+1], f
, max
);
4129 E
= DomainAddConstraints(D
, M
, MaxRays
);
4130 EVALUE_SET_DOMAIN(p
->arr
[2*(i
+j
)], E
);
4131 if (evalue_range_reduction_in_domain(&p
->arr
[2*(i
+j
)+1], E
))
4132 reduce_evalue(&p
->arr
[2*(i
+j
)+1]);
4133 value_decrement(max
, max
);
4135 value_clear(e
->x
.p
->arr
[2*i
].d
);
4150 void evalue_extract_affine(const evalue
*e
, Value
*coeff
, Value
*cst
, Value
*d
)
4152 value_set_si(*d
, 1);
4154 for ( ; value_zero_p(e
->d
); e
= &e
->x
.p
->arr
[0]) {
4156 assert(e
->x
.p
->type
== polynomial
);
4157 assert(e
->x
.p
->size
== 2);
4158 c
= &e
->x
.p
->arr
[1];
4159 value_multiply(coeff
[e
->x
.p
->pos
-1], *d
, c
->x
.n
);
4160 value_division(coeff
[e
->x
.p
->pos
-1], coeff
[e
->x
.p
->pos
-1], c
->d
);
4162 value_multiply(*cst
, *d
, e
->x
.n
);
4163 value_division(*cst
, *cst
, e
->d
);
4166 /* returns an evalue that corresponds to
4170 static evalue
*term(int param
, Value c
, Value den
)
4172 evalue
*EP
= ALLOC(evalue
);
4174 value_set_si(EP
->d
,0);
4175 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
4176 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
4177 evalue_set_reduce(&EP
->x
.p
->arr
[1], c
, den
);
4181 evalue
*affine2evalue(Value
*coeff
, Value denom
, int nvar
)
4184 evalue
*E
= ALLOC(evalue
);
4186 evalue_set_reduce(E
, coeff
[nvar
], denom
);
4187 for (i
= 0; i
< nvar
; ++i
) {
4189 if (value_zero_p(coeff
[i
]))
4191 t
= term(i
, coeff
[i
], denom
);
4198 void evalue_substitute(evalue
*e
, evalue
**subs
)
4204 if (value_notzero_p(e
->d
))
4208 assert(p
->type
!= partition
);
4210 for (i
= 0; i
< p
->size
; ++i
)
4211 evalue_substitute(&p
->arr
[i
], subs
);
4213 if (p
->type
== relation
) {
4214 /* For relation a ? b : c, compute (a' ? 1) * b' + (a' ? 0 : 1) * c' */
4218 value_set_si(v
->d
, 0);
4219 v
->x
.p
= new_enode(relation
, 3, 0);
4220 evalue_copy(&v
->x
.p
->arr
[0], &p
->arr
[0]);
4221 evalue_set_si(&v
->x
.p
->arr
[1], 0, 1);
4222 evalue_set_si(&v
->x
.p
->arr
[2], 1, 1);
4223 emul(v
, &p
->arr
[2]);
4228 value_set_si(v
->d
, 0);
4229 v
->x
.p
= new_enode(relation
, 2, 0);
4230 value_clear(v
->x
.p
->arr
[0].d
);
4231 v
->x
.p
->arr
[0] = p
->arr
[0];
4232 evalue_set_si(&v
->x
.p
->arr
[1], 1, 1);
4233 emul(v
, &p
->arr
[1]);
4236 eadd(&p
->arr
[2], &p
->arr
[1]);
4237 free_evalue_refs(&p
->arr
[2]);
4245 if (p
->type
== polynomial
)
4250 value_set_si(v
->d
, 0);
4251 v
->x
.p
= new_enode(p
->type
, 3, -1);
4252 value_clear(v
->x
.p
->arr
[0].d
);
4253 v
->x
.p
->arr
[0] = p
->arr
[0];
4254 evalue_set_si(&v
->x
.p
->arr
[1], 0, 1);
4255 evalue_set_si(&v
->x
.p
->arr
[2], 1, 1);
4258 offset
= type_offset(p
);
4260 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
4261 emul(v
, &p
->arr
[i
]);
4262 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
4263 free_evalue_refs(&(p
->arr
[i
]));
4266 if (p
->type
!= polynomial
)
4270 *e
= p
->arr
[offset
];
4274 /* evalue e is given in terms of "new" parameter; CP maps the new
4275 * parameters back to the old parameters.
4276 * Transforms e such that it refers back to the old parameters and
4277 * adds appropriate constraints to the domain.
4278 * In particular, if CP maps the new parameters onto an affine
4279 * subspace of the old parameters, then the corresponding equalities
4280 * are added to the domain.
4281 * Also, if any of the new parameters was a rational combination
4282 * of the old parameters $p' = (<a, p> + c)/m$, then modulo
4283 * constraints ${<a, p> + c)/m} = 0$ are added to ensure
4284 * the new evalue remains non-zero only for integer parameters
4285 * of the new parameters (which have been removed by the substitution).
4287 void evalue_backsubstitute(evalue
*e
, Matrix
*CP
, unsigned MaxRays
)
4294 unsigned nparam
= CP
->NbColumns
-1;
4298 if (EVALUE_IS_ZERO(*e
))
4301 assert(value_zero_p(e
->d
));
4303 assert(p
->type
== partition
);
4305 inv
= left_inverse(CP
, &eq
);
4306 subs
= ALLOCN(evalue
*, nparam
);
4307 for (i
= 0; i
< nparam
; ++i
)
4308 subs
[i
] = affine2evalue(inv
->p
[i
], inv
->p
[nparam
][inv
->NbColumns
-1],
4311 CEq
= Constraints2Polyhedron(eq
, MaxRays
);
4312 addeliminatedparams_partition(p
, inv
, CEq
, inv
->NbColumns
-1, MaxRays
);
4313 Polyhedron_Free(CEq
);
4315 for (i
= 0; i
< p
->size
/2; ++i
)
4316 evalue_substitute(&p
->arr
[2*i
+1], subs
);
4318 for (i
= 0; i
< nparam
; ++i
)
4319 evalue_free(subs
[i
]);
4323 for (i
= 0; i
< inv
->NbRows
-1; ++i
) {
4324 Vector_Gcd(inv
->p
[i
], inv
->NbColumns
, &gcd
);
4325 value_gcd(gcd
, gcd
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]);
4326 if (value_eq(gcd
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]))
4328 Vector_AntiScale(inv
->p
[i
], inv
->p
[i
], gcd
, inv
->NbColumns
);
4329 value_divexact(gcd
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1], gcd
);
4331 for (j
= 0; j
< p
->size
/2; ++j
) {
4332 evalue
*arg
= affine2evalue(inv
->p
[i
], gcd
, inv
->NbColumns
-1);
4337 value_set_si(rel
.d
, 0);
4338 rel
.x
.p
= new_enode(relation
, 2, 0);
4339 value_clear(rel
.x
.p
->arr
[1].d
);
4340 rel
.x
.p
->arr
[1] = p
->arr
[2*j
+1];
4341 ev
= &rel
.x
.p
->arr
[0];
4342 value_set_si(ev
->d
, 0);
4343 ev
->x
.p
= new_enode(fractional
, 3, -1);
4344 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
4345 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
4346 value_clear(ev
->x
.p
->arr
[0].d
);
4347 ev
->x
.p
->arr
[0] = *arg
;
4350 p
->arr
[2*j
+1] = rel
;
4361 * \sum_{i=0}^n c_i/d X^i
4363 * where d is the last element in the vector c.
4365 evalue
*evalue_polynomial(Vector
*c
, const evalue
* X
)
4367 unsigned dim
= c
->Size
-2;
4369 evalue
*EP
= ALLOC(evalue
);
4374 if (EVALUE_IS_ZERO(*X
) || dim
== 0) {
4375 evalue_set(EP
, c
->p
[0], c
->p
[dim
+1]);
4376 reduce_constant(EP
);
4380 evalue_set(EP
, c
->p
[dim
], c
->p
[dim
+1]);
4383 evalue_set(&EC
, c
->p
[dim
], c
->p
[dim
+1]);
4385 for (i
= dim
-1; i
>= 0; --i
) {
4387 value_assign(EC
.x
.n
, c
->p
[i
]);
4390 free_evalue_refs(&EC
);
4394 /* Create an evalue from an array of pairs of domains and evalues. */
4395 evalue
*evalue_from_section_array(struct evalue_section
*s
, int n
)
4400 res
= ALLOC(evalue
);
4404 evalue_set_si(res
, 0, 1);
4406 value_set_si(res
->d
, 0);
4407 res
->x
.p
= new_enode(partition
, 2*n
, s
[0].D
->Dimension
);
4408 for (i
= 0; i
< n
; ++i
) {
4409 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
], s
[i
].D
);
4410 value_clear(res
->x
.p
->arr
[2*i
+1].d
);
4411 res
->x
.p
->arr
[2*i
+1] = *s
[i
].E
;
4418 /* shift variables (>= first, 0-based) in polynomial n up (may be negative) */
4419 void evalue_shift_variables(evalue
*e
, int first
, int n
)
4422 if (value_notzero_p(e
->d
))
4424 assert(e
->x
.p
->type
== polynomial
||
4425 e
->x
.p
->type
== flooring
||
4426 e
->x
.p
->type
== fractional
);
4427 if (e
->x
.p
->type
== polynomial
&& e
->x
.p
->pos
>= first
+1) {
4428 assert(e
->x
.p
->pos
+ n
>= 1);
4431 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
4432 evalue_shift_variables(&e
->x
.p
->arr
[i
], first
, n
);
4435 static const evalue
*outer_floor(evalue
*e
, const evalue
*outer
)
4439 if (value_notzero_p(e
->d
))
4441 switch (e
->x
.p
->type
) {
4443 if (!outer
|| evalue_level_cmp(outer
, &e
->x
.p
->arr
[0]) > 0)
4444 return &e
->x
.p
->arr
[0];
4450 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
4451 outer
= outer_floor(&e
->x
.p
->arr
[i
], outer
);
4454 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
4455 outer
= outer_floor(&e
->x
.p
->arr
[2*i
+1], outer
);
4462 /* Find and return outermost floor argument or NULL if e has no floors */
4463 evalue
*evalue_outer_floor(evalue
*e
)
4465 const evalue
*floor
= outer_floor(e
, NULL
);
4466 return floor
? evalue_dup(floor
): NULL
;
4469 static void evalue_set_to_zero(evalue
*e
)
4471 if (EVALUE_IS_ZERO(*e
))
4473 if (value_zero_p(e
->d
)) {
4474 free_evalue_refs(e
);
4478 value_set_si(e
->d
, 1);
4479 value_set_si(e
->x
.n
, 0);
4482 /* Replace (outer) floor with argument "floor" by variable "var" (0-based)
4483 * and drop terms not containing the floor.
4484 * Returns true if e contains the floor.
4486 int evalue_replace_floor(evalue
*e
, const evalue
*floor
, int var
)
4492 if (value_notzero_p(e
->d
))
4494 switch (e
->x
.p
->type
) {
4496 if (!eequal(floor
, &e
->x
.p
->arr
[0]))
4498 e
->x
.p
->type
= polynomial
;
4499 e
->x
.p
->pos
= 1 + var
;
4501 free_evalue_refs(&e
->x
.p
->arr
[0]);
4502 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
4503 e
->x
.p
->arr
[i
] = e
->x
.p
->arr
[i
+1];
4504 evalue_set_to_zero(&e
->x
.p
->arr
[0]);
4509 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
) {
4510 int c
= evalue_replace_floor(&e
->x
.p
->arr
[i
], floor
, var
);
4513 evalue_set_to_zero(&e
->x
.p
->arr
[i
]);
4514 if (c
&& !reorder
&& evalue_level_cmp(&e
->x
.p
->arr
[i
], e
) < 0)
4517 evalue_reduce_size(e
);
4519 evalue_reorder_terms(e
);
4527 /* Replace (outer) floor with argument "floor" by variable zero */
4528 void evalue_drop_floor(evalue
*e
, const evalue
*floor
)
4533 if (value_notzero_p(e
->d
))
4535 switch (e
->x
.p
->type
) {
4537 if (!eequal(floor
, &e
->x
.p
->arr
[0]))
4540 free_evalue_refs(&p
->arr
[0]);
4541 for (i
= 2; i
< p
->size
; ++i
)
4542 free_evalue_refs(&p
->arr
[i
]);
4550 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
4551 evalue_drop_floor(&e
->x
.p
->arr
[i
], floor
);
4552 evalue_reduce_size(e
);