1 /***********************************************************************/
2 /* copyright 1997, Doran Wilde */
3 /* copyright 1997-2000, Vincent Loechner */
4 /* copyright 2003-2006, Sven Verdoolaege */
5 /* Permission is granted to copy, use, and distribute */
6 /* for any commercial or noncommercial purpose under the terms */
7 /* of the GNU General Public license, version 2, June 1991 */
8 /* (see file : LICENSE). */
9 /***********************************************************************/
15 #include <barvinok/evalue.h>
16 #include <barvinok/barvinok.h>
17 #include <barvinok/util.h>
19 #ifndef value_pmodulus
20 #define value_pmodulus(ref,val1,val2) (mpz_fdiv_r((ref),(val1),(val2)))
23 #define ALLOC(type) (type*)malloc(sizeof(type))
24 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
27 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
29 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
32 void evalue_set_si(evalue
*ev
, int n
, int d
) {
33 value_set_si(ev
->d
, d
);
35 value_set_si(ev
->x
.n
, n
);
38 void evalue_set(evalue
*ev
, Value n
, Value d
) {
39 value_assign(ev
->d
, d
);
41 value_assign(ev
->x
.n
, n
);
46 evalue
*EP
= ALLOC(evalue
);
48 evalue_set_si(EP
, 0, 1);
52 /* returns an evalue that corresponds to
56 evalue
*evalue_var(int var
)
58 evalue
*EP
= ALLOC(evalue
);
60 value_set_si(EP
->d
,0);
61 EP
->x
.p
= new_enode(polynomial
, 2, var
+ 1);
62 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
63 evalue_set_si(&EP
->x
.p
->arr
[1], 1, 1);
67 void aep_evalue(evalue
*e
, int *ref
) {
72 if (value_notzero_p(e
->d
))
73 return; /* a rational number, its already reduced */
75 return; /* hum... an overflow probably occured */
77 /* First check the components of p */
78 for (i
=0;i
<p
->size
;i
++)
79 aep_evalue(&p
->arr
[i
],ref
);
86 p
->pos
= ref
[p
->pos
-1]+1;
92 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
98 if (value_notzero_p(e
->d
))
99 return; /* a rational number, its already reduced */
101 return; /* hum... an overflow probably occured */
104 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
105 for(i
=0;i
<CT
->NbRows
-1;i
++)
106 for(j
=0;j
<CT
->NbColumns
;j
++)
107 if(value_notzero_p(CT
->p
[i
][j
])) {
112 /* Transform the references in e, using ref */
116 } /* addeliminatedparams_evalue */
118 static void addeliminatedparams_partition(enode
*p
, Matrix
*CT
, Polyhedron
*CEq
,
119 unsigned nparam
, unsigned MaxRays
)
122 assert(p
->type
== partition
);
125 for (i
= 0; i
< p
->size
/2; i
++) {
126 Polyhedron
*D
= EVALUE_DOMAIN(p
->arr
[2*i
]);
127 Polyhedron
*T
= DomainPreimage(D
, CT
, MaxRays
);
131 T
= DomainIntersection(D
, CEq
, MaxRays
);
134 EVALUE_SET_DOMAIN(p
->arr
[2*i
], T
);
138 void addeliminatedparams_enum(evalue
*e
, Matrix
*CT
, Polyhedron
*CEq
,
139 unsigned MaxRays
, unsigned nparam
)
144 if (CT
->NbRows
== CT
->NbColumns
)
147 if (EVALUE_IS_ZERO(*e
))
150 if (value_notzero_p(e
->d
)) {
153 value_set_si(res
.d
, 0);
154 res
.x
.p
= new_enode(partition
, 2, nparam
);
155 EVALUE_SET_DOMAIN(res
.x
.p
->arr
[0],
156 DomainConstraintSimplify(Polyhedron_Copy(CEq
), MaxRays
));
157 value_clear(res
.x
.p
->arr
[1].d
);
158 res
.x
.p
->arr
[1] = *e
;
166 addeliminatedparams_partition(p
, CT
, CEq
, nparam
, MaxRays
);
167 for (i
= 0; i
< p
->size
/2; i
++)
168 addeliminatedparams_evalue(&p
->arr
[2*i
+1], CT
);
171 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
179 assert(value_notzero_p(e1
->d
));
180 assert(value_notzero_p(e2
->d
));
181 value_multiply(m
, e1
->x
.n
, e2
->d
);
182 value_multiply(m2
, e2
->x
.n
, e1
->d
);
185 else if (value_gt(m
, m2
))
195 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
197 if (value_notzero_p(e1
->d
)) {
199 if (value_zero_p(e2
->d
))
201 r
= mod_rational_smaller(e1
, e2
);
202 return r
== -1 ? 0 : r
;
204 if (value_notzero_p(e2
->d
))
206 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
208 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
211 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
213 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
218 static int mod_term_smaller(const evalue
*e1
, const evalue
*e2
)
220 assert(value_zero_p(e1
->d
));
221 assert(value_zero_p(e2
->d
));
222 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
223 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
224 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
227 static void check_order(const evalue
*e
)
232 if (value_notzero_p(e
->d
))
235 switch (e
->x
.p
->type
) {
237 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
238 check_order(&e
->x
.p
->arr
[2*i
+1]);
241 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
243 if (value_notzero_p(a
->d
))
245 switch (a
->x
.p
->type
) {
247 assert(mod_term_smaller(&e
->x
.p
->arr
[0], &a
->x
.p
->arr
[0]));
256 for (i
= 0; i
< e
->x
.p
->size
; ++i
) {
258 if (value_notzero_p(a
->d
))
260 switch (a
->x
.p
->type
) {
262 assert(e
->x
.p
->pos
< a
->x
.p
->pos
);
273 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
275 if (value_notzero_p(a
->d
))
277 switch (a
->x
.p
->type
) {
288 /* Negative pos means inequality */
289 /* s is negative of substitution if m is not zero */
298 struct fixed_param
*fixed
;
303 static int relations_depth(evalue
*e
)
308 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
309 e
= &e
->x
.p
->arr
[1], ++d
);
313 static void poly_denom_not_constant(evalue
**pp
, Value
*d
)
318 while (value_zero_p(p
->d
)) {
319 assert(p
->x
.p
->type
== polynomial
);
320 assert(p
->x
.p
->size
== 2);
321 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
322 value_lcm(*d
, *d
, p
->x
.p
->arr
[1].d
);
328 static void poly_denom(evalue
*p
, Value
*d
)
330 poly_denom_not_constant(&p
, d
);
331 value_lcm(*d
, *d
, p
->d
);
334 static void realloc_substitution(struct subst
*s
, int d
)
336 struct fixed_param
*n
;
339 for (i
= 0; i
< s
->n
; ++i
)
346 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
352 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
355 /* May have been reduced already */
356 if (value_notzero_p(m
->d
))
359 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
360 assert(m
->x
.p
->size
== 3);
362 /* fractional was inverted during reduction
363 * invert it back and move constant in
365 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
366 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
367 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
368 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
369 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
370 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
371 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
372 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
373 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
374 value_set_si(m
->x
.p
->arr
[1].d
, 1);
377 /* Oops. Nested identical relations. */
378 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
381 if (s
->n
>= s
->max
) {
382 int d
= relations_depth(r
);
383 realloc_substitution(s
, d
);
387 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
388 assert(p
->x
.p
->size
== 2);
391 assert(value_pos_p(f
->x
.n
));
393 value_init(s
->fixed
[s
->n
].m
);
394 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
395 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
396 value_init(s
->fixed
[s
->n
].d
);
397 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
398 value_init(s
->fixed
[s
->n
].s
.d
);
399 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
405 static int type_offset(enode
*p
)
407 return p
->type
== fractional
? 1 :
408 p
->type
== flooring
? 1 : 0;
411 static void reorder_terms_about(enode
*p
, evalue
*v
)
414 int offset
= type_offset(p
);
416 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
418 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
419 free_evalue_refs(&(p
->arr
[i
]));
425 static void reorder_terms(evalue
*e
)
430 assert(value_zero_p(e
->d
));
432 assert(p
->type
== fractional
); /* for now */
435 value_set_si(f
.d
, 0);
436 f
.x
.p
= new_enode(fractional
, 3, -1);
437 value_clear(f
.x
.p
->arr
[0].d
);
438 f
.x
.p
->arr
[0] = p
->arr
[0];
439 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
440 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
441 reorder_terms_about(p
, &f
);
447 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
453 if (value_notzero_p(e
->d
)) {
455 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
456 return; /* a rational number, its already reduced */
460 return; /* hum... an overflow probably occured */
462 /* First reduce the components of p */
463 add
= p
->type
== relation
;
464 for (i
=0; i
<p
->size
; i
++) {
466 add
= add_modulo_substitution(s
, e
);
468 if (i
== 0 && p
->type
==fractional
)
469 _reduce_evalue(&p
->arr
[i
], s
, 1);
471 _reduce_evalue(&p
->arr
[i
], s
, fract
);
473 if (add
&& i
== p
->size
-1) {
475 value_clear(s
->fixed
[s
->n
].m
);
476 value_clear(s
->fixed
[s
->n
].d
);
477 free_evalue_refs(&s
->fixed
[s
->n
].s
);
478 } else if (add
&& i
== 1)
479 s
->fixed
[s
->n
-1].pos
*= -1;
482 if (p
->type
==periodic
) {
484 /* Try to reduce the period */
485 for (i
=1; i
<=(p
->size
)/2; i
++) {
486 if ((p
->size
% i
)==0) {
488 /* Can we reduce the size to i ? */
490 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
491 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
494 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
498 you_lose
: /* OK, lets not do it */
503 /* Try to reduce its strength */
506 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
510 else if (p
->type
==polynomial
) {
511 for (k
= 0; s
&& k
< s
->n
; ++k
) {
512 if (s
->fixed
[k
].pos
== p
->pos
) {
513 int divide
= value_notone_p(s
->fixed
[k
].d
);
516 if (value_notzero_p(s
->fixed
[k
].m
)) {
519 assert(p
->size
== 2);
520 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
522 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
529 value_assign(d
.d
, s
->fixed
[k
].d
);
531 if (value_notzero_p(s
->fixed
[k
].m
))
532 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
534 value_set_si(d
.x
.n
, 1);
537 for (i
=p
->size
-1;i
>=1;i
--) {
538 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
540 emul(&d
, &p
->arr
[i
]);
541 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
542 free_evalue_refs(&(p
->arr
[i
]));
545 _reduce_evalue(&p
->arr
[0], s
, fract
);
548 free_evalue_refs(&d
);
554 /* Try to reduce the degree */
555 for (i
=p
->size
-1;i
>=1;i
--) {
556 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
558 /* Zero coefficient */
559 free_evalue_refs(&(p
->arr
[i
]));
564 /* Try to reduce its strength */
567 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
571 else if (p
->type
==fractional
) {
575 if (value_notzero_p(p
->arr
[0].d
)) {
577 value_assign(v
.d
, p
->arr
[0].d
);
579 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
584 evalue
*pp
= &p
->arr
[0];
585 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
586 assert(pp
->x
.p
->size
== 2);
588 /* search for exact duplicate among the modulo inequalities */
590 f
= &pp
->x
.p
->arr
[1];
591 for (k
= 0; s
&& k
< s
->n
; ++k
) {
592 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
593 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
594 value_eq(s
->fixed
[k
].m
, f
->d
) &&
595 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
602 /* replace { E/m } by { (E-1)/m } + 1/m */
607 evalue_set_si(&extra
, 1, 1);
608 value_assign(extra
.d
, g
);
609 eadd(&extra
, &v
.x
.p
->arr
[1]);
610 free_evalue_refs(&extra
);
612 /* We've been going in circles; stop now */
613 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
614 free_evalue_refs(&v
);
616 evalue_set_si(&v
, 0, 1);
621 value_set_si(v
.d
, 0);
622 v
.x
.p
= new_enode(fractional
, 3, -1);
623 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
624 value_assign(v
.x
.p
->arr
[1].d
, g
);
625 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
626 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
629 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
632 value_division(f
->d
, g
, f
->d
);
633 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
634 value_assign(f
->d
, g
);
635 value_decrement(f
->x
.n
, f
->x
.n
);
636 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
638 value_gcd(g
, f
->d
, f
->x
.n
);
639 value_division(f
->d
, f
->d
, g
);
640 value_division(f
->x
.n
, f
->x
.n
, g
);
649 /* reduction may have made this fractional arg smaller */
650 i
= reorder
? p
->size
: 1;
651 for ( ; i
< p
->size
; ++i
)
652 if (value_zero_p(p
->arr
[i
].d
) &&
653 p
->arr
[i
].x
.p
->type
== fractional
&&
654 !mod_term_smaller(e
, &p
->arr
[i
]))
658 value_set_si(v
.d
, 0);
659 v
.x
.p
= new_enode(fractional
, 3, -1);
660 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
661 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
662 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
670 evalue
*pp
= &p
->arr
[0];
673 poly_denom_not_constant(&pp
, &m
);
674 mpz_fdiv_r(r
, m
, pp
->d
);
675 if (value_notzero_p(r
)) {
677 value_set_si(v
.d
, 0);
678 v
.x
.p
= new_enode(fractional
, 3, -1);
680 value_multiply(r
, m
, pp
->x
.n
);
681 value_multiply(v
.x
.p
->arr
[1].d
, m
, pp
->d
);
682 value_init(v
.x
.p
->arr
[1].x
.n
);
683 mpz_fdiv_r(v
.x
.p
->arr
[1].x
.n
, r
, pp
->d
);
684 mpz_fdiv_q(r
, r
, pp
->d
);
686 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
687 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
689 while (value_zero_p(pp
->d
))
690 pp
= &pp
->x
.p
->arr
[0];
692 value_assign(pp
->d
, m
);
693 value_assign(pp
->x
.n
, r
);
695 value_gcd(r
, pp
->d
, pp
->x
.n
);
696 value_division(pp
->d
, pp
->d
, r
);
697 value_division(pp
->x
.n
, pp
->x
.n
, r
);
710 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
711 pp
= &pp
->x
.p
->arr
[0]) {
712 f
= &pp
->x
.p
->arr
[1];
713 assert(value_pos_p(f
->d
));
714 mpz_mul_ui(twice
, f
->x
.n
, 2);
715 if (value_lt(twice
, f
->d
))
717 if (value_eq(twice
, f
->d
))
725 value_set_si(v
.d
, 0);
726 v
.x
.p
= new_enode(fractional
, 3, -1);
727 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
728 poly_denom(&p
->arr
[0], &twice
);
729 value_assign(v
.x
.p
->arr
[1].d
, twice
);
730 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
731 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
732 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
734 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
735 pp
= &pp
->x
.p
->arr
[0]) {
736 f
= &pp
->x
.p
->arr
[1];
737 value_oppose(f
->x
.n
, f
->x
.n
);
738 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
740 value_division(pp
->d
, twice
, pp
->d
);
741 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
742 value_assign(pp
->d
, twice
);
743 value_oppose(pp
->x
.n
, pp
->x
.n
);
744 value_decrement(pp
->x
.n
, pp
->x
.n
);
745 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
747 /* Maybe we should do this during reduction of
750 value_gcd(twice
, pp
->d
, pp
->x
.n
);
751 value_division(pp
->d
, pp
->d
, twice
);
752 value_division(pp
->x
.n
, pp
->x
.n
, twice
);
762 reorder_terms_about(p
, &v
);
763 _reduce_evalue(&p
->arr
[1], s
, fract
);
766 /* Try to reduce the degree */
767 for (i
=p
->size
-1;i
>=2;i
--) {
768 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
770 /* Zero coefficient */
771 free_evalue_refs(&(p
->arr
[i
]));
776 /* Try to reduce its strength */
779 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
780 free_evalue_refs(&(p
->arr
[0]));
784 else if (p
->type
== flooring
) {
785 /* Try to reduce the degree */
786 for (i
=p
->size
-1;i
>=2;i
--) {
787 if (!EVALUE_IS_ZERO(p
->arr
[i
]))
789 /* Zero coefficient */
790 free_evalue_refs(&(p
->arr
[i
]));
795 /* Try to reduce its strength */
798 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
799 free_evalue_refs(&(p
->arr
[0]));
803 else if (p
->type
== relation
) {
804 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
805 free_evalue_refs(&(p
->arr
[2]));
806 free_evalue_refs(&(p
->arr
[0]));
813 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
814 free_evalue_refs(&(p
->arr
[2]));
817 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
818 free_evalue_refs(&(p
->arr
[1]));
819 free_evalue_refs(&(p
->arr
[0]));
820 evalue_set_si(e
, 0, 1);
827 /* Relation was reduced by means of an identical
828 * inequality => remove
830 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
833 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
834 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
836 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
838 free_evalue_refs(&(p
->arr
[2]));
842 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
844 evalue_set_si(e
, 0, 1);
845 free_evalue_refs(&(p
->arr
[1]));
847 free_evalue_refs(&(p
->arr
[0]));
853 } /* reduce_evalue */
855 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
860 for (k
= 0; k
< dim
; ++k
)
861 if (value_notzero_p(row
[k
+1]))
864 Vector_Normalize_Positive(row
+1, dim
+1, k
);
865 assert(s
->n
< s
->max
);
866 value_init(s
->fixed
[s
->n
].d
);
867 value_init(s
->fixed
[s
->n
].m
);
868 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
869 s
->fixed
[s
->n
].pos
= k
+1;
870 value_set_si(s
->fixed
[s
->n
].m
, 0);
871 r
= &s
->fixed
[s
->n
].s
;
873 for (l
= k
+1; l
< dim
; ++l
)
874 if (value_notzero_p(row
[l
+1])) {
875 value_set_si(r
->d
, 0);
876 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
877 value_init(r
->x
.p
->arr
[1].x
.n
);
878 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
879 value_set_si(r
->x
.p
->arr
[1].d
, 1);
883 value_oppose(r
->x
.n
, row
[dim
+1]);
884 value_set_si(r
->d
, 1);
888 static void _reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
, struct subst
*s
)
891 Polyhedron
*orig
= D
;
896 D
= DomainConvex(D
, 0);
897 if (!D
->next
&& D
->NbEq
) {
901 realloc_substitution(s
, dim
);
903 int d
= relations_depth(e
);
905 NALLOC(s
->fixed
, s
->max
);
908 for (j
= 0; j
< D
->NbEq
; ++j
)
909 add_substitution(s
, D
->Constraint
[j
], dim
);
913 _reduce_evalue(e
, s
, 0);
916 for (j
= 0; j
< s
->n
; ++j
) {
917 value_clear(s
->fixed
[j
].d
);
918 value_clear(s
->fixed
[j
].m
);
919 free_evalue_refs(&s
->fixed
[j
].s
);
924 void reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
)
926 struct subst s
= { NULL
, 0, 0 };
928 if (EVALUE_IS_ZERO(*e
))
932 evalue_set_si(e
, 0, 1);
935 _reduce_evalue_in_domain(e
, D
, &s
);
940 void reduce_evalue (evalue
*e
) {
941 struct subst s
= { NULL
, 0, 0 };
943 if (value_notzero_p(e
->d
))
944 return; /* a rational number, its already reduced */
946 if (e
->x
.p
->type
== partition
) {
949 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
950 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
952 /* This shouldn't really happen;
953 * Empty domains should not be added.
955 POL_ENSURE_VERTICES(D
);
957 _reduce_evalue_in_domain(&e
->x
.p
->arr
[2*i
+1], D
, &s
);
959 if (emptyQ(D
) || EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
960 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
961 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
962 value_clear(e
->x
.p
->arr
[2*i
].d
);
964 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
965 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
969 if (e
->x
.p
->size
== 0) {
971 evalue_set_si(e
, 0, 1);
974 _reduce_evalue(e
, &s
, 0);
979 static void print_evalue_r(FILE *DST
, const evalue
*e
, const char *const *pname
)
981 if(value_notzero_p(e
->d
)) {
982 if(value_notone_p(e
->d
)) {
983 value_print(DST
,VALUE_FMT
,e
->x
.n
);
985 value_print(DST
,VALUE_FMT
,e
->d
);
988 value_print(DST
,VALUE_FMT
,e
->x
.n
);
992 print_enode(DST
,e
->x
.p
,pname
);
996 void print_evalue(FILE *DST
, const evalue
*e
, const char * const *pname
)
998 print_evalue_r(DST
, e
, pname
);
999 if (value_notzero_p(e
->d
))
1003 void print_enode(FILE *DST
, enode
*p
, const char *const *pname
)
1008 fprintf(DST
, "NULL");
1014 for (i
=0; i
<p
->size
; i
++) {
1015 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1019 fprintf(DST
, " }\n");
1023 for (i
=p
->size
-1; i
>=0; i
--) {
1024 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1025 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
1027 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
1029 fprintf(DST
, " )\n");
1033 for (i
=0; i
<p
->size
; i
++) {
1034 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1035 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
1037 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
1042 for (i
=p
->size
-1; i
>=1; i
--) {
1043 print_evalue_r(DST
, &p
->arr
[i
], pname
);
1045 fprintf(DST
, " * ");
1046 fprintf(DST
, p
->type
== flooring
? "[" : "{");
1047 print_evalue_r(DST
, &p
->arr
[0], pname
);
1048 fprintf(DST
, p
->type
== flooring
? "]" : "}");
1050 fprintf(DST
, "^%d + ", i
-1);
1052 fprintf(DST
, " + ");
1055 fprintf(DST
, " )\n");
1059 print_evalue_r(DST
, &p
->arr
[0], pname
);
1060 fprintf(DST
, "= 0 ] * \n");
1061 print_evalue_r(DST
, &p
->arr
[1], pname
);
1063 fprintf(DST
, " +\n [ ");
1064 print_evalue_r(DST
, &p
->arr
[0], pname
);
1065 fprintf(DST
, "!= 0 ] * \n");
1066 print_evalue_r(DST
, &p
->arr
[2], pname
);
1070 char **new_names
= NULL
;
1071 const char *const *names
= pname
;
1072 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
1073 if (!pname
|| p
->pos
< maxdim
) {
1074 new_names
= ALLOCN(char *, maxdim
);
1075 for (i
= 0; i
< p
->pos
; ++i
) {
1077 new_names
[i
] = (char *)pname
[i
];
1079 new_names
[i
] = ALLOCN(char, 10);
1080 snprintf(new_names
[i
], 10, "%c", 'P'+i
);
1083 for ( ; i
< maxdim
; ++i
) {
1084 new_names
[i
] = ALLOCN(char, 10);
1085 snprintf(new_names
[i
], 10, "_p%d", i
);
1087 names
= (const char**)new_names
;
1090 for (i
=0; i
<p
->size
/2; i
++) {
1091 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
1092 print_evalue_r(DST
, &p
->arr
[2*i
+1], names
);
1093 if (value_notzero_p(p
->arr
[2*i
+1].d
))
1097 if (!pname
|| p
->pos
< maxdim
) {
1098 for (i
= pname
? p
->pos
: 0; i
< maxdim
; ++i
)
1111 static void eadd_rev(const evalue
*e1
, evalue
*res
)
1115 evalue_copy(&ev
, e1
);
1117 free_evalue_refs(res
);
1121 static void eadd_rev_cst(const evalue
*e1
, evalue
*res
)
1125 evalue_copy(&ev
, e1
);
1126 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
1127 free_evalue_refs(res
);
1131 static int is_zero_on(evalue
*e
, Polyhedron
*D
)
1136 tmp
.x
.p
= new_enode(partition
, 2, D
->Dimension
);
1137 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Domain_Copy(D
));
1138 evalue_copy(&tmp
.x
.p
->arr
[1], e
);
1139 reduce_evalue(&tmp
);
1140 is_zero
= EVALUE_IS_ZERO(tmp
);
1141 free_evalue_refs(&tmp
);
1145 struct section
{ Polyhedron
* D
; evalue E
; };
1147 void eadd_partitions(const evalue
*e1
, evalue
*res
)
1152 s
= (struct section
*)
1153 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
1154 sizeof(struct section
));
1156 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1157 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1158 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1161 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1162 assert(res
->x
.p
->size
>= 2);
1163 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1164 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
1166 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
1168 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1177 /* See if we can extend one of the domains in res to cover fd */
1178 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1179 if (is_zero_on(&res
->x
.p
->arr
[2*i
+1], fd
))
1181 if (i
< res
->x
.p
->size
/2) {
1182 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
1183 DomainConcat(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
])));
1186 value_init(s
[n
].E
.d
);
1187 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
1191 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1192 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
1193 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1195 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1196 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1202 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
1203 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1205 value_init(s
[n
].E
.d
);
1206 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1207 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1208 if (!emptyQ(fd
) && is_zero_on(&e1
->x
.p
->arr
[2*j
+1], fd
)) {
1209 d
= DomainConcat(fd
, d
);
1210 fd
= Empty_Polyhedron(fd
->Dimension
);
1216 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
1220 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1223 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1224 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1225 value_clear(res
->x
.p
->arr
[2*i
].d
);
1230 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1231 for (j
= 0; j
< n
; ++j
) {
1232 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1233 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1234 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1235 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1241 static void explicit_complement(evalue
*res
)
1243 enode
*rel
= new_enode(relation
, 3, 0);
1245 value_clear(rel
->arr
[0].d
);
1246 rel
->arr
[0] = res
->x
.p
->arr
[0];
1247 value_clear(rel
->arr
[1].d
);
1248 rel
->arr
[1] = res
->x
.p
->arr
[1];
1249 value_set_si(rel
->arr
[2].d
, 1);
1250 value_init(rel
->arr
[2].x
.n
);
1251 value_set_si(rel
->arr
[2].x
.n
, 0);
1256 static void reduce_constant(evalue
*e
)
1261 value_gcd(g
, e
->x
.n
, e
->d
);
1262 if (value_notone_p(g
)) {
1263 value_division(e
->d
, e
->d
,g
);
1264 value_division(e
->x
.n
, e
->x
.n
,g
);
1269 void eadd(const evalue
*e1
, evalue
*res
)
1273 if (EVALUE_IS_ZERO(*e1
))
1276 if (EVALUE_IS_ZERO(*res
)) {
1277 if (value_notzero_p(e1
->d
)) {
1278 value_assign(res
->d
, e1
->d
);
1279 value_assign(res
->x
.n
, e1
->x
.n
);
1281 value_clear(res
->x
.n
);
1282 value_set_si(res
->d
, 0);
1283 res
->x
.p
= ecopy(e1
->x
.p
);
1288 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1289 /* Add two rational numbers */
1290 if (value_eq(e1
->d
, res
->d
))
1291 value_addto(res
->x
.n
, res
->x
.n
, e1
->x
.n
);
1293 value_multiply(res
->x
.n
, res
->x
.n
, e1
->d
);
1294 value_addmul(res
->x
.n
, e1
->x
.n
, res
->d
);
1295 value_multiply(res
->d
,e1
->d
,res
->d
);
1297 reduce_constant(res
);
1300 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
1301 switch (res
->x
.p
->type
) {
1303 /* Add the constant to the constant term of a polynomial*/
1304 eadd(e1
, &res
->x
.p
->arr
[0]);
1307 /* Add the constant to all elements of a periodic number */
1308 for (i
=0; i
<res
->x
.p
->size
; i
++) {
1309 eadd(e1
, &res
->x
.p
->arr
[i
]);
1313 fprintf(stderr
, "eadd: cannot add const with vector\n");
1317 eadd(e1
, &res
->x
.p
->arr
[1]);
1320 assert(EVALUE_IS_ZERO(*e1
));
1321 break; /* Do nothing */
1323 /* Create (zero) complement if needed */
1324 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
1325 explicit_complement(res
);
1326 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1327 eadd(e1
, &res
->x
.p
->arr
[i
]);
1333 /* add polynomial or periodic to constant
1334 * you have to exchange e1 and res, before doing addition */
1336 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1340 else { // ((e1->d==0) && (res->d==0))
1341 assert(!((e1
->x
.p
->type
== partition
) ^
1342 (res
->x
.p
->type
== partition
)));
1343 if (e1
->x
.p
->type
== partition
) {
1344 eadd_partitions(e1
, res
);
1347 if (e1
->x
.p
->type
== relation
&&
1348 (res
->x
.p
->type
!= relation
||
1349 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1353 if (res
->x
.p
->type
== relation
) {
1354 if (e1
->x
.p
->type
== relation
&&
1355 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1356 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1357 explicit_complement(res
);
1358 for (i
= 1; i
< e1
->x
.p
->size
; ++i
)
1359 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1362 if (res
->x
.p
->size
< 3)
1363 explicit_complement(res
);
1364 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1365 eadd(e1
, &res
->x
.p
->arr
[i
]);
1368 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
1369 /* adding to evalues of different type. two cases are possible
1370 * res is periodic and e1 is polynomial, you have to exchange
1371 * e1 and res then to add e1 to the constant term of res */
1372 if (e1
->x
.p
->type
== polynomial
) {
1373 eadd_rev_cst(e1
, res
);
1375 else if (res
->x
.p
->type
== polynomial
) {
1376 /* res is polynomial and e1 is periodic,
1377 add e1 to the constant term of res */
1379 eadd(e1
,&res
->x
.p
->arr
[0]);
1385 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
1386 ((res
->x
.p
->type
== fractional
||
1387 res
->x
.p
->type
== flooring
) &&
1388 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1389 /* adding evalues of different position (i.e function of different unknowns
1390 * to case are possible */
1392 switch (res
->x
.p
->type
) {
1395 if (mod_term_smaller(res
, e1
))
1396 eadd(e1
,&res
->x
.p
->arr
[1]);
1398 eadd_rev_cst(e1
, res
);
1400 case polynomial
: // res and e1 are polynomials
1401 // add e1 to the constant term of res
1403 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1404 eadd(e1
,&res
->x
.p
->arr
[0]);
1406 eadd_rev_cst(e1
, res
);
1407 // value_clear(g); value_clear(m1); value_clear(m2);
1409 case periodic
: // res and e1 are pointers to periodic numbers
1410 //add e1 to all elements of res
1412 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1413 for (i
=0;i
<res
->x
.p
->size
;i
++) {
1414 eadd(e1
,&res
->x
.p
->arr
[i
]);
1425 //same type , same pos and same size
1426 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1427 // add any element in e1 to the corresponding element in res
1428 i
= type_offset(res
->x
.p
);
1430 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1431 for (; i
<res
->x
.p
->size
; i
++) {
1432 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1437 /* Sizes are different */
1438 switch(res
->x
.p
->type
) {
1442 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1443 /* new enode and add res to that new node. If you do not do */
1444 /* that, you lose the the upper weight part of e1 ! */
1446 if(e1
->x
.p
->size
> res
->x
.p
->size
)
1449 i
= type_offset(res
->x
.p
);
1451 assert(eequal(&e1
->x
.p
->arr
[0],
1452 &res
->x
.p
->arr
[0]));
1453 for (; i
<e1
->x
.p
->size
; i
++) {
1454 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1461 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1464 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1465 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1466 to add periodicaly elements of e1 to elements of ne, and finaly to
1471 value_init(ex
); value_init(ey
);value_init(ep
);
1474 value_set_si(ex
,e1
->x
.p
->size
);
1475 value_set_si(ey
,res
->x
.p
->size
);
1476 value_assign (ep
,*Lcm(ex
,ey
));
1477 p
=(int)mpz_get_si(ep
);
1478 ne
= (evalue
*) malloc (sizeof(evalue
));
1480 value_set_si( ne
->d
,0);
1482 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1484 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1487 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1490 value_assign(res
->d
, ne
->d
);
1496 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
1505 static void emul_rev(const evalue
*e1
, evalue
*res
)
1509 evalue_copy(&ev
, e1
);
1511 free_evalue_refs(res
);
1515 static void emul_poly(const evalue
*e1
, evalue
*res
)
1517 int i
, j
, offset
= type_offset(res
->x
.p
);
1520 int size
= (e1
->x
.p
->size
+ res
->x
.p
->size
- offset
- 1);
1522 p
= new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1524 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1525 if (!EVALUE_IS_ZERO(e1
->x
.p
->arr
[i
]))
1528 /* special case pure power */
1529 if (i
== e1
->x
.p
->size
-1) {
1531 value_clear(p
->arr
[0].d
);
1532 p
->arr
[0] = res
->x
.p
->arr
[0];
1534 for (i
= offset
; i
< e1
->x
.p
->size
-1; ++i
)
1535 evalue_set_si(&p
->arr
[i
], 0, 1);
1536 for (i
= offset
; i
< res
->x
.p
->size
; ++i
) {
1537 value_clear(p
->arr
[i
+e1
->x
.p
->size
-offset
-1].d
);
1538 p
->arr
[i
+e1
->x
.p
->size
-offset
-1] = res
->x
.p
->arr
[i
];
1539 emul(&e1
->x
.p
->arr
[e1
->x
.p
->size
-1],
1540 &p
->arr
[i
+e1
->x
.p
->size
-offset
-1]);
1548 value_set_si(tmp
.d
,0);
1551 evalue_copy(&p
->arr
[0], &e1
->x
.p
->arr
[0]);
1552 for (i
= offset
; i
< e1
->x
.p
->size
; i
++) {
1553 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1554 emul(&res
->x
.p
->arr
[offset
], &tmp
.x
.p
->arr
[i
]);
1557 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1558 for (i
= offset
+1; i
<res
->x
.p
->size
; i
++)
1559 for (j
= offset
; j
<e1
->x
.p
->size
; j
++) {
1562 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1563 emul(&res
->x
.p
->arr
[i
], &ev
);
1564 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-offset
]);
1565 free_evalue_refs(&ev
);
1567 free_evalue_refs(res
);
1571 void emul_partitions(const evalue
*e1
, evalue
*res
)
1576 s
= (struct section
*)
1577 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1578 sizeof(struct section
));
1580 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1581 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1582 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1585 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1586 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1587 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1588 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1594 /* This code is only needed because the partitions
1595 are not true partitions.
1597 for (k
= 0; k
< n
; ++k
) {
1598 if (DomainIncludes(s
[k
].D
, d
))
1600 if (DomainIncludes(d
, s
[k
].D
)) {
1601 Domain_Free(s
[k
].D
);
1602 free_evalue_refs(&s
[k
].E
);
1613 value_init(s
[n
].E
.d
);
1614 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1615 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1619 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1620 value_clear(res
->x
.p
->arr
[2*i
].d
);
1621 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1626 evalue_set_si(res
, 0, 1);
1628 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1629 for (j
= 0; j
< n
; ++j
) {
1630 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1631 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1632 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1633 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1640 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1642 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1643 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1644 * evalues is not treated here */
1646 void emul(const evalue
*e1
, evalue
*res
)
1650 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1651 fprintf(stderr
, "emul: do not proced on evector type !\n");
1655 if (EVALUE_IS_ZERO(*res
))
1658 if (EVALUE_IS_ONE(*e1
))
1661 if (EVALUE_IS_ZERO(*e1
)) {
1662 if (value_notzero_p(res
->d
)) {
1663 value_assign(res
->d
, e1
->d
);
1664 value_assign(res
->x
.n
, e1
->x
.n
);
1666 free_evalue_refs(res
);
1668 evalue_set_si(res
, 0, 1);
1673 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1674 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1675 emul_partitions(e1
, res
);
1678 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1679 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1680 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1682 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1683 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1684 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1685 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3) {
1686 free_evalue_refs(&res
->x
.p
->arr
[2]);
1689 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1690 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1693 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1694 emul(e1
, &res
->x
.p
->arr
[i
]);
1696 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1697 switch(e1
->x
.p
->type
) {
1699 switch(res
->x
.p
->type
) {
1701 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1702 /* Product of two polynomials of the same variable */
1707 /* Product of two polynomials of different variables */
1709 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1710 for( i
=0; i
<res
->x
.p
->size
; i
++)
1711 emul(e1
, &res
->x
.p
->arr
[i
]);
1720 /* Product of a polynomial and a periodic or fractional */
1727 switch(res
->x
.p
->type
) {
1729 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1730 /* Product of two periodics of the same parameter and period */
1732 for(i
=0; i
<res
->x
.p
->size
;i
++)
1733 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1738 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1739 /* Product of two periodics of the same parameter and different periods */
1743 value_init(x
); value_init(y
);value_init(z
);
1746 value_set_si(x
,e1
->x
.p
->size
);
1747 value_set_si(y
,res
->x
.p
->size
);
1748 value_assign (z
,*Lcm(x
,y
));
1749 lcm
=(int)mpz_get_si(z
);
1750 newp
= (evalue
*) malloc (sizeof(evalue
));
1751 value_init(newp
->d
);
1752 value_set_si( newp
->d
,0);
1753 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1754 for(i
=0;i
<lcm
;i
++) {
1755 evalue_copy(&newp
->x
.p
->arr
[i
],
1756 &res
->x
.p
->arr
[i
%iy
]);
1759 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1761 value_assign(res
->d
,newp
->d
);
1764 value_clear(x
); value_clear(y
);value_clear(z
);
1768 /* Product of two periodics of different parameters */
1770 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1771 for(i
=0; i
<res
->x
.p
->size
; i
++)
1772 emul(e1
, &(res
->x
.p
->arr
[i
]));
1780 /* Product of a periodic and a polynomial */
1782 for(i
=0; i
<res
->x
.p
->size
; i
++)
1783 emul(e1
, &(res
->x
.p
->arr
[i
]));
1790 switch(res
->x
.p
->type
) {
1792 for(i
=0; i
<res
->x
.p
->size
; i
++)
1793 emul(e1
, &(res
->x
.p
->arr
[i
]));
1800 assert(e1
->x
.p
->type
== res
->x
.p
->type
);
1801 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1802 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1805 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1806 if (e1
->x
.p
->type
!= fractional
|| !value_two_p(d
.d
))
1811 value_set_si(d
.x
.n
, 1);
1812 /* { x }^2 == { x }/2 */
1813 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1814 assert(e1
->x
.p
->size
== 3);
1815 assert(res
->x
.p
->size
== 3);
1817 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1819 eadd(&res
->x
.p
->arr
[1], &tmp
);
1820 emul(&e1
->x
.p
->arr
[2], &tmp
);
1821 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[1]);
1822 emul(&e1
->x
.p
->arr
[1], &res
->x
.p
->arr
[2]);
1823 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1824 free_evalue_refs(&tmp
);
1829 if(mod_term_smaller(res
, e1
))
1830 for(i
=1; i
<res
->x
.p
->size
; i
++)
1831 emul(e1
, &(res
->x
.p
->arr
[i
]));
1846 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1847 /* Product of two rational numbers */
1848 value_multiply(res
->d
,e1
->d
,res
->d
);
1849 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1850 reduce_constant(res
);
1854 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1855 /* Product of an expression (polynomial or peririodic) and a rational number */
1861 /* Product of a rationel number and an expression (polynomial or peririodic) */
1863 i
= type_offset(res
->x
.p
);
1864 for (; i
<res
->x
.p
->size
; i
++)
1865 emul(e1
, &res
->x
.p
->arr
[i
]);
1875 /* Frees mask content ! */
1876 void emask(evalue
*mask
, evalue
*res
) {
1883 if (EVALUE_IS_ZERO(*res
)) {
1884 free_evalue_refs(mask
);
1888 assert(value_zero_p(mask
->d
));
1889 assert(mask
->x
.p
->type
== partition
);
1890 assert(value_zero_p(res
->d
));
1891 assert(res
->x
.p
->type
== partition
);
1892 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1893 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1894 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1895 pos
= res
->x
.p
->pos
;
1897 s
= (struct section
*)
1898 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1899 sizeof(struct section
));
1903 evalue_set_si(&mone
, -1, 1);
1906 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1907 assert(mask
->x
.p
->size
>= 2);
1908 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1909 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1911 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1913 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1922 value_init(s
[n
].E
.d
);
1923 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1927 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1928 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1931 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1932 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1933 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1934 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1936 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1937 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1943 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1944 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1946 value_init(s
[n
].E
.d
);
1947 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1948 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1954 /* Just ignore; this may have been previously masked off */
1956 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1960 free_evalue_refs(&mone
);
1961 free_evalue_refs(mask
);
1962 free_evalue_refs(res
);
1965 evalue_set_si(res
, 0, 1);
1967 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1968 for (j
= 0; j
< n
; ++j
) {
1969 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1970 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1971 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1978 void evalue_copy(evalue
*dst
, const evalue
*src
)
1980 value_assign(dst
->d
, src
->d
);
1981 if(value_notzero_p(src
->d
)) {
1982 value_init(dst
->x
.n
);
1983 value_assign(dst
->x
.n
, src
->x
.n
);
1985 dst
->x
.p
= ecopy(src
->x
.p
);
1988 evalue
*evalue_dup(const evalue
*e
)
1990 evalue
*res
= ALLOC(evalue
);
1992 evalue_copy(res
, e
);
1996 enode
*new_enode(enode_type type
,int size
,int pos
) {
2002 fprintf(stderr
, "Allocating enode of size 0 !\n" );
2005 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
2009 for(i
=0; i
<size
; i
++) {
2010 value_init(res
->arr
[i
].d
);
2011 value_set_si(res
->arr
[i
].d
,0);
2012 res
->arr
[i
].x
.p
= 0;
2017 enode
*ecopy(enode
*e
) {
2022 res
= new_enode(e
->type
,e
->size
,e
->pos
);
2023 for(i
=0;i
<e
->size
;++i
) {
2024 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
2025 if(value_zero_p(res
->arr
[i
].d
))
2026 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
2027 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
2028 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
2030 value_init(res
->arr
[i
].x
.n
);
2031 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
2037 int ecmp(const evalue
*e1
, const evalue
*e2
)
2043 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
2047 value_multiply(m
, e1
->x
.n
, e2
->d
);
2048 value_multiply(m2
, e2
->x
.n
, e1
->d
);
2050 if (value_lt(m
, m2
))
2052 else if (value_gt(m
, m2
))
2062 if (value_notzero_p(e1
->d
))
2064 if (value_notzero_p(e2
->d
))
2070 if (p1
->type
!= p2
->type
)
2071 return p1
->type
- p2
->type
;
2072 if (p1
->pos
!= p2
->pos
)
2073 return p1
->pos
- p2
->pos
;
2074 if (p1
->size
!= p2
->size
)
2075 return p1
->size
- p2
->size
;
2077 for (i
= p1
->size
-1; i
>= 0; --i
)
2078 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
2084 int eequal(const evalue
*e1
, const evalue
*e2
)
2089 if (value_ne(e1
->d
,e2
->d
))
2092 /* e1->d == e2->d */
2093 if (value_notzero_p(e1
->d
)) {
2094 if (value_ne(e1
->x
.n
,e2
->x
.n
))
2097 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2101 /* e1->d == e2->d == 0 */
2104 if (p1
->type
!= p2
->type
) return 0;
2105 if (p1
->size
!= p2
->size
) return 0;
2106 if (p1
->pos
!= p2
->pos
) return 0;
2107 for (i
=0; i
<p1
->size
; i
++)
2108 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
2113 void free_evalue_refs(evalue
*e
) {
2118 if (EVALUE_IS_DOMAIN(*e
)) {
2119 Domain_Free(EVALUE_DOMAIN(*e
));
2122 } else if (value_pos_p(e
->d
)) {
2124 /* 'e' stores a constant */
2126 value_clear(e
->x
.n
);
2129 assert(value_zero_p(e
->d
));
2132 if (!p
) return; /* null pointer */
2133 for (i
=0; i
<p
->size
; i
++) {
2134 free_evalue_refs(&(p
->arr
[i
]));
2138 } /* free_evalue_refs */
2140 void evalue_free(evalue
*e
)
2142 free_evalue_refs(e
);
2146 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
2147 Vector
* val
, evalue
*res
)
2149 unsigned nparam
= periods
->Size
;
2152 double d
= compute_evalue(e
, val
->p
);
2153 d
*= VALUE_TO_DOUBLE(m
);
2158 value_assign(res
->d
, m
);
2159 value_init(res
->x
.n
);
2160 value_set_double(res
->x
.n
, d
);
2161 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
2164 if (value_one_p(periods
->p
[p
]))
2165 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
2170 value_assign(tmp
, periods
->p
[p
]);
2171 value_set_si(res
->d
, 0);
2172 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
2174 value_decrement(tmp
, tmp
);
2175 value_assign(val
->p
[p
], tmp
);
2176 mod2table_r(e
, periods
, m
, p
+1, val
,
2177 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
2178 } while (value_pos_p(tmp
));
2184 static void rel2table(evalue
*e
, int zero
)
2186 if (value_pos_p(e
->d
)) {
2187 if (value_zero_p(e
->x
.n
) == zero
)
2188 value_set_si(e
->x
.n
, 1);
2190 value_set_si(e
->x
.n
, 0);
2191 value_set_si(e
->d
, 1);
2194 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
2195 rel2table(&e
->x
.p
->arr
[i
], zero
);
2199 void evalue_mod2table(evalue
*e
, int nparam
)
2204 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2207 for (i
=0; i
<p
->size
; i
++) {
2208 evalue_mod2table(&(p
->arr
[i
]), nparam
);
2210 if (p
->type
== relation
) {
2215 evalue_copy(©
, &p
->arr
[0]);
2217 rel2table(&p
->arr
[0], 1);
2218 emul(&p
->arr
[0], &p
->arr
[1]);
2220 rel2table(©
, 0);
2221 emul(©
, &p
->arr
[2]);
2222 eadd(&p
->arr
[2], &p
->arr
[1]);
2223 free_evalue_refs(&p
->arr
[2]);
2224 free_evalue_refs(©
);
2226 free_evalue_refs(&p
->arr
[0]);
2230 } else if (p
->type
== fractional
) {
2231 Vector
*periods
= Vector_Alloc(nparam
);
2232 Vector
*val
= Vector_Alloc(nparam
);
2238 value_set_si(tmp
, 1);
2239 Vector_Set(periods
->p
, 1, nparam
);
2240 Vector_Set(val
->p
, 0, nparam
);
2241 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
2244 assert(p
->type
== polynomial
);
2245 assert(p
->size
== 2);
2246 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
2247 value_lcm(tmp
, tmp
, p
->arr
[1].d
);
2249 value_lcm(tmp
, tmp
, ev
->d
);
2251 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
2254 evalue_set_si(&res
, 0, 1);
2255 /* Compute the polynomial using Horner's rule */
2256 for (i
=p
->size
-1;i
>1;i
--) {
2257 eadd(&p
->arr
[i
], &res
);
2260 eadd(&p
->arr
[1], &res
);
2262 free_evalue_refs(e
);
2263 free_evalue_refs(&EP
);
2268 Vector_Free(periods
);
2270 } /* evalue_mod2table */
2272 /********************************************************/
2273 /* function in domain */
2274 /* check if the parameters in list_args */
2275 /* verifies the constraints of Domain P */
2276 /********************************************************/
2277 int in_domain(Polyhedron
*P
, Value
*list_args
)
2280 Value v
; /* value of the constraint of a row when
2281 parameters are instantiated*/
2285 for (row
= 0; row
< P
->NbConstraints
; row
++) {
2286 Inner_Product(P
->Constraint
[row
]+1, list_args
, P
->Dimension
, &v
);
2287 value_addto(v
, v
, P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
2288 if (value_neg_p(v
) ||
2289 value_zero_p(P
->Constraint
[row
][0]) && value_notzero_p(v
)) {
2296 return in
|| (P
->next
&& in_domain(P
->next
, list_args
));
2299 /****************************************************/
2300 /* function compute enode */
2301 /* compute the value of enode p with parameters */
2302 /* list "list_args */
2303 /* compute the polynomial or the periodic */
2304 /****************************************************/
2306 static double compute_enode(enode
*p
, Value
*list_args
) {
2318 if (p
->type
== polynomial
) {
2320 value_assign(param
,list_args
[p
->pos
-1]);
2322 /* Compute the polynomial using Horner's rule */
2323 for (i
=p
->size
-1;i
>0;i
--) {
2324 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2325 res
*=VALUE_TO_DOUBLE(param
);
2327 res
+=compute_evalue(&p
->arr
[0],list_args
);
2329 else if (p
->type
== fractional
) {
2330 double d
= compute_evalue(&p
->arr
[0], list_args
);
2331 d
-= floor(d
+1e-10);
2333 /* Compute the polynomial using Horner's rule */
2334 for (i
=p
->size
-1;i
>1;i
--) {
2335 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2338 res
+=compute_evalue(&p
->arr
[1],list_args
);
2340 else if (p
->type
== flooring
) {
2341 double d
= compute_evalue(&p
->arr
[0], list_args
);
2344 /* Compute the polynomial using Horner's rule */
2345 for (i
=p
->size
-1;i
>1;i
--) {
2346 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2349 res
+=compute_evalue(&p
->arr
[1],list_args
);
2351 else if (p
->type
== periodic
) {
2352 value_assign(m
,list_args
[p
->pos
-1]);
2354 /* Choose the right element of the periodic */
2355 value_set_si(param
,p
->size
);
2356 value_pmodulus(m
,m
,param
);
2357 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2359 else if (p
->type
== relation
) {
2360 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2361 res
= compute_evalue(&p
->arr
[1], list_args
);
2362 else if (p
->size
> 2)
2363 res
= compute_evalue(&p
->arr
[2], list_args
);
2365 else if (p
->type
== partition
) {
2366 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2367 Value
*vals
= list_args
;
2370 for (i
= 0; i
< dim
; ++i
) {
2371 value_init(vals
[i
]);
2373 value_assign(vals
[i
], list_args
[i
]);
2376 for (i
= 0; i
< p
->size
/2; ++i
)
2377 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2378 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2382 for (i
= 0; i
< dim
; ++i
)
2383 value_clear(vals
[i
]);
2392 } /* compute_enode */
2394 /*************************************************/
2395 /* return the value of Ehrhart Polynomial */
2396 /* It returns a double, because since it is */
2397 /* a recursive function, some intermediate value */
2398 /* might not be integral */
2399 /*************************************************/
2401 double compute_evalue(const evalue
*e
, Value
*list_args
)
2405 if (value_notzero_p(e
->d
)) {
2406 if (value_notone_p(e
->d
))
2407 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2409 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2412 res
= compute_enode(e
->x
.p
,list_args
);
2414 } /* compute_evalue */
2417 /****************************************************/
2418 /* function compute_poly : */
2419 /* Check for the good validity domain */
2420 /* return the number of point in the Polyhedron */
2421 /* in allocated memory */
2422 /* Using the Ehrhart pseudo-polynomial */
2423 /****************************************************/
2424 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2427 /* double d; int i; */
2429 tmp
= (Value
*) malloc (sizeof(Value
));
2430 assert(tmp
!= NULL
);
2432 value_set_si(*tmp
,0);
2435 return(tmp
); /* no ehrhart polynomial */
2436 if(en
->ValidityDomain
) {
2437 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2438 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2443 return(tmp
); /* no Validity Domain */
2445 if(in_domain(en
->ValidityDomain
,list_args
)) {
2447 #ifdef EVAL_EHRHART_DEBUG
2448 Print_Domain(stdout
,en
->ValidityDomain
);
2449 print_evalue(stdout
,&en
->EP
);
2452 /* d = compute_evalue(&en->EP,list_args);
2454 printf("(double)%lf = %d\n", d, i ); */
2455 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2461 value_set_si(*tmp
,0);
2462 return(tmp
); /* no compatible domain with the arguments */
2463 } /* compute_poly */
2465 static evalue
*eval_polynomial(const enode
*p
, int offset
,
2466 evalue
*base
, Value
*values
)
2471 res
= evalue_zero();
2472 for (i
= p
->size
-1; i
> offset
; --i
) {
2473 c
= evalue_eval(&p
->arr
[i
], values
);
2478 c
= evalue_eval(&p
->arr
[offset
], values
);
2485 evalue
*evalue_eval(const evalue
*e
, Value
*values
)
2492 if (value_notzero_p(e
->d
)) {
2493 res
= ALLOC(evalue
);
2495 evalue_copy(res
, e
);
2498 switch (e
->x
.p
->type
) {
2500 value_init(param
.x
.n
);
2501 value_assign(param
.x
.n
, values
[e
->x
.p
->pos
-1]);
2502 value_init(param
.d
);
2503 value_set_si(param
.d
, 1);
2505 res
= eval_polynomial(e
->x
.p
, 0, ¶m
, values
);
2506 free_evalue_refs(¶m
);
2509 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2510 mpz_fdiv_r(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2512 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2513 evalue_free(param2
);
2516 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2517 mpz_fdiv_q(param2
->x
.n
, param2
->x
.n
, param2
->d
);
2518 value_set_si(param2
->d
, 1);
2520 res
= eval_polynomial(e
->x
.p
, 1, param2
, values
);
2521 evalue_free(param2
);
2524 param2
= evalue_eval(&e
->x
.p
->arr
[0], values
);
2525 if (value_zero_p(param2
->x
.n
))
2526 res
= evalue_eval(&e
->x
.p
->arr
[1], values
);
2527 else if (e
->x
.p
->size
> 2)
2528 res
= evalue_eval(&e
->x
.p
->arr
[2], values
);
2530 res
= evalue_zero();
2531 evalue_free(param2
);
2534 assert(e
->x
.p
->pos
== EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
);
2535 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2536 if (in_domain(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), values
)) {
2537 res
= evalue_eval(&e
->x
.p
->arr
[2*i
+1], values
);
2541 res
= evalue_zero();
2549 size_t value_size(Value v
) {
2550 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2551 * sizeof(v
[0]._mp_d
[0]);
2554 size_t domain_size(Polyhedron
*D
)
2557 size_t s
= sizeof(*D
);
2559 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2560 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2561 s
+= value_size(D
->Constraint
[i
][j
]);
2564 for (i = 0; i < D->NbRays; ++i)
2565 for (j = 0; j < D->Dimension+2; ++j)
2566 s += value_size(D->Ray[i][j]);
2569 return D
->next
? s
+domain_size(D
->next
) : s
;
2572 size_t enode_size(enode
*p
) {
2573 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2576 if (p
->type
== partition
)
2577 for (i
= 0; i
< p
->size
/2; ++i
) {
2578 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2579 s
+= evalue_size(&p
->arr
[2*i
+1]);
2582 for (i
= 0; i
< p
->size
; ++i
) {
2583 s
+= evalue_size(&p
->arr
[i
]);
2588 size_t evalue_size(evalue
*e
)
2590 size_t s
= sizeof(*e
);
2591 s
+= value_size(e
->d
);
2592 if (value_notzero_p(e
->d
))
2593 s
+= value_size(e
->x
.n
);
2595 s
+= enode_size(e
->x
.p
);
2599 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2601 evalue
*found
= NULL
;
2606 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2609 value_init(offset
.d
);
2610 value_init(offset
.x
.n
);
2611 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2612 value_lcm(offset
.d
, m
, offset
.d
);
2613 value_set_si(offset
.x
.n
, 1);
2616 evalue_copy(©
, cst
);
2619 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2621 if (eequal(base
, &e
->x
.p
->arr
[0]))
2622 found
= &e
->x
.p
->arr
[0];
2624 value_set_si(offset
.x
.n
, -2);
2627 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2629 if (eequal(base
, &e
->x
.p
->arr
[0]))
2632 free_evalue_refs(cst
);
2633 free_evalue_refs(&offset
);
2636 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2637 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2642 static evalue
*find_relation_pair(evalue
*e
)
2645 evalue
*found
= NULL
;
2647 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2650 if (e
->x
.p
->type
== fractional
) {
2655 poly_denom(&e
->x
.p
->arr
[0], &m
);
2657 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2658 cst
= &cst
->x
.p
->arr
[0])
2661 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2662 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2667 i
= e
->x
.p
->type
== relation
;
2668 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2669 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2674 void evalue_mod2relation(evalue
*e
) {
2677 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2680 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2681 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2682 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2683 value_clear(e
->x
.p
->arr
[2*i
].d
);
2684 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2686 if (2*i
< e
->x
.p
->size
) {
2687 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2688 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2693 if (e
->x
.p
->size
== 0) {
2695 evalue_set_si(e
, 0, 1);
2701 while ((d
= find_relation_pair(e
)) != NULL
) {
2705 value_init(split
.d
);
2706 value_set_si(split
.d
, 0);
2707 split
.x
.p
= new_enode(relation
, 3, 0);
2708 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2709 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2711 ev
= &split
.x
.p
->arr
[0];
2712 value_set_si(ev
->d
, 0);
2713 ev
->x
.p
= new_enode(fractional
, 3, -1);
2714 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2715 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2716 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2722 free_evalue_refs(&split
);
2726 static int evalue_comp(const void * a
, const void * b
)
2728 const evalue
*e1
= *(const evalue
**)a
;
2729 const evalue
*e2
= *(const evalue
**)b
;
2730 return ecmp(e1
, e2
);
2733 void evalue_combine(evalue
*e
)
2740 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2743 NALLOC(evs
, e
->x
.p
->size
/2);
2744 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2745 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2746 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2747 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2748 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2749 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2750 value_clear(p
->arr
[2*k
].d
);
2751 value_clear(p
->arr
[2*k
+1].d
);
2752 p
->arr
[2*k
] = *(evs
[i
]-1);
2753 p
->arr
[2*k
+1] = *(evs
[i
]);
2756 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2759 value_clear((evs
[i
]-1)->d
);
2763 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2764 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2765 free_evalue_refs(evs
[i
]);
2769 for (i
= 2*k
; i
< p
->size
; ++i
)
2770 value_clear(p
->arr
[i
].d
);
2777 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2779 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2781 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2784 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2785 Polyhedron
*D
, *N
, **P
;
2788 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2795 if (D
->NbEq
<= H
->NbEq
) {
2801 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2802 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2803 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2804 reduce_evalue(&tmp
);
2805 if (value_notzero_p(tmp
.d
) ||
2806 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2809 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2810 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2813 free_evalue_refs(&tmp
);
2819 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2821 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2823 value_clear(e
->x
.p
->arr
[2*i
].d
);
2824 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2826 if (2*i
< e
->x
.p
->size
) {
2827 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2828 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2835 H
= DomainConvex(D
, 0);
2836 E
= DomainDifference(H
, D
, 0);
2838 D
= DomainDifference(H
, E
, 0);
2841 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2845 /* Use smallest representative for coefficients in affine form in
2846 * argument of fractional.
2847 * Since any change will make the argument non-standard,
2848 * the containing evalue will have to be reduced again afterward.
2850 static void fractional_minimal_coefficients(enode
*p
)
2856 assert(p
->type
== fractional
);
2858 while (value_zero_p(pp
->d
)) {
2859 assert(pp
->x
.p
->type
== polynomial
);
2860 assert(pp
->x
.p
->size
== 2);
2861 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2862 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2863 if (value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2864 value_subtract(pp
->x
.p
->arr
[1].x
.n
,
2865 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2866 pp
= &pp
->x
.p
->arr
[0];
2872 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2877 unsigned dim
= D
->Dimension
;
2878 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2881 assert(p
->type
== fractional
|| p
->type
== flooring
);
2882 value_set_si(T
->p
[1][dim
], 1);
2883 evalue_extract_affine(&p
->arr
[0], T
->p
[0], &T
->p
[0][dim
], d
);
2884 I
= DomainImage(D
, T
, 0);
2885 H
= DomainConvex(I
, 0);
2895 static void replace_by_affine(evalue
*e
, Value offset
)
2902 value_init(inc
.x
.n
);
2903 value_set_si(inc
.d
, 1);
2904 value_oppose(inc
.x
.n
, offset
);
2905 eadd(&inc
, &p
->arr
[0]);
2906 reorder_terms_about(p
, &p
->arr
[0]); /* frees arr[0] */
2910 free_evalue_refs(&inc
);
2913 int evalue_range_reduction_in_domain(evalue
*e
, Polyhedron
*D
)
2922 if (value_notzero_p(e
->d
))
2927 if (p
->type
== relation
) {
2934 fractional_minimal_coefficients(p
->arr
[0].x
.p
);
2935 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
);
2936 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2937 equal
= value_eq(min
, max
);
2938 mpz_cdiv_q(min
, min
, d
);
2939 mpz_fdiv_q(max
, max
, d
);
2941 if (bounded
&& value_gt(min
, max
)) {
2947 evalue_set_si(e
, 0, 1);
2950 free_evalue_refs(&(p
->arr
[1]));
2951 free_evalue_refs(&(p
->arr
[0]));
2957 return r
? r
: evalue_range_reduction_in_domain(e
, D
);
2958 } else if (bounded
&& equal
) {
2961 free_evalue_refs(&(p
->arr
[2]));
2964 free_evalue_refs(&(p
->arr
[0]));
2970 return evalue_range_reduction_in_domain(e
, D
);
2971 } else if (bounded
&& value_eq(min
, max
)) {
2972 /* zero for a single value */
2974 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2975 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2976 value_multiply(min
, min
, d
);
2977 value_subtract(M
->p
[0][D
->Dimension
+1],
2978 M
->p
[0][D
->Dimension
+1], min
);
2979 E
= DomainAddConstraints(D
, M
, 0);
2985 r
= evalue_range_reduction_in_domain(&p
->arr
[1], E
);
2987 r
|= evalue_range_reduction_in_domain(&p
->arr
[2], D
);
2989 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2997 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
3000 i
= p
->type
== relation
? 1 :
3001 p
->type
== fractional
? 1 : 0;
3002 for (; i
<p
->size
; i
++)
3003 r
|= evalue_range_reduction_in_domain(&p
->arr
[i
], D
);
3005 if (p
->type
!= fractional
) {
3006 if (r
&& p
->type
== polynomial
) {
3009 value_set_si(f
.d
, 0);
3010 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3011 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3012 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3013 reorder_terms_about(p
, &f
);
3024 fractional_minimal_coefficients(p
);
3025 I
= polynomial_projection(p
, D
, &d
, NULL
);
3026 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
3027 mpz_fdiv_q(min
, min
, d
);
3028 mpz_fdiv_q(max
, max
, d
);
3029 value_subtract(d
, max
, min
);
3031 if (bounded
&& value_eq(min
, max
)) {
3032 replace_by_affine(e
, min
);
3034 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
3035 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
3036 * See pages 199-200 of PhD thesis.
3044 value_set_si(rem
.d
, 0);
3045 rem
.x
.p
= new_enode(fractional
, 3, -1);
3046 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
3047 value_clear(rem
.x
.p
->arr
[1].d
);
3048 value_clear(rem
.x
.p
->arr
[2].d
);
3049 rem
.x
.p
->arr
[1] = p
->arr
[1];
3050 rem
.x
.p
->arr
[2] = p
->arr
[2];
3051 for (i
= 3; i
< p
->size
; ++i
)
3052 p
->arr
[i
-2] = p
->arr
[i
];
3056 value_init(inc
.x
.n
);
3057 value_set_si(inc
.d
, 1);
3058 value_oppose(inc
.x
.n
, min
);
3061 evalue_copy(&t
, &p
->arr
[0]);
3065 value_set_si(f
.d
, 0);
3066 f
.x
.p
= new_enode(fractional
, 3, -1);
3067 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3068 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3069 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
3071 value_init(factor
.d
);
3072 evalue_set_si(&factor
, -1, 1);
3078 value_clear(f
.x
.p
->arr
[1].x
.n
);
3079 value_clear(f
.x
.p
->arr
[2].x
.n
);
3080 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3081 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3085 reorder_terms(&rem
);
3092 free_evalue_refs(&inc
);
3093 free_evalue_refs(&t
);
3094 free_evalue_refs(&f
);
3095 free_evalue_refs(&factor
);
3096 free_evalue_refs(&rem
);
3098 evalue_range_reduction_in_domain(e
, D
);
3102 _reduce_evalue(&p
->arr
[0], 0, 1);
3114 void evalue_range_reduction(evalue
*e
)
3117 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3120 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3121 if (evalue_range_reduction_in_domain(&e
->x
.p
->arr
[2*i
+1],
3122 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
3123 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3125 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
3126 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
3127 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3128 value_clear(e
->x
.p
->arr
[2*i
].d
);
3130 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
3131 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
3139 Enumeration
* partition2enumeration(evalue
*EP
)
3142 Enumeration
*en
, *res
= NULL
;
3144 if (EVALUE_IS_ZERO(*EP
)) {
3149 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
3150 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
3151 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
3154 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
3155 value_clear(EP
->x
.p
->arr
[2*i
].d
);
3156 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
3164 int evalue_frac2floor_in_domain3(evalue
*e
, Polyhedron
*D
, int shift
)
3173 if (value_notzero_p(e
->d
))
3178 i
= p
->type
== relation
? 1 :
3179 p
->type
== fractional
? 1 : 0;
3180 for (; i
<p
->size
; i
++)
3181 r
|= evalue_frac2floor_in_domain3(&p
->arr
[i
], D
, shift
);
3183 if (p
->type
!= fractional
) {
3184 if (r
&& p
->type
== polynomial
) {
3187 value_set_si(f
.d
, 0);
3188 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
3189 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
3190 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
3191 reorder_terms_about(p
, &f
);
3201 I
= polynomial_projection(p
, D
, &d
, NULL
);
3204 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3207 assert(I
->NbEq
== 0); /* Should have been reduced */
3210 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3211 if (value_pos_p(I
->Constraint
[i
][1]))
3214 if (i
< I
->NbConstraints
) {
3216 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3217 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3218 if (value_neg_p(min
)) {
3220 mpz_fdiv_q(min
, min
, d
);
3221 value_init(offset
.d
);
3222 value_set_si(offset
.d
, 1);
3223 value_init(offset
.x
.n
);
3224 value_oppose(offset
.x
.n
, min
);
3225 eadd(&offset
, &p
->arr
[0]);
3226 free_evalue_refs(&offset
);
3236 value_set_si(fl
.d
, 0);
3237 fl
.x
.p
= new_enode(flooring
, 3, -1);
3238 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
3239 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
3240 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
3242 eadd(&fl
, &p
->arr
[0]);
3243 reorder_terms_about(p
, &p
->arr
[0]);
3247 free_evalue_refs(&fl
);
3252 int evalue_frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
3254 return evalue_frac2floor_in_domain3(e
, D
, 1);
3257 void evalue_frac2floor2(evalue
*e
, int shift
)
3260 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
3262 if (evalue_frac2floor_in_domain3(e
, NULL
, 0))
3268 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3269 if (evalue_frac2floor_in_domain3(&e
->x
.p
->arr
[2*i
+1],
3270 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), shift
))
3271 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
3274 void evalue_frac2floor(evalue
*e
)
3276 evalue_frac2floor2(e
, 1);
3279 /* Add a new variable with lower bound 1 and upper bound specified
3280 * by row. If negative is true, then the new variable has upper
3281 * bound -1 and lower bound specified by row.
3283 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
3284 Vector
*row
, int negative
)
3288 int nparam
= D
->Dimension
- nvar
;
3291 nr
= D
->NbConstraints
+ 2;
3292 nc
= D
->Dimension
+ 2 + 1;
3293 C
= Matrix_Alloc(nr
, nc
);
3294 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
3295 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
3296 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3297 D
->Dimension
+ 1 - nvar
);
3302 nc
= C
->NbColumns
+ 1;
3303 C
= Matrix_Alloc(nr
, nc
);
3304 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
3305 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
3306 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3307 oldC
->NbColumns
- 1 - nvar
);
3310 value_set_si(C
->p
[nr
-2][0], 1);
3312 value_set_si(C
->p
[nr
-2][1 + nvar
], -1);
3314 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
3315 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
3317 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
3318 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
3324 static void floor2frac_r(evalue
*e
, int nvar
)
3331 if (value_notzero_p(e
->d
))
3336 assert(p
->type
== flooring
);
3337 for (i
= 1; i
< p
->size
; i
++)
3338 floor2frac_r(&p
->arr
[i
], nvar
);
3340 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
3341 assert(pp
->x
.p
->type
== polynomial
);
3342 pp
->x
.p
->pos
-= nvar
;
3346 value_set_si(f
.d
, 0);
3347 f
.x
.p
= new_enode(fractional
, 3, -1);
3348 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3349 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3350 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3352 eadd(&f
, &p
->arr
[0]);
3353 reorder_terms_about(p
, &p
->arr
[0]);
3357 free_evalue_refs(&f
);
3360 /* Convert flooring back to fractional and shift position
3361 * of the parameters by nvar
3363 static void floor2frac(evalue
*e
, int nvar
)
3365 floor2frac_r(e
, nvar
);
3369 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
3372 int nparam
= D
->Dimension
- nvar
;
3376 D
= Constraints2Polyhedron(C
, 0);
3380 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
3382 /* Double check that D was not unbounded. */
3383 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
3391 static evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
3392 int *signs
, Matrix
*C
, unsigned MaxRays
)
3398 evalue
*factor
= NULL
;
3402 if (EVALUE_IS_ZERO(*e
))
3406 Polyhedron
*DD
= Disjoint_Domain(D
, 0, MaxRays
);
3413 res
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3416 for (Q
= DD
; Q
; Q
= DD
) {
3422 t
= esum_over_domain(e
, nvar
, Q
, signs
, C
, MaxRays
);
3435 if (value_notzero_p(e
->d
)) {
3438 t
= esum_over_domain_cst(nvar
, D
, C
);
3440 if (!EVALUE_IS_ONE(*e
))
3446 switch (e
->x
.p
->type
) {
3448 evalue
*pp
= &e
->x
.p
->arr
[0];
3450 if (pp
->x
.p
->pos
> nvar
) {
3451 /* remainder is independent of the summated vars */
3457 floor2frac(&f
, nvar
);
3459 t
= esum_over_domain_cst(nvar
, D
, C
);
3463 free_evalue_refs(&f
);
3468 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3469 poly_denom(pp
, &row
->p
[1 + nvar
]);
3470 value_set_si(row
->p
[0], 1);
3471 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3472 pp
= &pp
->x
.p
->arr
[0]) {
3474 assert(pp
->x
.p
->type
== polynomial
);
3476 if (pos
>= 1 + nvar
)
3478 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3479 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3480 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3482 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3483 value_division(row
->p
[1 + D
->Dimension
+ 1],
3484 row
->p
[1 + D
->Dimension
+ 1],
3486 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3487 row
->p
[1 + D
->Dimension
+ 1],
3489 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3493 int pos
= e
->x
.p
->pos
;
3496 factor
= ALLOC(evalue
);
3497 value_init(factor
->d
);
3498 value_set_si(factor
->d
, 0);
3499 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3500 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3501 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3505 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3506 negative
= signs
[pos
-1] < 0;
3507 value_set_si(row
->p
[0], 1);
3509 value_set_si(row
->p
[pos
], -1);
3510 value_set_si(row
->p
[1 + nvar
], 1);
3512 value_set_si(row
->p
[pos
], 1);
3513 value_set_si(row
->p
[1 + nvar
], -1);
3521 offset
= type_offset(e
->x
.p
);
3523 res
= esum_over_domain(&e
->x
.p
->arr
[offset
], nvar
, D
, signs
, C
, MaxRays
);
3527 evalue_copy(&cum
, factor
);
3531 for (i
= 1; offset
+i
< e
->x
.p
->size
; ++i
) {
3535 C
= esum_add_constraint(nvar
, D
, C
, row
, negative
);
3541 Vector_Print(stderr, P_VALUE_FMT, row);
3543 Matrix_Print(stderr, P_VALUE_FMT, C);
3545 t
= esum_over_domain(&e
->x
.p
->arr
[offset
+i
], nvar
, D
, signs
, C
, MaxRays
);
3550 if (negative
&& (i
% 2))
3560 if (factor
&& offset
+i
+1 < e
->x
.p
->size
)
3567 free_evalue_refs(&cum
);
3568 evalue_free(factor
);
3579 static void domain_signs(Polyhedron
*D
, int *signs
)
3583 POL_ENSURE_VERTICES(D
);
3584 for (j
= 0; j
< D
->Dimension
; ++j
) {
3586 for (k
= 0; k
< D
->NbRays
; ++k
) {
3587 signs
[j
] = value_sign(D
->Ray
[k
][1+j
]);
3594 static void shift_floor_in_domain(evalue
*e
, Polyhedron
*D
)
3601 if (value_notzero_p(e
->d
))
3606 for (i
= type_offset(p
); i
< p
->size
; ++i
)
3607 shift_floor_in_domain(&p
->arr
[i
], D
);
3609 if (p
->type
!= flooring
)
3615 I
= polynomial_projection(p
, D
, &d
, NULL
);
3616 assert(I
->NbEq
== 0); /* Should have been reduced */
3618 for (i
= 0; i
< I
->NbConstraints
; ++i
)
3619 if (value_pos_p(I
->Constraint
[i
][1]))
3621 assert(i
< I
->NbConstraints
);
3622 if (i
< I
->NbConstraints
) {
3623 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
3624 mpz_fdiv_q(m
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
3625 if (value_neg_p(m
)) {
3626 /* replace [e] by [e-m]+m such that e-m >= 0 */
3631 value_set_si(f
.d
, 1);
3632 value_oppose(f
.x
.n
, m
);
3633 eadd(&f
, &p
->arr
[0]);
3636 value_set_si(f
.d
, 0);
3637 f
.x
.p
= new_enode(flooring
, 3, -1);
3638 value_clear(f
.x
.p
->arr
[0].d
);
3639 f
.x
.p
->arr
[0] = p
->arr
[0];
3640 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
3641 value_set_si(f
.x
.p
->arr
[1].d
, 1);
3642 value_init(f
.x
.p
->arr
[1].x
.n
);
3643 value_assign(f
.x
.p
->arr
[1].x
.n
, m
);
3644 reorder_terms_about(p
, &f
);
3655 /* Make arguments of all floors non-negative */
3656 static void shift_floor_arguments(evalue
*e
)
3660 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
3663 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3664 shift_floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
3665 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
3668 evalue
*evalue_sum(evalue
*e
, int nvar
, unsigned MaxRays
)
3672 evalue
*res
= ALLOC(evalue
);
3676 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3677 evalue_copy(res
, e
);
3681 evalue_split_domains_into_orthants(e
, MaxRays
);
3682 evalue_frac2floor2(e
, 0);
3683 evalue_set_si(res
, 0, 1);
3685 assert(value_zero_p(e
->d
));
3686 assert(e
->x
.p
->type
== partition
);
3687 shift_floor_arguments(e
);
3689 assert(e
->x
.p
->size
>= 2);
3690 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3692 signs
= alloca(sizeof(int) * dim
);
3694 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3696 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3697 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3698 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
, 0,
3709 evalue
*esum(evalue
*e
, int nvar
)
3711 return evalue_sum(e
, nvar
, 0);
3714 /* Initial silly implementation */
3715 void eor(evalue
*e1
, evalue
*res
)
3721 evalue_set_si(&mone
, -1, 1);
3723 evalue_copy(&E
, res
);
3729 free_evalue_refs(&E
);
3730 free_evalue_refs(&mone
);
3733 /* computes denominator of polynomial evalue
3734 * d should point to a value initialized to 1
3736 void evalue_denom(const evalue
*e
, Value
*d
)
3740 if (value_notzero_p(e
->d
)) {
3741 value_lcm(*d
, *d
, e
->d
);
3744 assert(e
->x
.p
->type
== polynomial
||
3745 e
->x
.p
->type
== fractional
||
3746 e
->x
.p
->type
== flooring
);
3747 offset
= type_offset(e
->x
.p
);
3748 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3749 evalue_denom(&e
->x
.p
->arr
[i
], d
);
3752 /* Divides the evalue e by the integer n */
3753 void evalue_div(evalue
*e
, Value n
)
3757 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3760 if (value_notzero_p(e
->d
)) {
3763 value_multiply(e
->d
, e
->d
, n
);
3764 value_gcd(gc
, e
->x
.n
, e
->d
);
3765 if (value_notone_p(gc
)) {
3766 value_division(e
->d
, e
->d
, gc
);
3767 value_division(e
->x
.n
, e
->x
.n
, gc
);
3772 if (e
->x
.p
->type
== partition
) {
3773 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3774 evalue_div(&e
->x
.p
->arr
[2*i
+1], n
);
3777 offset
= type_offset(e
->x
.p
);
3778 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3779 evalue_div(&e
->x
.p
->arr
[i
], n
);
3782 /* Multiplies the evalue e by the integer n */
3783 void evalue_mul(evalue
*e
, Value n
)
3787 if (value_one_p(n
) || EVALUE_IS_ZERO(*e
))
3790 if (value_notzero_p(e
->d
)) {
3793 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3794 value_gcd(gc
, e
->x
.n
, e
->d
);
3795 if (value_notone_p(gc
)) {
3796 value_division(e
->d
, e
->d
, gc
);
3797 value_division(e
->x
.n
, e
->x
.n
, gc
);
3802 if (e
->x
.p
->type
== partition
) {
3803 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3804 evalue_mul(&e
->x
.p
->arr
[2*i
+1], n
);
3807 offset
= type_offset(e
->x
.p
);
3808 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3809 evalue_mul(&e
->x
.p
->arr
[i
], n
);
3812 /* Multiplies the evalue e by the n/d */
3813 void evalue_mul_div(evalue
*e
, Value n
, Value d
)
3817 if ((value_one_p(n
) && value_one_p(d
)) || EVALUE_IS_ZERO(*e
))
3820 if (value_notzero_p(e
->d
)) {
3823 value_multiply(e
->x
.n
, e
->x
.n
, n
);
3824 value_multiply(e
->d
, e
->d
, d
);
3825 value_gcd(gc
, e
->x
.n
, e
->d
);
3826 if (value_notone_p(gc
)) {
3827 value_division(e
->d
, e
->d
, gc
);
3828 value_division(e
->x
.n
, e
->x
.n
, gc
);
3833 if (e
->x
.p
->type
== partition
) {
3834 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3835 evalue_mul_div(&e
->x
.p
->arr
[2*i
+1], n
, d
);
3838 offset
= type_offset(e
->x
.p
);
3839 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3840 evalue_mul_div(&e
->x
.p
->arr
[i
], n
, d
);
3843 void evalue_negate(evalue
*e
)
3847 if (value_notzero_p(e
->d
)) {
3848 value_oppose(e
->x
.n
, e
->x
.n
);
3851 if (e
->x
.p
->type
== partition
) {
3852 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3853 evalue_negate(&e
->x
.p
->arr
[2*i
+1]);
3856 offset
= type_offset(e
->x
.p
);
3857 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
3858 evalue_negate(&e
->x
.p
->arr
[i
]);
3861 void evalue_add_constant(evalue
*e
, const Value cst
)
3865 if (value_zero_p(e
->d
)) {
3866 if (e
->x
.p
->type
== partition
) {
3867 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
3868 evalue_add_constant(&e
->x
.p
->arr
[2*i
+1], cst
);
3871 if (e
->x
.p
->type
== relation
) {
3872 for (i
= 1; i
< e
->x
.p
->size
; ++i
)
3873 evalue_add_constant(&e
->x
.p
->arr
[i
], cst
);
3877 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
3878 } while (value_zero_p(e
->d
));
3880 value_addmul(e
->x
.n
, cst
, e
->d
);
3883 static void evalue_frac2polynomial_r(evalue
*e
, int *signs
, int sign
, int in_frac
)
3888 int sign_odd
= sign
;
3890 if (value_notzero_p(e
->d
)) {
3891 if (in_frac
&& sign
* value_sign(e
->x
.n
) < 0) {
3892 value_set_si(e
->x
.n
, 0);
3893 value_set_si(e
->d
, 1);
3898 if (e
->x
.p
->type
== relation
) {
3899 for (i
= e
->x
.p
->size
-1; i
>= 1; --i
)
3900 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
, sign
, in_frac
);
3904 if (e
->x
.p
->type
== polynomial
)
3905 sign_odd
*= signs
[e
->x
.p
->pos
-1];
3906 offset
= type_offset(e
->x
.p
);
3907 evalue_frac2polynomial_r(&e
->x
.p
->arr
[offset
], signs
, sign
, in_frac
);
3908 in_frac
|= e
->x
.p
->type
== fractional
;
3909 for (i
= e
->x
.p
->size
-1; i
> offset
; --i
)
3910 evalue_frac2polynomial_r(&e
->x
.p
->arr
[i
], signs
,
3911 (i
- offset
) % 2 ? sign_odd
: sign
, in_frac
);
3913 if (e
->x
.p
->type
!= fractional
)
3916 /* replace { a/m } by (m-1)/m if sign != 0
3917 * and by (m-1)/(2m) if sign == 0
3921 evalue_denom(&e
->x
.p
->arr
[0], &d
);
3922 free_evalue_refs(&e
->x
.p
->arr
[0]);
3923 value_init(e
->x
.p
->arr
[0].d
);
3924 value_init(e
->x
.p
->arr
[0].x
.n
);
3926 value_addto(e
->x
.p
->arr
[0].d
, d
, d
);
3928 value_assign(e
->x
.p
->arr
[0].d
, d
);
3929 value_decrement(e
->x
.p
->arr
[0].x
.n
, d
);
3933 reorder_terms_about(p
, &p
->arr
[0]);
3939 /* Approximate the evalue in fractional representation by a polynomial.
3940 * If sign > 0, the result is an upper bound;
3941 * if sign < 0, the result is a lower bound;
3942 * if sign = 0, the result is an intermediate approximation.
3944 void evalue_frac2polynomial(evalue
*e
, int sign
, unsigned MaxRays
)
3949 if (value_notzero_p(e
->d
))
3951 assert(e
->x
.p
->type
== partition
);
3952 /* make sure all variables in the domains have a fixed sign */
3954 evalue_split_domains_into_orthants(e
, MaxRays
);
3955 if (EVALUE_IS_ZERO(*e
))
3959 assert(e
->x
.p
->size
>= 2);
3960 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3962 signs
= alloca(sizeof(int) * dim
);
3965 for (i
= 0; i
< dim
; ++i
)
3967 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3969 domain_signs(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), signs
);
3970 evalue_frac2polynomial_r(&e
->x
.p
->arr
[2*i
+1], signs
, sign
, 0);
3974 /* Split the domains of e (which is assumed to be a partition)
3975 * such that each resulting domain lies entirely in one orthant.
3977 void evalue_split_domains_into_orthants(evalue
*e
, unsigned MaxRays
)
3980 assert(value_zero_p(e
->d
));
3981 assert(e
->x
.p
->type
== partition
);
3982 assert(e
->x
.p
->size
>= 2);
3983 dim
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
;
3985 for (i
= 0; i
< dim
; ++i
) {
3988 C
= Matrix_Alloc(1, 1 + dim
+ 1);
3989 value_set_si(C
->p
[0][0], 1);
3990 value_init(split
.d
);
3991 value_set_si(split
.d
, 0);
3992 split
.x
.p
= new_enode(partition
, 4, dim
);
3993 value_set_si(C
->p
[0][1+i
], 1);
3994 C2
= Matrix_Copy(C
);
3995 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[0], Constraints2Polyhedron(C2
, MaxRays
));
3997 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
3998 value_set_si(C
->p
[0][1+i
], -1);
3999 value_set_si(C
->p
[0][1+dim
], -1);
4000 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[2], Constraints2Polyhedron(C
, MaxRays
));
4001 evalue_set_si(&split
.x
.p
->arr
[3], 1, 1);
4003 free_evalue_refs(&split
);
4009 static evalue
*find_fractional_with_max_periods(evalue
*e
, Polyhedron
*D
,
4012 Value
*min
, Value
*max
)
4019 if (value_notzero_p(e
->d
))
4022 if (e
->x
.p
->type
== fractional
) {
4027 I
= polynomial_projection(e
->x
.p
, D
, &d
, &T
);
4028 bounded
= line_minmax(I
, min
, max
); /* frees I */
4032 value_set_si(mp
, max_periods
);
4033 mpz_fdiv_q(*min
, *min
, d
);
4034 mpz_fdiv_q(*max
, *max
, d
);
4035 value_assign(T
->p
[1][D
->Dimension
], d
);
4036 value_subtract(d
, *max
, *min
);
4037 if (value_ge(d
, mp
))
4040 f
= evalue_dup(&e
->x
.p
->arr
[0]);
4051 for (i
= type_offset(e
->x
.p
); i
< e
->x
.p
->size
; ++i
)
4052 if ((f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[i
], D
, max_periods
,
4059 static void replace_fract_by_affine(evalue
*e
, evalue
*f
, Value val
)
4063 if (value_notzero_p(e
->d
))
4066 offset
= type_offset(e
->x
.p
);
4067 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
4068 replace_fract_by_affine(&e
->x
.p
->arr
[i
], f
, val
);
4070 if (e
->x
.p
->type
!= fractional
)
4073 if (!eequal(&e
->x
.p
->arr
[0], f
))
4076 replace_by_affine(e
, val
);
4079 /* Look for fractional parts that can be removed by splitting the corresponding
4080 * domain into at most max_periods parts.
4081 * We use a very simply strategy that looks for the first fractional part
4082 * that satisfies the condition, performs the split and then continues
4083 * looking for other fractional parts in the split domains until no
4084 * such fractional part can be found anymore.
4086 void evalue_split_periods(evalue
*e
, int max_periods
, unsigned int MaxRays
)
4093 if (EVALUE_IS_ZERO(*e
))
4095 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
) {
4097 "WARNING: evalue_split_periods called on incorrect evalue type\n");
4105 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
4110 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
4112 f
= find_fractional_with_max_periods(&e
->x
.p
->arr
[2*i
+1], D
, max_periods
,
4117 M
= Matrix_Alloc(2, 2+D
->Dimension
);
4119 value_subtract(d
, max
, min
);
4120 n
= VALUE_TO_INT(d
)+1;
4122 value_set_si(M
->p
[0][0], 1);
4123 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
4124 value_multiply(d
, max
, T
->p
[1][D
->Dimension
]);
4125 value_subtract(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
], d
);
4126 value_set_si(d
, -1);
4127 value_set_si(M
->p
[1][0], 1);
4128 Vector_Scale(T
->p
[0], M
->p
[1]+1, d
, D
->Dimension
+1);
4129 value_addmul(M
->p
[1][1+D
->Dimension
], max
, T
->p
[1][D
->Dimension
]);
4130 value_addto(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4131 T
->p
[1][D
->Dimension
]);
4132 value_decrement(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
]);
4134 p
= new_enode(partition
, e
->x
.p
->size
+ (n
-1)*2, e
->x
.p
->pos
);
4135 for (j
= 0; j
< 2*i
; ++j
) {
4136 value_clear(p
->arr
[j
].d
);
4137 p
->arr
[j
] = e
->x
.p
->arr
[j
];
4139 for (j
= 2*i
+2; j
< e
->x
.p
->size
; ++j
) {
4140 value_clear(p
->arr
[j
+2*(n
-1)].d
);
4141 p
->arr
[j
+2*(n
-1)] = e
->x
.p
->arr
[j
];
4143 for (j
= n
-1; j
>= 0; --j
) {
4145 value_clear(p
->arr
[2*i
+1].d
);
4146 p
->arr
[2*i
+1] = e
->x
.p
->arr
[2*i
+1];
4148 evalue_copy(&p
->arr
[2*(i
+j
)+1], &e
->x
.p
->arr
[2*i
+1]);
4150 value_subtract(M
->p
[1][1+D
->Dimension
], M
->p
[1][1+D
->Dimension
],
4151 T
->p
[1][D
->Dimension
]);
4152 value_addto(M
->p
[0][1+D
->Dimension
], M
->p
[0][1+D
->Dimension
],
4153 T
->p
[1][D
->Dimension
]);
4155 replace_fract_by_affine(&p
->arr
[2*(i
+j
)+1], f
, max
);
4156 E
= DomainAddConstraints(D
, M
, MaxRays
);
4157 EVALUE_SET_DOMAIN(p
->arr
[2*(i
+j
)], E
);
4158 if (evalue_range_reduction_in_domain(&p
->arr
[2*(i
+j
)+1], E
))
4159 reduce_evalue(&p
->arr
[2*(i
+j
)+1]);
4160 value_decrement(max
, max
);
4162 value_clear(e
->x
.p
->arr
[2*i
].d
);
4177 void evalue_extract_affine(const evalue
*e
, Value
*coeff
, Value
*cst
, Value
*d
)
4179 value_set_si(*d
, 1);
4181 for ( ; value_zero_p(e
->d
); e
= &e
->x
.p
->arr
[0]) {
4183 assert(e
->x
.p
->type
== polynomial
);
4184 assert(e
->x
.p
->size
== 2);
4185 c
= &e
->x
.p
->arr
[1];
4186 value_multiply(coeff
[e
->x
.p
->pos
-1], *d
, c
->x
.n
);
4187 value_division(coeff
[e
->x
.p
->pos
-1], coeff
[e
->x
.p
->pos
-1], c
->d
);
4189 value_multiply(*cst
, *d
, e
->x
.n
);
4190 value_division(*cst
, *cst
, e
->d
);
4193 /* returns an evalue that corresponds to
4197 static evalue
*term(int param
, Value c
, Value den
)
4199 evalue
*EP
= ALLOC(evalue
);
4201 value_set_si(EP
->d
,0);
4202 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
4203 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
4204 value_init(EP
->x
.p
->arr
[1].x
.n
);
4205 value_assign(EP
->x
.p
->arr
[1].d
, den
);
4206 value_assign(EP
->x
.p
->arr
[1].x
.n
, c
);
4210 evalue
*affine2evalue(Value
*coeff
, Value denom
, int nvar
)
4213 evalue
*E
= ALLOC(evalue
);
4215 evalue_set(E
, coeff
[nvar
], denom
);
4216 for (i
= 0; i
< nvar
; ++i
) {
4218 if (value_zero_p(coeff
[i
]))
4220 t
= term(i
, coeff
[i
], denom
);
4227 void evalue_substitute(evalue
*e
, evalue
**subs
)
4233 if (value_notzero_p(e
->d
))
4237 assert(p
->type
!= partition
);
4239 for (i
= 0; i
< p
->size
; ++i
)
4240 evalue_substitute(&p
->arr
[i
], subs
);
4242 if (p
->type
== polynomial
)
4247 value_set_si(v
->d
, 0);
4248 v
->x
.p
= new_enode(p
->type
, 3, -1);
4249 value_clear(v
->x
.p
->arr
[0].d
);
4250 v
->x
.p
->arr
[0] = p
->arr
[0];
4251 evalue_set_si(&v
->x
.p
->arr
[1], 0, 1);
4252 evalue_set_si(&v
->x
.p
->arr
[2], 1, 1);
4255 offset
= type_offset(p
);
4257 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
4258 emul(v
, &p
->arr
[i
]);
4259 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
4260 free_evalue_refs(&(p
->arr
[i
]));
4263 if (p
->type
!= polynomial
)
4267 *e
= p
->arr
[offset
];
4271 /* evalue e is given in terms of "new" parameter; CP maps the new
4272 * parameters back to the old parameters.
4273 * Transforms e such that it refers back to the old parameters.
4275 void evalue_backsubstitute(evalue
*e
, Matrix
*CP
, unsigned MaxRays
)
4282 unsigned nparam
= CP
->NbColumns
-1;
4285 if (EVALUE_IS_ZERO(*e
))
4288 assert(value_zero_p(e
->d
));
4290 assert(p
->type
== partition
);
4292 inv
= left_inverse(CP
, &eq
);
4293 subs
= ALLOCN(evalue
*, nparam
);
4294 for (i
= 0; i
< nparam
; ++i
)
4295 subs
[i
] = affine2evalue(inv
->p
[i
], inv
->p
[nparam
][inv
->NbColumns
-1],
4298 CEq
= Constraints2Polyhedron(eq
, MaxRays
);
4299 addeliminatedparams_partition(p
, inv
, CEq
, inv
->NbColumns
-1, MaxRays
);
4300 Polyhedron_Free(CEq
);
4302 for (i
= 0; i
< p
->size
/2; ++i
)
4303 evalue_substitute(&p
->arr
[2*i
+1], subs
);
4305 for (i
= 0; i
< nparam
; ++i
)
4306 evalue_free(subs
[i
]);
4314 * \sum_{i=0}^n c_i/d X^i
4316 * where d is the last element in the vector c.
4318 evalue
*evalue_polynomial(Vector
*c
, const evalue
* X
)
4320 unsigned dim
= c
->Size
-2;
4322 evalue
*EP
= ALLOC(evalue
);
4327 if (EVALUE_IS_ZERO(*X
) || dim
== 0) {
4328 evalue_set(EP
, c
->p
[0], c
->p
[dim
+1]);
4329 reduce_constant(EP
);
4333 evalue_set(EP
, c
->p
[dim
], c
->p
[dim
+1]);
4336 evalue_set(&EC
, c
->p
[dim
], c
->p
[dim
+1]);
4338 for (i
= dim
-1; i
>= 0; --i
) {
4340 value_assign(EC
.x
.n
, c
->p
[i
]);
4343 free_evalue_refs(&EC
);
4347 /* Create an evalue from an array of pairs of domains and evalues. */
4348 evalue
*evalue_from_section_array(struct evalue_section
*s
, int n
)
4353 res
= ALLOC(evalue
);
4357 evalue_set_si(res
, 0, 1);
4359 value_set_si(res
->d
, 0);
4360 res
->x
.p
= new_enode(partition
, 2*n
, s
[0].D
->Dimension
);
4361 for (i
= 0; i
< n
; ++i
) {
4362 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
], s
[i
].D
);
4363 value_clear(res
->x
.p
->arr
[2*i
+1].d
);
4364 res
->x
.p
->arr
[2*i
+1] = *s
[i
].E
;