6 #include <barvinok/evalue.h>
7 #include <barvinok/barvinok.h>
8 #include <barvinok/util.h>
10 #ifndef value_pmodulus
11 #define value_pmodulus(ref,val1,val2) (mpz_fdiv_r((ref),(val1),(val2)))
14 #define ALLOC(type) (type*)malloc(sizeof(type))
17 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
19 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
22 void evalue_set_si(evalue
*ev
, int n
, int d
) {
23 value_set_si(ev
->d
, d
);
25 value_set_si(ev
->x
.n
, n
);
28 void evalue_set(evalue
*ev
, Value n
, Value d
) {
29 value_assign(ev
->d
, d
);
31 value_assign(ev
->x
.n
, n
);
36 evalue
*EP
= ALLOC(evalue
);
38 evalue_set_si(EP
, 0, 1);
42 void aep_evalue(evalue
*e
, int *ref
) {
47 if (value_notzero_p(e
->d
))
48 return; /* a rational number, its already reduced */
50 return; /* hum... an overflow probably occured */
52 /* First check the components of p */
53 for (i
=0;i
<p
->size
;i
++)
54 aep_evalue(&p
->arr
[i
],ref
);
61 p
->pos
= ref
[p
->pos
-1]+1;
67 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
73 if (value_notzero_p(e
->d
))
74 return; /* a rational number, its already reduced */
76 return; /* hum... an overflow probably occured */
79 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
80 for(i
=0;i
<CT
->NbRows
-1;i
++)
81 for(j
=0;j
<CT
->NbColumns
;j
++)
82 if(value_notzero_p(CT
->p
[i
][j
])) {
87 /* Transform the references in e, using ref */
91 } /* addeliminatedparams_evalue */
93 void addeliminatedparams_enum(evalue
*e
, Matrix
*CT
, Polyhedron
*CEq
,
94 unsigned MaxRays
, unsigned nparam
)
99 if (CT
->NbRows
== CT
->NbColumns
)
102 if (EVALUE_IS_ZERO(*e
))
105 if (value_notzero_p(e
->d
)) {
108 value_set_si(res
.d
, 0);
109 res
.x
.p
= new_enode(partition
, 2, nparam
);
110 EVALUE_SET_DOMAIN(res
.x
.p
->arr
[0],
111 DomainConstraintSimplify(Polyhedron_Copy(CEq
), MaxRays
));
112 value_clear(res
.x
.p
->arr
[1].d
);
113 res
.x
.p
->arr
[1] = *e
;
120 assert(p
->type
== partition
);
123 for (i
=0; i
<p
->size
/2; i
++) {
124 Polyhedron
*D
= EVALUE_DOMAIN(p
->arr
[2*i
]);
125 Polyhedron
*T
= DomainPreimage(D
, CT
, MaxRays
);
128 T
= DomainIntersection(D
, CEq
, MaxRays
);
130 EVALUE_SET_DOMAIN(p
->arr
[2*i
], T
);
131 addeliminatedparams_evalue(&p
->arr
[2*i
+1], CT
);
135 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
143 assert(value_notzero_p(e1
->d
));
144 assert(value_notzero_p(e2
->d
));
145 value_multiply(m
, e1
->x
.n
, e2
->d
);
146 value_multiply(m2
, e2
->x
.n
, e1
->d
);
149 else if (value_gt(m
, m2
))
159 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
161 if (value_notzero_p(e1
->d
)) {
163 if (value_zero_p(e2
->d
))
165 r
= mod_rational_smaller(e1
, e2
);
166 return r
== -1 ? 0 : r
;
168 if (value_notzero_p(e2
->d
))
170 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
172 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
175 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
177 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
182 static int mod_term_smaller(evalue
*e1
, evalue
*e2
)
184 assert(value_zero_p(e1
->d
));
185 assert(value_zero_p(e2
->d
));
186 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
187 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
188 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
191 /* Negative pos means inequality */
192 /* s is negative of substitution if m is not zero */
201 struct fixed_param
*fixed
;
206 static int relations_depth(evalue
*e
)
211 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
212 e
= &e
->x
.p
->arr
[1], ++d
);
216 static void poly_denom_not_constant(evalue
**pp
, Value
*d
)
221 while (value_zero_p(p
->d
)) {
222 assert(p
->x
.p
->type
== polynomial
);
223 assert(p
->x
.p
->size
== 2);
224 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
225 value_lcm(*d
, p
->x
.p
->arr
[1].d
, d
);
231 static void poly_denom(evalue
*p
, Value
*d
)
233 poly_denom_not_constant(&p
, d
);
234 value_lcm(*d
, p
->d
, d
);
237 static void realloc_substitution(struct subst
*s
, int d
)
239 struct fixed_param
*n
;
242 for (i
= 0; i
< s
->n
; ++i
)
249 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
255 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
258 /* May have been reduced already */
259 if (value_notzero_p(m
->d
))
262 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
263 assert(m
->x
.p
->size
== 3);
265 /* fractional was inverted during reduction
266 * invert it back and move constant in
268 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
269 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
270 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
271 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
272 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
273 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
274 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
275 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
276 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
277 value_set_si(m
->x
.p
->arr
[1].d
, 1);
280 /* Oops. Nested identical relations. */
281 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
284 if (s
->n
>= s
->max
) {
285 int d
= relations_depth(r
);
286 realloc_substitution(s
, d
);
290 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
291 assert(p
->x
.p
->size
== 2);
294 assert(value_pos_p(f
->x
.n
));
296 value_init(s
->fixed
[s
->n
].m
);
297 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
298 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
299 value_init(s
->fixed
[s
->n
].d
);
300 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
301 value_init(s
->fixed
[s
->n
].s
.d
);
302 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
308 static int type_offset(enode
*p
)
310 return p
->type
== fractional
? 1 :
311 p
->type
== flooring
? 1 : 0;
314 static void reorder_terms(enode
*p
, evalue
*v
)
317 int offset
= type_offset(p
);
319 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
321 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
322 free_evalue_refs(&(p
->arr
[i
]));
328 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
334 if (value_notzero_p(e
->d
)) {
336 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
337 return; /* a rational number, its already reduced */
341 return; /* hum... an overflow probably occured */
343 /* First reduce the components of p */
344 add
= p
->type
== relation
;
345 for (i
=0; i
<p
->size
; i
++) {
347 add
= add_modulo_substitution(s
, e
);
349 if (i
== 0 && p
->type
==fractional
)
350 _reduce_evalue(&p
->arr
[i
], s
, 1);
352 _reduce_evalue(&p
->arr
[i
], s
, fract
);
354 if (add
&& i
== p
->size
-1) {
356 value_clear(s
->fixed
[s
->n
].m
);
357 value_clear(s
->fixed
[s
->n
].d
);
358 free_evalue_refs(&s
->fixed
[s
->n
].s
);
359 } else if (add
&& i
== 1)
360 s
->fixed
[s
->n
-1].pos
*= -1;
363 if (p
->type
==periodic
) {
365 /* Try to reduce the period */
366 for (i
=1; i
<=(p
->size
)/2; i
++) {
367 if ((p
->size
% i
)==0) {
369 /* Can we reduce the size to i ? */
371 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
372 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
375 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
379 you_lose
: /* OK, lets not do it */
384 /* Try to reduce its strength */
387 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
391 else if (p
->type
==polynomial
) {
392 for (k
= 0; s
&& k
< s
->n
; ++k
) {
393 if (s
->fixed
[k
].pos
== p
->pos
) {
394 int divide
= value_notone_p(s
->fixed
[k
].d
);
397 if (value_notzero_p(s
->fixed
[k
].m
)) {
400 assert(p
->size
== 2);
401 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
403 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
410 value_assign(d
.d
, s
->fixed
[k
].d
);
412 if (value_notzero_p(s
->fixed
[k
].m
))
413 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
415 value_set_si(d
.x
.n
, 1);
418 for (i
=p
->size
-1;i
>=1;i
--) {
419 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
421 emul(&d
, &p
->arr
[i
]);
422 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
423 free_evalue_refs(&(p
->arr
[i
]));
426 _reduce_evalue(&p
->arr
[0], s
, fract
);
429 free_evalue_refs(&d
);
435 /* Try to reduce the degree */
436 for (i
=p
->size
-1;i
>=1;i
--) {
437 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
439 /* Zero coefficient */
440 free_evalue_refs(&(p
->arr
[i
]));
445 /* Try to reduce its strength */
448 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
452 else if (p
->type
==fractional
) {
456 if (value_notzero_p(p
->arr
[0].d
)) {
458 value_assign(v
.d
, p
->arr
[0].d
);
460 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
465 evalue
*pp
= &p
->arr
[0];
466 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
467 assert(pp
->x
.p
->size
== 2);
469 /* search for exact duplicate among the modulo inequalities */
471 f
= &pp
->x
.p
->arr
[1];
472 for (k
= 0; s
&& k
< s
->n
; ++k
) {
473 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
474 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
475 value_eq(s
->fixed
[k
].m
, f
->d
) &&
476 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
483 /* replace { E/m } by { (E-1)/m } + 1/m */
488 evalue_set_si(&extra
, 1, 1);
489 value_assign(extra
.d
, g
);
490 eadd(&extra
, &v
.x
.p
->arr
[1]);
491 free_evalue_refs(&extra
);
493 /* We've been going in circles; stop now */
494 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
495 free_evalue_refs(&v
);
497 evalue_set_si(&v
, 0, 1);
502 value_set_si(v
.d
, 0);
503 v
.x
.p
= new_enode(fractional
, 3, -1);
504 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
505 value_assign(v
.x
.p
->arr
[1].d
, g
);
506 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
507 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
510 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
513 value_division(f
->d
, g
, f
->d
);
514 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
515 value_assign(f
->d
, g
);
516 value_decrement(f
->x
.n
, f
->x
.n
);
517 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
519 Gcd(f
->d
, f
->x
.n
, &g
);
520 value_division(f
->d
, f
->d
, g
);
521 value_division(f
->x
.n
, f
->x
.n
, g
);
530 /* reduction may have made this fractional arg smaller */
531 i
= reorder
? p
->size
: 1;
532 for ( ; i
< p
->size
; ++i
)
533 if (value_zero_p(p
->arr
[i
].d
) &&
534 p
->arr
[i
].x
.p
->type
== fractional
&&
535 !mod_term_smaller(e
, &p
->arr
[i
]))
539 value_set_si(v
.d
, 0);
540 v
.x
.p
= new_enode(fractional
, 3, -1);
541 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
542 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
543 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
551 evalue
*pp
= &p
->arr
[0];
554 poly_denom_not_constant(&pp
, &m
);
555 mpz_fdiv_r(r
, m
, pp
->d
);
556 if (value_notzero_p(r
)) {
558 value_set_si(v
.d
, 0);
559 v
.x
.p
= new_enode(fractional
, 3, -1);
561 value_multiply(r
, m
, pp
->x
.n
);
562 value_multiply(v
.x
.p
->arr
[1].d
, m
, pp
->d
);
563 value_init(v
.x
.p
->arr
[1].x
.n
);
564 mpz_fdiv_r(v
.x
.p
->arr
[1].x
.n
, r
, pp
->d
);
565 mpz_fdiv_q(r
, r
, pp
->d
);
567 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
568 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
570 while (value_zero_p(pp
->d
))
571 pp
= &pp
->x
.p
->arr
[0];
573 value_assign(pp
->d
, m
);
574 value_assign(pp
->x
.n
, r
);
576 Gcd(pp
->d
, pp
->x
.n
, &r
);
577 value_division(pp
->d
, pp
->d
, r
);
578 value_division(pp
->x
.n
, pp
->x
.n
, r
);
591 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
592 pp
= &pp
->x
.p
->arr
[0]) {
593 f
= &pp
->x
.p
->arr
[1];
594 assert(value_pos_p(f
->d
));
595 mpz_mul_ui(twice
, f
->x
.n
, 2);
596 if (value_lt(twice
, f
->d
))
598 if (value_eq(twice
, f
->d
))
606 value_set_si(v
.d
, 0);
607 v
.x
.p
= new_enode(fractional
, 3, -1);
608 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
609 poly_denom(&p
->arr
[0], &twice
);
610 value_assign(v
.x
.p
->arr
[1].d
, twice
);
611 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
612 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
613 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
615 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
616 pp
= &pp
->x
.p
->arr
[0]) {
617 f
= &pp
->x
.p
->arr
[1];
618 value_oppose(f
->x
.n
, f
->x
.n
);
619 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
621 value_division(pp
->d
, twice
, pp
->d
);
622 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
623 value_assign(pp
->d
, twice
);
624 value_oppose(pp
->x
.n
, pp
->x
.n
);
625 value_decrement(pp
->x
.n
, pp
->x
.n
);
626 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
628 /* Maybe we should do this during reduction of
631 Gcd(pp
->d
, pp
->x
.n
, &twice
);
632 value_division(pp
->d
, pp
->d
, twice
);
633 value_division(pp
->x
.n
, pp
->x
.n
, twice
);
643 reorder_terms(p
, &v
);
644 _reduce_evalue(&p
->arr
[1], s
, fract
);
647 /* Try to reduce the degree */
648 for (i
=p
->size
-1;i
>=2;i
--) {
649 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
651 /* Zero coefficient */
652 free_evalue_refs(&(p
->arr
[i
]));
657 /* Try to reduce its strength */
660 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
661 free_evalue_refs(&(p
->arr
[0]));
665 else if (p
->type
== flooring
) {
666 /* Try to reduce the degree */
667 for (i
=p
->size
-1;i
>=2;i
--) {
668 if (!EVALUE_IS_ZERO(p
->arr
[i
]))
670 /* Zero coefficient */
671 free_evalue_refs(&(p
->arr
[i
]));
676 /* Try to reduce its strength */
679 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
680 free_evalue_refs(&(p
->arr
[0]));
684 else if (p
->type
== relation
) {
685 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
686 free_evalue_refs(&(p
->arr
[2]));
687 free_evalue_refs(&(p
->arr
[0]));
694 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
695 free_evalue_refs(&(p
->arr
[2]));
698 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
699 free_evalue_refs(&(p
->arr
[1]));
700 free_evalue_refs(&(p
->arr
[0]));
701 evalue_set_si(e
, 0, 1);
708 /* Relation was reduced by means of an identical
709 * inequality => remove
711 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
714 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
715 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
717 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
719 free_evalue_refs(&(p
->arr
[2]));
723 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
725 evalue_set_si(e
, 0, 1);
726 free_evalue_refs(&(p
->arr
[1]));
728 free_evalue_refs(&(p
->arr
[0]));
734 } /* reduce_evalue */
736 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
741 for (k
= 0; k
< dim
; ++k
)
742 if (value_notzero_p(row
[k
+1]))
745 Vector_Normalize_Positive(row
+1, dim
+1, k
);
746 assert(s
->n
< s
->max
);
747 value_init(s
->fixed
[s
->n
].d
);
748 value_init(s
->fixed
[s
->n
].m
);
749 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
750 s
->fixed
[s
->n
].pos
= k
+1;
751 value_set_si(s
->fixed
[s
->n
].m
, 0);
752 r
= &s
->fixed
[s
->n
].s
;
754 for (l
= k
+1; l
< dim
; ++l
)
755 if (value_notzero_p(row
[l
+1])) {
756 value_set_si(r
->d
, 0);
757 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
758 value_init(r
->x
.p
->arr
[1].x
.n
);
759 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
760 value_set_si(r
->x
.p
->arr
[1].d
, 1);
764 value_oppose(r
->x
.n
, row
[dim
+1]);
765 value_set_si(r
->d
, 1);
769 static void _reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
, struct subst
*s
)
772 Polyhedron
*orig
= D
;
777 D
= DomainConvex(D
, 0);
778 if (!D
->next
&& D
->NbEq
) {
782 realloc_substitution(s
, dim
);
784 int d
= relations_depth(e
);
786 NALLOC(s
->fixed
, s
->max
);
789 for (j
= 0; j
< D
->NbEq
; ++j
)
790 add_substitution(s
, D
->Constraint
[j
], dim
);
794 _reduce_evalue(e
, s
, 0);
797 for (j
= 0; j
< s
->n
; ++j
) {
798 value_clear(s
->fixed
[j
].d
);
799 value_clear(s
->fixed
[j
].m
);
800 free_evalue_refs(&s
->fixed
[j
].s
);
805 void reduce_evalue_in_domain(evalue
*e
, Polyhedron
*D
)
807 struct subst s
= { NULL
, 0, 0 };
808 _reduce_evalue_in_domain(e
, D
, &s
);
813 void reduce_evalue (evalue
*e
) {
814 struct subst s
= { NULL
, 0, 0 };
816 if (value_notzero_p(e
->d
))
817 return; /* a rational number, its already reduced */
819 if (e
->x
.p
->type
== partition
) {
822 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
823 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
825 /* This shouldn't really happen;
826 * Empty domains should not be added.
828 POL_ENSURE_VERTICES(D
);
830 _reduce_evalue_in_domain(&e
->x
.p
->arr
[2*i
+1], D
, &s
);
832 if (emptyQ(D
) || EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
833 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
834 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
835 value_clear(e
->x
.p
->arr
[2*i
].d
);
837 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
838 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
842 if (e
->x
.p
->size
== 0) {
844 evalue_set_si(e
, 0, 1);
847 _reduce_evalue(e
, &s
, 0);
852 void print_evalue(FILE *DST
,evalue
*e
,char **pname
) {
854 if(value_notzero_p(e
->d
)) {
855 if(value_notone_p(e
->d
)) {
856 value_print(DST
,VALUE_FMT
,e
->x
.n
);
858 value_print(DST
,VALUE_FMT
,e
->d
);
861 value_print(DST
,VALUE_FMT
,e
->x
.n
);
865 print_enode(DST
,e
->x
.p
,pname
);
869 void print_enode(FILE *DST
,enode
*p
,char **pname
) {
874 fprintf(DST
, "NULL");
880 for (i
=0; i
<p
->size
; i
++) {
881 print_evalue(DST
, &p
->arr
[i
], pname
);
885 fprintf(DST
, " }\n");
889 for (i
=p
->size
-1; i
>=0; i
--) {
890 print_evalue(DST
, &p
->arr
[i
], pname
);
891 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
893 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
895 fprintf(DST
, " )\n");
899 for (i
=0; i
<p
->size
; i
++) {
900 print_evalue(DST
, &p
->arr
[i
], pname
);
901 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
903 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
908 for (i
=p
->size
-1; i
>=1; i
--) {
909 print_evalue(DST
, &p
->arr
[i
], pname
);
912 fprintf(DST
, p
->type
== flooring
? "[" : "{");
913 print_evalue(DST
, &p
->arr
[0], pname
);
914 fprintf(DST
, p
->type
== flooring
? "]" : "}");
916 fprintf(DST
, "^%d + ", i
-1);
921 fprintf(DST
, " )\n");
925 print_evalue(DST
, &p
->arr
[0], pname
);
926 fprintf(DST
, "= 0 ] * \n");
927 print_evalue(DST
, &p
->arr
[1], pname
);
929 fprintf(DST
, " +\n [ ");
930 print_evalue(DST
, &p
->arr
[0], pname
);
931 fprintf(DST
, "!= 0 ] * \n");
932 print_evalue(DST
, &p
->arr
[2], pname
);
936 char **names
= pname
;
937 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
938 if (!pname
|| p
->pos
< maxdim
) {
939 NALLOC(names
, maxdim
);
940 for (i
= 0; i
< p
->pos
; ++i
) {
944 NALLOC(names
[i
], 10);
945 snprintf(names
[i
], 10, "%c", 'P'+i
);
948 for ( ; i
< maxdim
; ++i
) {
949 NALLOC(names
[i
], 10);
950 snprintf(names
[i
], 10, "_p%d", i
);
954 for (i
=0; i
<p
->size
/2; i
++) {
955 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
956 print_evalue(DST
, &p
->arr
[2*i
+1], names
);
959 if (!pname
|| p
->pos
< maxdim
) {
960 for (i
= pname
? p
->pos
: 0; i
< maxdim
; ++i
)
973 static void eadd_rev(evalue
*e1
, evalue
*res
)
977 evalue_copy(&ev
, e1
);
979 free_evalue_refs(res
);
983 static void eadd_rev_cst (evalue
*e1
, evalue
*res
)
987 evalue_copy(&ev
, e1
);
988 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
989 free_evalue_refs(res
);
993 static int is_zero_on(evalue
*e
, Polyhedron
*D
)
998 tmp
.x
.p
= new_enode(partition
, 2, D
->Dimension
);
999 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Domain_Copy(D
));
1000 evalue_copy(&tmp
.x
.p
->arr
[1], e
);
1001 reduce_evalue(&tmp
);
1002 is_zero
= EVALUE_IS_ZERO(tmp
);
1003 free_evalue_refs(&tmp
);
1007 struct section
{ Polyhedron
* D
; evalue E
; };
1009 void eadd_partitions (evalue
*e1
,evalue
*res
)
1014 s
= (struct section
*)
1015 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
1016 sizeof(struct section
));
1018 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1019 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1020 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1023 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1024 assert(res
->x
.p
->size
>= 2);
1025 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1026 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
1028 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
1030 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1039 /* See if we can extend one of the domains in res to cover fd */
1040 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1041 if (is_zero_on(&res
->x
.p
->arr
[2*i
+1], fd
))
1043 if (i
< res
->x
.p
->size
/2) {
1044 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
1045 DomainConcat(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
])));
1048 value_init(s
[n
].E
.d
);
1049 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
1053 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1054 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
1055 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1057 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1058 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1064 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
1065 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1067 value_init(s
[n
].E
.d
);
1068 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1069 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1070 if (!emptyQ(fd
) && is_zero_on(&e1
->x
.p
->arr
[2*j
+1], fd
)) {
1071 d
= DomainConcat(fd
, d
);
1072 fd
= Empty_Polyhedron(fd
->Dimension
);
1078 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
1082 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1085 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1086 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1087 value_clear(res
->x
.p
->arr
[2*i
].d
);
1092 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1093 for (j
= 0; j
< n
; ++j
) {
1094 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1095 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1096 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1097 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1103 static void explicit_complement(evalue
*res
)
1105 enode
*rel
= new_enode(relation
, 3, 0);
1107 value_clear(rel
->arr
[0].d
);
1108 rel
->arr
[0] = res
->x
.p
->arr
[0];
1109 value_clear(rel
->arr
[1].d
);
1110 rel
->arr
[1] = res
->x
.p
->arr
[1];
1111 value_set_si(rel
->arr
[2].d
, 1);
1112 value_init(rel
->arr
[2].x
.n
);
1113 value_set_si(rel
->arr
[2].x
.n
, 0);
1118 void eadd(evalue
*e1
,evalue
*res
) {
1121 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1122 /* Add two rational numbers */
1128 value_multiply(m1
,e1
->x
.n
,res
->d
);
1129 value_multiply(m2
,res
->x
.n
,e1
->d
);
1130 value_addto(res
->x
.n
,m1
,m2
);
1131 value_multiply(res
->d
,e1
->d
,res
->d
);
1132 Gcd(res
->x
.n
,res
->d
,&g
);
1133 if (value_notone_p(g
)) {
1134 value_division(res
->d
,res
->d
,g
);
1135 value_division(res
->x
.n
,res
->x
.n
,g
);
1137 value_clear(g
); value_clear(m1
); value_clear(m2
);
1140 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
1141 switch (res
->x
.p
->type
) {
1143 /* Add the constant to the constant term of a polynomial*/
1144 eadd(e1
, &res
->x
.p
->arr
[0]);
1147 /* Add the constant to all elements of a periodic number */
1148 for (i
=0; i
<res
->x
.p
->size
; i
++) {
1149 eadd(e1
, &res
->x
.p
->arr
[i
]);
1153 fprintf(stderr
, "eadd: cannot add const with vector\n");
1157 eadd(e1
, &res
->x
.p
->arr
[1]);
1160 assert(EVALUE_IS_ZERO(*e1
));
1161 break; /* Do nothing */
1163 /* Create (zero) complement if needed */
1164 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
1165 explicit_complement(res
);
1166 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1167 eadd(e1
, &res
->x
.p
->arr
[i
]);
1173 /* add polynomial or periodic to constant
1174 * you have to exchange e1 and res, before doing addition */
1176 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1180 else { // ((e1->d==0) && (res->d==0))
1181 assert(!((e1
->x
.p
->type
== partition
) ^
1182 (res
->x
.p
->type
== partition
)));
1183 if (e1
->x
.p
->type
== partition
) {
1184 eadd_partitions(e1
, res
);
1187 if (e1
->x
.p
->type
== relation
&&
1188 (res
->x
.p
->type
!= relation
||
1189 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1193 if (res
->x
.p
->type
== relation
) {
1194 if (e1
->x
.p
->type
== relation
&&
1195 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1196 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1197 explicit_complement(res
);
1198 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1199 explicit_complement(e1
);
1200 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1201 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1204 if (res
->x
.p
->size
< 3)
1205 explicit_complement(res
);
1206 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1207 eadd(e1
, &res
->x
.p
->arr
[i
]);
1210 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
1211 /* adding to evalues of different type. two cases are possible
1212 * res is periodic and e1 is polynomial, you have to exchange
1213 * e1 and res then to add e1 to the constant term of res */
1214 if (e1
->x
.p
->type
== polynomial
) {
1215 eadd_rev_cst(e1
, res
);
1217 else if (res
->x
.p
->type
== polynomial
) {
1218 /* res is polynomial and e1 is periodic,
1219 add e1 to the constant term of res */
1221 eadd(e1
,&res
->x
.p
->arr
[0]);
1227 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
1228 ((res
->x
.p
->type
== fractional
||
1229 res
->x
.p
->type
== flooring
) &&
1230 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1231 /* adding evalues of different position (i.e function of different unknowns
1232 * to case are possible */
1234 switch (res
->x
.p
->type
) {
1237 if(mod_term_smaller(res
, e1
))
1238 eadd(e1
,&res
->x
.p
->arr
[1]);
1240 eadd_rev_cst(e1
, res
);
1242 case polynomial
: // res and e1 are polynomials
1243 // add e1 to the constant term of res
1245 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1246 eadd(e1
,&res
->x
.p
->arr
[0]);
1248 eadd_rev_cst(e1
, res
);
1249 // value_clear(g); value_clear(m1); value_clear(m2);
1251 case periodic
: // res and e1 are pointers to periodic numbers
1252 //add e1 to all elements of res
1254 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1255 for (i
=0;i
<res
->x
.p
->size
;i
++) {
1256 eadd(e1
,&res
->x
.p
->arr
[i
]);
1267 //same type , same pos and same size
1268 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1269 // add any element in e1 to the corresponding element in res
1270 i
= type_offset(res
->x
.p
);
1272 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1273 for (; i
<res
->x
.p
->size
; i
++) {
1274 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1279 /* Sizes are different */
1280 switch(res
->x
.p
->type
) {
1284 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1285 /* new enode and add res to that new node. If you do not do */
1286 /* that, you lose the the upper weight part of e1 ! */
1288 if(e1
->x
.p
->size
> res
->x
.p
->size
)
1291 i
= type_offset(res
->x
.p
);
1293 assert(eequal(&e1
->x
.p
->arr
[0],
1294 &res
->x
.p
->arr
[0]));
1295 for (; i
<e1
->x
.p
->size
; i
++) {
1296 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1303 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1306 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1307 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1308 to add periodicaly elements of e1 to elements of ne, and finaly to
1313 value_init(ex
); value_init(ey
);value_init(ep
);
1316 value_set_si(ex
,e1
->x
.p
->size
);
1317 value_set_si(ey
,res
->x
.p
->size
);
1318 value_assign (ep
,*Lcm(ex
,ey
));
1319 p
=(int)mpz_get_si(ep
);
1320 ne
= (evalue
*) malloc (sizeof(evalue
));
1322 value_set_si( ne
->d
,0);
1324 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1326 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1329 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1332 value_assign(res
->d
, ne
->d
);
1338 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
1347 static void emul_rev (evalue
*e1
, evalue
*res
)
1351 evalue_copy(&ev
, e1
);
1353 free_evalue_refs(res
);
1357 static void emul_poly (evalue
*e1
, evalue
*res
)
1359 int i
, j
, o
= type_offset(res
->x
.p
);
1361 int size
=(e1
->x
.p
->size
+ res
->x
.p
->size
- o
- 1);
1363 value_set_si(tmp
.d
,0);
1364 tmp
.x
.p
=new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1366 evalue_copy(&tmp
.x
.p
->arr
[0], &e1
->x
.p
->arr
[0]);
1367 for (i
=o
; i
< e1
->x
.p
->size
; i
++) {
1368 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1369 emul(&res
->x
.p
->arr
[o
], &tmp
.x
.p
->arr
[i
]);
1372 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1373 for (i
=o
+1; i
<res
->x
.p
->size
; i
++)
1374 for (j
=o
; j
<e1
->x
.p
->size
; j
++) {
1377 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1378 emul(&res
->x
.p
->arr
[i
], &ev
);
1379 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-o
]);
1380 free_evalue_refs(&ev
);
1382 free_evalue_refs(res
);
1386 void emul_partitions (evalue
*e1
,evalue
*res
)
1391 s
= (struct section
*)
1392 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1393 sizeof(struct section
));
1395 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1396 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1397 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1400 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1401 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1402 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1403 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1409 /* This code is only needed because the partitions
1410 are not true partitions.
1412 for (k
= 0; k
< n
; ++k
) {
1413 if (DomainIncludes(s
[k
].D
, d
))
1415 if (DomainIncludes(d
, s
[k
].D
)) {
1416 Domain_Free(s
[k
].D
);
1417 free_evalue_refs(&s
[k
].E
);
1428 value_init(s
[n
].E
.d
);
1429 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1430 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1434 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1435 value_clear(res
->x
.p
->arr
[2*i
].d
);
1436 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1441 evalue_set_si(res
, 0, 1);
1443 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1444 for (j
= 0; j
< n
; ++j
) {
1445 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1446 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1447 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1448 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1455 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1457 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1458 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1459 * evalues is not treated here */
1461 void emul (evalue
*e1
, evalue
*res
){
1464 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1465 fprintf(stderr
, "emul: do not proced on evector type !\n");
1469 if (EVALUE_IS_ZERO(*res
))
1472 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1473 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1474 emul_partitions(e1
, res
);
1477 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1478 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1479 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1481 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1482 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1483 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1484 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1485 explicit_complement(res
);
1486 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1487 explicit_complement(e1
);
1488 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1489 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1492 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1493 emul(e1
, &res
->x
.p
->arr
[i
]);
1495 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1496 switch(e1
->x
.p
->type
) {
1498 switch(res
->x
.p
->type
) {
1500 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1501 /* Product of two polynomials of the same variable */
1506 /* Product of two polynomials of different variables */
1508 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1509 for( i
=0; i
<res
->x
.p
->size
; i
++)
1510 emul(e1
, &res
->x
.p
->arr
[i
]);
1519 /* Product of a polynomial and a periodic or fractional */
1526 switch(res
->x
.p
->type
) {
1528 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1529 /* Product of two periodics of the same parameter and period */
1531 for(i
=0; i
<res
->x
.p
->size
;i
++)
1532 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1537 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1538 /* Product of two periodics of the same parameter and different periods */
1542 value_init(x
); value_init(y
);value_init(z
);
1545 value_set_si(x
,e1
->x
.p
->size
);
1546 value_set_si(y
,res
->x
.p
->size
);
1547 value_assign (z
,*Lcm(x
,y
));
1548 lcm
=(int)mpz_get_si(z
);
1549 newp
= (evalue
*) malloc (sizeof(evalue
));
1550 value_init(newp
->d
);
1551 value_set_si( newp
->d
,0);
1552 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1553 for(i
=0;i
<lcm
;i
++) {
1554 evalue_copy(&newp
->x
.p
->arr
[i
],
1555 &res
->x
.p
->arr
[i
%iy
]);
1558 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1560 value_assign(res
->d
,newp
->d
);
1563 value_clear(x
); value_clear(y
);value_clear(z
);
1567 /* Product of two periodics of different parameters */
1569 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1570 for(i
=0; i
<res
->x
.p
->size
; i
++)
1571 emul(e1
, &(res
->x
.p
->arr
[i
]));
1579 /* Product of a periodic and a polynomial */
1581 for(i
=0; i
<res
->x
.p
->size
; i
++)
1582 emul(e1
, &(res
->x
.p
->arr
[i
]));
1589 switch(res
->x
.p
->type
) {
1591 for(i
=0; i
<res
->x
.p
->size
; i
++)
1592 emul(e1
, &(res
->x
.p
->arr
[i
]));
1599 assert(e1
->x
.p
->type
== res
->x
.p
->type
);
1600 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1601 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1604 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1605 if (e1
->x
.p
->type
!= fractional
|| !value_two_p(d
.d
))
1610 value_set_si(d
.x
.n
, 1);
1611 /* { x }^2 == { x }/2 */
1612 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1613 assert(e1
->x
.p
->size
== 3);
1614 assert(res
->x
.p
->size
== 3);
1616 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1618 eadd(&res
->x
.p
->arr
[1], &tmp
);
1619 emul(&e1
->x
.p
->arr
[2], &tmp
);
1620 emul(&e1
->x
.p
->arr
[1], res
);
1621 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1622 free_evalue_refs(&tmp
);
1627 if(mod_term_smaller(res
, e1
))
1628 for(i
=1; i
<res
->x
.p
->size
; i
++)
1629 emul(e1
, &(res
->x
.p
->arr
[i
]));
1644 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1645 /* Product of two rational numbers */
1649 value_multiply(res
->d
,e1
->d
,res
->d
);
1650 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1651 Gcd(res
->x
.n
, res
->d
,&g
);
1652 if (value_notone_p(g
)) {
1653 value_division(res
->d
,res
->d
,g
);
1654 value_division(res
->x
.n
,res
->x
.n
,g
);
1660 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1661 /* Product of an expression (polynomial or peririodic) and a rational number */
1667 /* Product of a rationel number and an expression (polynomial or peririodic) */
1669 i
= type_offset(res
->x
.p
);
1670 for (; i
<res
->x
.p
->size
; i
++)
1671 emul(e1
, &res
->x
.p
->arr
[i
]);
1681 /* Frees mask content ! */
1682 void emask(evalue
*mask
, evalue
*res
) {
1689 if (EVALUE_IS_ZERO(*res
)) {
1690 free_evalue_refs(mask
);
1694 assert(value_zero_p(mask
->d
));
1695 assert(mask
->x
.p
->type
== partition
);
1696 assert(value_zero_p(res
->d
));
1697 assert(res
->x
.p
->type
== partition
);
1698 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1699 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1700 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1701 pos
= res
->x
.p
->pos
;
1703 s
= (struct section
*)
1704 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1705 sizeof(struct section
));
1709 evalue_set_si(&mone
, -1, 1);
1712 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1713 assert(mask
->x
.p
->size
>= 2);
1714 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1715 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1717 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1719 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1728 value_init(s
[n
].E
.d
);
1729 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1733 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1734 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1737 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1738 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1739 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1740 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1742 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1743 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1749 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1750 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1752 value_init(s
[n
].E
.d
);
1753 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1754 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1760 /* Just ignore; this may have been previously masked off */
1762 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1766 free_evalue_refs(&mone
);
1767 free_evalue_refs(mask
);
1768 free_evalue_refs(res
);
1771 evalue_set_si(res
, 0, 1);
1773 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1774 for (j
= 0; j
< n
; ++j
) {
1775 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1776 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1777 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1784 void evalue_copy(evalue
*dst
, evalue
*src
)
1786 value_assign(dst
->d
, src
->d
);
1787 if(value_notzero_p(src
->d
)) {
1788 value_init(dst
->x
.n
);
1789 value_assign(dst
->x
.n
, src
->x
.n
);
1791 dst
->x
.p
= ecopy(src
->x
.p
);
1794 enode
*new_enode(enode_type type
,int size
,int pos
) {
1800 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1803 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1807 for(i
=0; i
<size
; i
++) {
1808 value_init(res
->arr
[i
].d
);
1809 value_set_si(res
->arr
[i
].d
,0);
1810 res
->arr
[i
].x
.p
= 0;
1815 enode
*ecopy(enode
*e
) {
1820 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1821 for(i
=0;i
<e
->size
;++i
) {
1822 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1823 if(value_zero_p(res
->arr
[i
].d
))
1824 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1825 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1826 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1828 value_init(res
->arr
[i
].x
.n
);
1829 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1835 int ecmp(const evalue
*e1
, const evalue
*e2
)
1841 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
1845 value_multiply(m
, e1
->x
.n
, e2
->d
);
1846 value_multiply(m2
, e2
->x
.n
, e1
->d
);
1848 if (value_lt(m
, m2
))
1850 else if (value_gt(m
, m2
))
1860 if (value_notzero_p(e1
->d
))
1862 if (value_notzero_p(e2
->d
))
1868 if (p1
->type
!= p2
->type
)
1869 return p1
->type
- p2
->type
;
1870 if (p1
->pos
!= p2
->pos
)
1871 return p1
->pos
- p2
->pos
;
1872 if (p1
->size
!= p2
->size
)
1873 return p1
->size
- p2
->size
;
1875 for (i
= p1
->size
-1; i
>= 0; --i
)
1876 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
1882 int eequal(evalue
*e1
,evalue
*e2
) {
1887 if (value_ne(e1
->d
,e2
->d
))
1890 /* e1->d == e2->d */
1891 if (value_notzero_p(e1
->d
)) {
1892 if (value_ne(e1
->x
.n
,e2
->x
.n
))
1895 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1899 /* e1->d == e2->d == 0 */
1902 if (p1
->type
!= p2
->type
) return 0;
1903 if (p1
->size
!= p2
->size
) return 0;
1904 if (p1
->pos
!= p2
->pos
) return 0;
1905 for (i
=0; i
<p1
->size
; i
++)
1906 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
1911 void free_evalue_refs(evalue
*e
) {
1916 if (EVALUE_IS_DOMAIN(*e
)) {
1917 Domain_Free(EVALUE_DOMAIN(*e
));
1920 } else if (value_pos_p(e
->d
)) {
1922 /* 'e' stores a constant */
1924 value_clear(e
->x
.n
);
1927 assert(value_zero_p(e
->d
));
1930 if (!p
) return; /* null pointer */
1931 for (i
=0; i
<p
->size
; i
++) {
1932 free_evalue_refs(&(p
->arr
[i
]));
1936 } /* free_evalue_refs */
1938 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
1939 Vector
* val
, evalue
*res
)
1941 unsigned nparam
= periods
->Size
;
1944 double d
= compute_evalue(e
, val
->p
);
1945 d
*= VALUE_TO_DOUBLE(m
);
1950 value_assign(res
->d
, m
);
1951 value_init(res
->x
.n
);
1952 value_set_double(res
->x
.n
, d
);
1953 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
1956 if (value_one_p(periods
->p
[p
]))
1957 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
1962 value_assign(tmp
, periods
->p
[p
]);
1963 value_set_si(res
->d
, 0);
1964 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
1966 value_decrement(tmp
, tmp
);
1967 value_assign(val
->p
[p
], tmp
);
1968 mod2table_r(e
, periods
, m
, p
+1, val
,
1969 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
1970 } while (value_pos_p(tmp
));
1976 static void rel2table(evalue
*e
, int zero
)
1978 if (value_pos_p(e
->d
)) {
1979 if (value_zero_p(e
->x
.n
) == zero
)
1980 value_set_si(e
->x
.n
, 1);
1982 value_set_si(e
->x
.n
, 0);
1983 value_set_si(e
->d
, 1);
1986 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
1987 rel2table(&e
->x
.p
->arr
[i
], zero
);
1991 void evalue_mod2table(evalue
*e
, int nparam
)
1996 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
1999 for (i
=0; i
<p
->size
; i
++) {
2000 evalue_mod2table(&(p
->arr
[i
]), nparam
);
2002 if (p
->type
== relation
) {
2007 evalue_copy(©
, &p
->arr
[0]);
2009 rel2table(&p
->arr
[0], 1);
2010 emul(&p
->arr
[0], &p
->arr
[1]);
2012 rel2table(©
, 0);
2013 emul(©
, &p
->arr
[2]);
2014 eadd(&p
->arr
[2], &p
->arr
[1]);
2015 free_evalue_refs(&p
->arr
[2]);
2016 free_evalue_refs(©
);
2018 free_evalue_refs(&p
->arr
[0]);
2022 } else if (p
->type
== fractional
) {
2023 Vector
*periods
= Vector_Alloc(nparam
);
2024 Vector
*val
= Vector_Alloc(nparam
);
2030 value_set_si(tmp
, 1);
2031 Vector_Set(periods
->p
, 1, nparam
);
2032 Vector_Set(val
->p
, 0, nparam
);
2033 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
2036 assert(p
->type
== polynomial
);
2037 assert(p
->size
== 2);
2038 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
2039 value_lcm(tmp
, p
->arr
[1].d
, &tmp
);
2041 value_lcm(tmp
, ev
->d
, &tmp
);
2043 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
2046 evalue_set_si(&res
, 0, 1);
2047 /* Compute the polynomial using Horner's rule */
2048 for (i
=p
->size
-1;i
>1;i
--) {
2049 eadd(&p
->arr
[i
], &res
);
2052 eadd(&p
->arr
[1], &res
);
2054 free_evalue_refs(e
);
2055 free_evalue_refs(&EP
);
2060 Vector_Free(periods
);
2062 } /* evalue_mod2table */
2064 /********************************************************/
2065 /* function in domain */
2066 /* check if the parameters in list_args */
2067 /* verifies the constraints of Domain P */
2068 /********************************************************/
2069 int in_domain(Polyhedron
*P
, Value
*list_args
) {
2072 Value v
; /* value of the constraint of a row when
2073 parameters are instanciated*/
2079 /*P->Constraint constraint matrice of polyhedron P*/
2080 for(row
=0;row
<P
->NbConstraints
;row
++) {
2081 value_assign(v
,P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
2082 for(col
=1;col
<P
->Dimension
+1;col
++) {
2083 value_multiply(tmp
,P
->Constraint
[row
][col
],list_args
[col
-1]);
2084 value_addto(v
,v
,tmp
);
2086 if (value_notzero_p(P
->Constraint
[row
][0])) {
2088 /*if v is not >=0 then this constraint is not respected */
2089 if (value_neg_p(v
)) {
2093 return P
->next
? in_domain(P
->next
, list_args
) : 0;
2098 /*if v is not = 0 then this constraint is not respected */
2099 if (value_notzero_p(v
))
2104 /*if not return before this point => all
2105 the constraints are respected */
2111 /****************************************************/
2112 /* function compute enode */
2113 /* compute the value of enode p with parameters */
2114 /* list "list_args */
2115 /* compute the polynomial or the periodic */
2116 /****************************************************/
2118 static double compute_enode(enode
*p
, Value
*list_args
) {
2130 if (p
->type
== polynomial
) {
2132 value_assign(param
,list_args
[p
->pos
-1]);
2134 /* Compute the polynomial using Horner's rule */
2135 for (i
=p
->size
-1;i
>0;i
--) {
2136 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2137 res
*=VALUE_TO_DOUBLE(param
);
2139 res
+=compute_evalue(&p
->arr
[0],list_args
);
2141 else if (p
->type
== fractional
) {
2142 double d
= compute_evalue(&p
->arr
[0], list_args
);
2143 d
-= floor(d
+1e-10);
2145 /* Compute the polynomial using Horner's rule */
2146 for (i
=p
->size
-1;i
>1;i
--) {
2147 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2150 res
+=compute_evalue(&p
->arr
[1],list_args
);
2152 else if (p
->type
== flooring
) {
2153 double d
= compute_evalue(&p
->arr
[0], list_args
);
2156 /* Compute the polynomial using Horner's rule */
2157 for (i
=p
->size
-1;i
>1;i
--) {
2158 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2161 res
+=compute_evalue(&p
->arr
[1],list_args
);
2163 else if (p
->type
== periodic
) {
2164 value_assign(m
,list_args
[p
->pos
-1]);
2166 /* Choose the right element of the periodic */
2167 value_set_si(param
,p
->size
);
2168 value_pmodulus(m
,m
,param
);
2169 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2171 else if (p
->type
== relation
) {
2172 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2173 res
= compute_evalue(&p
->arr
[1], list_args
);
2174 else if (p
->size
> 2)
2175 res
= compute_evalue(&p
->arr
[2], list_args
);
2177 else if (p
->type
== partition
) {
2178 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2179 Value
*vals
= list_args
;
2182 for (i
= 0; i
< dim
; ++i
) {
2183 value_init(vals
[i
]);
2185 value_assign(vals
[i
], list_args
[i
]);
2188 for (i
= 0; i
< p
->size
/2; ++i
)
2189 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2190 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2194 for (i
= 0; i
< dim
; ++i
)
2195 value_clear(vals
[i
]);
2204 } /* compute_enode */
2206 /*************************************************/
2207 /* return the value of Ehrhart Polynomial */
2208 /* It returns a double, because since it is */
2209 /* a recursive function, some intermediate value */
2210 /* might not be integral */
2211 /*************************************************/
2213 double compute_evalue(evalue
*e
,Value
*list_args
) {
2217 if (value_notzero_p(e
->d
)) {
2218 if (value_notone_p(e
->d
))
2219 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2221 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2224 res
= compute_enode(e
->x
.p
,list_args
);
2226 } /* compute_evalue */
2229 /****************************************************/
2230 /* function compute_poly : */
2231 /* Check for the good validity domain */
2232 /* return the number of point in the Polyhedron */
2233 /* in allocated memory */
2234 /* Using the Ehrhart pseudo-polynomial */
2235 /****************************************************/
2236 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2239 /* double d; int i; */
2241 tmp
= (Value
*) malloc (sizeof(Value
));
2242 assert(tmp
!= NULL
);
2244 value_set_si(*tmp
,0);
2247 return(tmp
); /* no ehrhart polynomial */
2248 if(en
->ValidityDomain
) {
2249 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2250 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2255 return(tmp
); /* no Validity Domain */
2257 if(in_domain(en
->ValidityDomain
,list_args
)) {
2259 #ifdef EVAL_EHRHART_DEBUG
2260 Print_Domain(stdout
,en
->ValidityDomain
);
2261 print_evalue(stdout
,&en
->EP
);
2264 /* d = compute_evalue(&en->EP,list_args);
2266 printf("(double)%lf = %d\n", d, i ); */
2267 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2273 value_set_si(*tmp
,0);
2274 return(tmp
); /* no compatible domain with the arguments */
2275 } /* compute_poly */
2277 size_t value_size(Value v
) {
2278 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2279 * sizeof(v
[0]._mp_d
[0]);
2282 size_t domain_size(Polyhedron
*D
)
2285 size_t s
= sizeof(*D
);
2287 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2288 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2289 s
+= value_size(D
->Constraint
[i
][j
]);
2292 for (i = 0; i < D->NbRays; ++i)
2293 for (j = 0; j < D->Dimension+2; ++j)
2294 s += value_size(D->Ray[i][j]);
2297 return D
->next
? s
+domain_size(D
->next
) : s
;
2300 size_t enode_size(enode
*p
) {
2301 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2304 if (p
->type
== partition
)
2305 for (i
= 0; i
< p
->size
/2; ++i
) {
2306 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2307 s
+= evalue_size(&p
->arr
[2*i
+1]);
2310 for (i
= 0; i
< p
->size
; ++i
) {
2311 s
+= evalue_size(&p
->arr
[i
]);
2316 size_t evalue_size(evalue
*e
)
2318 size_t s
= sizeof(*e
);
2319 s
+= value_size(e
->d
);
2320 if (value_notzero_p(e
->d
))
2321 s
+= value_size(e
->x
.n
);
2323 s
+= enode_size(e
->x
.p
);
2327 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2329 evalue
*found
= NULL
;
2334 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2337 value_init(offset
.d
);
2338 value_init(offset
.x
.n
);
2339 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2340 value_lcm(m
, offset
.d
, &offset
.d
);
2341 value_set_si(offset
.x
.n
, 1);
2344 evalue_copy(©
, cst
);
2347 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2349 if (eequal(base
, &e
->x
.p
->arr
[0]))
2350 found
= &e
->x
.p
->arr
[0];
2352 value_set_si(offset
.x
.n
, -2);
2355 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2357 if (eequal(base
, &e
->x
.p
->arr
[0]))
2360 free_evalue_refs(cst
);
2361 free_evalue_refs(&offset
);
2364 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2365 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2370 static evalue
*find_relation_pair(evalue
*e
)
2373 evalue
*found
= NULL
;
2375 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2378 if (e
->x
.p
->type
== fractional
) {
2383 poly_denom(&e
->x
.p
->arr
[0], &m
);
2385 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2386 cst
= &cst
->x
.p
->arr
[0])
2389 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2390 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2395 i
= e
->x
.p
->type
== relation
;
2396 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2397 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2402 void evalue_mod2relation(evalue
*e
) {
2405 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2408 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2409 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2410 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2411 value_clear(e
->x
.p
->arr
[2*i
].d
);
2412 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2414 if (2*i
< e
->x
.p
->size
) {
2415 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2416 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2421 if (e
->x
.p
->size
== 0) {
2423 evalue_set_si(e
, 0, 1);
2429 while ((d
= find_relation_pair(e
)) != NULL
) {
2433 value_init(split
.d
);
2434 value_set_si(split
.d
, 0);
2435 split
.x
.p
= new_enode(relation
, 3, 0);
2436 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2437 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2439 ev
= &split
.x
.p
->arr
[0];
2440 value_set_si(ev
->d
, 0);
2441 ev
->x
.p
= new_enode(fractional
, 3, -1);
2442 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2443 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2444 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2450 free_evalue_refs(&split
);
2454 static int evalue_comp(const void * a
, const void * b
)
2456 const evalue
*e1
= *(const evalue
**)a
;
2457 const evalue
*e2
= *(const evalue
**)b
;
2458 return ecmp(e1
, e2
);
2461 void evalue_combine(evalue
*e
)
2468 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2471 NALLOC(evs
, e
->x
.p
->size
/2);
2472 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2473 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2474 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2475 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2476 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2477 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2478 value_clear(p
->arr
[2*k
].d
);
2479 value_clear(p
->arr
[2*k
+1].d
);
2480 p
->arr
[2*k
] = *(evs
[i
]-1);
2481 p
->arr
[2*k
+1] = *(evs
[i
]);
2484 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2487 value_clear((evs
[i
]-1)->d
);
2491 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2492 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2493 free_evalue_refs(evs
[i
]);
2497 for (i
= 2*k
; i
< p
->size
; ++i
)
2498 value_clear(p
->arr
[i
].d
);
2505 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2507 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2509 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2512 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2513 Polyhedron
*D
, *N
, **P
;
2516 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2523 if (D
->NbEq
<= H
->NbEq
) {
2529 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2530 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2531 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2532 reduce_evalue(&tmp
);
2533 if (value_notzero_p(tmp
.d
) ||
2534 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2537 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2538 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2541 free_evalue_refs(&tmp
);
2547 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2549 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2551 value_clear(e
->x
.p
->arr
[2*i
].d
);
2552 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2554 if (2*i
< e
->x
.p
->size
) {
2555 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2556 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2563 H
= DomainConvex(D
, 0);
2564 E
= DomainDifference(H
, D
, 0);
2566 D
= DomainDifference(H
, E
, 0);
2569 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2573 /* May change coefficients to become non-standard if fiddle is set
2574 * => reduce p afterwards to correct
2576 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2577 Matrix
**R
, int fiddle
)
2581 unsigned dim
= D
->Dimension
;
2582 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2587 assert(p
->type
== fractional
);
2589 value_set_si(T
->p
[1][dim
], 1);
2591 while (value_zero_p(pp
->d
)) {
2592 assert(pp
->x
.p
->type
== polynomial
);
2593 assert(pp
->x
.p
->size
== 2);
2594 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2595 value_division(T
->p
[0][pp
->x
.p
->pos
-1], *d
, pp
->x
.p
->arr
[1].d
);
2596 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2597 if (fiddle
&& value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2598 value_subtract(pp
->x
.p
->arr
[1].x
.n
,
2599 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2600 value_multiply(T
->p
[0][pp
->x
.p
->pos
-1],
2601 T
->p
[0][pp
->x
.p
->pos
-1], pp
->x
.p
->arr
[1].x
.n
);
2602 pp
= &pp
->x
.p
->arr
[0];
2604 value_division(T
->p
[0][dim
], *d
, pp
->d
);
2605 value_multiply(T
->p
[0][dim
], T
->p
[0][dim
], pp
->x
.n
);
2606 I
= DomainImage(D
, T
, 0);
2607 H
= DomainConvex(I
, 0);
2616 int evalue_range_reduction_in_domain(evalue
*e
, Polyhedron
*D
)
2626 if (value_notzero_p(e
->d
))
2631 if (p
->type
== relation
) {
2637 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
, 1);
2638 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2639 equal
= value_eq(min
, max
);
2640 mpz_cdiv_q(min
, min
, d
);
2641 mpz_fdiv_q(max
, max
, d
);
2643 if (bounded
&& value_gt(min
, max
)) {
2649 evalue_set_si(e
, 0, 1);
2652 free_evalue_refs(&(p
->arr
[1]));
2653 free_evalue_refs(&(p
->arr
[0]));
2659 return r
? r
: evalue_range_reduction_in_domain(e
, D
);
2660 } else if (bounded
&& equal
) {
2663 free_evalue_refs(&(p
->arr
[2]));
2666 free_evalue_refs(&(p
->arr
[0]));
2672 return evalue_range_reduction_in_domain(e
, D
);
2673 } else if (bounded
&& value_eq(min
, max
)) {
2674 /* zero for a single value */
2676 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2677 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2678 value_multiply(min
, min
, d
);
2679 value_subtract(M
->p
[0][D
->Dimension
+1],
2680 M
->p
[0][D
->Dimension
+1], min
);
2681 E
= DomainAddConstraints(D
, M
, 0);
2687 r
= evalue_range_reduction_in_domain(&p
->arr
[1], E
);
2689 r
|= evalue_range_reduction_in_domain(&p
->arr
[2], D
);
2691 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2699 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2702 i
= p
->type
== relation
? 1 :
2703 p
->type
== fractional
? 1 : 0;
2704 for (; i
<p
->size
; i
++)
2705 r
|= evalue_range_reduction_in_domain(&p
->arr
[i
], D
);
2707 if (p
->type
!= fractional
) {
2708 if (r
&& p
->type
== polynomial
) {
2711 value_set_si(f
.d
, 0);
2712 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2713 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2714 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2715 reorder_terms(p
, &f
);
2726 I
= polynomial_projection(p
, D
, &d
, &T
, 1);
2728 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2729 mpz_fdiv_q(min
, min
, d
);
2730 mpz_fdiv_q(max
, max
, d
);
2731 value_subtract(d
, max
, min
);
2733 if (bounded
&& value_eq(min
, max
)) {
2736 value_init(inc
.x
.n
);
2737 value_set_si(inc
.d
, 1);
2738 value_oppose(inc
.x
.n
, min
);
2739 eadd(&inc
, &p
->arr
[0]);
2740 reorder_terms(p
, &p
->arr
[0]); /* frees arr[0] */
2744 free_evalue_refs(&inc
);
2746 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
2753 value_set_si(rem
.d
, 0);
2754 rem
.x
.p
= new_enode(fractional
, 3, -1);
2755 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
2756 value_clear(rem
.x
.p
->arr
[1].d
);
2757 value_clear(rem
.x
.p
->arr
[2].d
);
2758 rem
.x
.p
->arr
[1] = p
->arr
[1];
2759 rem
.x
.p
->arr
[2] = p
->arr
[2];
2760 for (i
= 3; i
< p
->size
; ++i
)
2761 p
->arr
[i
-2] = p
->arr
[i
];
2765 value_init(inc
.x
.n
);
2766 value_set_si(inc
.d
, 1);
2767 value_oppose(inc
.x
.n
, min
);
2770 evalue_copy(&t
, &p
->arr
[0]);
2774 value_set_si(f
.d
, 0);
2775 f
.x
.p
= new_enode(fractional
, 3, -1);
2776 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
2777 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2778 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
2780 value_init(factor
.d
);
2781 evalue_set_si(&factor
, -1, 1);
2787 value_clear(f
.x
.p
->arr
[1].x
.n
);
2788 value_clear(f
.x
.p
->arr
[2].x
.n
);
2789 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2790 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
2796 free_evalue_refs(&inc
);
2797 free_evalue_refs(&t
);
2798 free_evalue_refs(&f
);
2799 free_evalue_refs(&factor
);
2800 free_evalue_refs(&rem
);
2802 evalue_range_reduction_in_domain(e
, D
);
2806 _reduce_evalue(&p
->arr
[0], 0, 1);
2810 value_set_si(f
.d
, 0);
2811 f
.x
.p
= new_enode(fractional
, 3, -1);
2812 value_clear(f
.x
.p
->arr
[0].d
);
2813 f
.x
.p
->arr
[0] = p
->arr
[0];
2814 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2815 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
2816 reorder_terms(p
, &f
);
2830 void evalue_range_reduction(evalue
*e
)
2833 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2836 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2837 if (evalue_range_reduction_in_domain(&e
->x
.p
->arr
[2*i
+1],
2838 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
2839 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
2841 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2842 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2843 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
2844 value_clear(e
->x
.p
->arr
[2*i
].d
);
2846 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2847 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2855 Enumeration
* partition2enumeration(evalue
*EP
)
2858 Enumeration
*en
, *res
= NULL
;
2860 if (EVALUE_IS_ZERO(*EP
)) {
2865 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
2866 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
2867 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
2870 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
2871 value_clear(EP
->x
.p
->arr
[2*i
].d
);
2872 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
2880 int evalue_frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
2890 if (value_notzero_p(e
->d
))
2895 i
= p
->type
== relation
? 1 :
2896 p
->type
== fractional
? 1 : 0;
2897 for (; i
<p
->size
; i
++)
2898 r
|= evalue_frac2floor_in_domain(&p
->arr
[i
], D
);
2900 if (p
->type
!= fractional
) {
2901 if (r
&& p
->type
== polynomial
) {
2904 value_set_si(f
.d
, 0);
2905 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2906 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2907 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2908 reorder_terms(p
, &f
);
2917 I
= polynomial_projection(p
, D
, &d
, &T
, 0);
2920 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2923 assert(I
->NbEq
== 0); /* Should have been reduced */
2926 for (i
= 0; i
< I
->NbConstraints
; ++i
)
2927 if (value_pos_p(I
->Constraint
[i
][1]))
2930 if (i
< I
->NbConstraints
) {
2932 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
2933 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
2934 if (value_neg_p(min
)) {
2936 mpz_fdiv_q(min
, min
, d
);
2937 value_init(offset
.d
);
2938 value_set_si(offset
.d
, 1);
2939 value_init(offset
.x
.n
);
2940 value_oppose(offset
.x
.n
, min
);
2941 eadd(&offset
, &p
->arr
[0]);
2942 free_evalue_refs(&offset
);
2952 value_set_si(fl
.d
, 0);
2953 fl
.x
.p
= new_enode(flooring
, 3, -1);
2954 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
2955 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
2956 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
2958 eadd(&fl
, &p
->arr
[0]);
2959 reorder_terms(p
, &p
->arr
[0]);
2963 free_evalue_refs(&fl
);
2968 void evalue_frac2floor(evalue
*e
)
2971 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2974 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2975 if (evalue_frac2floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
2976 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
])))
2977 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
2980 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
2985 int nparam
= D
->Dimension
- nvar
;
2988 nr
= D
->NbConstraints
+ 2;
2989 nc
= D
->Dimension
+ 2 + 1;
2990 C
= Matrix_Alloc(nr
, nc
);
2991 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
2992 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
2993 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
2994 D
->Dimension
+ 1 - nvar
);
2999 nc
= C
->NbColumns
+ 1;
3000 C
= Matrix_Alloc(nr
, nc
);
3001 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
3002 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
3003 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
3004 oldC
->NbColumns
- 1 - nvar
);
3007 value_set_si(C
->p
[nr
-2][0], 1);
3008 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
3009 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
3011 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
3012 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
3018 static void floor2frac_r(evalue
*e
, int nvar
)
3025 if (value_notzero_p(e
->d
))
3030 assert(p
->type
== flooring
);
3031 for (i
= 1; i
< p
->size
; i
++)
3032 floor2frac_r(&p
->arr
[i
], nvar
);
3034 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
3035 assert(pp
->x
.p
->type
== polynomial
);
3036 pp
->x
.p
->pos
-= nvar
;
3040 value_set_si(f
.d
, 0);
3041 f
.x
.p
= new_enode(fractional
, 3, -1);
3042 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
3043 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
3044 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
3046 eadd(&f
, &p
->arr
[0]);
3047 reorder_terms(p
, &p
->arr
[0]);
3050 free_evalue_refs(&f
);
3053 /* Convert flooring back to fractional and shift position
3054 * of the parameters by nvar
3056 static void floor2frac(evalue
*e
, int nvar
)
3058 floor2frac_r(e
, nvar
);
3062 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
3065 int nparam
= D
->Dimension
- nvar
;
3069 D
= Constraints2Polyhedron(C
, 0);
3073 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
3075 /* Double check that D was not unbounded. */
3076 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
3084 evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
3091 evalue
*factor
= NULL
;
3094 if (EVALUE_IS_ZERO(*e
))
3098 Polyhedron
*DD
= Disjoint_Domain(D
, 0, 0);
3105 res
= esum_over_domain(e
, nvar
, Q
, C
);
3108 for (Q
= DD
; Q
; Q
= DD
) {
3114 t
= esum_over_domain(e
, nvar
, Q
, C
);
3121 free_evalue_refs(t
);
3128 if (value_notzero_p(e
->d
)) {
3131 t
= esum_over_domain_cst(nvar
, D
, C
);
3133 if (!EVALUE_IS_ONE(*e
))
3139 switch (e
->x
.p
->type
) {
3141 evalue
*pp
= &e
->x
.p
->arr
[0];
3143 if (pp
->x
.p
->pos
> nvar
) {
3144 /* remainder is independent of the summated vars */
3150 floor2frac(&f
, nvar
);
3152 t
= esum_over_domain_cst(nvar
, D
, C
);
3156 free_evalue_refs(&f
);
3161 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3162 poly_denom(pp
, &row
->p
[1 + nvar
]);
3163 value_set_si(row
->p
[0], 1);
3164 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3165 pp
= &pp
->x
.p
->arr
[0]) {
3167 assert(pp
->x
.p
->type
== polynomial
);
3169 if (pos
>= 1 + nvar
)
3171 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3172 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3173 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3175 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3176 value_division(row
->p
[1 + D
->Dimension
+ 1],
3177 row
->p
[1 + D
->Dimension
+ 1],
3179 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3180 row
->p
[1 + D
->Dimension
+ 1],
3182 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3186 int pos
= e
->x
.p
->pos
;
3189 factor
= ALLOC(evalue
);
3190 value_init(factor
->d
);
3191 value_set_si(factor
->d
, 0);
3192 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3193 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3194 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3198 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3199 for (i
= 0; i
< D
->NbRays
; ++i
)
3200 if (value_notzero_p(D
->Ray
[i
][pos
]))
3202 assert(i
< D
->NbRays
);
3203 if (value_neg_p(D
->Ray
[i
][pos
])) {
3204 factor
= ALLOC(evalue
);
3205 value_init(factor
->d
);
3206 evalue_set_si(factor
, -1, 1);
3208 value_set_si(row
->p
[0], 1);
3209 value_set_si(row
->p
[pos
], 1);
3210 value_set_si(row
->p
[1 + nvar
], -1);
3217 i
= type_offset(e
->x
.p
);
3219 res
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
3224 evalue_copy(&cum
, factor
);
3228 for (; i
< e
->x
.p
->size
; ++i
) {
3232 C
= esum_add_constraint(nvar
, D
, C
, row
);
3238 Vector_Print(stderr, P_VALUE_FMT, row);
3240 Matrix_Print(stderr, P_VALUE_FMT, C);
3242 t
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
3251 free_evalue_refs(t
);
3254 if (factor
&& i
+1 < e
->x
.p
->size
)
3261 free_evalue_refs(factor
);
3262 free_evalue_refs(&cum
);
3274 evalue
*esum(evalue
*e
, int nvar
)
3277 evalue
*res
= ALLOC(evalue
);
3281 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3282 evalue_copy(res
, e
);
3286 evalue_set_si(res
, 0, 1);
3288 assert(value_zero_p(e
->d
));
3289 assert(e
->x
.p
->type
== partition
);
3291 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3293 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3294 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
3296 free_evalue_refs(t
);
3305 /* Initial silly implementation */
3306 void eor(evalue
*e1
, evalue
*res
)
3312 evalue_set_si(&mone
, -1, 1);
3314 evalue_copy(&E
, res
);
3320 free_evalue_refs(&E
);
3321 free_evalue_refs(&mone
);