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);
1174 fd
= DomainConstraintSimplify(fd
, 0);
1179 /* See if we can extend one of the domains in res to cover fd */
1180 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1181 if (is_zero_on(&res
->x
.p
->arr
[2*i
+1], fd
))
1183 if (i
< res
->x
.p
->size
/2) {
1184 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
1185 DomainConcat(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
])));
1188 value_init(s
[n
].E
.d
);
1189 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
1193 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1194 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
1195 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1197 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1198 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1199 d
= DomainConstraintSimplify(d
, 0);
1205 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
1206 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1208 value_init(s
[n
].E
.d
);
1209 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1210 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1211 fd
= DomainConstraintSimplify(fd
, 0);
1212 if (!emptyQ(fd
) && is_zero_on(&e1
->x
.p
->arr
[2*j
+1], fd
)) {
1213 d
= DomainConcat(fd
, d
);
1214 fd
= Empty_Polyhedron(fd
->Dimension
);
1220 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
1224 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1227 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1228 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1229 value_clear(res
->x
.p
->arr
[2*i
].d
);
1234 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1235 for (j
= 0; j
< n
; ++j
) {
1236 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1237 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1238 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1244 static void explicit_complement(evalue
*res
)
1246 enode
*rel
= new_enode(relation
, 3, 0);
1248 value_clear(rel
->arr
[0].d
);
1249 rel
->arr
[0] = res
->x
.p
->arr
[0];
1250 value_clear(rel
->arr
[1].d
);
1251 rel
->arr
[1] = res
->x
.p
->arr
[1];
1252 value_set_si(rel
->arr
[2].d
, 1);
1253 value_init(rel
->arr
[2].x
.n
);
1254 value_set_si(rel
->arr
[2].x
.n
, 0);
1259 static void reduce_constant(evalue
*e
)
1264 value_gcd(g
, e
->x
.n
, e
->d
);
1265 if (value_notone_p(g
)) {
1266 value_division(e
->d
, e
->d
,g
);
1267 value_division(e
->x
.n
, e
->x
.n
,g
);
1272 void eadd(const evalue
*e1
, evalue
*res
)
1276 if (EVALUE_IS_ZERO(*e1
))
1279 if (EVALUE_IS_ZERO(*res
)) {
1280 if (value_notzero_p(e1
->d
)) {
1281 value_assign(res
->d
, e1
->d
);
1282 value_assign(res
->x
.n
, e1
->x
.n
);
1284 value_clear(res
->x
.n
);
1285 value_set_si(res
->d
, 0);
1286 res
->x
.p
= ecopy(e1
->x
.p
);
1291 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1292 /* Add two rational numbers */
1293 if (value_eq(e1
->d
, res
->d
))
1294 value_addto(res
->x
.n
, res
->x
.n
, e1
->x
.n
);
1296 value_multiply(res
->x
.n
, res
->x
.n
, e1
->d
);
1297 value_addmul(res
->x
.n
, e1
->x
.n
, res
->d
);
1298 value_multiply(res
->d
,e1
->d
,res
->d
);
1300 reduce_constant(res
);
1303 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
1304 switch (res
->x
.p
->type
) {
1306 /* Add the constant to the constant term of a polynomial*/
1307 eadd(e1
, &res
->x
.p
->arr
[0]);
1310 /* Add the constant to all elements of a periodic number */
1311 for (i
=0; i
<res
->x
.p
->size
; i
++) {
1312 eadd(e1
, &res
->x
.p
->arr
[i
]);
1316 fprintf(stderr
, "eadd: cannot add const with vector\n");
1320 eadd(e1
, &res
->x
.p
->arr
[1]);
1323 assert(EVALUE_IS_ZERO(*e1
));
1324 break; /* Do nothing */
1326 /* Create (zero) complement if needed */
1327 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
1328 explicit_complement(res
);
1329 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1330 eadd(e1
, &res
->x
.p
->arr
[i
]);
1336 /* add polynomial or periodic to constant
1337 * you have to exchange e1 and res, before doing addition */
1339 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1343 else { // ((e1->d==0) && (res->d==0))
1344 assert(!((e1
->x
.p
->type
== partition
) ^
1345 (res
->x
.p
->type
== partition
)));
1346 if (e1
->x
.p
->type
== partition
) {
1347 eadd_partitions(e1
, res
);
1350 if (e1
->x
.p
->type
== relation
&&
1351 (res
->x
.p
->type
!= relation
||
1352 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1356 if (res
->x
.p
->type
== relation
) {
1357 if (e1
->x
.p
->type
== relation
&&
1358 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1359 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1360 explicit_complement(res
);
1361 for (i
= 1; i
< e1
->x
.p
->size
; ++i
)
1362 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1365 if (res
->x
.p
->size
< 3)
1366 explicit_complement(res
);
1367 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1368 eadd(e1
, &res
->x
.p
->arr
[i
]);
1371 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
1372 /* adding to evalues of different type. two cases are possible
1373 * res is periodic and e1 is polynomial, you have to exchange
1374 * e1 and res then to add e1 to the constant term of res */
1375 if (e1
->x
.p
->type
== polynomial
) {
1376 eadd_rev_cst(e1
, res
);
1378 else if (res
->x
.p
->type
== polynomial
) {
1379 /* res is polynomial and e1 is periodic,
1380 add e1 to the constant term of res */
1382 eadd(e1
,&res
->x
.p
->arr
[0]);
1388 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
1389 ((res
->x
.p
->type
== fractional
||
1390 res
->x
.p
->type
== flooring
) &&
1391 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1392 /* adding evalues of different position (i.e function of different unknowns
1393 * to case are possible */
1395 switch (res
->x
.p
->type
) {
1398 if (mod_term_smaller(res
, e1
))
1399 eadd(e1
,&res
->x
.p
->arr
[1]);
1401 eadd_rev_cst(e1
, res
);
1403 case polynomial
: // res and e1 are polynomials
1404 // add e1 to the constant term of res
1406 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1407 eadd(e1
,&res
->x
.p
->arr
[0]);
1409 eadd_rev_cst(e1
, res
);
1410 // value_clear(g); value_clear(m1); value_clear(m2);
1412 case periodic
: // res and e1 are pointers to periodic numbers
1413 //add e1 to all elements of res
1415 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1416 for (i
=0;i
<res
->x
.p
->size
;i
++) {
1417 eadd(e1
,&res
->x
.p
->arr
[i
]);
1428 //same type , same pos and same size
1429 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1430 // add any element in e1 to the corresponding element in res
1431 i
= type_offset(res
->x
.p
);
1433 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1434 for (; i
<res
->x
.p
->size
; i
++) {
1435 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1440 /* Sizes are different */
1441 switch(res
->x
.p
->type
) {
1445 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1446 /* new enode and add res to that new node. If you do not do */
1447 /* that, you lose the the upper weight part of e1 ! */
1449 if(e1
->x
.p
->size
> res
->x
.p
->size
)
1452 i
= type_offset(res
->x
.p
);
1454 assert(eequal(&e1
->x
.p
->arr
[0],
1455 &res
->x
.p
->arr
[0]));
1456 for (; i
<e1
->x
.p
->size
; i
++) {
1457 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1464 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1467 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1468 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1469 to add periodicaly elements of e1 to elements of ne, and finaly to
1474 value_init(ex
); value_init(ey
);value_init(ep
);
1477 value_set_si(ex
,e1
->x
.p
->size
);
1478 value_set_si(ey
,res
->x
.p
->size
);
1479 value_assign (ep
,*Lcm(ex
,ey
));
1480 p
=(int)mpz_get_si(ep
);
1481 ne
= (evalue
*) malloc (sizeof(evalue
));
1483 value_set_si( ne
->d
,0);
1485 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1487 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1490 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1493 value_assign(res
->d
, ne
->d
);
1499 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
1508 static void emul_rev(const evalue
*e1
, evalue
*res
)
1512 evalue_copy(&ev
, e1
);
1514 free_evalue_refs(res
);
1518 static void emul_poly(const evalue
*e1
, evalue
*res
)
1520 int i
, j
, offset
= type_offset(res
->x
.p
);
1523 int size
= (e1
->x
.p
->size
+ res
->x
.p
->size
- offset
- 1);
1525 p
= new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1527 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1528 if (!EVALUE_IS_ZERO(e1
->x
.p
->arr
[i
]))
1531 /* special case pure power */
1532 if (i
== e1
->x
.p
->size
-1) {
1534 value_clear(p
->arr
[0].d
);
1535 p
->arr
[0] = res
->x
.p
->arr
[0];
1537 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1538 evalue_set_si(&p
->arr
[i
], 0, 1);
1539 for (i
= offset
; i
< res
->x
.p
->size
; ++i
) {
1540 value_clear(p
->arr
[i
+e1
->x
.p
->size
-offset
-1].d
);
1541 p
->arr
[i
+e1
->x
.p
->size
-offset
-1] = res
->x
.p
->arr
[i
];
1542 emul(&e1
->x
.p
->arr
[e1
->x
.p
->size
-1],
1543 &p
->arr
[i
+e1
->x
.p
->size
-offset
-1]);
1551 value_set_si(tmp
.d
,0);
1554 evalue_copy(&p
->arr
[0], &e1
->x
.p
->arr
[0]);
1555 for (i
= offset
; i
< e1
->x
.p
->size
; i
++) {
1556 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1557 emul(&res
->x
.p
->arr
[offset
], &tmp
.x
.p
->arr
[i
]);
1560 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1561 for (i
= offset
+1; i
<res
->x
.p
->size
; i
++)
1562 for (j
= offset
; j
<e1
->x
.p
->size
; j
++) {
1565 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1566 emul(&res
->x
.p
->arr
[i
], &ev
);
1567 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-offset
]);
1568 free_evalue_refs(&ev
);
1570 free_evalue_refs(res
);
1574 void emul_partitions(const evalue
*e1
, evalue
*res
)
1579 s
= (struct section
*)
1580 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1581 sizeof(struct section
));
1583 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1584 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1585 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1588 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1589 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1590 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1591 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1592 d
= DomainConstraintSimplify(d
, 0);
1598 /* This code is only needed because the partitions
1599 are not true partitions.
1601 for (k
= 0; k
< n
; ++k
) {
1602 if (DomainIncludes(s
[k
].D
, d
))
1604 if (DomainIncludes(d
, s
[k
].D
)) {
1605 Domain_Free(s
[k
].D
);
1606 free_evalue_refs(&s
[k
].E
);
1617 value_init(s
[n
].E
.d
);
1618 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1619 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1623 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1624 value_clear(res
->x
.p
->arr
[2*i
].d
);
1625 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1630 evalue_set_si(res
, 0, 1);
1632 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1633 for (j
= 0; j
< n
; ++j
) {
1634 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1635 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1636 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1643 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1645 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1646 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1647 * evalues is not treated here */
1649 void emul(const evalue
*e1
, evalue
*res
)
1653 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1654 fprintf(stderr
, "emul: do not proced on evector type !\n");
1658 if (EVALUE_IS_ZERO(*res
))
1661 if (EVALUE_IS_ONE(*e1
))
1664 if (EVALUE_IS_ZERO(*e1
)) {
1665 if (value_notzero_p(res
->d
)) {
1666 value_assign(res
->d
, e1
->d
);
1667 value_assign(res
->x
.n
, e1
->x
.n
);
1669 free_evalue_refs(res
);
1671 evalue_set_si(res
, 0, 1);
1676 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1677 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1678 emul_partitions(e1
, res
);
1681 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1682 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1683 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1685 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1686 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1687 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1688 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3) {
1689 free_evalue_refs(&res
->x
.p
->arr
[2]);
1692 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1693 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1696 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1697 emul(e1
, &res
->x
.p
->arr
[i
]);
1699 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1700 switch(e1
->x
.p
->type
) {
1702 switch(res
->x
.p
->type
) {
1704 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1705 /* Product of two polynomials of the same variable */
1710 /* Product of two polynomials of different variables */
1712 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1713 for( i
=0; i
<res
->x
.p
->size
; i
++)
1714 emul(e1
, &res
->x
.p
->arr
[i
]);
1723 /* Product of a polynomial and a periodic or fractional */
1730 switch(res
->x
.p
->type
) {
1732 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1733 /* Product of two periodics of the same parameter and period */
1735 for(i
=0; i
<res
->x
.p
->size
;i
++)
1736 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1741 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1742 /* Product of two periodics of the same parameter and different periods */
1746 value_init(x
); value_init(y
);value_init(z
);
1749 value_set_si(x
,e1
->x
.p
->size
);
1750 value_set_si(y
,res
->x
.p
->size
);
1751 value_assign (z
,*Lcm(x
,y
));
1752 lcm
=(int)mpz_get_si(z
);
1753 newp
= (evalue
*) malloc (sizeof(evalue
));
1754 value_init(newp
->d
);
1755 value_set_si( newp
->d
,0);
1756 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1757 for(i
=0;i
<lcm
;i
++) {
1758 evalue_copy(&newp
->x
.p
->arr
[i
],
1759 &res
->x
.p
->arr
[i
%iy
]);
1762 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1764 value_assign(res
->d
,newp
->d
);
1767 value_clear(x
); value_clear(y
);value_clear(z
);
1771 /* Product of two periodics of different parameters */
1773 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1774 for(i
=0; i
<res
->x
.p
->size
; i
++)
1775 emul(e1
, &(res
->x
.p
->arr
[i
]));
1783 /* Product of a periodic and a polynomial */
1785 for(i
=0; i
<res
->x
.p
->size
; i
++)
1786 emul(e1
, &(res
->x
.p
->arr
[i
]));
1793 switch(res
->x
.p
->type
) {
1795 for(i
=0; i
<res
->x
.p
->size
; i
++)
1796 emul(e1
, &(res
->x
.p
->arr
[i
]));
1803 assert(e1
->x
.p
->type
== res
->x
.p
->type
);
1804 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1805 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1808 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1809 if (e1
->x
.p
->type
!= fractional
|| !value_two_p(d
.d
))
1814 value_set_si(d
.x
.n
, 1);
1815 /* { x }^2 == { x }/2 */
1816 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1817 assert(e1
->x
.p
->size
== 3);
1818 assert(res
->x
.p
->size
== 3);
1820 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1822 eadd(&res
->x
.p
->arr
[1], &tmp
);
1823 emul(&e1
->x
.p
->arr
[2], &tmp
);
1824 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[1]);
1825 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[2]);
1826 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1827 free_evalue_refs(&tmp
);
1832 if(mod_term_smaller(res
, e1
))
1833 for(i
=1; i
<res
->x
.p
->size
; i
++)
1834 emul(e1
, &(res
->x
.p
->arr
[i
]));
1849 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1850 /* Product of two rational numbers */
1851 value_multiply(res
->d
,e1
->d
,res
->d
);
1852 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1853 reduce_constant(res
);
1857 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1858 /* Product of an expression (polynomial or peririodic) and a rational number */
1864 /* Product of a rationel number and an expression (polynomial or peririodic) */
1866 i
= type_offset(res
->x
.p
);
1867 for (; i
<res
->x
.p
->size
; i
++)
1868 emul(e1
, &res
->x
.p
->arr
[i
]);
1878 /* Frees mask content ! */
1879 void emask(evalue
*mask
, evalue
*res
) {
1886 if (EVALUE_IS_ZERO(*res
)) {
1887 free_evalue_refs(mask
);
1891 assert(value_zero_p(mask
->d
));
1892 assert(mask
->x
.p
->type
== partition
);
1893 assert(value_zero_p(res
->d
));
1894 assert(res
->x
.p
->type
== partition
);
1895 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1896 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1897 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1898 pos
= res
->x
.p
->pos
;
1900 s
= (struct section
*)
1901 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1902 sizeof(struct section
));
1906 evalue_set_si(&mone
, -1, 1);
1909 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1910 assert(mask
->x
.p
->size
>= 2);
1911 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1912 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1914 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1916 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1925 value_init(s
[n
].E
.d
);
1926 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1930 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1931 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1934 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1935 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1936 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1937 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1939 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1940 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1946 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1947 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1949 value_init(s
[n
].E
.d
);
1950 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1951 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1957 /* Just ignore; this may have been previously masked off */
1959 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1963 free_evalue_refs(&mone
);
1964 free_evalue_refs(mask
);
1965 free_evalue_refs(res
);
1968 evalue_set_si(res
, 0, 1);
1970 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1971 for (j
= 0; j
< n
; ++j
) {
1972 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1973 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1974 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1981 void evalue_copy(evalue
*dst
, const evalue
*src
)
1983 value_assign(dst
->d
, src
->d
);
1984 if(value_notzero_p(src
->d
)) {
1985 value_init(dst
->x
.n
);
1986 value_assign(dst
->x
.n
, src
->x
.n
);
1988 dst
->x
.p
= ecopy(src
->x
.p
);
1991 evalue
*evalue_dup(const evalue
*e
)
1993 evalue
*res
= ALLOC(evalue
);
1995 evalue_copy(res
, e
);
1999 enode
*new_enode(enode_type type
,int size
,int pos
) {
2005 fprintf(stderr
, "Allocating enode of size 0 !\n" );
2008 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
2012 for(i
=0; i
<size
; i
++) {
2013 value_init(res
->arr
[i
].d
);
2014 value_set_si(res
->arr
[i
].d
,0);
2015 res
->arr
[i
].x
.p
= 0;
2020 enode
*ecopy(enode
*e
) {
2025 res
= new_enode(e
->type
,e
->size
,e
->pos
);
2026 for(i
=0;i
<e
->size
;++i
) {
2027 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
2028 if(value_zero_p(res
->arr
[i
].d
))
2029 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
2030 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
2031 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
2033 value_init(res
->arr
[i
].x
.n
);
2034 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
2040 int ecmp(const evalue
*e1
, const evalue
*e2
)
2046 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
2050 value_multiply(m
, e1
->x
.n
, e2
->d
);
2051 value_multiply(m2
, e2
->x
.n
, e1
->d
);
2053 if (value_lt(m
, m2
))
2055 else if (value_gt(m
, m2
))
2065 if (value_notzero_p(e1
->d
))
2067 if (value_notzero_p(e2
->d
))
2073 if (p1
->type
!= p2
->type
)
2074 return p1
->type
- p2
->type
;
2075 if (p1
->pos
!= p2
->pos
)
2076 return p1
->pos
- p2
->pos
;
2077 if (p1
->size
!= p2
->size
)
2078 return p1
->size
- p2
->size
;
2080 for (i
= p1
->size
-1; i
>= 0; --i
)
2081 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
2087 int eequal(const evalue
*e1
, const evalue
*e2
)
2092 if (value_ne(e1
->d
,e2
->d
))
2095 /* e1->d == e2->d */
2096 if (value_notzero_p(e1
->d
)) {
2097 if (value_ne(e1
->x
.n
,e2
->x
.n
))
2100 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2104 /* e1->d == e2->d == 0 */
2107 if (p1
->type
!= p2
->type
) return 0;
2108 if (p1
->size
!= p2
->size
) return 0;
2109 if (p1
->pos
!= p2
->pos
) return 0;
2110 for (i
=0; i
<p1
->size
; i
++)
2111 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
2116 void free_evalue_refs(evalue
*e
) {
2121 if (EVALUE_IS_DOMAIN(*e
)) {
2122 Domain_Free(EVALUE_DOMAIN(*e
));
2125 } else if (value_pos_p(e
->d
)) {
2127 /* 'e' stores a constant */
2129 value_clear(e
->x
.n
);
2132 assert(value_zero_p(e
->d
));
2135 if (!p
) return; /* null pointer */
2136 for (i
=0; i
<p
->size
; i
++) {
2137 free_evalue_refs(&(p
->arr
[i
]));
2141 } /* free_evalue_refs */
2143 void evalue_free(evalue
*e
)
2145 free_evalue_refs(e
);
2149 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
2150 Vector
* val
, evalue
*res
)
2152 unsigned nparam
= periods
->Size
;
2155 double d
= compute_evalue(e
, val
->p
);
2156 d
*= VALUE_TO_DOUBLE(m
);
2161 value_assign(res
->d
, m
);
2162 value_init(res
->x
.n
);
2163 value_set_double(res
->x
.n
, d
);
2164 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
2167 if (value_one_p(periods
->p
[p
]))
2168 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
2173 value_assign(tmp
, periods
->p
[p
]);
2174 value_set_si(res
->d
, 0);
2175 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
2177 value_decrement(tmp
, tmp
);
2178 value_assign(val
->p
[p
], tmp
);
2179 mod2table_r(e
, periods
, m
, p
+1, val
,
2180 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
2181 } while (value_pos_p(tmp
));
2187 static void rel2table(evalue
*e
, int zero
)
2189 if (value_pos_p(e
->d
)) {
2190 if (value_zero_p(e
->x
.n
) == zero
)
2191 value_set_si(e
->x
.n
, 1);
2193 value_set_si(e
->x
.n
, 0);
2194 value_set_si(e
->d
, 1);
2197 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
2198 rel2table(&e
->x
.p
->arr
[i
], zero
);
2202 void evalue_mod2table(evalue
*e
, int nparam
)
2207 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2210 for (i
=0; i
<p
->size
; i
++) {
2211 evalue_mod2table(&(p
->arr
[i
]), nparam
);
2213 if (p
->type
== relation
) {
2218 evalue_copy(©
, &p
->arr
[0]);
2220 rel2table(&p
->arr
[0], 1);
2221 emul(&p
->arr
[0], &p
->arr
[1]);
2223 rel2table(©
, 0);
2224 emul(©
, &p
->arr
[2]);
2225 eadd(&p
->arr
[2], &p
->arr
[1]);
2226 free_evalue_refs(&p
->arr
[2]);
2227 free_evalue_refs(©
);
2229 free_evalue_refs(&p
->arr
[0]);
2233 } else if (p
->type
== fractional
) {
2234 Vector
*periods
= Vector_Alloc(nparam
);
2235 Vector
*val
= Vector_Alloc(nparam
);
2241 value_set_si(tmp
, 1);
2242 Vector_Set(periods
->p
, 1, nparam
);
2243 Vector_Set(val
->p
, 0, nparam
);
2244 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
2247 assert(p
->type
== polynomial
);
2248 assert(p
->size
== 2);
2249 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
2250 value_lcm(tmp
, tmp
, p
->arr
[1].d
);
2252 value_lcm(tmp
, tmp
, ev
->d
);
2254 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
2257 evalue_set_si(&res
, 0, 1);
2258 /* Compute the polynomial using Horner's rule */
2259 for (i
=p
->size
-1;i
>1;i
--) {
2260 eadd(&p
->arr
[i
], &res
);
2263 eadd(&p
->arr
[1], &res
);
2265 free_evalue_refs(e
);
2266 free_evalue_refs(&EP
);
2271 Vector_Free(periods
);
2273 } /* evalue_mod2table */
2275 /********************************************************/
2276 /* function in domain */
2277 /* check if the parameters in list_args */
2278 /* verifies the constraints of Domain P */
2279 /********************************************************/
2280 int in_domain(Polyhedron
*P
, Value
*list_args
)
2283 Value v
; /* value of the constraint of a row when
2284 parameters are instantiated*/
2288 for (row
= 0; row
< P
->NbConstraints
; row
++) {
2289 Inner_Product(P
->Constraint
[row
]+1, list_args
, P
->Dimension
, &v
);
2290 value_addto(v
, v
, P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
2291 if (value_neg_p(v
) ||
2292 value_zero_p(P
->Constraint
[row
][0]) && value_notzero_p(v
)) {
2299 return in
|| (P
->next
&& in_domain(P
->next
, list_args
));
2302 /****************************************************/
2303 /* function compute enode */
2304 /* compute the value of enode p with parameters */
2305 /* list "list_args */
2306 /* compute the polynomial or the periodic */
2307 /****************************************************/
2309 static double compute_enode(enode
*p
, Value
*list_args
) {
2321 if (p
->type
== polynomial
) {
2323 value_assign(param
,list_args
[p
->pos
-1]);
2325 /* Compute the polynomial using Horner's rule */
2326 for (i
=p
->size
-1;i
>0;i
--) {
2327 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2328 res
*=VALUE_TO_DOUBLE(param
);
2330 res
+=compute_evalue(&p
->arr
[0],list_args
);
2332 else if (p
->type
== fractional
) {
2333 double d
= compute_evalue(&p
->arr
[0], list_args
);
2334 d
-= floor(d
+1e-10);
2336 /* Compute the polynomial using Horner's rule */
2337 for (i
=p
->size
-1;i
>1;i
--) {
2338 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2341 res
+=compute_evalue(&p
->arr
[1],list_args
);
2343 else if (p
->type
== flooring
) {
2344 double d
= compute_evalue(&p
->arr
[0], list_args
);
2347 /* Compute the polynomial using Horner's rule */
2348 for (i
=p
->size
-1;i
>1;i
--) {
2349 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2352 res
+=compute_evalue(&p
->arr
[1],list_args
);
2354 else if (p
->type
== periodic
) {
2355 value_assign(m
,list_args
[p
->pos
-1]);
2357 /* Choose the right element of the periodic */
2358 value_set_si(param
,p
->size
);
2359 value_pmodulus(m
,m
,param
);
2360 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2362 else if (p
->type
== relation
) {
2363 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2364 res
= compute_evalue(&p
->arr
[1], list_args
);
2365 else if (p
->size
> 2)
2366 res
= compute_evalue(&p
->arr
[2], list_args
);
2368 else if (p
->type
== partition
) {
2369 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2370 Value
*vals
= list_args
;
2373 for (i
= 0; i
< dim
; ++i
) {
2374 value_init(vals
[i
]);
2376 value_assign(vals
[i
], list_args
[i
]);
2379 for (i
= 0; i
< p
->size
/2; ++i
)
2380 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2381 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2385 for (i
= 0; i
< dim
; ++i
)
2386 value_clear(vals
[i
]);
2395 } /* compute_enode */
2397 /*************************************************/
2398 /* return the value of Ehrhart Polynomial */
2399 /* It returns a double, because since it is */
2400 /* a recursive function, some intermediate value */
2401 /* might not be integral */
2402 /*************************************************/
2404 double compute_evalue(const evalue
*e
, Value
*list_args
)
2408 if (value_notzero_p(e
->d
)) {
2409 if (value_notone_p(e
->d
))
2410 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2412 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2415 res
= compute_enode(e
->x
.p
,list_args
);
2417 } /* compute_evalue */
2420 /****************************************************/
2421 /* function compute_poly : */
2422 /* Check for the good validity domain */
2423 /* return the number of point in the Polyhedron */
2424 /* in allocated memory */
2425 /* Using the Ehrhart pseudo-polynomial */
2426 /****************************************************/
2427 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2430 /* double d; int i; */
2432 tmp
= (Value
*) malloc (sizeof(Value
));
2433 assert(tmp
!= NULL
);
2435 value_set_si(*tmp
,0);
2438 return(tmp
); /* no ehrhart polynomial */
2439 if(en
->ValidityDomain
) {
2440 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2441 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2446 return(tmp
); /* no Validity Domain */
2448 if(in_domain(en
->ValidityDomain
,list_args
)) {
2450 #ifdef EVAL_EHRHART_DEBUG
2451 Print_Domain(stdout
,en
->ValidityDomain
);
2452 print_evalue(stdout
,&en
->EP
);
2455 /* d = compute_evalue(&en->EP,list_args);
2457 printf("(double)%lf = %d\n", d, i ); */
2458 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2464 value_set_si(*tmp
,0);
2465 return(tmp
); /* no compatible domain with the arguments */
2466 } /* compute_poly */
2468 static evalue
*eval_polynomial(const enode
*p
, int offset
,
2469 evalue
*base
, Value
*values
)
2474 res
= evalue_zero();
2475 for (i
= p
->size
-1; i
> offset
; --i
) {
2476 c
= evalue_eval(&p
->arr
[i
], values
);
2481 c
= evalue_eval(&p
->arr
[offset
], values
);
2488 evalue
*evalue_eval(const evalue
*e
, Value
*values
)
2495 if (value_notzero_p(e
->d
)) {
2496 res
= ALLOC(evalue
);
2498 evalue_copy(res
, e
);
2501 switch (e
->x
.p
->type
) {
2503 value_init(param
.x
.n
);
2504 value_assign(param
.x
.n
, values
[e
->x
.p
->pos
-1]);
2505 value_init(param
.d
);
2506 value_set_si(param
.d
, 1);
2508 res
= eval_polynomial(e
->x
.p
, 0, ¶m
, values
);
2509 free_evalue_refs(¶m
);
2512 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2513 mpz_fdiv_r(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2515 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2516 evalue_free(param2
);
2519 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2520 mpz_fdiv_q(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2521 value_set_si(param2
->d
, 1);
2523 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2524 evalue_free(param2
);
2527 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2528 if (value_zero_p(param2
->x
.n
))
2529 res
= evalue_eval(&e
->x
.p
->arr
[1], values
);
2530 else if (e
->x
.p
->size
> 2)
2531 res
= evalue_eval(&e
->x
.p
->arr
[2], values
);
2533 res
= evalue_zero();
2534 evalue_free(param2
);
2537 assert(e
->x
.p
->pos
== EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
);
2538 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2539 if (in_domain(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), values
)) {
2540 res
= evalue_eval(&e
->x
.p
->arr
[2*i
+1], values
);
2544 res
= evalue_zero();
2552 size_t value_size(Value v
) {
2553 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2554 * sizeof(v
[0]._mp_d
[0]);
2557 size_t domain_size(Polyhedron
*D
)
2560 size_t s
= sizeof(*D
);
2562 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2563 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2564 s
+= value_size(D
->Constraint
[i
][j
]);
2567 for (i = 0; i < D->NbRays; ++i)
2568 for (j = 0; j < D->Dimension+2; ++j)
2569 s += value_size(D->Ray[i][j]);
2572 return D
->next
? s
+domain_size(D
->next
) : s
;
2575 size_t enode_size(enode
*p
) {
2576 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2579 if (p
->type
== partition
)
2580 for (i
= 0; i
< p
->size
/2; ++i
) {
2581 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2582 s
+= evalue_size(&p
->arr
[2*i
+1]);
2585 for (i
= 0; i
< p
->size
; ++i
) {
2586 s
+= evalue_size(&p
->arr
[i
]);
2591 size_t evalue_size(evalue
*e
)
2593 size_t s
= sizeof(*e
);
2594 s
+= value_size(e
->d
);
2595 if (value_notzero_p(e
->d
))
2596 s
+= value_size(e
->x
.n
);
2598 s
+= enode_size(e
->x
.p
);
2602 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2604 evalue
*found
= NULL
;
2609 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2612 value_init(offset
.d
);
2613 value_init(offset
.x
.n
);
2614 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2615 value_lcm(offset
.d
, m
, offset
.d
);
2616 value_set_si(offset
.x
.n
, 1);
2619 evalue_copy(©
, cst
);
2622 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2624 if (eequal(base
, &e
->x
.p
->arr
[0]))
2625 found
= &e
->x
.p
->arr
[0];
2627 value_set_si(offset
.x
.n
, -2);
2630 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2632 if (eequal(base
, &e
->x
.p
->arr
[0]))
2635 free_evalue_refs(cst
);
2636 free_evalue_refs(&offset
);
2639 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2640 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2645 static evalue
*find_relation_pair(evalue
*e
)
2648 evalue
*found
= NULL
;
2650 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2653 if (e
->x
.p
->type
== fractional
) {
2658 poly_denom(&e
->x
.p
->arr
[0], &m
);
2660 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2661 cst
= &cst
->x
.p
->arr
[0])
2664 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2665 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2670 i
= e
->x
.p
->type
== relation
;
2671 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2672 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2677 void evalue_mod2relation(evalue
*e
) {
2680 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2683 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2684 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2685 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2686 value_clear(e
->x
.p
->arr
[2*i
].d
);
2687 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2689 if (2*i
< e
->x
.p
->size
) {
2690 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2691 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2696 if (e
->x
.p
->size
== 0) {
2698 evalue_set_si(e
, 0, 1);
2704 while ((d
= find_relation_pair(e
)) != NULL
) {
2708 value_init(split
.d
);
2709 value_set_si(split
.d
, 0);
2710 split
.x
.p
= new_enode(relation
, 3, 0);
2711 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2712 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2714 ev
= &split
.x
.p
->arr
[0];
2715 value_set_si(ev
->d
, 0);
2716 ev
->x
.p
= new_enode(fractional
, 3, -1);
2717 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2718 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2719 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2725 free_evalue_refs(&split
);
2729 static int evalue_comp(const void * a
, const void * b
)
2731 const evalue
*e1
= *(const evalue
**)a
;
2732 const evalue
*e2
= *(const evalue
**)b
;
2733 return ecmp(e1
, e2
);
2736 void evalue_combine(evalue
*e
)
2743 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2746 NALLOC(evs
, e
->x
.p
->size
/2);
2747 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2748 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2749 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2750 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2751 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2752 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2753 value_clear(p
->arr
[2*k
].d
);
2754 value_clear(p
->arr
[2*k
+1].d
);
2755 p
->arr
[2*k
] = *(evs
[i
]-1);
2756 p
->arr
[2*k
+1] = *(evs
[i
]);
2759 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2762 value_clear((evs
[i
]-1)->d
);
2766 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2767 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2768 free_evalue_refs(evs
[i
]);
2772 for (i
= 2*k
; i
< p
->size
; ++i
)
2773 value_clear(p
->arr
[i
].d
);
2780 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2782 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2784 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2787 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2788 Polyhedron
*D
, *N
, **P
;
2791 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2798 if (D
->NbEq
<= H
->NbEq
) {
2804 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2805 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2806 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2807 reduce_evalue(&tmp
);
2808 if (value_notzero_p(tmp
.d
) ||
2809 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2812 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2813 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2816 free_evalue_refs(&tmp
);
2822 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2824 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2826 value_clear(e
->x
.p
->arr
[2*i
].d
);
2827 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2829 if (2*i
< e
->x
.p
->size
) {
2830 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2831 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2838 H
= DomainConvex(D
, 0);
2839 E
= DomainDifference(H
, D
, 0);
2841 D
= DomainDifference(H
, E
, 0);
2844 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2848 /* Use smallest representative for coefficients in affine form in
2849 * argument of fractional.
2850 * Since any change will make the argument non-standard,
2851 * the containing evalue will have to be reduced again afterward.
2853 static void fractional_minimal_coefficients(enode
*p
)
2859 assert(p
->type
== fractional
);
2861 while (value_zero_p(pp
->d
)) {
2862 assert(pp
->x
.p
->type
== polynomial
);
2863 assert(pp
->x
.p
->size
== 2);
2864 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2865 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2866 if (value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2867 value_subtract(pp
->x
.p
->arr
[1].x
.n
,
2868 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2869 pp
= &pp
->x
.p
->arr
[0];
2875 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2880 unsigned dim
= D
->Dimension
;
2881 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2884 assert(p
->type
== fractional
|| p
->type
== flooring
);
2885 value_set_si(T
->p
[1][dim
], 1);
2886 evalue_extract_affine(&p
->arr
[0], T
->p
[0], &T
->p
[0][dim
], d
);
2887 I
= DomainImage(D
, T
, 0);
2888 H
= DomainConvex(I
, 0);
2898 static void replace_by_affine(evalue
*e
, Value offset
)
2905 value_init(inc
.x
.n
);
2906 value_set_si(inc
.d
, 1);
2907 value_oppose(inc
.x
.n
, offset
);
2908 eadd(&inc
, &p
->arr
[0]);
2909 reorder_terms_about(p
, &p
->arr
[0]); /* frees arr[0] */
2913 free_evalue_refs(&inc
);
2916 int evalue_range_reduction_in_domain(evalue
*e
, Polyhedron
*D
)
2925 if (value_notzero_p(e
->d
))
2930 if (p
->type
== relation
) {
2937 fractional_minimal_coefficients(p
->arr
[0].x
.p
);
2938 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
);
2939 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2940 equal
= value_eq(min
, max
);
2941 mpz_cdiv_q(min
, min
, d
);
2942 mpz_fdiv_q(max
, max
, d
);
2944 if (bounded
&& value_gt(min
, max
)) {
2950 evalue_set_si(e
, 0, 1);
2953 free_evalue_refs(&(p
->arr
[1]));
2954 free_evalue_refs(&(p
->arr
[0]));
2960 return r
? r
: evalue_range_reduction_in_domain(e
, D
);
2961 } else if (bounded
&& equal
) {
2964 free_evalue_refs(&(p
->arr
[2]));
2967 free_evalue_refs(&(p
->arr
[0]));
2973 return evalue_range_reduction_in_domain(e
, D
);
2974 } else if (bounded
&& value_eq(min
, max
)) {
2975 /* zero for a single value */
2977 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2978 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2979 value_multiply(min
, min
, d
);
2980 value_subtract(M
->p
[0][D
->Dimension
+1],
2981 M
->p
[0][D
->Dimension
+1], min
);
2982 E
= DomainAddConstraints(D
, M
, 0);
2988 r
= evalue_range_reduction_in_domain(&p
->arr
[1], E
);
2990 r
|= evalue_range_reduction_in_domain(&p
->arr
[2], D
);
2992 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
3000 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
3003 i
= p
->type
== relation
? 1 :
3004 p
->type
== fractional
? 1 : 0;
3005 for (; i
<p
->size
; i
++)
3006 r
|= evalue_range_reduction_in_domain(&p
->arr
[i
], D
);
3008 if (p
->type
!= fractional
) {
3009 if (r
&& p
->type
== polynomial
) {
3012 value_set_si(f
.d
, 0);
3013 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3014 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3015 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3016 reorder_terms_about(p
, &f
);
3027 fractional_minimal_coefficients(p
);
3028 I
= polynomial_projection(p
, D
, &d
, NULL
);
3029 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
3030 mpz_fdiv_q(min
, min
, d
);
3031 mpz_fdiv_q(max
, max
, d
);
3032 value_subtract(d
, max
, min
);
3034 if (bounded
&& value_eq(min
, max
)) {
3035 replace_by_affine(e
, min
);
3037 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
3038 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
3039 * See pages 199-200 of PhD thesis.
3047 value_set_si(rem
.d
, 0);
3048 rem
.x
.p
= new_enode(fractional
, 3, -1);
3049 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
3050 value_clear(rem
.x
.p
->arr
[1].d
);
3051 value_clear(rem
.x
.p
->arr
[2].d
);
3052 rem
.x
.p
->arr
[1] = p
->arr
[1];
3053 rem
.x
.p
->arr
[2] = p
->arr
[2];
3054 for (i
= 3; i
< p
->size
; ++i
)
3055 p
->arr
[i
-2] = p
->arr
[i
];
3059 value_init(inc
.x
.n
);
3060 value_set_si(inc
.d
, 1);
3061 value_oppose(inc
.x
.n
, min
);
3064 evalue_copy(&t
, &p
->arr
[0]);
3068 value_set_si(f
.d
, 0);
3069 f
.x
.p
= new_enode(fractional
, 3, -1);
3070 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3071 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3072 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
3074 value_init(factor
.d
);
3075 evalue_set_si(&factor
, -1, 1);
3081 value_clear(f
.x
.p
->arr
[1].x
.n
);
3082 value_clear(f
.x
.p
->arr
[2].x
.n
);
3083 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3084 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3088 reorder_terms(&rem
);
3095 free_evalue_refs(&inc
);
3096 free_evalue_refs(&t
);
3097 free_evalue_refs(&f
);
3098 free_evalue_refs(&factor
);
3099 free_evalue_refs(&rem
);
3101 evalue_range_reduction_in_domain(e
, D
);
3105 _reduce_evalue(&p
->arr
[0], 0, 1);
3117 void evalue_range_reduction(evalue
*e
)
3120 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3123 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3124 if (evalue_range_reduction_in_domain(&e
->x
.p
->arr
[2*i
+1],
3125 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
3126 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3128 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
3129 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
3130 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3131 value_clear(e
->x
.p
->arr
[2*i
].d
);
3133 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
3134 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
3142 Enumeration
* partition2enumeration(evalue
*EP
)
3145 Enumeration
*en
, *res
= NULL
;
3147 if (EVALUE_IS_ZERO(*EP
)) {
3152 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
3153 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
3154 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
3157 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
3158 value_clear(EP
->x
.p
->arr
[2*i
].d
);
3159 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
3167 int evalue_frac2floor_in_domain3(evalue
*e
, Polyhedron
*D
, int shift
)
3176 if (value_notzero_p(e
->d
))
3181 i
= p
->type
== relation
? 1 :
3182 p
->type
== fractional
? 1 : 0;
3183 for (; i
<p
->size
; i
++)
3184 r
|= evalue_frac2floor_in_domain3(&p
->arr
[i
], D
, shift
);
3186 if (p
->type
!= fractional
) {
3187 if (r
&& p
->type
== polynomial
) {
3190 value_set_si(f
.d
, 0);
3191 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3192 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3193 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3194 reorder_terms_about(p
, &f
);
3204 I
= polynomial_projection(p
, D
, &d
, NULL
);
3207 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3210 assert(I
->NbEq
== 0); /* Should have been reduced */
3213 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3214 if (value_pos_p(I
->Constraint
[i
][1]))
3217 if (i
< I
->NbConstraints
) {
3219 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3220 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3221 if (value_neg_p(min
)) {
3223 mpz_fdiv_q(min
, min
, d
);
3224 value_init(offset
.d
);
3225 value_set_si(offset
.d
, 1);
3226 value_init(offset
.x
.n
);
3227 value_oppose(offset
.x
.n
, min
);
3228 eadd(&offset
, &p
->arr
[0]);
3229 free_evalue_refs(&offset
);
3239 value_set_si(fl
.d
, 0);
3240 fl
.x
.p
= new_enode(flooring
, 3, -1);
3241 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
3242 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
3243 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
3245 eadd(&fl
, &p
->arr
[0]);
3246 reorder_terms_about(p
, &p
->arr
[0]);
3250 free_evalue_refs(&fl
);
3255 int evalue_frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
3257 return evalue_frac2floor_in_domain3(e
, D
, 1);
3260 void evalue_frac2floor2(evalue
*e
, int shift
)
3263 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
3265 if (evalue_frac2floor_in_domain3(e
, NULL
, 0))
3271 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3272 if (evalue_frac2floor_in_domain3(&e
->x
.p
->arr
[2*i
+1],
3273 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), shift
))
3274 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3277 void evalue_frac2floor(evalue
*e
)
3279 evalue_frac2floor2(e
, 1);
3282 /* Add a new variable with lower bound 1 and upper bound specified
3283 * by row. If negative is true, then the new variable has upper
3284 * bound -1 and lower bound specified by row.
3286 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
3287 Vector
*row
, int negative
)
3291 int nparam
= D
->Dimension
- nvar
;
3294 nr
= D
->NbConstraints
+ 2;
3295 nc
= D
->Dimension
+ 2 + 1;
3296 C
= Matrix_Alloc(nr
, nc
);
3297 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
3298 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
3299 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3300 D
->Dimension
+ 1 - nvar
);
3305 nc
= C
->NbColumns
+ 1;
3306 C
= Matrix_Alloc(nr
, nc
);
3307 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
3308 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
3309 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3310 oldC
->NbColumns
- 1 - nvar
);
3313 value_set_si(C
->p
[nr
-2][0], 1);
3315 value_set_si(C
->p
[nr
-2][1 + nvar
], -1);
3317 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
3318 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
3320 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
3321 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
3327 static void floor2frac_r(evalue
*e
, int nvar
)
3334 if (value_notzero_p(e
->d
))
3339 assert(p
->type
== flooring
);
3340 for (i
= 1; i
< p
->size
; i
++)
3341 floor2frac_r(&p
->arr
[i
], nvar
);
3343 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
3344 assert(pp
->x
.p
->type
== polynomial
);
3345 pp
->x
.p
->pos
-= nvar
;
3349 value_set_si(f
.d
, 0);
3350 f
.x
.p
= new_enode(fractional
, 3, -1);
3351 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3352 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3353 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3355 eadd(&f
, &p
->arr
[0]);
3356 reorder_terms_about(p
, &p
->arr
[0]);
3360 free_evalue_refs(&f
);
3363 /* Convert flooring back to fractional and shift position
3364 * of the parameters by nvar
3366 static void floor2frac(evalue
*e
, int nvar
)
3368 floor2frac_r(e
, nvar
);
3372 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
3375 int nparam
= D
->Dimension
- nvar
;
3379 D
= Constraints2Polyhedron(C
, 0);
3383 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
3385 /* Double check that D was not unbounded. */
3386 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
3394 static evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
3395 int *signs
, Matrix
*C
, unsigned MaxRays
)
3401 evalue
*factor
= NULL
;
3405 if (EVALUE_IS_ZERO(*e
))
3409 Polyhedron
*DD
= Disjoint_Domain(D
, 0, MaxRays
);
3416 res
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3419 for (Q
= DD
; Q
; Q
= DD
) {
3425 t
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3438 if (value_notzero_p(e
->d
)) {
3441 t
= esum_over_domain_cst(nvar
, D
, C
);
3443 if (!EVALUE_IS_ONE(*e
))
3449 switch (e
->x
.p
->type
) {
3451 evalue
*pp
= &e
->x
.p
->arr
[0];
3453 if (pp
->x
.p
->pos
> nvar
) {
3454 /* remainder is independent of the summated vars */
3460 floor2frac(&f
, nvar
);
3462 t
= esum_over_domain_cst(nvar
, D
, C
);
3466 free_evalue_refs(&f
);
3471 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3472 poly_denom(pp
, &row
->p
[1 + nvar
]);
3473 value_set_si(row
->p
[0], 1);
3474 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3475 pp
= &pp
->x
.p
->arr
[0]) {
3477 assert(pp
->x
.p
->type
== polynomial
);
3479 if (pos
>= 1 + nvar
)
3481 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3482 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3483 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3485 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3486 value_division(row
->p
[1 + D
->Dimension
+ 1],
3487 row
->p
[1 + D
->Dimension
+ 1],
3489 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3490 row
->p
[1 + D
->Dimension
+ 1],
3492 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3496 int pos
= e
->x
.p
->pos
;
3499 factor
= ALLOC(evalue
);
3500 value_init(factor
->d
);
3501 value_set_si(factor
->d
, 0);
3502 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3503 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3504 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3508 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3509 negative
= signs
[pos
-1] < 0;
3510 value_set_si(row
->p
[0], 1);
3512 value_set_si(row
->p
[pos
], -1);
3513 value_set_si(row
->p
[1 + nvar
], 1);
3515 value_set_si(row
->p
[pos
], 1);
3516 value_set_si(row
->p
[1 + nvar
], -1);
3524 offset
= type_offset(e
->x
.p
);
3526 res
= esum_over_domain(&e
->x
.p
->arr
[offset
], nvar
, D
, signs
, C
, MaxRays
);
3530 evalue_copy(&cum
, factor
);
3534 for (i
= 1; offset
+i
< e
->x
.p
->size
; ++i
) {
3538 C
= esum_add_constraint(nvar
, D
, C
, row
, negative
);
3544 Vector_Print(stderr, P_VALUE_FMT, row);
3546 Matrix_Print(stderr, P_VALUE_FMT, C);
3548 t
= esum_over_domain(&e
->x
.p
->arr
[offset
+i
], nvar
, D
, signs
, C
, MaxRays
);
3553 if (negative
&& (i
% 2))
3563 if (factor
&& offset
+i
+1 < e
->x
.p
->size
)
3570 free_evalue_refs(&cum
);
3571 evalue_free(factor
);
3582 static void domain_signs(Polyhedron
*D
, int *signs
)
3586 POL_ENSURE_VERTICES(D
);
3587 for (j
= 0; j
< D
->Dimension
; ++j
) {
3589 for (k
= 0; k
< D
->NbRays
; ++k
) {
3590 signs
[j
] = value_sign(D
->Ray
[k
][1+j
]);
3597 static void shift_floor_in_domain(evalue
*e
, Polyhedron
*D
)
3604 if (value_notzero_p(e
->d
))
3609 for (i
= type_offset(p
); i
< p
->size
; ++i
)
3610 shift_floor_in_domain(&p
->arr
[i
], D
);
3612 if (p
->type
!= flooring
)
3618 I
= polynomial_projection(p
, D
, &d
, NULL
);
3619 assert(I
->NbEq
== 0); /* Should have been reduced */
3621 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3622 if (value_pos_p(I
->Constraint
[i
][1]))
3624 assert(i
< I
->NbConstraints
);
3625 if (i
< I
->NbConstraints
) {
3626 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3627 mpz_fdiv_q(m
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3628 if (value_neg_p(m
)) {
3629 /* replace [e] by [e-m]+m such that e-m >= 0 */
3634 value_set_si(f
.d
, 1);
3635 value_oppose(f
.x
.n
, m
);
3636 eadd(&f
, &p
->arr
[0]);
3639 value_set_si(f
.d
, 0);
3640 f
.x
.p
= new_enode(flooring
, 3, -1);
3641 value_clear(f
.x
.p
->arr
[0].d
);
3642 f
.x
.p
->arr
[0] = p
->arr
[0];
3643 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
3644 value_set_si(f
.x
.p
->arr
[1].d
, 1);
3645 value_init(f
.x
.p
->arr
[1].x
.n
);
3646 value_assign(f
.x
.p
->arr
[1].x
.n
, m
);
3647 reorder_terms_about(p
, &f
);
3658 /* Make arguments of all floors non-negative */
3659 static void shift_floor_arguments(evalue
*e
)
3663 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3666 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3667 shift_floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
3668 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3671 evalue
*evalue_sum(evalue
*e
, int nvar
, unsigned MaxRays
)
3675 evalue
*res
= ALLOC(evalue
);
3679 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3680 evalue_copy(res
, e
);
3684 evalue_split_domains_into_orthants(e
, MaxRays
);
3685 evalue_frac2floor2(e
, 0);
3686 evalue_set_si(res
, 0, 1);
3688 assert(value_zero_p(e
->d
));
3689 assert(e
->x
.p
->type
== partition
);
3690 shift_floor_arguments(e
);
3692 assert(e
->x
.p
->size
>= 2);
3693 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3695 signs
= alloca(sizeof(int) * dim
);
3697 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3699 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3700 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3701 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
, 0,
3712 evalue
*esum(evalue
*e
, int nvar
)
3714 return evalue_sum(e
, nvar
, 0);
3717 /* Initial silly implementation */
3718 void eor(evalue
*e1
, evalue
*res
)
3724 evalue_set_si(&mone
, -1, 1);
3726 evalue_copy(&E
, res
);
3732 free_evalue_refs(&E
);
3733 free_evalue_refs(&mone
);
3736 /* computes denominator of polynomial evalue
3737 * d should point to a value initialized to 1
3739 void evalue_denom(const evalue
*e
, Value
*d
)
3743 if (value_notzero_p(e
->d
)) {
3744 value_lcm(*d
, *d
, e
->d
);
3747 assert(e
->x
.p
->type
== polynomial
||
3748 e
->x
.p
->type
== fractional
||
3749 e
->x
.p
->type
== flooring
);
3750 offset
= type_offset(e
->x
.p
);
3751 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3752 evalue_denom(&e
->x
.p
->arr
[i
], d
);
3755 /* Divides the evalue e by the integer n */
3756 void evalue_div(evalue
*e
, Value n
)
3760 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3763 if (value_notzero_p(e
->d
)) {
3766 value_multiply(e
->d
, e
->d
, n
);
3767 value_gcd(gc
, e
->x
.n
, e
->d
);
3768 if (value_notone_p(gc
)) {
3769 value_division(e
->d
, e
->d
, gc
);
3770 value_division(e
->x
.n
, e
->x
.n
, gc
);
3775 if (e
->x
.p
->type
== partition
) {
3776 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3777 evalue_div(&e
->x
.p
->arr
[2*i
+1], n
);
3780 offset
= type_offset(e
->x
.p
);
3781 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3782 evalue_div(&e
->x
.p
->arr
[i
], n
);
3785 /* Multiplies the evalue e by the integer n */
3786 void evalue_mul(evalue
*e
, Value n
)
3790 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3793 if (value_notzero_p(e
->d
)) {
3796 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3797 value_gcd(gc
, e
->x
.n
, e
->d
);
3798 if (value_notone_p(gc
)) {
3799 value_division(e
->d
, e
->d
, gc
);
3800 value_division(e
->x
.n
, e
->x
.n
, gc
);
3805 if (e
->x
.p
->type
== partition
) {
3806 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3807 evalue_mul(&e
->x
.p
->arr
[2*i
+1], n
);
3810 offset
= type_offset(e
->x
.p
);
3811 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3812 evalue_mul(&e
->x
.p
->arr
[i
], n
);
3815 /* Multiplies the evalue e by the n/d */
3816 void evalue_mul_div(evalue
*e
, Value n
, Value d
)
3820 if ((value_one_p(n
) && value_one_p(d
)) || EVALUE_IS_ZERO(*e
))
3823 if (value_notzero_p(e
->d
)) {
3826 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3827 value_multiply(e
->d
, e
->d
, d
);
3828 value_gcd(gc
, e
->x
.n
, e
->d
);
3829 if (value_notone_p(gc
)) {
3830 value_division(e
->d
, e
->d
, gc
);
3831 value_division(e
->x
.n
, e
->x
.n
, gc
);
3836 if (e
->x
.p
->type
== partition
) {
3837 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3838 evalue_mul_div(&e
->x
.p
->arr
[2*i
+1], n
, d
);
3841 offset
= type_offset(e
->x
.p
);
3842 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3843 evalue_mul_div(&e
->x
.p
->arr
[i
], n
, d
);
3846 void evalue_negate(evalue
*e
)
3850 if (value_notzero_p(e
->d
)) {
3851 value_oppose(e
->x
.n
, e
->x
.n
);
3854 if (e
->x
.p
->type
== partition
) {
3855 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3856 evalue_negate(&e
->x
.p
->arr
[2*i
+1]);
3859 offset
= type_offset(e
->x
.p
);
3860 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3861 evalue_negate(&e
->x
.p
->arr
[i
]);
3864 void evalue_add_constant(evalue
*e
, const Value cst
)
3868 if (value_zero_p(e
->d
)) {
3869 if (e
->x
.p
->type
== partition
) {
3870 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3871 evalue_add_constant(&e
->x
.p
->arr
[2*i
+1], cst
);
3874 if (e
->x
.p
->type
== relation
) {
3875 for (i
= 1; i
< e
->x
.p
->size
; ++i
)
3876 evalue_add_constant(&e
->x
.p
->arr
[i
], cst
);
3880 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
3881 } while (value_zero_p(e
->d
));
3883 value_addmul(e
->x
.n
, cst
, e
->d
);
3886 static void evalue_frac2polynomial_r(evalue
*e
, int *signs
, int sign
, int in_frac
)
3891 int sign_odd
= sign
;
3893 if (value_notzero_p(e
->d
)) {
3894 if (in_frac
&& sign
* value_sign(e
->x
.n
) < 0) {
3895 value_set_si(e
->x
.n
, 0);
3896 value_set_si(e
->d
, 1);
3901 if (e
->x
.p
->type
== relation
) {
3902 for (i
= e
->x
.p
->size
-1; i
>= 1; --i
)
3903 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
, sign
, in_frac
);
3907 if (e
->x
.p
->type
== polynomial
)
3908 sign_odd
*= signs
[e
->x
.p
->pos
-1];
3909 offset
= type_offset(e
->x
.p
);
3910 evalue_frac2polynomial_r(&e
->x
.p
->arr
[offset
], signs
, sign
, in_frac
);
3911 in_frac
|= e
->x
.p
->type
== fractional
;
3912 for (i
= e
->x
.p
->size
-1; i
> offset
; --i
)
3913 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
,
3914 (i
- offset
) % 2 ? sign_odd
: sign
, in_frac
);
3916 if (e
->x
.p
->type
!= fractional
)
3919 /* replace { a/m } by (m-1)/m if sign != 0
3920 * and by (m-1)/(2m) if sign == 0
3924 evalue_denom(&e
->x
.p
->arr
[0], &d
);
3925 free_evalue_refs(&e
->x
.p
->arr
[0]);
3926 value_init(e
->x
.p
->arr
[0].d
);
3927 value_init(e
->x
.p
->arr
[0].x
.n
);
3929 value_addto(e
->x
.p
->arr
[0].d
, d
, d
);
3931 value_assign(e
->x
.p
->arr
[0].d
, d
);
3932 value_decrement(e
->x
.p
->arr
[0].x
.n
, d
);
3936 reorder_terms_about(p
, &p
->arr
[0]);
3942 /* Approximate the evalue in fractional representation by a polynomial.
3943 * If sign > 0, the result is an upper bound;
3944 * if sign < 0, the result is a lower bound;
3945 * if sign = 0, the result is an intermediate approximation.
3947 void evalue_frac2polynomial(evalue
*e
, int sign
, unsigned MaxRays
)
3952 if (value_notzero_p(e
->d
))
3954 assert(e
->x
.p
->type
== partition
);
3955 /* make sure all variables in the domains have a fixed sign */
3957 evalue_split_domains_into_orthants(e
, MaxRays
);
3958 if (EVALUE_IS_ZERO(*e
))
3962 assert(e
->x
.p
->size
>= 2);
3963 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3965 signs
= alloca(sizeof(int) * dim
);
3968 for (i
= 0; i
< dim
; ++i
)
3970 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3972 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3973 evalue_frac2polynomial_r(&e
->x
.p
->arr
[2*i
+1], signs
, sign
, 0);
3977 /* Split the domains of e (which is assumed to be a partition)
3978 * such that each resulting domain lies entirely in one orthant.
3980 void evalue_split_domains_into_orthants(evalue
*e
, unsigned MaxRays
)
3983 assert(value_zero_p(e
->d
));
3984 assert(e
->x
.p
->type
== partition
);
3985 assert(e
->x
.p
->size
>= 2);
3986 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3988 for (i
= 0; i
< dim
; ++i
) {
3991 C
= Matrix_Alloc(1, 1 + dim
+ 1);
3992 value_set_si(C
->p
[0][0], 1);
3993 value_init(split
.d
);
3994 value_set_si(split
.d
, 0);
3995 split
.x
.p
= new_enode(partition
, 4, dim
);
3996 value_set_si(C
->p
[0][1+i
], 1);
3997 C2
= Matrix_Copy(C
);
3998 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[0], Constraints2Polyhedron(C2
, MaxRays
));
4000 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
4001 value_set_si(C
->p
[0][1+i
], -1);
4002 value_set_si(C
->p
[0][1+dim
], -1);
4003 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[2], Constraints2Polyhedron(C
, MaxRays
));
4004 evalue_set_si(&split
.x
.p
->arr
[3], 1, 1);
4006 free_evalue_refs(&split
);
4012 static evalue
*find_fractional_with_max_periods(evalue
*e
, Polyhedron
*D
,
4015 Value
*min
, Value
*max
)
4022 if (value_notzero_p(e
->d
))
4025 if (e
->x
.p
->type
== fractional
) {
4030 I
= polynomial_projection(e
->x
.p
, D
, &d
, &T
);
4031 bounded
= line_minmax(I
, min
, max
); /* frees I */
4035 value_set_si(mp
, max_periods
);
4036 mpz_fdiv_q(*min
, *min
, d
);
4037 mpz_fdiv_q(*max
, *max
, d
);
4038 value_assign(T
->p
[1][D
->Dimension
], d
);
4039 value_subtract(d
, *max
, *min
);
4040 if (value_ge(d
, mp
))
4043 f
= evalue_dup(&e
->x
.p
->arr
[0]);
4054 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
4055 if ((f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[i
], D
, max_periods
,
4062 static void replace_fract_by_affine(evalue
*e
, evalue
*f
, Value val
)
4066 if (value_notzero_p(e
->d
))
4069 offset
= type_offset(e
->x
.p
);
4070 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
4071 replace_fract_by_affine(&e
->x
.p
->arr
[i
], f
, val
);
4073 if (e
->x
.p
->type
!= fractional
)
4076 if (!eequal(&e
->x
.p
->arr
[0], f
))
4079 replace_by_affine(e
, val
);
4082 /* Look for fractional parts that can be removed by splitting the corresponding
4083 * domain into at most max_periods parts.
4084 * We use a very simply strategy that looks for the first fractional part
4085 * that satisfies the condition, performs the split and then continues
4086 * looking for other fractional parts in the split domains until no
4087 * such fractional part can be found anymore.
4089 void evalue_split_periods(evalue
*e
, int max_periods
, unsigned int MaxRays
)
4096 if (EVALUE_IS_ZERO(*e
))
4098 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
4100 "WARNING: evalue_split_periods called on incorrect evalue type\n");
4108 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
4113 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
4115 f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[2*i
+1], D
, max_periods
,
4120 M
= Matrix_Alloc(2, 2+D
->Dimension
);
4122 value_subtract(d
, max
, min
);
4123 n
= VALUE_TO_INT(d
)+1;
4125 value_set_si(M
->p
[0][0], 1);
4126 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
4127 value_multiply(d
, max
, T
->p
[1][D
->Dimension
]);
4128 value_subtract(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
], d
);
4129 value_set_si(d
, -1);
4130 value_set_si(M
->p
[1][0], 1);
4131 Vector_Scale(T
->p
[0], M
->p
[1]+1, d
, D
->Dimension
+1);
4132 value_addmul(M
->p
[1][1+D
->Dimension
], max
, T
->p
[1][D
->Dimension
]);
4133 value_addto(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4134 T
->p
[1][D
->Dimension
]);
4135 value_decrement(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
]);
4137 p
= new_enode(partition
, e
->x
.p
->size
+ (n
-1)*2, e
->x
.p
->pos
);
4138 for (j
= 0; j
< 2*i
; ++j
) {
4139 value_clear(p
->arr
[j
].d
);
4140 p
->arr
[j
] = e
->x
.p
->arr
[j
];
4142 for (j
= 2*i
+2; j
< e
->x
.p
->size
; ++j
) {
4143 value_clear(p
->arr
[j
+2*(n
-1)].d
);
4144 p
->arr
[j
+2*(n
-1)] = e
->x
.p
->arr
[j
];
4146 for (j
= n
-1; j
>= 0; --j
) {
4148 value_clear(p
->arr
[2*i
+1].d
);
4149 p
->arr
[2*i
+1] = e
->x
.p
->arr
[2*i
+1];
4151 evalue_copy(&p
->arr
[2*(i
+j
)+1], &e
->x
.p
->arr
[2*i
+1]);
4153 value_subtract(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4154 T
->p
[1][D
->Dimension
]);
4155 value_addto(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
],
4156 T
->p
[1][D
->Dimension
]);
4158 replace_fract_by_affine(&p
->arr
[2*(i
+j
)+1], f
, max
);
4159 E
= DomainAddConstraints(D
, M
, MaxRays
);
4160 EVALUE_SET_DOMAIN(p
->arr
[2*(i
+j
)], E
);
4161 if (evalue_range_reduction_in_domain(&p
->arr
[2*(i
+j
)+1], E
))
4162 reduce_evalue(&p
->arr
[2*(i
+j
)+1]);
4163 value_decrement(max
, max
);
4165 value_clear(e
->x
.p
->arr
[2*i
].d
);
4180 void evalue_extract_affine(const evalue
*e
, Value
*coeff
, Value
*cst
, Value
*d
)
4182 value_set_si(*d
, 1);
4184 for ( ; value_zero_p(e
->d
); e
= &e
->x
.p
->arr
[0]) {
4186 assert(e
->x
.p
->type
== polynomial
);
4187 assert(e
->x
.p
->size
== 2);
4188 c
= &e
->x
.p
->arr
[1];
4189 value_multiply(coeff
[e
->x
.p
->pos
-1], *d
, c
->x
.n
);
4190 value_division(coeff
[e
->x
.p
->pos
-1], coeff
[e
->x
.p
->pos
-1], c
->d
);
4192 value_multiply(*cst
, *d
, e
->x
.n
);
4193 value_division(*cst
, *cst
, e
->d
);
4196 /* returns an evalue that corresponds to
4200 static evalue
*term(int param
, Value c
, Value den
)
4202 evalue
*EP
= ALLOC(evalue
);
4204 value_set_si(EP
->d
,0);
4205 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
4206 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
4207 value_init(EP
->x
.p
->arr
[1].x
.n
);
4208 value_assign(EP
->x
.p
->arr
[1].d
, den
);
4209 value_assign(EP
->x
.p
->arr
[1].x
.n
, c
);
4213 evalue
*affine2evalue(Value
*coeff
, Value denom
, int nvar
)
4216 evalue
*E
= ALLOC(evalue
);
4218 evalue_set(E
, coeff
[nvar
], denom
);
4219 for (i
= 0; i
< nvar
; ++i
) {
4221 if (value_zero_p(coeff
[i
]))
4223 t
= term(i
, coeff
[i
], denom
);
4230 void evalue_substitute(evalue
*e
, evalue
**subs
)
4236 if (value_notzero_p(e
->d
))
4240 assert(p
->type
!= partition
);
4242 for (i
= 0; i
< p
->size
; ++i
)
4243 evalue_substitute(&p
->arr
[i
], subs
);
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
;