1 /***********************************************************************/
2 /* copyright 1997, Doran Wilde */
3 /* copyright 1997-2000, Vincent Loechner */
4 /* copyright 2003-2006, Sven Verdoolaege */
5 /* Permission is granted to copy, use, and distribute */
6 /* for any commercial or noncommercial purpose under the terms */
7 /* of the GNU General Public license, version 2, June 1991 */
8 /* (see file : LICENSE). */
9 /***********************************************************************/
16 #include <barvinok/evalue.h>
17 #include <barvinok/barvinok.h>
18 #include <barvinok/util.h>
20 #ifndef value_pmodulus
21 #define value_pmodulus(ref,val1,val2) (mpz_fdiv_r((ref),(val1),(val2)))
24 #define ALLOC(type) (type*)malloc(sizeof(type))
25 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
28 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
30 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
33 void evalue_set_si(evalue
*ev
, int n
, int d
) {
34 value_set_si(ev
->d
, d
);
36 value_set_si(ev
->x
.n
, n
);
39 void evalue_set(evalue
*ev
, Value n
, Value d
) {
40 value_assign(ev
->d
, d
);
42 value_assign(ev
->x
.n
, n
);
47 evalue
*EP
= ALLOC(evalue
);
49 evalue_set_si(EP
, 0, 1);
53 /* returns an evalue that corresponds to
57 evalue
*evalue_var(int var
)
59 evalue
*EP
= ALLOC(evalue
);
61 value_set_si(EP
->d
,0);
62 EP
->x
.p
= new_enode(polynomial
, 2, var
+ 1);
63 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
64 evalue_set_si(&EP
->x
.p
->arr
[1], 1, 1);
68 void aep_evalue(evalue
*e
, int *ref
) {
73 if (value_notzero_p(e
->d
))
74 return; /* a rational number, its already reduced */
76 return; /* hum... an overflow probably occured */
78 /* First check the components of p */
79 for (i
=0;i
<p
->size
;i
++)
80 aep_evalue(&p
->arr
[i
],ref
);
87 p
->pos
= ref
[p
->pos
-1]+1;
93 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
99 if (value_notzero_p(e
->d
))
100 return; /* a rational number, its already reduced */
102 return; /* hum... an overflow probably occured */
105 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
106 for(i
=0;i
<CT
->NbRows
-1;i
++)
107 for(j
=0;j
<CT
->NbColumns
;j
++)
108 if(value_notzero_p(CT
->p
[i
][j
])) {
113 /* Transform the references in e, using ref */
117 } /* addeliminatedparams_evalue */
119 static void addeliminatedparams_partition(enode
*p
, Matrix
*CT
, Polyhedron
*CEq
,
120 unsigned nparam
, unsigned MaxRays
)
123 assert(p
->type
== partition
);
126 for (i
= 0; i
< p
->size
/2; i
++) {
127 Polyhedron
*D
= EVALUE_DOMAIN(p
->arr
[2*i
]);
128 Polyhedron
*T
= DomainPreimage(D
, CT
, MaxRays
);
132 T
= DomainIntersection(D
, CEq
, MaxRays
);
135 EVALUE_SET_DOMAIN(p
->arr
[2*i
], T
);
139 void addeliminatedparams_enum(evalue
*e
, Matrix
*CT
, Polyhedron
*CEq
,
140 unsigned MaxRays
, unsigned nparam
)
145 if (CT
->NbRows
== CT
->NbColumns
)
148 if (EVALUE_IS_ZERO(*e
))
151 if (value_notzero_p(e
->d
)) {
154 value_set_si(res
.d
, 0);
155 res
.x
.p
= new_enode(partition
, 2, nparam
);
156 EVALUE_SET_DOMAIN(res
.x
.p
->arr
[0],
157 DomainConstraintSimplify(Polyhedron_Copy(CEq
), MaxRays
));
158 value_clear(res
.x
.p
->arr
[1].d
);
159 res
.x
.p
->arr
[1] = *e
;
167 addeliminatedparams_partition(p
, CT
, CEq
, nparam
, MaxRays
);
168 for (i
= 0; i
< p
->size
/2; i
++)
169 addeliminatedparams_evalue(&p
->arr
[2*i
+1], CT
);
172 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
180 assert(value_notzero_p(e1
->d
));
181 assert(value_notzero_p(e2
->d
));
182 value_multiply(m
, e1
->x
.n
, e2
->d
);
183 value_multiply(m2
, e2
->x
.n
, e1
->d
);
186 else if (value_gt(m
, m2
))
196 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
198 if (value_notzero_p(e1
->d
)) {
200 if (value_zero_p(e2
->d
))
202 r
= mod_rational_smaller(e1
, e2
);
203 return r
== -1 ? 0 : r
;
205 if (value_notzero_p(e2
->d
))
207 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
209 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
212 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
214 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
219 static int mod_term_smaller(const evalue
*e1
, const evalue
*e2
)
221 assert(value_zero_p(e1
->d
));
222 assert(value_zero_p(e2
->d
));
223 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
224 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
225 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
228 static void check_order(const evalue
*e
)
233 if (value_notzero_p(e
->d
))
236 switch (e
->x
.p
->type
) {
238 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
239 check_order(&e
->x
.p
->arr
[2*i
+1]);
242 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
244 if (value_notzero_p(a
->d
))
246 switch (a
->x
.p
->type
) {
248 assert(mod_term_smaller(&e
->x
.p
->arr
[0], &a
->x
.p
->arr
[0]));
257 for (i
= 0; i
< e
->x
.p
->size
; ++i
) {
259 if (value_notzero_p(a
->d
))
261 switch (a
->x
.p
->type
) {
263 assert(e
->x
.p
->pos
< a
->x
.p
->pos
);
274 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
276 if (value_notzero_p(a
->d
))
278 switch (a
->x
.p
->type
) {
289 /* Negative pos means inequality */
290 /* s is negative of substitution if m is not zero */
299 struct fixed_param
*fixed
;
304 static int relations_depth(evalue
*e
)
309 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
310 e
= &e
->x
.p
->arr
[1], ++d
);
314 static void poly_denom_not_constant(evalue
**pp
, Value
*d
)
319 while (value_zero_p(p
->d
)) {
320 assert(p
->x
.p
->type
== polynomial
);
321 assert(p
->x
.p
->size
== 2);
322 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
323 value_lcm(*d
, *d
, p
->x
.p
->arr
[1].d
);
329 static void poly_denom(evalue
*p
, Value
*d
)
331 poly_denom_not_constant(&p
, d
);
332 value_lcm(*d
, *d
, p
->d
);
335 static void realloc_substitution(struct subst
*s
, int d
)
337 struct fixed_param
*n
;
340 for (i
= 0; i
< s
->n
; ++i
)
347 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
353 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
356 /* May have been reduced already */
357 if (value_notzero_p(m
->d
))
360 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
361 assert(m
->x
.p
->size
== 3);
363 /* fractional was inverted during reduction
364 * invert it back and move constant in
366 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
367 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
368 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
369 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
370 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
371 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
372 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
373 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
374 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
375 value_set_si(m
->x
.p
->arr
[1].d
, 1);
378 /* Oops. Nested identical relations. */
379 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
382 if (s
->n
>= s
->max
) {
383 int d
= relations_depth(r
);
384 realloc_substitution(s
, d
);
388 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
389 assert(p
->x
.p
->size
== 2);
392 assert(value_pos_p(f
->x
.n
));
394 value_init(s
->fixed
[s
->n
].m
);
395 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
396 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
397 value_init(s
->fixed
[s
->n
].d
);
398 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
399 value_init(s
->fixed
[s
->n
].s
.d
);
400 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
406 static int type_offset(enode
*p
)
408 return p
->type
== fractional
? 1 :
409 p
->type
== flooring
? 1 : 0;
412 static void reorder_terms_about(enode
*p
, evalue
*v
)
415 int offset
= type_offset(p
);
417 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
419 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
420 free_evalue_refs(&(p
->arr
[i
]));
426 static void reorder_terms(evalue
*e
)
431 assert(value_zero_p(e
->d
));
433 assert(p
->type
== fractional
); /* for now */
436 value_set_si(f
.d
, 0);
437 f
.x
.p
= new_enode(fractional
, 3, -1);
438 value_clear(f
.x
.p
->arr
[0].d
);
439 f
.x
.p
->arr
[0] = p
->arr
[0];
440 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
441 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
442 reorder_terms_about(p
, &f
);
448 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
454 if (value_notzero_p(e
->d
)) {
456 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
457 return; /* a rational number, its already reduced */
461 return; /* hum... an overflow probably occured */
463 /* First reduce the components of p */
464 add
= p
->type
== relation
;
465 for (i
=0; i
<p
->size
; i
++) {
467 add
= add_modulo_substitution(s
, e
);
469 if (i
== 0 && p
->type
==fractional
)
470 _reduce_evalue(&p
->arr
[i
], s
, 1);
472 _reduce_evalue(&p
->arr
[i
], s
, fract
);
474 if (add
&& i
== p
->size
-1) {
476 value_clear(s
->fixed
[s
->n
].m
);
477 value_clear(s
->fixed
[s
->n
].d
);
478 free_evalue_refs(&s
->fixed
[s
->n
].s
);
479 } else if (add
&& i
== 1)
480 s
->fixed
[s
->n
-1].pos
*= -1;
483 if (p
->type
==periodic
) {
485 /* Try to reduce the period */
486 for (i
=1; i
<=(p
->size
)/2; i
++) {
487 if ((p
->size
% i
)==0) {
489 /* Can we reduce the size to i ? */
491 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
492 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
495 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
499 you_lose
: /* OK, lets not do it */
504 /* Try to reduce its strength */
507 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
511 else if (p
->type
==polynomial
) {
512 for (k
= 0; s
&& k
< s
->n
; ++k
) {
513 if (s
->fixed
[k
].pos
== p
->pos
) {
514 int divide
= value_notone_p(s
->fixed
[k
].d
);
517 if (value_notzero_p(s
->fixed
[k
].m
)) {
520 assert(p
->size
== 2);
521 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
523 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
530 value_assign(d
.d
, s
->fixed
[k
].d
);
532 if (value_notzero_p(s
->fixed
[k
].m
))
533 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
535 value_set_si(d
.x
.n
, 1);
538 for (i
=p
->size
-1;i
>=1;i
--) {
539 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
541 emul(&d
, &p
->arr
[i
]);
542 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
543 free_evalue_refs(&(p
->arr
[i
]));
546 _reduce_evalue(&p
->arr
[0], s
, fract
);
549 free_evalue_refs(&d
);
555 /* Try to reduce the degree */
556 for (i
=p
->size
-1;i
>=1;i
--) {
557 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
559 /* Zero coefficient */
560 free_evalue_refs(&(p
->arr
[i
]));
565 /* Try to reduce its strength */
568 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
572 else if (p
->type
==fractional
) {
576 if (value_notzero_p(p
->arr
[0].d
)) {
578 value_assign(v
.d
, p
->arr
[0].d
);
580 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
585 evalue
*pp
= &p
->arr
[0];
586 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
587 assert(pp
->x
.p
->size
== 2);
589 /* search for exact duplicate among the modulo inequalities */
591 f
= &pp
->x
.p
->arr
[1];
592 for (k
= 0; s
&& k
< s
->n
; ++k
) {
593 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
594 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
595 value_eq(s
->fixed
[k
].m
, f
->d
) &&
596 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
603 /* replace { E/m } by { (E-1)/m } + 1/m */
608 evalue_set_si(&extra
, 1, 1);
609 value_assign(extra
.d
, g
);
610 eadd(&extra
, &v
.x
.p
->arr
[1]);
611 free_evalue_refs(&extra
);
613 /* We've been going in circles; stop now */
614 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
615 free_evalue_refs(&v
);
617 evalue_set_si(&v
, 0, 1);
622 value_set_si(v
.d
, 0);
623 v
.x
.p
= new_enode(fractional
, 3, -1);
624 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
625 value_assign(v
.x
.p
->arr
[1].d
, g
);
626 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
627 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
630 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
633 value_division(f
->d
, g
, f
->d
);
634 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
635 value_assign(f
->d
, g
);
636 value_decrement(f
->x
.n
, f
->x
.n
);
637 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
639 value_gcd(g
, f
->d
, f
->x
.n
);
640 value_division(f
->d
, f
->d
, g
);
641 value_division(f
->x
.n
, f
->x
.n
, g
);
650 /* reduction may have made this fractional arg smaller */
651 i
= reorder
? p
->size
: 1;
652 for ( ; i
< p
->size
; ++i
)
653 if (value_zero_p(p
->arr
[i
].d
) &&
654 p
->arr
[i
].x
.p
->type
== fractional
&&
655 !mod_term_smaller(e
, &p
->arr
[i
]))
659 value_set_si(v
.d
, 0);
660 v
.x
.p
= new_enode(fractional
, 3, -1);
661 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
662 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
663 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
671 evalue
*pp
= &p
->arr
[0];
674 poly_denom_not_constant(&pp
, &m
);
675 mpz_fdiv_r(r
, m
, pp
->d
);
676 if (value_notzero_p(r
)) {
678 value_set_si(v
.d
, 0);
679 v
.x
.p
= new_enode(fractional
, 3, -1);
681 value_multiply(r
, m
, pp
->x
.n
);
682 value_multiply(v
.x
.p
->arr
[1].d
, m
, pp
->d
);
683 value_init(v
.x
.p
->arr
[1].x
.n
);
684 mpz_fdiv_r(v
.x
.p
->arr
[1].x
.n
, r
, pp
->d
);
685 mpz_fdiv_q(r
, r
, pp
->d
);
687 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
688 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
690 while (value_zero_p(pp
->d
))
691 pp
= &pp
->x
.p
->arr
[0];
693 value_assign(pp
->d
, m
);
694 value_assign(pp
->x
.n
, r
);
696 value_gcd(r
, pp
->d
, pp
->x
.n
);
697 value_division(pp
->d
, pp
->d
, r
);
698 value_division(pp
->x
.n
, pp
->x
.n
, r
);
711 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
712 pp
= &pp
->x
.p
->arr
[0]) {
713 f
= &pp
->x
.p
->arr
[1];
714 assert(value_pos_p(f
->d
));
715 mpz_mul_ui(twice
, f
->x
.n
, 2);
716 if (value_lt(twice
, f
->d
))
718 if (value_eq(twice
, f
->d
))
726 value_set_si(v
.d
, 0);
727 v
.x
.p
= new_enode(fractional
, 3, -1);
728 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
729 poly_denom(&p
->arr
[0], &twice
);
730 value_assign(v
.x
.p
->arr
[1].d
, twice
);
731 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
732 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
733 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
735 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
736 pp
= &pp
->x
.p
->arr
[0]) {
737 f
= &pp
->x
.p
->arr
[1];
738 value_oppose(f
->x
.n
, f
->x
.n
);
739 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
741 value_division(pp
->d
, twice
, pp
->d
);
742 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
743 value_assign(pp
->d
, twice
);
744 value_oppose(pp
->x
.n
, pp
->x
.n
);
745 value_decrement(pp
->x
.n
, pp
->x
.n
);
746 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
748 /* Maybe we should do this during reduction of
751 value_gcd(twice
, pp
->d
, pp
->x
.n
);
752 value_division(pp
->d
, pp
->d
, twice
);
753 value_division(pp
->x
.n
, pp
->x
.n
, twice
);
763 reorder_terms_about(p
, &v
);
764 _reduce_evalue(&p
->arr
[1], s
, fract
);
767 /* Try to reduce the degree */
768 for (i
=p
->size
-1;i
>=2;i
--) {
769 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
771 /* Zero coefficient */
772 free_evalue_refs(&(p
->arr
[i
]));
777 /* Try to reduce its strength */
780 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
781 free_evalue_refs(&(p
->arr
[0]));
785 else if (p
->type
== flooring
) {
786 /* Try to reduce the degree */
787 for (i
=p
->size
-1;i
>=2;i
--) {
788 if (!EVALUE_IS_ZERO(p
->arr
[i
]))
790 /* Zero coefficient */
791 free_evalue_refs(&(p
->arr
[i
]));
796 /* Try to reduce its strength */
799 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
800 free_evalue_refs(&(p
->arr
[0]));
804 else if (p
->type
== relation
) {
805 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
806 free_evalue_refs(&(p
->arr
[2]));
807 free_evalue_refs(&(p
->arr
[0]));
814 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
815 free_evalue_refs(&(p
->arr
[2]));
818 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
819 free_evalue_refs(&(p
->arr
[1]));
820 free_evalue_refs(&(p
->arr
[0]));
821 evalue_set_si(e
, 0, 1);
828 /* Relation was reduced by means of an identical
829 * inequality => remove
831 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
834 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
835 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
837 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
839 free_evalue_refs(&(p
->arr
[2]));
843 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
845 evalue_set_si(e
, 0, 1);
846 free_evalue_refs(&(p
->arr
[1]));
848 free_evalue_refs(&(p
->arr
[0]));
854 } /* reduce_evalue */
856 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
861 for (k
= 0; k
< dim
; ++k
)
862 if (value_notzero_p(row
[k
+1]))
865 Vector_Normalize_Positive(row
+1, dim
+1, k
);
866 assert(s
->n
< s
->max
);
867 value_init(s
->fixed
[s
->n
].d
);
868 value_init(s
->fixed
[s
->n
].m
);
869 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
870 s
->fixed
[s
->n
].pos
= k
+1;
871 value_set_si(s
->fixed
[s
->n
].m
, 0);
872 r
= &s
->fixed
[s
->n
].s
;
874 for (l
= k
+1; l
< dim
; ++l
)
875 if (value_notzero_p(row
[l
+1])) {
876 value_set_si(r
->d
, 0);
877 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
878 value_init(r
->x
.p
->arr
[1].x
.n
);
879 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
880 value_set_si(r
->x
.p
->arr
[1].d
, 1);
884 value_oppose(r
->x
.n
, row
[dim
+1]);
885 value_set_si(r
->d
, 1);
889 static void _reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
, struct subst
*s
)
892 Polyhedron
*orig
= D
;
897 D
= DomainConvex(D
, 0);
898 if (!D
->next
&& D
->NbEq
) {
902 realloc_substitution(s
, dim
);
904 int d
= relations_depth(e
);
906 NALLOC(s
->fixed
, s
->max
);
909 for (j
= 0; j
< D
->NbEq
; ++j
)
910 add_substitution(s
, D
->Constraint
[j
], dim
);
914 _reduce_evalue(e
, s
, 0);
917 for (j
= 0; j
< s
->n
; ++j
) {
918 value_clear(s
->fixed
[j
].d
);
919 value_clear(s
->fixed
[j
].m
);
920 free_evalue_refs(&s
->fixed
[j
].s
);
925 void reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
)
927 struct subst s
= { NULL
, 0, 0 };
929 if (EVALUE_IS_ZERO(*e
))
933 evalue_set_si(e
, 0, 1);
936 _reduce_evalue_in_domain(e
, D
, &s
);
941 void reduce_evalue (evalue
*e
) {
942 struct subst s
= { NULL
, 0, 0 };
944 if (value_notzero_p(e
->d
))
945 return; /* a rational number, its already reduced */
947 if (e
->x
.p
->type
== partition
) {
950 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
951 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
953 /* This shouldn't really happen;
954 * Empty domains should not be added.
956 POL_ENSURE_VERTICES(D
);
958 _reduce_evalue_in_domain(&e
->x
.p
->arr
[2*i
+1], D
, &s
);
960 if (emptyQ(D
) || EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
961 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
962 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
963 value_clear(e
->x
.p
->arr
[2*i
].d
);
965 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
966 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
970 if (e
->x
.p
->size
== 0) {
972 evalue_set_si(e
, 0, 1);
975 _reduce_evalue(e
, &s
, 0);
980 static void print_evalue_r(FILE *DST
, const evalue
*e
, const char *const *pname
)
982 if(value_notzero_p(e
->d
)) {
983 if(value_notone_p(e
->d
)) {
984 value_print(DST
,VALUE_FMT
,e
->x
.n
);
986 value_print(DST
,VALUE_FMT
,e
->d
);
989 value_print(DST
,VALUE_FMT
,e
->x
.n
);
993 print_enode(DST
,e
->x
.p
,pname
);
997 void print_evalue(FILE *DST
, const evalue
*e
, const char * const *pname
)
999 print_evalue_r(DST
, e
, pname
);
1000 if (value_notzero_p(e
->d
))
1004 void print_enode(FILE *DST
, enode
*p
, const char *const *pname
)
1009 fprintf(DST
, "NULL");
1015 for (i
=0; i
<p
->size
; i
++) {
1016 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1020 fprintf(DST
, " }\n");
1024 for (i
=p
->size
-1; i
>=0; i
--) {
1025 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1026 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
1028 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
1030 fprintf(DST
, " )\n");
1034 for (i
=0; i
<p
->size
; i
++) {
1035 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1036 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
1038 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
1043 for (i
=p
->size
-1; i
>=1; i
--) {
1044 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1046 fprintf(DST
, " * ");
1047 fprintf(DST
, p
->type
== flooring
? "[" : "{");
1048 print_evalue_r(DST
, &p
->arr
[0], pname
);
1049 fprintf(DST
, p
->type
== flooring
? "]" : "}");
1051 fprintf(DST
, "^%d + ", i
-1);
1053 fprintf(DST
, " + ");
1056 fprintf(DST
, " )\n");
1060 print_evalue_r(DST
, &p
->arr
[0], pname
);
1061 fprintf(DST
, "= 0 ] * \n");
1062 print_evalue_r(DST
, &p
->arr
[1], pname
);
1064 fprintf(DST
, " +\n [ ");
1065 print_evalue_r(DST
, &p
->arr
[0], pname
);
1066 fprintf(DST
, "!= 0 ] * \n");
1067 print_evalue_r(DST
, &p
->arr
[2], pname
);
1071 char **new_names
= NULL
;
1072 const char *const *names
= pname
;
1073 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
1074 if (!pname
|| p
->pos
< maxdim
) {
1075 new_names
= ALLOCN(char *, maxdim
);
1076 for (i
= 0; i
< p
->pos
; ++i
) {
1078 new_names
[i
] = (char *)pname
[i
];
1080 new_names
[i
] = ALLOCN(char, 10);
1081 snprintf(new_names
[i
], 10, "%c", 'P'+i
);
1084 for ( ; i
< maxdim
; ++i
) {
1085 new_names
[i
] = ALLOCN(char, 10);
1086 snprintf(new_names
[i
], 10, "_p%d", i
);
1088 names
= (const char**)new_names
;
1091 for (i
=0; i
<p
->size
/2; i
++) {
1092 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
1093 print_evalue_r(DST
, &p
->arr
[2*i
+1], names
);
1094 if (value_notzero_p(p
->arr
[2*i
+1].d
))
1098 if (!pname
|| p
->pos
< maxdim
) {
1099 for (i
= pname
? p
->pos
: 0; i
< maxdim
; ++i
)
1112 static void eadd_rev(const evalue
*e1
, evalue
*res
)
1116 evalue_copy(&ev
, e1
);
1118 free_evalue_refs(res
);
1122 static void eadd_rev_cst(const evalue
*e1
, evalue
*res
)
1126 evalue_copy(&ev
, e1
);
1127 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
1128 free_evalue_refs(res
);
1132 static int is_zero_on(evalue
*e
, Polyhedron
*D
)
1137 tmp
.x
.p
= new_enode(partition
, 2, D
->Dimension
);
1138 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Domain_Copy(D
));
1139 evalue_copy(&tmp
.x
.p
->arr
[1], e
);
1140 reduce_evalue(&tmp
);
1141 is_zero
= EVALUE_IS_ZERO(tmp
);
1142 free_evalue_refs(&tmp
);
1146 struct section
{ Polyhedron
* D
; evalue E
; };
1148 void eadd_partitions(const evalue
*e1
, evalue
*res
)
1153 s
= (struct section
*)
1154 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
1155 sizeof(struct section
));
1157 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1158 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1159 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1162 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1163 assert(res
->x
.p
->size
>= 2);
1164 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1165 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
1167 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
1169 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1178 /* See if we can extend one of the domains in res to cover fd */
1179 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1180 if (is_zero_on(&res
->x
.p
->arr
[2*i
+1], fd
))
1182 if (i
< res
->x
.p
->size
/2) {
1183 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
1184 DomainConcat(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
])));
1187 value_init(s
[n
].E
.d
);
1188 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
1192 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1193 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
1194 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1196 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1197 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1203 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
1204 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1206 value_init(s
[n
].E
.d
);
1207 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1208 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1209 if (!emptyQ(fd
) && is_zero_on(&e1
->x
.p
->arr
[2*j
+1], fd
)) {
1210 d
= DomainConcat(fd
, d
);
1211 fd
= Empty_Polyhedron(fd
->Dimension
);
1217 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
1221 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1224 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1225 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1226 value_clear(res
->x
.p
->arr
[2*i
].d
);
1231 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1232 for (j
= 0; j
< n
; ++j
) {
1233 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1234 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1235 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1236 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1242 static void explicit_complement(evalue
*res
)
1244 enode
*rel
= new_enode(relation
, 3, 0);
1246 value_clear(rel
->arr
[0].d
);
1247 rel
->arr
[0] = res
->x
.p
->arr
[0];
1248 value_clear(rel
->arr
[1].d
);
1249 rel
->arr
[1] = res
->x
.p
->arr
[1];
1250 value_set_si(rel
->arr
[2].d
, 1);
1251 value_init(rel
->arr
[2].x
.n
);
1252 value_set_si(rel
->arr
[2].x
.n
, 0);
1257 static void reduce_constant(evalue
*e
)
1262 value_gcd(g
, e
->x
.n
, e
->d
);
1263 if (value_notone_p(g
)) {
1264 value_division(e
->d
, e
->d
,g
);
1265 value_division(e
->x
.n
, e
->x
.n
,g
);
1270 void eadd(const evalue
*e1
, evalue
*res
)
1274 if (EVALUE_IS_ZERO(*e1
))
1277 if (EVALUE_IS_ZERO(*res
)) {
1278 if (value_notzero_p(e1
->d
)) {
1279 value_assign(res
->d
, e1
->d
);
1280 value_assign(res
->x
.n
, e1
->x
.n
);
1282 value_clear(res
->x
.n
);
1283 value_set_si(res
->d
, 0);
1284 res
->x
.p
= ecopy(e1
->x
.p
);
1289 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1290 /* Add two rational numbers */
1291 if (value_eq(e1
->d
, res
->d
))
1292 value_addto(res
->x
.n
, res
->x
.n
, e1
->x
.n
);
1294 value_multiply(res
->x
.n
, res
->x
.n
, e1
->d
);
1295 value_addmul(res
->x
.n
, e1
->x
.n
, res
->d
);
1296 value_multiply(res
->d
,e1
->d
,res
->d
);
1298 reduce_constant(res
);
1301 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
1302 switch (res
->x
.p
->type
) {
1304 /* Add the constant to the constant term of a polynomial*/
1305 eadd(e1
, &res
->x
.p
->arr
[0]);
1308 /* Add the constant to all elements of a periodic number */
1309 for (i
=0; i
<res
->x
.p
->size
; i
++) {
1310 eadd(e1
, &res
->x
.p
->arr
[i
]);
1314 fprintf(stderr
, "eadd: cannot add const with vector\n");
1318 eadd(e1
, &res
->x
.p
->arr
[1]);
1321 assert(EVALUE_IS_ZERO(*e1
));
1322 break; /* Do nothing */
1324 /* Create (zero) complement if needed */
1325 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
1326 explicit_complement(res
);
1327 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1328 eadd(e1
, &res
->x
.p
->arr
[i
]);
1334 /* add polynomial or periodic to constant
1335 * you have to exchange e1 and res, before doing addition */
1337 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1341 else { // ((e1->d==0) && (res->d==0))
1342 assert(!((e1
->x
.p
->type
== partition
) ^
1343 (res
->x
.p
->type
== partition
)));
1344 if (e1
->x
.p
->type
== partition
) {
1345 eadd_partitions(e1
, res
);
1348 if (e1
->x
.p
->type
== relation
&&
1349 (res
->x
.p
->type
!= relation
||
1350 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1354 if (res
->x
.p
->type
== relation
) {
1355 if (e1
->x
.p
->type
== relation
&&
1356 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1357 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1358 explicit_complement(res
);
1359 for (i
= 1; i
< e1
->x
.p
->size
; ++i
)
1360 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1363 if (res
->x
.p
->size
< 3)
1364 explicit_complement(res
);
1365 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1366 eadd(e1
, &res
->x
.p
->arr
[i
]);
1369 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
1370 /* adding to evalues of different type. two cases are possible
1371 * res is periodic and e1 is polynomial, you have to exchange
1372 * e1 and res then to add e1 to the constant term of res */
1373 if (e1
->x
.p
->type
== polynomial
) {
1374 eadd_rev_cst(e1
, res
);
1376 else if (res
->x
.p
->type
== polynomial
) {
1377 /* res is polynomial and e1 is periodic,
1378 add e1 to the constant term of res */
1380 eadd(e1
,&res
->x
.p
->arr
[0]);
1386 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
1387 ((res
->x
.p
->type
== fractional
||
1388 res
->x
.p
->type
== flooring
) &&
1389 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1390 /* adding evalues of different position (i.e function of different unknowns
1391 * to case are possible */
1393 switch (res
->x
.p
->type
) {
1396 if (mod_term_smaller(res
, e1
))
1397 eadd(e1
,&res
->x
.p
->arr
[1]);
1399 eadd_rev_cst(e1
, res
);
1401 case polynomial
: // res and e1 are polynomials
1402 // add e1 to the constant term of res
1404 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1405 eadd(e1
,&res
->x
.p
->arr
[0]);
1407 eadd_rev_cst(e1
, res
);
1408 // value_clear(g); value_clear(m1); value_clear(m2);
1410 case periodic
: // res and e1 are pointers to periodic numbers
1411 //add e1 to all elements of res
1413 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1414 for (i
=0;i
<res
->x
.p
->size
;i
++) {
1415 eadd(e1
,&res
->x
.p
->arr
[i
]);
1426 //same type , same pos and same size
1427 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1428 // add any element in e1 to the corresponding element in res
1429 i
= type_offset(res
->x
.p
);
1431 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1432 for (; i
<res
->x
.p
->size
; i
++) {
1433 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1438 /* Sizes are different */
1439 switch(res
->x
.p
->type
) {
1443 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1444 /* new enode and add res to that new node. If you do not do */
1445 /* that, you lose the the upper weight part of e1 ! */
1447 if(e1
->x
.p
->size
> res
->x
.p
->size
)
1450 i
= type_offset(res
->x
.p
);
1452 assert(eequal(&e1
->x
.p
->arr
[0],
1453 &res
->x
.p
->arr
[0]));
1454 for (; i
<e1
->x
.p
->size
; i
++) {
1455 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1462 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1465 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1466 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1467 to add periodicaly elements of e1 to elements of ne, and finaly to
1472 value_init(ex
); value_init(ey
);value_init(ep
);
1475 value_set_si(ex
,e1
->x
.p
->size
);
1476 value_set_si(ey
,res
->x
.p
->size
);
1477 value_assign (ep
,*Lcm(ex
,ey
));
1478 p
=(int)mpz_get_si(ep
);
1479 ne
= (evalue
*) malloc (sizeof(evalue
));
1481 value_set_si( ne
->d
,0);
1483 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1485 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1488 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1491 value_assign(res
->d
, ne
->d
);
1497 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
1506 static void emul_rev(const evalue
*e1
, evalue
*res
)
1510 evalue_copy(&ev
, e1
);
1512 free_evalue_refs(res
);
1516 static void emul_poly(const evalue
*e1
, evalue
*res
)
1518 int i
, j
, offset
= type_offset(res
->x
.p
);
1521 int size
= (e1
->x
.p
->size
+ res
->x
.p
->size
- offset
- 1);
1523 p
= new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1525 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1526 if (!EVALUE_IS_ZERO(e1
->x
.p
->arr
[i
]))
1529 /* special case pure power */
1530 if (i
== e1
->x
.p
->size
-1) {
1532 value_clear(p
->arr
[0].d
);
1533 p
->arr
[0] = res
->x
.p
->arr
[0];
1535 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1536 evalue_set_si(&p
->arr
[i
], 0, 1);
1537 for (i
= offset
; i
< res
->x
.p
->size
; ++i
) {
1538 value_clear(p
->arr
[i
+e1
->x
.p
->size
-offset
-1].d
);
1539 p
->arr
[i
+e1
->x
.p
->size
-offset
-1] = res
->x
.p
->arr
[i
];
1540 emul(&e1
->x
.p
->arr
[e1
->x
.p
->size
-1],
1541 &p
->arr
[i
+e1
->x
.p
->size
-offset
-1]);
1549 value_set_si(tmp
.d
,0);
1552 evalue_copy(&p
->arr
[0], &e1
->x
.p
->arr
[0]);
1553 for (i
= offset
; i
< e1
->x
.p
->size
; i
++) {
1554 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1555 emul(&res
->x
.p
->arr
[offset
], &tmp
.x
.p
->arr
[i
]);
1558 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1559 for (i
= offset
+1; i
<res
->x
.p
->size
; i
++)
1560 for (j
= offset
; j
<e1
->x
.p
->size
; j
++) {
1563 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1564 emul(&res
->x
.p
->arr
[i
], &ev
);
1565 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-offset
]);
1566 free_evalue_refs(&ev
);
1568 free_evalue_refs(res
);
1572 void emul_partitions(const evalue
*e1
, evalue
*res
)
1577 s
= (struct section
*)
1578 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1579 sizeof(struct section
));
1581 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1582 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1583 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1586 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1587 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1588 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1589 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1595 /* This code is only needed because the partitions
1596 are not true partitions.
1598 for (k
= 0; k
< n
; ++k
) {
1599 if (DomainIncludes(s
[k
].D
, d
))
1601 if (DomainIncludes(d
, s
[k
].D
)) {
1602 Domain_Free(s
[k
].D
);
1603 free_evalue_refs(&s
[k
].E
);
1614 value_init(s
[n
].E
.d
);
1615 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1616 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1620 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1621 value_clear(res
->x
.p
->arr
[2*i
].d
);
1622 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1627 evalue_set_si(res
, 0, 1);
1629 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1630 for (j
= 0; j
< n
; ++j
) {
1631 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1632 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1633 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1634 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1641 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1643 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1644 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1645 * evalues is not treated here */
1647 void emul(const evalue
*e1
, evalue
*res
)
1651 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1652 fprintf(stderr
, "emul: do not proced on evector type !\n");
1656 if (EVALUE_IS_ZERO(*res
))
1659 if (EVALUE_IS_ONE(*e1
))
1662 if (EVALUE_IS_ZERO(*e1
)) {
1663 if (value_notzero_p(res
->d
)) {
1664 value_assign(res
->d
, e1
->d
);
1665 value_assign(res
->x
.n
, e1
->x
.n
);
1667 free_evalue_refs(res
);
1669 evalue_set_si(res
, 0, 1);
1674 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1675 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1676 emul_partitions(e1
, res
);
1679 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1680 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1681 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1683 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1684 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1685 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1686 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3) {
1687 free_evalue_refs(&res
->x
.p
->arr
[2]);
1690 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1691 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1694 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1695 emul(e1
, &res
->x
.p
->arr
[i
]);
1697 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1698 switch(e1
->x
.p
->type
) {
1700 switch(res
->x
.p
->type
) {
1702 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1703 /* Product of two polynomials of the same variable */
1708 /* Product of two polynomials of different variables */
1710 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1711 for( i
=0; i
<res
->x
.p
->size
; i
++)
1712 emul(e1
, &res
->x
.p
->arr
[i
]);
1721 /* Product of a polynomial and a periodic or fractional */
1728 switch(res
->x
.p
->type
) {
1730 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1731 /* Product of two periodics of the same parameter and period */
1733 for(i
=0; i
<res
->x
.p
->size
;i
++)
1734 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1739 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1740 /* Product of two periodics of the same parameter and different periods */
1744 value_init(x
); value_init(y
);value_init(z
);
1747 value_set_si(x
,e1
->x
.p
->size
);
1748 value_set_si(y
,res
->x
.p
->size
);
1749 value_assign (z
,*Lcm(x
,y
));
1750 lcm
=(int)mpz_get_si(z
);
1751 newp
= (evalue
*) malloc (sizeof(evalue
));
1752 value_init(newp
->d
);
1753 value_set_si( newp
->d
,0);
1754 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1755 for(i
=0;i
<lcm
;i
++) {
1756 evalue_copy(&newp
->x
.p
->arr
[i
],
1757 &res
->x
.p
->arr
[i
%iy
]);
1760 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1762 value_assign(res
->d
,newp
->d
);
1765 value_clear(x
); value_clear(y
);value_clear(z
);
1769 /* Product of two periodics of different parameters */
1771 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1772 for(i
=0; i
<res
->x
.p
->size
; i
++)
1773 emul(e1
, &(res
->x
.p
->arr
[i
]));
1781 /* Product of a periodic and a polynomial */
1783 for(i
=0; i
<res
->x
.p
->size
; i
++)
1784 emul(e1
, &(res
->x
.p
->arr
[i
]));
1791 switch(res
->x
.p
->type
) {
1793 for(i
=0; i
<res
->x
.p
->size
; i
++)
1794 emul(e1
, &(res
->x
.p
->arr
[i
]));
1801 assert(e1
->x
.p
->type
== res
->x
.p
->type
);
1802 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1803 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1806 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1807 if (e1
->x
.p
->type
!= fractional
|| !value_two_p(d
.d
))
1812 value_set_si(d
.x
.n
, 1);
1813 /* { x }^2 == { x }/2 */
1814 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1815 assert(e1
->x
.p
->size
== 3);
1816 assert(res
->x
.p
->size
== 3);
1818 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1820 eadd(&res
->x
.p
->arr
[1], &tmp
);
1821 emul(&e1
->x
.p
->arr
[2], &tmp
);
1822 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[1]);
1823 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[2]);
1824 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1825 free_evalue_refs(&tmp
);
1830 if(mod_term_smaller(res
, e1
))
1831 for(i
=1; i
<res
->x
.p
->size
; i
++)
1832 emul(e1
, &(res
->x
.p
->arr
[i
]));
1847 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1848 /* Product of two rational numbers */
1849 value_multiply(res
->d
,e1
->d
,res
->d
);
1850 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1851 reduce_constant(res
);
1855 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1856 /* Product of an expression (polynomial or peririodic) and a rational number */
1862 /* Product of a rationel number and an expression (polynomial or peririodic) */
1864 i
= type_offset(res
->x
.p
);
1865 for (; i
<res
->x
.p
->size
; i
++)
1866 emul(e1
, &res
->x
.p
->arr
[i
]);
1876 /* Frees mask content ! */
1877 void emask(evalue
*mask
, evalue
*res
) {
1884 if (EVALUE_IS_ZERO(*res
)) {
1885 free_evalue_refs(mask
);
1889 assert(value_zero_p(mask
->d
));
1890 assert(mask
->x
.p
->type
== partition
);
1891 assert(value_zero_p(res
->d
));
1892 assert(res
->x
.p
->type
== partition
);
1893 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1894 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1895 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1896 pos
= res
->x
.p
->pos
;
1898 s
= (struct section
*)
1899 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1900 sizeof(struct section
));
1904 evalue_set_si(&mone
, -1, 1);
1907 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1908 assert(mask
->x
.p
->size
>= 2);
1909 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1910 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1912 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1914 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1923 value_init(s
[n
].E
.d
);
1924 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1928 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1929 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1932 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1933 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1934 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1935 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1937 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1938 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1944 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1945 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1947 value_init(s
[n
].E
.d
);
1948 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1949 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1955 /* Just ignore; this may have been previously masked off */
1957 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1961 free_evalue_refs(&mone
);
1962 free_evalue_refs(mask
);
1963 free_evalue_refs(res
);
1966 evalue_set_si(res
, 0, 1);
1968 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1969 for (j
= 0; j
< n
; ++j
) {
1970 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1971 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1972 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1979 void evalue_copy(evalue
*dst
, const evalue
*src
)
1981 value_assign(dst
->d
, src
->d
);
1982 if(value_notzero_p(src
->d
)) {
1983 value_init(dst
->x
.n
);
1984 value_assign(dst
->x
.n
, src
->x
.n
);
1986 dst
->x
.p
= ecopy(src
->x
.p
);
1989 evalue
*evalue_dup(const evalue
*e
)
1991 evalue
*res
= ALLOC(evalue
);
1993 evalue_copy(res
, e
);
1997 enode
*new_enode(enode_type type
,int size
,int pos
) {
2003 fprintf(stderr
, "Allocating enode of size 0 !\n" );
2006 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
2010 for(i
=0; i
<size
; i
++) {
2011 value_init(res
->arr
[i
].d
);
2012 value_set_si(res
->arr
[i
].d
,0);
2013 res
->arr
[i
].x
.p
= 0;
2018 enode
*ecopy(enode
*e
) {
2023 res
= new_enode(e
->type
,e
->size
,e
->pos
);
2024 for(i
=0;i
<e
->size
;++i
) {
2025 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
2026 if(value_zero_p(res
->arr
[i
].d
))
2027 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
2028 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
2029 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
2031 value_init(res
->arr
[i
].x
.n
);
2032 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
2038 int ecmp(const evalue
*e1
, const evalue
*e2
)
2044 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
2048 value_multiply(m
, e1
->x
.n
, e2
->d
);
2049 value_multiply(m2
, e2
->x
.n
, e1
->d
);
2051 if (value_lt(m
, m2
))
2053 else if (value_gt(m
, m2
))
2063 if (value_notzero_p(e1
->d
))
2065 if (value_notzero_p(e2
->d
))
2071 if (p1
->type
!= p2
->type
)
2072 return p1
->type
- p2
->type
;
2073 if (p1
->pos
!= p2
->pos
)
2074 return p1
->pos
- p2
->pos
;
2075 if (p1
->size
!= p2
->size
)
2076 return p1
->size
- p2
->size
;
2078 for (i
= p1
->size
-1; i
>= 0; --i
)
2079 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
2085 int eequal(const evalue
*e1
, const evalue
*e2
)
2090 if (value_ne(e1
->d
,e2
->d
))
2093 /* e1->d == e2->d */
2094 if (value_notzero_p(e1
->d
)) {
2095 if (value_ne(e1
->x
.n
,e2
->x
.n
))
2098 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2102 /* e1->d == e2->d == 0 */
2105 if (p1
->type
!= p2
->type
) return 0;
2106 if (p1
->size
!= p2
->size
) return 0;
2107 if (p1
->pos
!= p2
->pos
) return 0;
2108 for (i
=0; i
<p1
->size
; i
++)
2109 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
2114 void free_evalue_refs(evalue
*e
) {
2119 if (EVALUE_IS_DOMAIN(*e
)) {
2120 Domain_Free(EVALUE_DOMAIN(*e
));
2123 } else if (value_pos_p(e
->d
)) {
2125 /* 'e' stores a constant */
2127 value_clear(e
->x
.n
);
2130 assert(value_zero_p(e
->d
));
2133 if (!p
) return; /* null pointer */
2134 for (i
=0; i
<p
->size
; i
++) {
2135 free_evalue_refs(&(p
->arr
[i
]));
2139 } /* free_evalue_refs */
2141 void evalue_free(evalue
*e
)
2143 free_evalue_refs(e
);
2147 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
2148 Vector
* val
, evalue
*res
)
2150 unsigned nparam
= periods
->Size
;
2153 double d
= compute_evalue(e
, val
->p
);
2154 d
*= VALUE_TO_DOUBLE(m
);
2159 value_assign(res
->d
, m
);
2160 value_init(res
->x
.n
);
2161 value_set_double(res
->x
.n
, d
);
2162 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
2165 if (value_one_p(periods
->p
[p
]))
2166 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
2171 value_assign(tmp
, periods
->p
[p
]);
2172 value_set_si(res
->d
, 0);
2173 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
2175 value_decrement(tmp
, tmp
);
2176 value_assign(val
->p
[p
], tmp
);
2177 mod2table_r(e
, periods
, m
, p
+1, val
,
2178 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
2179 } while (value_pos_p(tmp
));
2185 static void rel2table(evalue
*e
, int zero
)
2187 if (value_pos_p(e
->d
)) {
2188 if (value_zero_p(e
->x
.n
) == zero
)
2189 value_set_si(e
->x
.n
, 1);
2191 value_set_si(e
->x
.n
, 0);
2192 value_set_si(e
->d
, 1);
2195 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
2196 rel2table(&e
->x
.p
->arr
[i
], zero
);
2200 void evalue_mod2table(evalue
*e
, int nparam
)
2205 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2208 for (i
=0; i
<p
->size
; i
++) {
2209 evalue_mod2table(&(p
->arr
[i
]), nparam
);
2211 if (p
->type
== relation
) {
2216 evalue_copy(©
, &p
->arr
[0]);
2218 rel2table(&p
->arr
[0], 1);
2219 emul(&p
->arr
[0], &p
->arr
[1]);
2221 rel2table(©
, 0);
2222 emul(©
, &p
->arr
[2]);
2223 eadd(&p
->arr
[2], &p
->arr
[1]);
2224 free_evalue_refs(&p
->arr
[2]);
2225 free_evalue_refs(©
);
2227 free_evalue_refs(&p
->arr
[0]);
2231 } else if (p
->type
== fractional
) {
2232 Vector
*periods
= Vector_Alloc(nparam
);
2233 Vector
*val
= Vector_Alloc(nparam
);
2239 value_set_si(tmp
, 1);
2240 Vector_Set(periods
->p
, 1, nparam
);
2241 Vector_Set(val
->p
, 0, nparam
);
2242 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
2245 assert(p
->type
== polynomial
);
2246 assert(p
->size
== 2);
2247 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
2248 value_lcm(tmp
, tmp
, p
->arr
[1].d
);
2250 value_lcm(tmp
, tmp
, ev
->d
);
2252 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
2255 evalue_set_si(&res
, 0, 1);
2256 /* Compute the polynomial using Horner's rule */
2257 for (i
=p
->size
-1;i
>1;i
--) {
2258 eadd(&p
->arr
[i
], &res
);
2261 eadd(&p
->arr
[1], &res
);
2263 free_evalue_refs(e
);
2264 free_evalue_refs(&EP
);
2269 Vector_Free(periods
);
2271 } /* evalue_mod2table */
2273 /********************************************************/
2274 /* function in domain */
2275 /* check if the parameters in list_args */
2276 /* verifies the constraints of Domain P */
2277 /********************************************************/
2278 int in_domain(Polyhedron
*P
, Value
*list_args
)
2281 Value v
; /* value of the constraint of a row when
2282 parameters are instantiated*/
2286 for (row
= 0; row
< P
->NbConstraints
; row
++) {
2287 Inner_Product(P
->Constraint
[row
]+1, list_args
, P
->Dimension
, &v
);
2288 value_addto(v
, v
, P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
2289 if (value_neg_p(v
) ||
2290 value_zero_p(P
->Constraint
[row
][0]) && value_notzero_p(v
)) {
2297 return in
|| (P
->next
&& in_domain(P
->next
, list_args
));
2300 /****************************************************/
2301 /* function compute enode */
2302 /* compute the value of enode p with parameters */
2303 /* list "list_args */
2304 /* compute the polynomial or the periodic */
2305 /****************************************************/
2307 static double compute_enode(enode
*p
, Value
*list_args
) {
2319 if (p
->type
== polynomial
) {
2321 value_assign(param
,list_args
[p
->pos
-1]);
2323 /* Compute the polynomial using Horner's rule */
2324 for (i
=p
->size
-1;i
>0;i
--) {
2325 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2326 res
*=VALUE_TO_DOUBLE(param
);
2328 res
+=compute_evalue(&p
->arr
[0],list_args
);
2330 else if (p
->type
== fractional
) {
2331 double d
= compute_evalue(&p
->arr
[0], list_args
);
2332 d
-= floor(d
+1e-10);
2334 /* Compute the polynomial using Horner's rule */
2335 for (i
=p
->size
-1;i
>1;i
--) {
2336 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2339 res
+=compute_evalue(&p
->arr
[1],list_args
);
2341 else if (p
->type
== flooring
) {
2342 double d
= compute_evalue(&p
->arr
[0], list_args
);
2345 /* Compute the polynomial using Horner's rule */
2346 for (i
=p
->size
-1;i
>1;i
--) {
2347 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2350 res
+=compute_evalue(&p
->arr
[1],list_args
);
2352 else if (p
->type
== periodic
) {
2353 value_assign(m
,list_args
[p
->pos
-1]);
2355 /* Choose the right element of the periodic */
2356 value_set_si(param
,p
->size
);
2357 value_pmodulus(m
,m
,param
);
2358 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2360 else if (p
->type
== relation
) {
2361 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2362 res
= compute_evalue(&p
->arr
[1], list_args
);
2363 else if (p
->size
> 2)
2364 res
= compute_evalue(&p
->arr
[2], list_args
);
2366 else if (p
->type
== partition
) {
2367 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2368 Value
*vals
= list_args
;
2371 for (i
= 0; i
< dim
; ++i
) {
2372 value_init(vals
[i
]);
2374 value_assign(vals
[i
], list_args
[i
]);
2377 for (i
= 0; i
< p
->size
/2; ++i
)
2378 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2379 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2383 for (i
= 0; i
< dim
; ++i
)
2384 value_clear(vals
[i
]);
2393 } /* compute_enode */
2395 /*************************************************/
2396 /* return the value of Ehrhart Polynomial */
2397 /* It returns a double, because since it is */
2398 /* a recursive function, some intermediate value */
2399 /* might not be integral */
2400 /*************************************************/
2402 double compute_evalue(const evalue
*e
, Value
*list_args
)
2406 if (value_notzero_p(e
->d
)) {
2407 if (value_notone_p(e
->d
))
2408 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2410 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2413 res
= compute_enode(e
->x
.p
,list_args
);
2415 } /* compute_evalue */
2418 /****************************************************/
2419 /* function compute_poly : */
2420 /* Check for the good validity domain */
2421 /* return the number of point in the Polyhedron */
2422 /* in allocated memory */
2423 /* Using the Ehrhart pseudo-polynomial */
2424 /****************************************************/
2425 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2428 /* double d; int i; */
2430 tmp
= (Value
*) malloc (sizeof(Value
));
2431 assert(tmp
!= NULL
);
2433 value_set_si(*tmp
,0);
2436 return(tmp
); /* no ehrhart polynomial */
2437 if(en
->ValidityDomain
) {
2438 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2439 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2444 return(tmp
); /* no Validity Domain */
2446 if(in_domain(en
->ValidityDomain
,list_args
)) {
2448 #ifdef EVAL_EHRHART_DEBUG
2449 Print_Domain(stdout
,en
->ValidityDomain
);
2450 print_evalue(stdout
,&en
->EP
);
2453 /* d = compute_evalue(&en->EP,list_args);
2455 printf("(double)%lf = %d\n", d, i ); */
2456 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2462 value_set_si(*tmp
,0);
2463 return(tmp
); /* no compatible domain with the arguments */
2464 } /* compute_poly */
2466 static evalue
*eval_polynomial(const enode
*p
, int offset
,
2467 evalue
*base
, Value
*values
)
2472 res
= evalue_zero();
2473 for (i
= p
->size
-1; i
> offset
; --i
) {
2474 c
= evalue_eval(&p
->arr
[i
], values
);
2479 c
= evalue_eval(&p
->arr
[offset
], values
);
2486 evalue
*evalue_eval(const evalue
*e
, Value
*values
)
2493 if (value_notzero_p(e
->d
)) {
2494 res
= ALLOC(evalue
);
2496 evalue_copy(res
, e
);
2499 switch (e
->x
.p
->type
) {
2501 value_init(param
.x
.n
);
2502 value_assign(param
.x
.n
, values
[e
->x
.p
->pos
-1]);
2503 value_init(param
.d
);
2504 value_set_si(param
.d
, 1);
2506 res
= eval_polynomial(e
->x
.p
, 0, ¶m
, values
);
2507 free_evalue_refs(¶m
);
2510 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2511 mpz_fdiv_r(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2513 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2514 evalue_free(param2
);
2517 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2518 mpz_fdiv_q(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2519 value_set_si(param2
->d
, 1);
2521 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2522 evalue_free(param2
);
2525 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2526 if (value_zero_p(param2
->x
.n
))
2527 res
= evalue_eval(&e
->x
.p
->arr
[1], values
);
2528 else if (e
->x
.p
->size
> 2)
2529 res
= evalue_eval(&e
->x
.p
->arr
[2], values
);
2531 res
= evalue_zero();
2532 evalue_free(param2
);
2535 assert(e
->x
.p
->pos
== EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
);
2536 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2537 if (in_domain(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), values
)) {
2538 res
= evalue_eval(&e
->x
.p
->arr
[2*i
+1], values
);
2542 res
= evalue_zero();
2550 size_t value_size(Value v
) {
2551 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2552 * sizeof(v
[0]._mp_d
[0]);
2555 size_t domain_size(Polyhedron
*D
)
2558 size_t s
= sizeof(*D
);
2560 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2561 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2562 s
+= value_size(D
->Constraint
[i
][j
]);
2565 for (i = 0; i < D->NbRays; ++i)
2566 for (j = 0; j < D->Dimension+2; ++j)
2567 s += value_size(D->Ray[i][j]);
2570 return D
->next
? s
+domain_size(D
->next
) : s
;
2573 size_t enode_size(enode
*p
) {
2574 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2577 if (p
->type
== partition
)
2578 for (i
= 0; i
< p
->size
/2; ++i
) {
2579 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2580 s
+= evalue_size(&p
->arr
[2*i
+1]);
2583 for (i
= 0; i
< p
->size
; ++i
) {
2584 s
+= evalue_size(&p
->arr
[i
]);
2589 size_t evalue_size(evalue
*e
)
2591 size_t s
= sizeof(*e
);
2592 s
+= value_size(e
->d
);
2593 if (value_notzero_p(e
->d
))
2594 s
+= value_size(e
->x
.n
);
2596 s
+= enode_size(e
->x
.p
);
2600 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2602 evalue
*found
= NULL
;
2607 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2610 value_init(offset
.d
);
2611 value_init(offset
.x
.n
);
2612 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2613 value_lcm(offset
.d
, m
, offset
.d
);
2614 value_set_si(offset
.x
.n
, 1);
2617 evalue_copy(©
, cst
);
2620 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2622 if (eequal(base
, &e
->x
.p
->arr
[0]))
2623 found
= &e
->x
.p
->arr
[0];
2625 value_set_si(offset
.x
.n
, -2);
2628 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2630 if (eequal(base
, &e
->x
.p
->arr
[0]))
2633 free_evalue_refs(cst
);
2634 free_evalue_refs(&offset
);
2637 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2638 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2643 static evalue
*find_relation_pair(evalue
*e
)
2646 evalue
*found
= NULL
;
2648 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2651 if (e
->x
.p
->type
== fractional
) {
2656 poly_denom(&e
->x
.p
->arr
[0], &m
);
2658 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2659 cst
= &cst
->x
.p
->arr
[0])
2662 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2663 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2668 i
= e
->x
.p
->type
== relation
;
2669 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2670 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2675 void evalue_mod2relation(evalue
*e
) {
2678 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2681 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2682 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2683 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2684 value_clear(e
->x
.p
->arr
[2*i
].d
);
2685 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2687 if (2*i
< e
->x
.p
->size
) {
2688 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2689 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2694 if (e
->x
.p
->size
== 0) {
2696 evalue_set_si(e
, 0, 1);
2702 while ((d
= find_relation_pair(e
)) != NULL
) {
2706 value_init(split
.d
);
2707 value_set_si(split
.d
, 0);
2708 split
.x
.p
= new_enode(relation
, 3, 0);
2709 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2710 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2712 ev
= &split
.x
.p
->arr
[0];
2713 value_set_si(ev
->d
, 0);
2714 ev
->x
.p
= new_enode(fractional
, 3, -1);
2715 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2716 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2717 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2723 free_evalue_refs(&split
);
2727 static int evalue_comp(const void * a
, const void * b
)
2729 const evalue
*e1
= *(const evalue
**)a
;
2730 const evalue
*e2
= *(const evalue
**)b
;
2731 return ecmp(e1
, e2
);
2734 void evalue_combine(evalue
*e
)
2741 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2744 NALLOC(evs
, e
->x
.p
->size
/2);
2745 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2746 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2747 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2748 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2749 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2750 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2751 value_clear(p
->arr
[2*k
].d
);
2752 value_clear(p
->arr
[2*k
+1].d
);
2753 p
->arr
[2*k
] = *(evs
[i
]-1);
2754 p
->arr
[2*k
+1] = *(evs
[i
]);
2757 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2760 value_clear((evs
[i
]-1)->d
);
2764 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2765 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2766 free_evalue_refs(evs
[i
]);
2770 for (i
= 2*k
; i
< p
->size
; ++i
)
2771 value_clear(p
->arr
[i
].d
);
2778 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2780 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2782 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2785 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2786 Polyhedron
*D
, *N
, **P
;
2789 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2796 if (D
->NbEq
<= H
->NbEq
) {
2802 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2803 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2804 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2805 reduce_evalue(&tmp
);
2806 if (value_notzero_p(tmp
.d
) ||
2807 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2810 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2811 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2814 free_evalue_refs(&tmp
);
2820 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2822 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2824 value_clear(e
->x
.p
->arr
[2*i
].d
);
2825 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2827 if (2*i
< e
->x
.p
->size
) {
2828 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2829 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2836 H
= DomainConvex(D
, 0);
2837 E
= DomainDifference(H
, D
, 0);
2839 D
= DomainDifference(H
, E
, 0);
2842 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2846 /* Use smallest representative for coefficients in affine form in
2847 * argument of fractional.
2848 * Since any change will make the argument non-standard,
2849 * the containing evalue will have to be reduced again afterward.
2851 static void fractional_minimal_coefficients(enode
*p
)
2857 assert(p
->type
== fractional
);
2859 while (value_zero_p(pp
->d
)) {
2860 assert(pp
->x
.p
->type
== polynomial
);
2861 assert(pp
->x
.p
->size
== 2);
2862 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2863 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2864 if (value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2865 value_subtract(pp
->x
.p
->arr
[1].x
.n
,
2866 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2867 pp
= &pp
->x
.p
->arr
[0];
2873 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2878 unsigned dim
= D
->Dimension
;
2879 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2882 assert(p
->type
== fractional
|| p
->type
== flooring
);
2883 value_set_si(T
->p
[1][dim
], 1);
2884 evalue_extract_affine(&p
->arr
[0], T
->p
[0], &T
->p
[0][dim
], d
);
2885 I
= DomainImage(D
, T
, 0);
2886 H
= DomainConvex(I
, 0);
2896 static void replace_by_affine(evalue
*e
, Value offset
)
2903 value_init(inc
.x
.n
);
2904 value_set_si(inc
.d
, 1);
2905 value_oppose(inc
.x
.n
, offset
);
2906 eadd(&inc
, &p
->arr
[0]);
2907 reorder_terms_about(p
, &p
->arr
[0]); /* frees arr[0] */
2911 free_evalue_refs(&inc
);
2914 int evalue_range_reduction_in_domain(evalue
*e
, Polyhedron
*D
)
2923 if (value_notzero_p(e
->d
))
2928 if (p
->type
== relation
) {
2935 fractional_minimal_coefficients(p
->arr
[0].x
.p
);
2936 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
);
2937 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2938 equal
= value_eq(min
, max
);
2939 mpz_cdiv_q(min
, min
, d
);
2940 mpz_fdiv_q(max
, max
, d
);
2942 if (bounded
&& value_gt(min
, max
)) {
2948 evalue_set_si(e
, 0, 1);
2951 free_evalue_refs(&(p
->arr
[1]));
2952 free_evalue_refs(&(p
->arr
[0]));
2958 return r
? r
: evalue_range_reduction_in_domain(e
, D
);
2959 } else if (bounded
&& equal
) {
2962 free_evalue_refs(&(p
->arr
[2]));
2965 free_evalue_refs(&(p
->arr
[0]));
2971 return evalue_range_reduction_in_domain(e
, D
);
2972 } else if (bounded
&& value_eq(min
, max
)) {
2973 /* zero for a single value */
2975 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2976 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2977 value_multiply(min
, min
, d
);
2978 value_subtract(M
->p
[0][D
->Dimension
+1],
2979 M
->p
[0][D
->Dimension
+1], min
);
2980 E
= DomainAddConstraints(D
, M
, 0);
2986 r
= evalue_range_reduction_in_domain(&p
->arr
[1], E
);
2988 r
|= evalue_range_reduction_in_domain(&p
->arr
[2], D
);
2990 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2998 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
3001 i
= p
->type
== relation
? 1 :
3002 p
->type
== fractional
? 1 : 0;
3003 for (; i
<p
->size
; i
++)
3004 r
|= evalue_range_reduction_in_domain(&p
->arr
[i
], D
);
3006 if (p
->type
!= fractional
) {
3007 if (r
&& p
->type
== polynomial
) {
3010 value_set_si(f
.d
, 0);
3011 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3012 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3013 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3014 reorder_terms_about(p
, &f
);
3025 fractional_minimal_coefficients(p
);
3026 I
= polynomial_projection(p
, D
, &d
, NULL
);
3027 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
3028 mpz_fdiv_q(min
, min
, d
);
3029 mpz_fdiv_q(max
, max
, d
);
3030 value_subtract(d
, max
, min
);
3032 if (bounded
&& value_eq(min
, max
)) {
3033 replace_by_affine(e
, min
);
3035 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
3036 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
3037 * See pages 199-200 of PhD thesis.
3045 value_set_si(rem
.d
, 0);
3046 rem
.x
.p
= new_enode(fractional
, 3, -1);
3047 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
3048 value_clear(rem
.x
.p
->arr
[1].d
);
3049 value_clear(rem
.x
.p
->arr
[2].d
);
3050 rem
.x
.p
->arr
[1] = p
->arr
[1];
3051 rem
.x
.p
->arr
[2] = p
->arr
[2];
3052 for (i
= 3; i
< p
->size
; ++i
)
3053 p
->arr
[i
-2] = p
->arr
[i
];
3057 value_init(inc
.x
.n
);
3058 value_set_si(inc
.d
, 1);
3059 value_oppose(inc
.x
.n
, min
);
3062 evalue_copy(&t
, &p
->arr
[0]);
3066 value_set_si(f
.d
, 0);
3067 f
.x
.p
= new_enode(fractional
, 3, -1);
3068 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3069 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3070 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
3072 value_init(factor
.d
);
3073 evalue_set_si(&factor
, -1, 1);
3079 value_clear(f
.x
.p
->arr
[1].x
.n
);
3080 value_clear(f
.x
.p
->arr
[2].x
.n
);
3081 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3082 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3086 reorder_terms(&rem
);
3093 free_evalue_refs(&inc
);
3094 free_evalue_refs(&t
);
3095 free_evalue_refs(&f
);
3096 free_evalue_refs(&factor
);
3097 free_evalue_refs(&rem
);
3099 evalue_range_reduction_in_domain(e
, D
);
3103 _reduce_evalue(&p
->arr
[0], 0, 1);
3115 void evalue_range_reduction(evalue
*e
)
3118 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3121 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3122 if (evalue_range_reduction_in_domain(&e
->x
.p
->arr
[2*i
+1],
3123 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
3124 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3126 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
3127 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
3128 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3129 value_clear(e
->x
.p
->arr
[2*i
].d
);
3131 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
3132 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
3140 Enumeration
* partition2enumeration(evalue
*EP
)
3143 Enumeration
*en
, *res
= NULL
;
3145 if (EVALUE_IS_ZERO(*EP
)) {
3150 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
3151 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
3152 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
3155 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
3156 value_clear(EP
->x
.p
->arr
[2*i
].d
);
3157 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
3165 int evalue_frac2floor_in_domain3(evalue
*e
, Polyhedron
*D
, int shift
)
3174 if (value_notzero_p(e
->d
))
3179 i
= p
->type
== relation
? 1 :
3180 p
->type
== fractional
? 1 : 0;
3181 for (; i
<p
->size
; i
++)
3182 r
|= evalue_frac2floor_in_domain3(&p
->arr
[i
], D
, shift
);
3184 if (p
->type
!= fractional
) {
3185 if (r
&& p
->type
== polynomial
) {
3188 value_set_si(f
.d
, 0);
3189 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3190 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3191 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3192 reorder_terms_about(p
, &f
);
3202 I
= polynomial_projection(p
, D
, &d
, NULL
);
3205 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3208 assert(I
->NbEq
== 0); /* Should have been reduced */
3211 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3212 if (value_pos_p(I
->Constraint
[i
][1]))
3215 if (i
< I
->NbConstraints
) {
3217 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3218 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3219 if (value_neg_p(min
)) {
3221 mpz_fdiv_q(min
, min
, d
);
3222 value_init(offset
.d
);
3223 value_set_si(offset
.d
, 1);
3224 value_init(offset
.x
.n
);
3225 value_oppose(offset
.x
.n
, min
);
3226 eadd(&offset
, &p
->arr
[0]);
3227 free_evalue_refs(&offset
);
3237 value_set_si(fl
.d
, 0);
3238 fl
.x
.p
= new_enode(flooring
, 3, -1);
3239 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
3240 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
3241 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
3243 eadd(&fl
, &p
->arr
[0]);
3244 reorder_terms_about(p
, &p
->arr
[0]);
3248 free_evalue_refs(&fl
);
3253 int evalue_frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
3255 return evalue_frac2floor_in_domain3(e
, D
, 1);
3258 void evalue_frac2floor2(evalue
*e
, int shift
)
3261 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
3263 if (evalue_frac2floor_in_domain3(e
, NULL
, 0))
3269 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3270 if (evalue_frac2floor_in_domain3(&e
->x
.p
->arr
[2*i
+1],
3271 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), shift
))
3272 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3275 void evalue_frac2floor(evalue
*e
)
3277 evalue_frac2floor2(e
, 1);
3280 /* Add a new variable with lower bound 1 and upper bound specified
3281 * by row. If negative is true, then the new variable has upper
3282 * bound -1 and lower bound specified by row.
3284 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
3285 Vector
*row
, int negative
)
3289 int nparam
= D
->Dimension
- nvar
;
3292 nr
= D
->NbConstraints
+ 2;
3293 nc
= D
->Dimension
+ 2 + 1;
3294 C
= Matrix_Alloc(nr
, nc
);
3295 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
3296 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
3297 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3298 D
->Dimension
+ 1 - nvar
);
3303 nc
= C
->NbColumns
+ 1;
3304 C
= Matrix_Alloc(nr
, nc
);
3305 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
3306 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
3307 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3308 oldC
->NbColumns
- 1 - nvar
);
3311 value_set_si(C
->p
[nr
-2][0], 1);
3313 value_set_si(C
->p
[nr
-2][1 + nvar
], -1);
3315 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
3316 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
3318 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
3319 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
3325 static void floor2frac_r(evalue
*e
, int nvar
)
3332 if (value_notzero_p(e
->d
))
3337 assert(p
->type
== flooring
);
3338 for (i
= 1; i
< p
->size
; i
++)
3339 floor2frac_r(&p
->arr
[i
], nvar
);
3341 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
3342 assert(pp
->x
.p
->type
== polynomial
);
3343 pp
->x
.p
->pos
-= nvar
;
3347 value_set_si(f
.d
, 0);
3348 f
.x
.p
= new_enode(fractional
, 3, -1);
3349 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3350 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3351 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3353 eadd(&f
, &p
->arr
[0]);
3354 reorder_terms_about(p
, &p
->arr
[0]);
3358 free_evalue_refs(&f
);
3361 /* Convert flooring back to fractional and shift position
3362 * of the parameters by nvar
3364 static void floor2frac(evalue
*e
, int nvar
)
3366 floor2frac_r(e
, nvar
);
3370 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
3373 int nparam
= D
->Dimension
- nvar
;
3377 D
= Constraints2Polyhedron(C
, 0);
3381 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
3383 /* Double check that D was not unbounded. */
3384 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
3392 static evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
3393 int *signs
, Matrix
*C
, unsigned MaxRays
)
3399 evalue
*factor
= NULL
;
3403 if (EVALUE_IS_ZERO(*e
))
3407 Polyhedron
*DD
= Disjoint_Domain(D
, 0, MaxRays
);
3414 res
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3417 for (Q
= DD
; Q
; Q
= DD
) {
3423 t
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3436 if (value_notzero_p(e
->d
)) {
3439 t
= esum_over_domain_cst(nvar
, D
, C
);
3441 if (!EVALUE_IS_ONE(*e
))
3447 switch (e
->x
.p
->type
) {
3449 evalue
*pp
= &e
->x
.p
->arr
[0];
3451 if (pp
->x
.p
->pos
> nvar
) {
3452 /* remainder is independent of the summated vars */
3458 floor2frac(&f
, nvar
);
3460 t
= esum_over_domain_cst(nvar
, D
, C
);
3464 free_evalue_refs(&f
);
3469 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3470 poly_denom(pp
, &row
->p
[1 + nvar
]);
3471 value_set_si(row
->p
[0], 1);
3472 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3473 pp
= &pp
->x
.p
->arr
[0]) {
3475 assert(pp
->x
.p
->type
== polynomial
);
3477 if (pos
>= 1 + nvar
)
3479 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3480 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3481 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3483 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3484 value_division(row
->p
[1 + D
->Dimension
+ 1],
3485 row
->p
[1 + D
->Dimension
+ 1],
3487 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3488 row
->p
[1 + D
->Dimension
+ 1],
3490 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3494 int pos
= e
->x
.p
->pos
;
3497 factor
= ALLOC(evalue
);
3498 value_init(factor
->d
);
3499 value_set_si(factor
->d
, 0);
3500 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3501 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3502 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3506 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3507 negative
= signs
[pos
-1] < 0;
3508 value_set_si(row
->p
[0], 1);
3510 value_set_si(row
->p
[pos
], -1);
3511 value_set_si(row
->p
[1 + nvar
], 1);
3513 value_set_si(row
->p
[pos
], 1);
3514 value_set_si(row
->p
[1 + nvar
], -1);
3522 offset
= type_offset(e
->x
.p
);
3524 res
= esum_over_domain(&e
->x
.p
->arr
[offset
], nvar
, D
, signs
, C
, MaxRays
);
3528 evalue_copy(&cum
, factor
);
3532 for (i
= 1; offset
+i
< e
->x
.p
->size
; ++i
) {
3536 C
= esum_add_constraint(nvar
, D
, C
, row
, negative
);
3542 Vector_Print(stderr, P_VALUE_FMT, row);
3544 Matrix_Print(stderr, P_VALUE_FMT, C);
3546 t
= esum_over_domain(&e
->x
.p
->arr
[offset
+i
], nvar
, D
, signs
, C
, MaxRays
);
3551 if (negative
&& (i
% 2))
3561 if (factor
&& offset
+i
+1 < e
->x
.p
->size
)
3568 free_evalue_refs(&cum
);
3569 evalue_free(factor
);
3580 static void domain_signs(Polyhedron
*D
, int *signs
)
3584 POL_ENSURE_VERTICES(D
);
3585 for (j
= 0; j
< D
->Dimension
; ++j
) {
3587 for (k
= 0; k
< D
->NbRays
; ++k
) {
3588 signs
[j
] = value_sign(D
->Ray
[k
][1+j
]);
3595 static void shift_floor_in_domain(evalue
*e
, Polyhedron
*D
)
3602 if (value_notzero_p(e
->d
))
3607 for (i
= type_offset(p
); i
< p
->size
; ++i
)
3608 shift_floor_in_domain(&p
->arr
[i
], D
);
3610 if (p
->type
!= flooring
)
3616 I
= polynomial_projection(p
, D
, &d
, NULL
);
3617 assert(I
->NbEq
== 0); /* Should have been reduced */
3619 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3620 if (value_pos_p(I
->Constraint
[i
][1]))
3622 assert(i
< I
->NbConstraints
);
3623 if (i
< I
->NbConstraints
) {
3624 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3625 mpz_fdiv_q(m
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3626 if (value_neg_p(m
)) {
3627 /* replace [e] by [e-m]+m such that e-m >= 0 */
3632 value_set_si(f
.d
, 1);
3633 value_oppose(f
.x
.n
, m
);
3634 eadd(&f
, &p
->arr
[0]);
3637 value_set_si(f
.d
, 0);
3638 f
.x
.p
= new_enode(flooring
, 3, -1);
3639 value_clear(f
.x
.p
->arr
[0].d
);
3640 f
.x
.p
->arr
[0] = p
->arr
[0];
3641 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
3642 value_set_si(f
.x
.p
->arr
[1].d
, 1);
3643 value_init(f
.x
.p
->arr
[1].x
.n
);
3644 value_assign(f
.x
.p
->arr
[1].x
.n
, m
);
3645 reorder_terms_about(p
, &f
);
3656 /* Make arguments of all floors non-negative */
3657 static void shift_floor_arguments(evalue
*e
)
3661 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3664 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3665 shift_floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
3666 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3669 evalue
*evalue_sum(evalue
*e
, int nvar
, unsigned MaxRays
)
3673 evalue
*res
= ALLOC(evalue
);
3677 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3678 evalue_copy(res
, e
);
3682 evalue_split_domains_into_orthants(e
, MaxRays
);
3683 evalue_frac2floor2(e
, 0);
3684 evalue_set_si(res
, 0, 1);
3686 assert(value_zero_p(e
->d
));
3687 assert(e
->x
.p
->type
== partition
);
3688 shift_floor_arguments(e
);
3690 assert(e
->x
.p
->size
>= 2);
3691 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3693 signs
= alloca(sizeof(int) * dim
);
3695 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3697 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3698 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3699 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
, 0,
3710 evalue
*esum(evalue
*e
, int nvar
)
3712 return evalue_sum(e
, nvar
, 0);
3715 /* Initial silly implementation */
3716 void eor(evalue
*e1
, evalue
*res
)
3722 evalue_set_si(&mone
, -1, 1);
3724 evalue_copy(&E
, res
);
3730 free_evalue_refs(&E
);
3731 free_evalue_refs(&mone
);
3734 /* computes denominator of polynomial evalue
3735 * d should point to a value initialized to 1
3737 void evalue_denom(const evalue
*e
, Value
*d
)
3741 if (value_notzero_p(e
->d
)) {
3742 value_lcm(*d
, *d
, e
->d
);
3745 assert(e
->x
.p
->type
== polynomial
||
3746 e
->x
.p
->type
== fractional
||
3747 e
->x
.p
->type
== flooring
);
3748 offset
= type_offset(e
->x
.p
);
3749 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3750 evalue_denom(&e
->x
.p
->arr
[i
], d
);
3753 /* Divides the evalue e by the integer n */
3754 void evalue_div(evalue
*e
, Value n
)
3758 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3761 if (value_notzero_p(e
->d
)) {
3764 value_multiply(e
->d
, e
->d
, n
);
3765 value_gcd(gc
, e
->x
.n
, e
->d
);
3766 if (value_notone_p(gc
)) {
3767 value_division(e
->d
, e
->d
, gc
);
3768 value_division(e
->x
.n
, e
->x
.n
, gc
);
3773 if (e
->x
.p
->type
== partition
) {
3774 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3775 evalue_div(&e
->x
.p
->arr
[2*i
+1], n
);
3778 offset
= type_offset(e
->x
.p
);
3779 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3780 evalue_div(&e
->x
.p
->arr
[i
], n
);
3783 /* Multiplies the evalue e by the integer n */
3784 void evalue_mul(evalue
*e
, Value n
)
3788 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3791 if (value_notzero_p(e
->d
)) {
3794 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3795 value_gcd(gc
, e
->x
.n
, e
->d
);
3796 if (value_notone_p(gc
)) {
3797 value_division(e
->d
, e
->d
, gc
);
3798 value_division(e
->x
.n
, e
->x
.n
, gc
);
3803 if (e
->x
.p
->type
== partition
) {
3804 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3805 evalue_mul(&e
->x
.p
->arr
[2*i
+1], n
);
3808 offset
= type_offset(e
->x
.p
);
3809 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3810 evalue_mul(&e
->x
.p
->arr
[i
], n
);
3813 /* Multiplies the evalue e by the n/d */
3814 void evalue_mul_div(evalue
*e
, Value n
, Value d
)
3818 if ((value_one_p(n
) && value_one_p(d
)) || EVALUE_IS_ZERO(*e
))
3821 if (value_notzero_p(e
->d
)) {
3824 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3825 value_multiply(e
->d
, e
->d
, d
);
3826 value_gcd(gc
, e
->x
.n
, e
->d
);
3827 if (value_notone_p(gc
)) {
3828 value_division(e
->d
, e
->d
, gc
);
3829 value_division(e
->x
.n
, e
->x
.n
, gc
);
3834 if (e
->x
.p
->type
== partition
) {
3835 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3836 evalue_mul_div(&e
->x
.p
->arr
[2*i
+1], n
, d
);
3839 offset
= type_offset(e
->x
.p
);
3840 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3841 evalue_mul_div(&e
->x
.p
->arr
[i
], n
, d
);
3844 void evalue_negate(evalue
*e
)
3848 if (value_notzero_p(e
->d
)) {
3849 value_oppose(e
->x
.n
, e
->x
.n
);
3852 if (e
->x
.p
->type
== partition
) {
3853 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3854 evalue_negate(&e
->x
.p
->arr
[2*i
+1]);
3857 offset
= type_offset(e
->x
.p
);
3858 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3859 evalue_negate(&e
->x
.p
->arr
[i
]);
3862 void evalue_add_constant(evalue
*e
, const Value cst
)
3866 if (value_zero_p(e
->d
)) {
3867 if (e
->x
.p
->type
== partition
) {
3868 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3869 evalue_add_constant(&e
->x
.p
->arr
[2*i
+1], cst
);
3872 if (e
->x
.p
->type
== relation
) {
3873 for (i
= 1; i
< e
->x
.p
->size
; ++i
)
3874 evalue_add_constant(&e
->x
.p
->arr
[i
], cst
);
3878 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
3879 } while (value_zero_p(e
->d
));
3881 value_addmul(e
->x
.n
, cst
, e
->d
);
3884 static void evalue_frac2polynomial_r(evalue
*e
, int *signs
, int sign
, int in_frac
)
3889 int sign_odd
= sign
;
3891 if (value_notzero_p(e
->d
)) {
3892 if (in_frac
&& sign
* value_sign(e
->x
.n
) < 0) {
3893 value_set_si(e
->x
.n
, 0);
3894 value_set_si(e
->d
, 1);
3899 if (e
->x
.p
->type
== relation
) {
3900 for (i
= e
->x
.p
->size
-1; i
>= 1; --i
)
3901 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
, sign
, in_frac
);
3905 if (e
->x
.p
->type
== polynomial
)
3906 sign_odd
*= signs
[e
->x
.p
->pos
-1];
3907 offset
= type_offset(e
->x
.p
);
3908 evalue_frac2polynomial_r(&e
->x
.p
->arr
[offset
], signs
, sign
, in_frac
);
3909 in_frac
|= e
->x
.p
->type
== fractional
;
3910 for (i
= e
->x
.p
->size
-1; i
> offset
; --i
)
3911 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
,
3912 (i
- offset
) % 2 ? sign_odd
: sign
, in_frac
);
3914 if (e
->x
.p
->type
!= fractional
)
3917 /* replace { a/m } by (m-1)/m if sign != 0
3918 * and by (m-1)/(2m) if sign == 0
3922 evalue_denom(&e
->x
.p
->arr
[0], &d
);
3923 free_evalue_refs(&e
->x
.p
->arr
[0]);
3924 value_init(e
->x
.p
->arr
[0].d
);
3925 value_init(e
->x
.p
->arr
[0].x
.n
);
3927 value_addto(e
->x
.p
->arr
[0].d
, d
, d
);
3929 value_assign(e
->x
.p
->arr
[0].d
, d
);
3930 value_decrement(e
->x
.p
->arr
[0].x
.n
, d
);
3934 reorder_terms_about(p
, &p
->arr
[0]);
3940 /* Approximate the evalue in fractional representation by a polynomial.
3941 * If sign > 0, the result is an upper bound;
3942 * if sign < 0, the result is a lower bound;
3943 * if sign = 0, the result is an intermediate approximation.
3945 void evalue_frac2polynomial(evalue
*e
, int sign
, unsigned MaxRays
)
3950 if (value_notzero_p(e
->d
))
3952 assert(e
->x
.p
->type
== partition
);
3953 /* make sure all variables in the domains have a fixed sign */
3955 evalue_split_domains_into_orthants(e
, MaxRays
);
3956 if (EVALUE_IS_ZERO(*e
))
3960 assert(e
->x
.p
->size
>= 2);
3961 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3963 signs
= alloca(sizeof(int) * dim
);
3966 for (i
= 0; i
< dim
; ++i
)
3968 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3970 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3971 evalue_frac2polynomial_r(&e
->x
.p
->arr
[2*i
+1], signs
, sign
, 0);
3975 /* Split the domains of e (which is assumed to be a partition)
3976 * such that each resulting domain lies entirely in one orthant.
3978 void evalue_split_domains_into_orthants(evalue
*e
, unsigned MaxRays
)
3981 assert(value_zero_p(e
->d
));
3982 assert(e
->x
.p
->type
== partition
);
3983 assert(e
->x
.p
->size
>= 2);
3984 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3986 for (i
= 0; i
< dim
; ++i
) {
3989 C
= Matrix_Alloc(1, 1 + dim
+ 1);
3990 value_set_si(C
->p
[0][0], 1);
3991 value_init(split
.d
);
3992 value_set_si(split
.d
, 0);
3993 split
.x
.p
= new_enode(partition
, 4, dim
);
3994 value_set_si(C
->p
[0][1+i
], 1);
3995 C2
= Matrix_Copy(C
);
3996 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[0], Constraints2Polyhedron(C2
, MaxRays
));
3998 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
3999 value_set_si(C
->p
[0][1+i
], -1);
4000 value_set_si(C
->p
[0][1+dim
], -1);
4001 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[2], Constraints2Polyhedron(C
, MaxRays
));
4002 evalue_set_si(&split
.x
.p
->arr
[3], 1, 1);
4004 free_evalue_refs(&split
);
4010 static evalue
*find_fractional_with_max_periods(evalue
*e
, Polyhedron
*D
,
4013 Value
*min
, Value
*max
)
4020 if (value_notzero_p(e
->d
))
4023 if (e
->x
.p
->type
== fractional
) {
4028 I
= polynomial_projection(e
->x
.p
, D
, &d
, &T
);
4029 bounded
= line_minmax(I
, min
, max
); /* frees I */
4033 value_set_si(mp
, max_periods
);
4034 mpz_fdiv_q(*min
, *min
, d
);
4035 mpz_fdiv_q(*max
, *max
, d
);
4036 value_assign(T
->p
[1][D
->Dimension
], d
);
4037 value_subtract(d
, *max
, *min
);
4038 if (value_ge(d
, mp
))
4041 f
= evalue_dup(&e
->x
.p
->arr
[0]);
4052 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
4053 if ((f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[i
], D
, max_periods
,
4060 static void replace_fract_by_affine(evalue
*e
, evalue
*f
, Value val
)
4064 if (value_notzero_p(e
->d
))
4067 offset
= type_offset(e
->x
.p
);
4068 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
4069 replace_fract_by_affine(&e
->x
.p
->arr
[i
], f
, val
);
4071 if (e
->x
.p
->type
!= fractional
)
4074 if (!eequal(&e
->x
.p
->arr
[0], f
))
4077 replace_by_affine(e
, val
);
4080 /* Look for fractional parts that can be removed by splitting the corresponding
4081 * domain into at most max_periods parts.
4082 * We use a very simply strategy that looks for the first fractional part
4083 * that satisfies the condition, performs the split and then continues
4084 * looking for other fractional parts in the split domains until no
4085 * such fractional part can be found anymore.
4087 void evalue_split_periods(evalue
*e
, int max_periods
, unsigned int MaxRays
)
4094 if (EVALUE_IS_ZERO(*e
))
4096 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
4098 "WARNING: evalue_split_periods called on incorrect evalue type\n");
4106 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
4111 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
4113 f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[2*i
+1], D
, max_periods
,
4118 M
= Matrix_Alloc(2, 2+D
->Dimension
);
4120 value_subtract(d
, max
, min
);
4121 n
= VALUE_TO_INT(d
)+1;
4123 value_set_si(M
->p
[0][0], 1);
4124 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
4125 value_multiply(d
, max
, T
->p
[1][D
->Dimension
]);
4126 value_subtract(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
], d
);
4127 value_set_si(d
, -1);
4128 value_set_si(M
->p
[1][0], 1);
4129 Vector_Scale(T
->p
[0], M
->p
[1]+1, d
, D
->Dimension
+1);
4130 value_addmul(M
->p
[1][1+D
->Dimension
], max
, T
->p
[1][D
->Dimension
]);
4131 value_addto(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4132 T
->p
[1][D
->Dimension
]);
4133 value_decrement(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
]);
4135 p
= new_enode(partition
, e
->x
.p
->size
+ (n
-1)*2, e
->x
.p
->pos
);
4136 for (j
= 0; j
< 2*i
; ++j
) {
4137 value_clear(p
->arr
[j
].d
);
4138 p
->arr
[j
] = e
->x
.p
->arr
[j
];
4140 for (j
= 2*i
+2; j
< e
->x
.p
->size
; ++j
) {
4141 value_clear(p
->arr
[j
+2*(n
-1)].d
);
4142 p
->arr
[j
+2*(n
-1)] = e
->x
.p
->arr
[j
];
4144 for (j
= n
-1; j
>= 0; --j
) {
4146 value_clear(p
->arr
[2*i
+1].d
);
4147 p
->arr
[2*i
+1] = e
->x
.p
->arr
[2*i
+1];
4149 evalue_copy(&p
->arr
[2*(i
+j
)+1], &e
->x
.p
->arr
[2*i
+1]);
4151 value_subtract(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4152 T
->p
[1][D
->Dimension
]);
4153 value_addto(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
],
4154 T
->p
[1][D
->Dimension
]);
4156 replace_fract_by_affine(&p
->arr
[2*(i
+j
)+1], f
, max
);
4157 E
= DomainAddConstraints(D
, M
, MaxRays
);
4158 EVALUE_SET_DOMAIN(p
->arr
[2*(i
+j
)], E
);
4159 if (evalue_range_reduction_in_domain(&p
->arr
[2*(i
+j
)+1], E
))
4160 reduce_evalue(&p
->arr
[2*(i
+j
)+1]);
4161 value_decrement(max
, max
);
4163 value_clear(e
->x
.p
->arr
[2*i
].d
);
4178 void evalue_extract_affine(const evalue
*e
, Value
*coeff
, Value
*cst
, Value
*d
)
4180 value_set_si(*d
, 1);
4182 for ( ; value_zero_p(e
->d
); e
= &e
->x
.p
->arr
[0]) {
4184 assert(e
->x
.p
->type
== polynomial
);
4185 assert(e
->x
.p
->size
== 2);
4186 c
= &e
->x
.p
->arr
[1];
4187 value_multiply(coeff
[e
->x
.p
->pos
-1], *d
, c
->x
.n
);
4188 value_division(coeff
[e
->x
.p
->pos
-1], coeff
[e
->x
.p
->pos
-1], c
->d
);
4190 value_multiply(*cst
, *d
, e
->x
.n
);
4191 value_division(*cst
, *cst
, e
->d
);
4194 /* returns an evalue that corresponds to
4198 static evalue
*term(int param
, Value c
, Value den
)
4200 evalue
*EP
= ALLOC(evalue
);
4202 value_set_si(EP
->d
,0);
4203 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
4204 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
4205 value_init(EP
->x
.p
->arr
[1].x
.n
);
4206 value_assign(EP
->x
.p
->arr
[1].d
, den
);
4207 value_assign(EP
->x
.p
->arr
[1].x
.n
, c
);
4211 evalue
*affine2evalue(Value
*coeff
, Value denom
, int nvar
)
4214 evalue
*E
= ALLOC(evalue
);
4216 evalue_set(E
, coeff
[nvar
], denom
);
4217 for (i
= 0; i
< nvar
; ++i
) {
4219 if (value_zero_p(coeff
[i
]))
4221 t
= term(i
, coeff
[i
], denom
);
4228 void evalue_substitute(evalue
*e
, evalue
**subs
)
4234 if (value_notzero_p(e
->d
))
4238 assert(p
->type
!= partition
);
4240 for (i
= 0; i
< p
->size
; ++i
)
4241 evalue_substitute(&p
->arr
[i
], subs
);
4243 if (p
->type
== polynomial
)
4248 value_set_si(v
->d
, 0);
4249 v
->x
.p
= new_enode(p
->type
, 3, -1);
4250 value_clear(v
->x
.p
->arr
[0].d
);
4251 v
->x
.p
->arr
[0] = p
->arr
[0];
4252 evalue_set_si(&v
->x
.p
->arr
[1], 0, 1);
4253 evalue_set_si(&v
->x
.p
->arr
[2], 1, 1);
4256 offset
= type_offset(p
);
4258 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
4259 emul(v
, &p
->arr
[i
]);
4260 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
4261 free_evalue_refs(&(p
->arr
[i
]));
4264 if (p
->type
!= polynomial
)
4268 *e
= p
->arr
[offset
];
4272 /* evalue e is given in terms of "new" parameter; CP maps the new
4273 * parameters back to the old parameters.
4274 * Transforms e such that it refers back to the old parameters.
4276 void evalue_backsubstitute(evalue
*e
, Matrix
*CP
, unsigned MaxRays
)
4283 unsigned nparam
= CP
->NbColumns
-1;
4286 if (EVALUE_IS_ZERO(*e
))
4289 assert(value_zero_p(e
->d
));
4291 assert(p
->type
== partition
);
4293 inv
= left_inverse(CP
, &eq
);
4294 subs
= ALLOCN(evalue
*, nparam
);
4295 for (i
= 0; i
< nparam
; ++i
)
4296 subs
[i
] = affine2evalue(inv
->p
[i
], inv
->p
[nparam
][inv
->NbColumns
-1],
4299 CEq
= Constraints2Polyhedron(eq
, MaxRays
);
4300 addeliminatedparams_partition(p
, inv
, CEq
, inv
->NbColumns
-1, MaxRays
);
4301 Polyhedron_Free(CEq
);
4303 for (i
= 0; i
< p
->size
/2; ++i
)
4304 evalue_substitute(&p
->arr
[2*i
+1], subs
);
4306 for (i
= 0; i
< nparam
; ++i
)
4307 evalue_free(subs
[i
]);
4315 * \sum_{i=0}^n c_i/d X^i
4317 * where d is the last element in the vector c.
4319 evalue
*evalue_polynomial(Vector
*c
, const evalue
* X
)
4321 unsigned dim
= c
->Size
-2;
4323 evalue
*EP
= ALLOC(evalue
);
4328 if (EVALUE_IS_ZERO(*X
) || dim
== 0) {
4329 evalue_set(EP
, c
->p
[0], c
->p
[dim
+1]);
4330 reduce_constant(EP
);
4334 evalue_set(EP
, c
->p
[dim
], c
->p
[dim
+1]);
4337 evalue_set(&EC
, c
->p
[dim
], c
->p
[dim
+1]);
4339 for (i
= dim
-1; i
>= 0; --i
) {
4341 value_assign(EC
.x
.n
, c
->p
[i
]);
4344 free_evalue_refs(&EC
);
4348 /* Create an evalue from an array of pairs of domains and evalues. */
4349 evalue
*evalue_from_section_array(struct evalue_section
*s
, int n
)
4354 res
= ALLOC(evalue
);
4358 evalue_set_si(res
, 0, 1);
4360 value_set_si(res
->d
, 0);
4361 res
->x
.p
= new_enode(partition
, 2*n
, s
[0].D
->Dimension
);
4362 for (i
= 0; i
< n
; ++i
) {
4363 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
], s
[i
].D
);
4364 value_clear(res
->x
.p
->arr
[2*i
+1].d
);
4365 res
->x
.p
->arr
[2*i
+1] = *s
[i
].E
;