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)))
15 #define ALLOC(p) p = (typeof(p))malloc(sizeof(*p))
16 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
18 #define ALLOC(p) p = (void *)malloc(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
);
34 void aep_evalue(evalue
*e
, int *ref
) {
39 if (value_notzero_p(e
->d
))
40 return; /* a rational number, its already reduced */
42 return; /* hum... an overflow probably occured */
44 /* First check the components of p */
45 for (i
=0;i
<p
->size
;i
++)
46 aep_evalue(&p
->arr
[i
],ref
);
53 p
->pos
= ref
[p
->pos
-1]+1;
59 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
65 if (value_notzero_p(e
->d
))
66 return; /* a rational number, its already reduced */
68 return; /* hum... an overflow probably occured */
71 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
72 for(i
=0;i
<CT
->NbRows
-1;i
++)
73 for(j
=0;j
<CT
->NbColumns
;j
++)
74 if(value_notzero_p(CT
->p
[i
][j
])) {
79 /* Transform the references in e, using ref */
83 } /* addeliminatedparams_evalue */
85 void addeliminatedparams_enum(evalue
*e
, Matrix
*CT
, Polyhedron
*CEq
,
86 unsigned MaxRays
, unsigned nparam
)
91 if (CT
->NbRows
== CT
->NbColumns
)
94 if (EVALUE_IS_ZERO(*e
))
97 if (value_notzero_p(e
->d
)) {
100 value_set_si(res
.d
, 0);
101 res
.x
.p
= new_enode(partition
, 2, nparam
);
102 EVALUE_SET_DOMAIN(res
.x
.p
->arr
[0],
103 DomainConstraintSimplify(Polyhedron_Copy(CEq
), MaxRays
));
104 value_clear(res
.x
.p
->arr
[1].d
);
105 res
.x
.p
->arr
[1] = *e
;
112 assert(p
->type
== partition
);
115 for (i
=0; i
<p
->size
/2; i
++) {
116 Polyhedron
*D
= EVALUE_DOMAIN(p
->arr
[2*i
]);
117 Polyhedron
*T
= DomainPreimage(D
, CT
, MaxRays
);
120 T
= DomainIntersection(D
, CEq
, MaxRays
);
122 EVALUE_SET_DOMAIN(p
->arr
[2*i
], T
);
123 addeliminatedparams_evalue(&p
->arr
[2*i
+1], CT
);
127 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
135 assert(value_notzero_p(e1
->d
));
136 assert(value_notzero_p(e2
->d
));
137 value_multiply(m
, e1
->x
.n
, e2
->d
);
138 value_multiply(m2
, e2
->x
.n
, e1
->d
);
141 else if (value_gt(m
, m2
))
151 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
153 if (value_notzero_p(e1
->d
)) {
155 if (value_zero_p(e2
->d
))
157 r
= mod_rational_smaller(e1
, e2
);
158 return r
== -1 ? 0 : r
;
160 if (value_notzero_p(e2
->d
))
162 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
164 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
167 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
169 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
174 static int mod_term_smaller(evalue
*e1
, evalue
*e2
)
176 assert(value_zero_p(e1
->d
));
177 assert(value_zero_p(e2
->d
));
178 assert(e1
->x
.p
->type
== fractional
|| e1
->x
.p
->type
== flooring
);
179 assert(e2
->x
.p
->type
== fractional
|| e2
->x
.p
->type
== flooring
);
180 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
183 /* Negative pos means inequality */
184 /* s is negative of substitution if m is not zero */
193 struct fixed_param
*fixed
;
198 static int relations_depth(evalue
*e
)
203 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
204 e
= &e
->x
.p
->arr
[1], ++d
);
208 static void poly_denom_not_constant(evalue
**pp
, Value
*d
)
213 while (value_zero_p(p
->d
)) {
214 assert(p
->x
.p
->type
== polynomial
);
215 assert(p
->x
.p
->size
== 2);
216 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
217 value_lcm(*d
, p
->x
.p
->arr
[1].d
, d
);
223 static void poly_denom(evalue
*p
, Value
*d
)
225 poly_denom_not_constant(&p
, d
);
226 value_lcm(*d
, p
->d
, d
);
229 #define EVALUE_IS_ONE(ev) (value_one_p((ev).d) && value_one_p((ev).x.n))
231 static void realloc_substitution(struct subst
*s
, int d
)
233 struct fixed_param
*n
;
236 for (i
= 0; i
< s
->n
; ++i
)
243 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
249 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
252 /* May have been reduced already */
253 if (value_notzero_p(m
->d
))
256 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
257 assert(m
->x
.p
->size
== 3);
259 /* fractional was inverted during reduction
260 * invert it back and move constant in
262 if (!EVALUE_IS_ONE(m
->x
.p
->arr
[2])) {
263 assert(value_pos_p(m
->x
.p
->arr
[2].d
));
264 assert(value_mone_p(m
->x
.p
->arr
[2].x
.n
));
265 value_set_si(m
->x
.p
->arr
[2].x
.n
, 1);
266 value_increment(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].x
.n
);
267 assert(value_eq(m
->x
.p
->arr
[1].x
.n
, m
->x
.p
->arr
[1].d
));
268 value_set_si(m
->x
.p
->arr
[1].x
.n
, 1);
269 eadd(&m
->x
.p
->arr
[1], &m
->x
.p
->arr
[0]);
270 value_set_si(m
->x
.p
->arr
[1].x
.n
, 0);
271 value_set_si(m
->x
.p
->arr
[1].d
, 1);
274 /* Oops. Nested identical relations. */
275 if (!EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
278 if (s
->n
>= s
->max
) {
279 int d
= relations_depth(r
);
280 realloc_substitution(s
, d
);
284 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
285 assert(p
->x
.p
->size
== 2);
288 assert(value_pos_p(f
->x
.n
));
290 value_init(s
->fixed
[s
->n
].m
);
291 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
292 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
293 value_init(s
->fixed
[s
->n
].d
);
294 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
295 value_init(s
->fixed
[s
->n
].s
.d
);
296 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
302 static int type_offset(enode
*p
)
304 return p
->type
== fractional
? 1 :
305 p
->type
== flooring
? 1 : 0;
308 static void reorder_terms(enode
*p
, evalue
*v
)
311 int offset
= type_offset(p
);
313 for (i
= p
->size
-1; i
>= offset
+1; i
--) {
315 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
316 free_evalue_refs(&(p
->arr
[i
]));
322 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
328 if (value_notzero_p(e
->d
)) {
330 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
331 return; /* a rational number, its already reduced */
335 return; /* hum... an overflow probably occured */
337 /* First reduce the components of p */
338 add
= p
->type
== relation
;
339 for (i
=0; i
<p
->size
; i
++) {
341 add
= add_modulo_substitution(s
, e
);
343 if (i
== 0 && p
->type
==fractional
)
344 _reduce_evalue(&p
->arr
[i
], s
, 1);
346 _reduce_evalue(&p
->arr
[i
], s
, fract
);
348 if (add
&& i
== p
->size
-1) {
350 value_clear(s
->fixed
[s
->n
].m
);
351 value_clear(s
->fixed
[s
->n
].d
);
352 free_evalue_refs(&s
->fixed
[s
->n
].s
);
353 } else if (add
&& i
== 1)
354 s
->fixed
[s
->n
-1].pos
*= -1;
357 if (p
->type
==periodic
) {
359 /* Try to reduce the period */
360 for (i
=1; i
<=(p
->size
)/2; i
++) {
361 if ((p
->size
% i
)==0) {
363 /* Can we reduce the size to i ? */
365 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
366 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
369 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
373 you_lose
: /* OK, lets not do it */
378 /* Try to reduce its strength */
381 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
385 else if (p
->type
==polynomial
) {
386 for (k
= 0; s
&& k
< s
->n
; ++k
) {
387 if (s
->fixed
[k
].pos
== p
->pos
) {
388 int divide
= value_notone_p(s
->fixed
[k
].d
);
391 if (value_notzero_p(s
->fixed
[k
].m
)) {
394 assert(p
->size
== 2);
395 if (divide
&& value_ne(s
->fixed
[k
].d
, p
->arr
[1].x
.n
))
397 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
404 value_assign(d
.d
, s
->fixed
[k
].d
);
406 if (value_notzero_p(s
->fixed
[k
].m
))
407 value_oppose(d
.x
.n
, s
->fixed
[k
].m
);
409 value_set_si(d
.x
.n
, 1);
412 for (i
=p
->size
-1;i
>=1;i
--) {
413 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
415 emul(&d
, &p
->arr
[i
]);
416 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
417 free_evalue_refs(&(p
->arr
[i
]));
420 _reduce_evalue(&p
->arr
[0], s
, fract
);
423 free_evalue_refs(&d
);
429 /* Try to reduce the degree */
430 for (i
=p
->size
-1;i
>=1;i
--) {
431 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
433 /* Zero coefficient */
434 free_evalue_refs(&(p
->arr
[i
]));
439 /* Try to reduce its strength */
442 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
446 else if (p
->type
==fractional
) {
450 if (value_notzero_p(p
->arr
[0].d
)) {
452 value_assign(v
.d
, p
->arr
[0].d
);
454 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
459 evalue
*pp
= &p
->arr
[0];
460 assert(value_zero_p(pp
->d
) && pp
->x
.p
->type
== polynomial
);
461 assert(pp
->x
.p
->size
== 2);
463 /* search for exact duplicate among the modulo inequalities */
465 f
= &pp
->x
.p
->arr
[1];
466 for (k
= 0; s
&& k
< s
->n
; ++k
) {
467 if (-s
->fixed
[k
].pos
== pp
->x
.p
->pos
&&
468 value_eq(s
->fixed
[k
].d
, f
->x
.n
) &&
469 value_eq(s
->fixed
[k
].m
, f
->d
) &&
470 eequal(&s
->fixed
[k
].s
, &pp
->x
.p
->arr
[0]))
477 /* replace { E/m } by { (E-1)/m } + 1/m */
482 evalue_set_si(&extra
, 1, 1);
483 value_assign(extra
.d
, g
);
484 eadd(&extra
, &v
.x
.p
->arr
[1]);
485 free_evalue_refs(&extra
);
487 /* We've been going in circles; stop now */
488 if (value_ge(v
.x
.p
->arr
[1].x
.n
, v
.x
.p
->arr
[1].d
)) {
489 free_evalue_refs(&v
);
491 evalue_set_si(&v
, 0, 1);
496 value_set_si(v
.d
, 0);
497 v
.x
.p
= new_enode(fractional
, 3, -1);
498 evalue_set_si(&v
.x
.p
->arr
[1], 1, 1);
499 value_assign(v
.x
.p
->arr
[1].d
, g
);
500 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
501 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
504 for (f
= &v
.x
.p
->arr
[0]; value_zero_p(f
->d
);
507 value_division(f
->d
, g
, f
->d
);
508 value_multiply(f
->x
.n
, f
->x
.n
, f
->d
);
509 value_assign(f
->d
, g
);
510 value_decrement(f
->x
.n
, f
->x
.n
);
511 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
513 Gcd(f
->d
, f
->x
.n
, &g
);
514 value_division(f
->d
, f
->d
, g
);
515 value_division(f
->x
.n
, f
->x
.n
, g
);
524 /* reduction may have made this fractional arg smaller */
525 i
= reorder
? p
->size
: 1;
526 for ( ; i
< p
->size
; ++i
)
527 if (value_zero_p(p
->arr
[i
].d
) &&
528 p
->arr
[i
].x
.p
->type
== fractional
&&
529 !mod_term_smaller(e
, &p
->arr
[i
]))
533 value_set_si(v
.d
, 0);
534 v
.x
.p
= new_enode(fractional
, 3, -1);
535 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
536 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
537 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
545 evalue
*pp
= &p
->arr
[0];
548 poly_denom_not_constant(&pp
, &m
);
549 mpz_fdiv_r(r
, m
, pp
->d
);
550 if (value_notzero_p(r
)) {
552 value_set_si(v
.d
, 0);
553 v
.x
.p
= new_enode(fractional
, 3, -1);
555 value_multiply(r
, m
, pp
->x
.n
);
556 value_multiply(v
.x
.p
->arr
[1].d
, m
, pp
->d
);
557 value_init(v
.x
.p
->arr
[1].x
.n
);
558 mpz_fdiv_r(v
.x
.p
->arr
[1].x
.n
, r
, pp
->d
);
559 mpz_fdiv_q(r
, r
, pp
->d
);
561 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
562 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
564 while (value_zero_p(pp
->d
))
565 pp
= &pp
->x
.p
->arr
[0];
567 value_assign(pp
->d
, m
);
568 value_assign(pp
->x
.n
, r
);
570 Gcd(pp
->d
, pp
->x
.n
, &r
);
571 value_division(pp
->d
, pp
->d
, r
);
572 value_division(pp
->x
.n
, pp
->x
.n
, r
);
585 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
);
586 pp
= &pp
->x
.p
->arr
[0]) {
587 f
= &pp
->x
.p
->arr
[1];
588 assert(value_pos_p(f
->d
));
589 mpz_mul_ui(twice
, f
->x
.n
, 2);
590 if (value_lt(twice
, f
->d
))
592 if (value_eq(twice
, f
->d
))
600 value_set_si(v
.d
, 0);
601 v
.x
.p
= new_enode(fractional
, 3, -1);
602 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
603 poly_denom(&p
->arr
[0], &twice
);
604 value_assign(v
.x
.p
->arr
[1].d
, twice
);
605 value_decrement(v
.x
.p
->arr
[1].x
.n
, twice
);
606 evalue_set_si(&v
.x
.p
->arr
[2], -1, 1);
607 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
609 for (pp
= &v
.x
.p
->arr
[0]; value_zero_p(pp
->d
);
610 pp
= &pp
->x
.p
->arr
[0]) {
611 f
= &pp
->x
.p
->arr
[1];
612 value_oppose(f
->x
.n
, f
->x
.n
);
613 mpz_fdiv_r(f
->x
.n
, f
->x
.n
, f
->d
);
615 value_division(pp
->d
, twice
, pp
->d
);
616 value_multiply(pp
->x
.n
, pp
->x
.n
, pp
->d
);
617 value_assign(pp
->d
, twice
);
618 value_oppose(pp
->x
.n
, pp
->x
.n
);
619 value_decrement(pp
->x
.n
, pp
->x
.n
);
620 mpz_fdiv_r(pp
->x
.n
, pp
->x
.n
, pp
->d
);
622 /* Maybe we should do this during reduction of
625 Gcd(pp
->d
, pp
->x
.n
, &twice
);
626 value_division(pp
->d
, pp
->d
, twice
);
627 value_division(pp
->x
.n
, pp
->x
.n
, twice
);
637 reorder_terms(p
, &v
);
638 _reduce_evalue(&p
->arr
[1], s
, fract
);
641 /* Try to reduce the degree */
642 for (i
=p
->size
-1;i
>=2;i
--) {
643 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
645 /* Zero coefficient */
646 free_evalue_refs(&(p
->arr
[i
]));
651 /* Try to reduce its strength */
654 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
655 free_evalue_refs(&(p
->arr
[0]));
659 else if (p
->type
== flooring
) {
660 /* Try to reduce the degree */
661 for (i
=p
->size
-1;i
>=2;i
--) {
662 if (!EVALUE_IS_ZERO(p
->arr
[i
]))
664 /* Zero coefficient */
665 free_evalue_refs(&(p
->arr
[i
]));
670 /* Try to reduce its strength */
673 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
674 free_evalue_refs(&(p
->arr
[0]));
678 else if (p
->type
== relation
) {
679 if (p
->size
== 3 && eequal(&p
->arr
[1], &p
->arr
[2])) {
680 free_evalue_refs(&(p
->arr
[2]));
681 free_evalue_refs(&(p
->arr
[0]));
688 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
689 free_evalue_refs(&(p
->arr
[2]));
692 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
693 free_evalue_refs(&(p
->arr
[1]));
694 free_evalue_refs(&(p
->arr
[0]));
695 evalue_set_si(e
, 0, 1);
702 /* Relation was reduced by means of an identical
703 * inequality => remove
705 if (value_zero_p(m
->d
) && !EVALUE_IS_ZERO(m
->x
.p
->arr
[1]))
708 if (reduced
|| value_notzero_p(p
->arr
[0].d
)) {
709 if (!reduced
&& value_zero_p(p
->arr
[0].x
.n
)) {
711 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
713 free_evalue_refs(&(p
->arr
[2]));
717 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
719 evalue_set_si(e
, 0, 1);
720 free_evalue_refs(&(p
->arr
[1]));
722 free_evalue_refs(&(p
->arr
[0]));
728 } /* reduce_evalue */
730 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
735 for (k
= 0; k
< dim
; ++k
)
736 if (value_notzero_p(row
[k
+1]))
739 Vector_Normalize_Positive(row
+1, dim
+1, k
);
740 assert(s
->n
< s
->max
);
741 value_init(s
->fixed
[s
->n
].d
);
742 value_init(s
->fixed
[s
->n
].m
);
743 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
744 s
->fixed
[s
->n
].pos
= k
+1;
745 value_set_si(s
->fixed
[s
->n
].m
, 0);
746 r
= &s
->fixed
[s
->n
].s
;
748 for (l
= k
+1; l
< dim
; ++l
)
749 if (value_notzero_p(row
[l
+1])) {
750 value_set_si(r
->d
, 0);
751 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
752 value_init(r
->x
.p
->arr
[1].x
.n
);
753 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
754 value_set_si(r
->x
.p
->arr
[1].d
, 1);
758 value_oppose(r
->x
.n
, row
[dim
+1]);
759 value_set_si(r
->d
, 1);
763 void reduce_evalue (evalue
*e
) {
764 struct subst s
= { NULL
, 0, 0 };
766 if (value_notzero_p(e
->d
))
767 return; /* a rational number, its already reduced */
769 if (e
->x
.p
->type
== partition
) {
772 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
773 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
775 /* This shouldn't really happen;
776 * Empty domains should not be added.
783 D
= DomainConvex(D
, 0);
784 if (!D
->next
&& D
->NbEq
) {
788 realloc_substitution(&s
, dim
);
790 int d
= relations_depth(&e
->x
.p
->arr
[2*i
+1]);
792 NALLOC(s
.fixed
, s
.max
);
795 for (j
= 0; j
< D
->NbEq
; ++j
)
796 add_substitution(&s
, D
->Constraint
[j
], dim
);
798 if (D
!= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))
800 _reduce_evalue(&e
->x
.p
->arr
[2*i
+1], &s
, 0);
801 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
803 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
804 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
805 value_clear(e
->x
.p
->arr
[2*i
].d
);
807 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
808 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
813 for (j
= 0; j
< s
.n
; ++j
) {
814 value_clear(s
.fixed
[j
].d
);
815 value_clear(s
.fixed
[j
].m
);
816 free_evalue_refs(&s
.fixed
[j
].s
);
820 if (e
->x
.p
->size
== 0) {
822 evalue_set_si(e
, 0, 1);
825 _reduce_evalue(e
, &s
, 0);
830 void print_evalue(FILE *DST
,evalue
*e
,char **pname
) {
832 if(value_notzero_p(e
->d
)) {
833 if(value_notone_p(e
->d
)) {
834 value_print(DST
,VALUE_FMT
,e
->x
.n
);
836 value_print(DST
,VALUE_FMT
,e
->d
);
839 value_print(DST
,VALUE_FMT
,e
->x
.n
);
843 print_enode(DST
,e
->x
.p
,pname
);
847 void print_enode(FILE *DST
,enode
*p
,char **pname
) {
852 fprintf(DST
, "NULL");
858 for (i
=0; i
<p
->size
; i
++) {
859 print_evalue(DST
, &p
->arr
[i
], pname
);
863 fprintf(DST
, " }\n");
867 for (i
=p
->size
-1; i
>=0; i
--) {
868 print_evalue(DST
, &p
->arr
[i
], pname
);
869 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
871 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
873 fprintf(DST
, " )\n");
877 for (i
=0; i
<p
->size
; i
++) {
878 print_evalue(DST
, &p
->arr
[i
], pname
);
879 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
881 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
886 for (i
=p
->size
-1; i
>=1; i
--) {
887 print_evalue(DST
, &p
->arr
[i
], pname
);
890 fprintf(DST
, p
->type
== flooring
? "[" : "{");
891 print_evalue(DST
, &p
->arr
[0], pname
);
892 fprintf(DST
, p
->type
== flooring
? "]" : "}");
894 fprintf(DST
, "^%d + ", i
-1);
899 fprintf(DST
, " )\n");
903 print_evalue(DST
, &p
->arr
[0], pname
);
904 fprintf(DST
, "= 0 ] * \n");
905 print_evalue(DST
, &p
->arr
[1], pname
);
907 fprintf(DST
, " +\n [ ");
908 print_evalue(DST
, &p
->arr
[0], pname
);
909 fprintf(DST
, "!= 0 ] * \n");
910 print_evalue(DST
, &p
->arr
[2], pname
);
914 char **names
= pname
;
915 int maxdim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
916 if (p
->pos
< maxdim
) {
917 NALLOC(names
, maxdim
);
918 for (i
= 0; i
< p
->pos
; ++i
)
920 for ( ; i
< maxdim
; ++i
) {
921 NALLOC(names
[i
], 10);
922 snprintf(names
[i
], 10, "_p%d", i
);
926 for (i
=0; i
<p
->size
/2; i
++) {
927 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), names
);
928 print_evalue(DST
, &p
->arr
[2*i
+1], names
);
931 if (p
->pos
< maxdim
) {
932 for (i
= p
->pos
; i
< maxdim
; ++i
)
945 static void eadd_rev(evalue
*e1
, evalue
*res
)
949 evalue_copy(&ev
, e1
);
951 free_evalue_refs(res
);
955 static void eadd_rev_cst (evalue
*e1
, evalue
*res
)
959 evalue_copy(&ev
, e1
);
960 eadd(res
, &ev
.x
.p
->arr
[type_offset(ev
.x
.p
)]);
961 free_evalue_refs(res
);
965 struct section
{ Polyhedron
* D
; evalue E
; };
967 void eadd_partitions (evalue
*e1
,evalue
*res
)
972 s
= (struct section
*)
973 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
974 sizeof(struct section
));
976 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
977 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
978 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
981 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
982 assert(res
->x
.p
->size
>= 2);
983 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
984 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
986 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
988 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
997 value_init(s
[n
].E
.d
);
998 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
1002 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1003 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
1004 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1006 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1007 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1013 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
1014 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1016 value_init(s
[n
].E
.d
);
1017 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1018 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1023 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
1027 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1030 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
1031 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1032 value_clear(res
->x
.p
->arr
[2*i
].d
);
1037 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1038 for (j
= 0; j
< n
; ++j
) {
1039 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1040 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1041 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1042 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1048 static void explicit_complement(evalue
*res
)
1050 enode
*rel
= new_enode(relation
, 3, 0);
1052 value_clear(rel
->arr
[0].d
);
1053 rel
->arr
[0] = res
->x
.p
->arr
[0];
1054 value_clear(rel
->arr
[1].d
);
1055 rel
->arr
[1] = res
->x
.p
->arr
[1];
1056 value_set_si(rel
->arr
[2].d
, 1);
1057 value_init(rel
->arr
[2].x
.n
);
1058 value_set_si(rel
->arr
[2].x
.n
, 0);
1063 void eadd(evalue
*e1
,evalue
*res
) {
1066 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1067 /* Add two rational numbers */
1073 value_multiply(m1
,e1
->x
.n
,res
->d
);
1074 value_multiply(m2
,res
->x
.n
,e1
->d
);
1075 value_addto(res
->x
.n
,m1
,m2
);
1076 value_multiply(res
->d
,e1
->d
,res
->d
);
1077 Gcd(res
->x
.n
,res
->d
,&g
);
1078 if (value_notone_p(g
)) {
1079 value_division(res
->d
,res
->d
,g
);
1080 value_division(res
->x
.n
,res
->x
.n
,g
);
1082 value_clear(g
); value_clear(m1
); value_clear(m2
);
1085 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
1086 switch (res
->x
.p
->type
) {
1088 /* Add the constant to the constant term of a polynomial*/
1089 eadd(e1
, &res
->x
.p
->arr
[0]);
1092 /* Add the constant to all elements of a periodic number */
1093 for (i
=0; i
<res
->x
.p
->size
; i
++) {
1094 eadd(e1
, &res
->x
.p
->arr
[i
]);
1098 fprintf(stderr
, "eadd: cannot add const with vector\n");
1102 eadd(e1
, &res
->x
.p
->arr
[1]);
1105 assert(EVALUE_IS_ZERO(*e1
));
1106 break; /* Do nothing */
1108 /* Create (zero) complement if needed */
1109 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
1110 explicit_complement(res
);
1111 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1112 eadd(e1
, &res
->x
.p
->arr
[i
]);
1118 /* add polynomial or periodic to constant
1119 * you have to exchange e1 and res, before doing addition */
1121 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
1125 else { // ((e1->d==0) && (res->d==0))
1126 assert(!((e1
->x
.p
->type
== partition
) ^
1127 (res
->x
.p
->type
== partition
)));
1128 if (e1
->x
.p
->type
== partition
) {
1129 eadd_partitions(e1
, res
);
1132 if (e1
->x
.p
->type
== relation
&&
1133 (res
->x
.p
->type
!= relation
||
1134 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1138 if (res
->x
.p
->type
== relation
) {
1139 if (e1
->x
.p
->type
== relation
&&
1140 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1141 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1142 explicit_complement(res
);
1143 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1144 explicit_complement(e1
);
1145 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1146 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1149 if (res
->x
.p
->size
< 3)
1150 explicit_complement(res
);
1151 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1152 eadd(e1
, &res
->x
.p
->arr
[i
]);
1155 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
1156 /* adding to evalues of different type. two cases are possible
1157 * res is periodic and e1 is polynomial, you have to exchange
1158 * e1 and res then to add e1 to the constant term of res */
1159 if (e1
->x
.p
->type
== polynomial
) {
1160 eadd_rev_cst(e1
, res
);
1162 else if (res
->x
.p
->type
== polynomial
) {
1163 /* res is polynomial and e1 is periodic,
1164 add e1 to the constant term of res */
1166 eadd(e1
,&res
->x
.p
->arr
[0]);
1172 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
1173 ((res
->x
.p
->type
== fractional
||
1174 res
->x
.p
->type
== flooring
) &&
1175 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
1176 /* adding evalues of different position (i.e function of different unknowns
1177 * to case are possible */
1179 switch (res
->x
.p
->type
) {
1182 if(mod_term_smaller(res
, e1
))
1183 eadd(e1
,&res
->x
.p
->arr
[1]);
1185 eadd_rev_cst(e1
, res
);
1187 case polynomial
: // res and e1 are polynomials
1188 // add e1 to the constant term of res
1190 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1191 eadd(e1
,&res
->x
.p
->arr
[0]);
1193 eadd_rev_cst(e1
, res
);
1194 // value_clear(g); value_clear(m1); value_clear(m2);
1196 case periodic
: // res and e1 are pointers to periodic numbers
1197 //add e1 to all elements of res
1199 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1200 for (i
=0;i
<res
->x
.p
->size
;i
++) {
1201 eadd(e1
,&res
->x
.p
->arr
[i
]);
1212 //same type , same pos and same size
1213 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
1214 // add any element in e1 to the corresponding element in res
1215 i
= type_offset(res
->x
.p
);
1217 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
1218 for (; i
<res
->x
.p
->size
; i
++) {
1219 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1224 /* Sizes are different */
1225 switch(res
->x
.p
->type
) {
1229 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1230 /* new enode and add res to that new node. If you do not do */
1231 /* that, you lose the the upper weight part of e1 ! */
1233 if(e1
->x
.p
->size
> res
->x
.p
->size
)
1236 i
= type_offset(res
->x
.p
);
1238 assert(eequal(&e1
->x
.p
->arr
[0],
1239 &res
->x
.p
->arr
[0]));
1240 for (; i
<e1
->x
.p
->size
; i
++) {
1241 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1248 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1251 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1252 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1253 to add periodicaly elements of e1 to elements of ne, and finaly to
1258 value_init(ex
); value_init(ey
);value_init(ep
);
1261 value_set_si(ex
,e1
->x
.p
->size
);
1262 value_set_si(ey
,res
->x
.p
->size
);
1263 value_assign (ep
,*Lcm(ex
,ey
));
1264 p
=(int)mpz_get_si(ep
);
1265 ne
= (evalue
*) malloc (sizeof(evalue
));
1267 value_set_si( ne
->d
,0);
1269 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
1271 evalue_copy(&ne
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
%y
]);
1274 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
1277 value_assign(res
->d
, ne
->d
);
1283 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
1292 static void emul_rev (evalue
*e1
, evalue
*res
)
1296 evalue_copy(&ev
, e1
);
1298 free_evalue_refs(res
);
1302 static void emul_poly (evalue
*e1
, evalue
*res
)
1304 int i
, j
, o
= type_offset(res
->x
.p
);
1306 int size
=(e1
->x
.p
->size
+ res
->x
.p
->size
- o
- 1);
1308 value_set_si(tmp
.d
,0);
1309 tmp
.x
.p
=new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
1311 evalue_copy(&tmp
.x
.p
->arr
[0], &e1
->x
.p
->arr
[0]);
1312 for (i
=o
; i
< e1
->x
.p
->size
; i
++) {
1313 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
1314 emul(&res
->x
.p
->arr
[o
], &tmp
.x
.p
->arr
[i
]);
1317 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
1318 for (i
=o
+1; i
<res
->x
.p
->size
; i
++)
1319 for (j
=o
; j
<e1
->x
.p
->size
; j
++) {
1322 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
1323 emul(&res
->x
.p
->arr
[i
], &ev
);
1324 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-o
]);
1325 free_evalue_refs(&ev
);
1327 free_evalue_refs(res
);
1331 void emul_partitions (evalue
*e1
,evalue
*res
)
1336 s
= (struct section
*)
1337 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1338 sizeof(struct section
));
1340 assert(e1
->x
.p
->pos
== res
->x
.p
->pos
);
1341 assert(e1
->x
.p
->pos
== EVALUE_DOMAIN(e1
->x
.p
->arr
[0])->Dimension
);
1342 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1345 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1346 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1347 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1348 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1354 /* This code is only needed because the partitions
1355 are not true partitions.
1357 for (k
= 0; k
< n
; ++k
) {
1358 if (DomainIncludes(s
[k
].D
, d
))
1360 if (DomainIncludes(d
, s
[k
].D
)) {
1361 Domain_Free(s
[k
].D
);
1362 free_evalue_refs(&s
[k
].E
);
1373 value_init(s
[n
].E
.d
);
1374 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1375 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1379 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1380 value_clear(res
->x
.p
->arr
[2*i
].d
);
1381 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1386 evalue_set_si(res
, 0, 1);
1388 res
->x
.p
= new_enode(partition
, 2*n
, e1
->x
.p
->pos
);
1389 for (j
= 0; j
< n
; ++j
) {
1390 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1391 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1392 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1393 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1400 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1402 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1403 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1404 * evalues is not treated here */
1406 void emul (evalue
*e1
, evalue
*res
){
1409 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1410 fprintf(stderr
, "emul: do not proced on evector type !\n");
1414 if (EVALUE_IS_ZERO(*res
))
1417 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1418 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1419 emul_partitions(e1
, res
);
1422 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1423 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1424 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1426 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1427 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1428 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1429 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1430 explicit_complement(res
);
1431 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1432 explicit_complement(e1
);
1433 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1434 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1437 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1438 emul(e1
, &res
->x
.p
->arr
[i
]);
1440 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1441 switch(e1
->x
.p
->type
) {
1443 switch(res
->x
.p
->type
) {
1445 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1446 /* Product of two polynomials of the same variable */
1451 /* Product of two polynomials of different variables */
1453 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1454 for( i
=0; i
<res
->x
.p
->size
; i
++)
1455 emul(e1
, &res
->x
.p
->arr
[i
]);
1464 /* Product of a polynomial and a periodic or fractional */
1471 switch(res
->x
.p
->type
) {
1473 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1474 /* Product of two periodics of the same parameter and period */
1476 for(i
=0; i
<res
->x
.p
->size
;i
++)
1477 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1482 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1483 /* Product of two periodics of the same parameter and different periods */
1487 value_init(x
); value_init(y
);value_init(z
);
1490 value_set_si(x
,e1
->x
.p
->size
);
1491 value_set_si(y
,res
->x
.p
->size
);
1492 value_assign (z
,*Lcm(x
,y
));
1493 lcm
=(int)mpz_get_si(z
);
1494 newp
= (evalue
*) malloc (sizeof(evalue
));
1495 value_init(newp
->d
);
1496 value_set_si( newp
->d
,0);
1497 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1498 for(i
=0;i
<lcm
;i
++) {
1499 evalue_copy(&newp
->x
.p
->arr
[i
],
1500 &res
->x
.p
->arr
[i
%iy
]);
1503 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1505 value_assign(res
->d
,newp
->d
);
1508 value_clear(x
); value_clear(y
);value_clear(z
);
1512 /* Product of two periodics of different parameters */
1514 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1515 for(i
=0; i
<res
->x
.p
->size
; i
++)
1516 emul(e1
, &(res
->x
.p
->arr
[i
]));
1524 /* Product of a periodic and a polynomial */
1526 for(i
=0; i
<res
->x
.p
->size
; i
++)
1527 emul(e1
, &(res
->x
.p
->arr
[i
]));
1534 switch(res
->x
.p
->type
) {
1536 for(i
=0; i
<res
->x
.p
->size
; i
++)
1537 emul(e1
, &(res
->x
.p
->arr
[i
]));
1544 assert(e1
->x
.p
->type
== res
->x
.p
->type
);
1545 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1546 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1549 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1550 if (e1
->x
.p
->type
!= fractional
|| !value_two_p(d
.d
))
1555 value_set_si(d
.x
.n
, 1);
1556 /* { x }^2 == { x }/2 */
1557 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1558 assert(e1
->x
.p
->size
== 3);
1559 assert(res
->x
.p
->size
== 3);
1561 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1563 eadd(&res
->x
.p
->arr
[1], &tmp
);
1564 emul(&e1
->x
.p
->arr
[2], &tmp
);
1565 emul(&e1
->x
.p
->arr
[1], res
);
1566 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1567 free_evalue_refs(&tmp
);
1572 if(mod_term_smaller(res
, e1
))
1573 for(i
=1; i
<res
->x
.p
->size
; i
++)
1574 emul(e1
, &(res
->x
.p
->arr
[i
]));
1589 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1590 /* Product of two rational numbers */
1594 value_multiply(res
->d
,e1
->d
,res
->d
);
1595 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1596 Gcd(res
->x
.n
, res
->d
,&g
);
1597 if (value_notone_p(g
)) {
1598 value_division(res
->d
,res
->d
,g
);
1599 value_division(res
->x
.n
,res
->x
.n
,g
);
1605 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1606 /* Product of an expression (polynomial or peririodic) and a rational number */
1612 /* Product of a rationel number and an expression (polynomial or peririodic) */
1614 i
= type_offset(res
->x
.p
);
1615 for (; i
<res
->x
.p
->size
; i
++)
1616 emul(e1
, &res
->x
.p
->arr
[i
]);
1626 /* Frees mask content ! */
1627 void emask(evalue
*mask
, evalue
*res
) {
1634 if (EVALUE_IS_ZERO(*res
)) {
1635 free_evalue_refs(mask
);
1639 assert(value_zero_p(mask
->d
));
1640 assert(mask
->x
.p
->type
== partition
);
1641 assert(value_zero_p(res
->d
));
1642 assert(res
->x
.p
->type
== partition
);
1643 assert(mask
->x
.p
->pos
== res
->x
.p
->pos
);
1644 assert(res
->x
.p
->pos
== EVALUE_DOMAIN(res
->x
.p
->arr
[0])->Dimension
);
1645 assert(mask
->x
.p
->pos
== EVALUE_DOMAIN(mask
->x
.p
->arr
[0])->Dimension
);
1646 pos
= res
->x
.p
->pos
;
1648 s
= (struct section
*)
1649 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1650 sizeof(struct section
));
1654 evalue_set_si(&mone
, -1, 1);
1657 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1658 assert(mask
->x
.p
->size
>= 2);
1659 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1660 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1662 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1664 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1673 value_init(s
[n
].E
.d
);
1674 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1678 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1679 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1682 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1683 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1684 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1685 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1687 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1688 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1694 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1695 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1697 value_init(s
[n
].E
.d
);
1698 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1699 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1705 /* Just ignore; this may have been previously masked off */
1707 if (fd
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1711 free_evalue_refs(&mone
);
1712 free_evalue_refs(mask
);
1713 free_evalue_refs(res
);
1716 evalue_set_si(res
, 0, 1);
1718 res
->x
.p
= new_enode(partition
, 2*n
, pos
);
1719 for (j
= 0; j
< n
; ++j
) {
1720 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1721 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1722 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1729 void evalue_copy(evalue
*dst
, evalue
*src
)
1731 value_assign(dst
->d
, src
->d
);
1732 if(value_notzero_p(src
->d
)) {
1733 value_init(dst
->x
.n
);
1734 value_assign(dst
->x
.n
, src
->x
.n
);
1736 dst
->x
.p
= ecopy(src
->x
.p
);
1739 enode
*new_enode(enode_type type
,int size
,int pos
) {
1745 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1748 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1752 for(i
=0; i
<size
; i
++) {
1753 value_init(res
->arr
[i
].d
);
1754 value_set_si(res
->arr
[i
].d
,0);
1755 res
->arr
[i
].x
.p
= 0;
1760 enode
*ecopy(enode
*e
) {
1765 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1766 for(i
=0;i
<e
->size
;++i
) {
1767 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1768 if(value_zero_p(res
->arr
[i
].d
))
1769 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1770 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1771 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1773 value_init(res
->arr
[i
].x
.n
);
1774 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1780 int ecmp(const evalue
*e1
, const evalue
*e2
)
1786 if (value_notzero_p(e1
->d
) && value_notzero_p(e2
->d
)) {
1790 value_multiply(m
, e1
->x
.n
, e2
->d
);
1791 value_multiply(m2
, e2
->x
.n
, e1
->d
);
1793 if (value_lt(m
, m2
))
1795 else if (value_gt(m
, m2
))
1805 if (value_notzero_p(e1
->d
))
1807 if (value_notzero_p(e2
->d
))
1813 if (p1
->type
!= p2
->type
)
1814 return p1
->type
- p2
->type
;
1815 if (p1
->pos
!= p2
->pos
)
1816 return p1
->pos
- p2
->pos
;
1817 if (p1
->size
!= p2
->size
)
1818 return p1
->size
- p2
->size
;
1820 for (i
= p1
->size
-1; i
>= 0; --i
)
1821 if ((r
= ecmp(&p1
->arr
[i
], &p2
->arr
[i
])) != 0)
1827 int eequal(evalue
*e1
,evalue
*e2
) {
1832 if (value_ne(e1
->d
,e2
->d
))
1835 /* e1->d == e2->d */
1836 if (value_notzero_p(e1
->d
)) {
1837 if (value_ne(e1
->x
.n
,e2
->x
.n
))
1840 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1844 /* e1->d == e2->d == 0 */
1847 if (p1
->type
!= p2
->type
) return 0;
1848 if (p1
->size
!= p2
->size
) return 0;
1849 if (p1
->pos
!= p2
->pos
) return 0;
1850 for (i
=0; i
<p1
->size
; i
++)
1851 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
1856 void free_evalue_refs(evalue
*e
) {
1861 if (EVALUE_IS_DOMAIN(*e
)) {
1862 Domain_Free(EVALUE_DOMAIN(*e
));
1865 } else if (value_pos_p(e
->d
)) {
1867 /* 'e' stores a constant */
1869 value_clear(e
->x
.n
);
1872 assert(value_zero_p(e
->d
));
1875 if (!p
) return; /* null pointer */
1876 for (i
=0; i
<p
->size
; i
++) {
1877 free_evalue_refs(&(p
->arr
[i
]));
1881 } /* free_evalue_refs */
1883 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
1884 Vector
* val
, evalue
*res
)
1886 unsigned nparam
= periods
->Size
;
1889 double d
= compute_evalue(e
, val
->p
);
1890 d
*= VALUE_TO_DOUBLE(m
);
1895 value_assign(res
->d
, m
);
1896 value_init(res
->x
.n
);
1897 value_set_double(res
->x
.n
, d
);
1898 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
1901 if (value_one_p(periods
->p
[p
]))
1902 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
1907 value_assign(tmp
, periods
->p
[p
]);
1908 value_set_si(res
->d
, 0);
1909 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
1911 value_decrement(tmp
, tmp
);
1912 value_assign(val
->p
[p
], tmp
);
1913 mod2table_r(e
, periods
, m
, p
+1, val
,
1914 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
1915 } while (value_pos_p(tmp
));
1921 static void rel2table(evalue
*e
, int zero
)
1923 if (value_pos_p(e
->d
)) {
1924 if (value_zero_p(e
->x
.n
) == zero
)
1925 value_set_si(e
->x
.n
, 1);
1927 value_set_si(e
->x
.n
, 0);
1928 value_set_si(e
->d
, 1);
1931 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
1932 rel2table(&e
->x
.p
->arr
[i
], zero
);
1936 void evalue_mod2table(evalue
*e
, int nparam
)
1941 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
1944 for (i
=0; i
<p
->size
; i
++) {
1945 evalue_mod2table(&(p
->arr
[i
]), nparam
);
1947 if (p
->type
== relation
) {
1952 evalue_copy(©
, &p
->arr
[0]);
1954 rel2table(&p
->arr
[0], 1);
1955 emul(&p
->arr
[0], &p
->arr
[1]);
1957 rel2table(©
, 0);
1958 emul(©
, &p
->arr
[2]);
1959 eadd(&p
->arr
[2], &p
->arr
[1]);
1960 free_evalue_refs(&p
->arr
[2]);
1961 free_evalue_refs(©
);
1963 free_evalue_refs(&p
->arr
[0]);
1967 } else if (p
->type
== fractional
) {
1968 Vector
*periods
= Vector_Alloc(nparam
);
1969 Vector
*val
= Vector_Alloc(nparam
);
1975 value_set_si(tmp
, 1);
1976 Vector_Set(periods
->p
, 1, nparam
);
1977 Vector_Set(val
->p
, 0, nparam
);
1978 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
1981 assert(p
->type
== polynomial
);
1982 assert(p
->size
== 2);
1983 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
1984 value_lcm(tmp
, p
->arr
[1].d
, &tmp
);
1986 value_lcm(tmp
, ev
->d
, &tmp
);
1988 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
1991 evalue_set_si(&res
, 0, 1);
1992 /* Compute the polynomial using Horner's rule */
1993 for (i
=p
->size
-1;i
>1;i
--) {
1994 eadd(&p
->arr
[i
], &res
);
1997 eadd(&p
->arr
[1], &res
);
1999 free_evalue_refs(e
);
2000 free_evalue_refs(&EP
);
2005 Vector_Free(periods
);
2007 } /* evalue_mod2table */
2009 /********************************************************/
2010 /* function in domain */
2011 /* check if the parameters in list_args */
2012 /* verifies the constraints of Domain P */
2013 /********************************************************/
2014 int in_domain(Polyhedron
*P
, Value
*list_args
) {
2017 Value v
; /* value of the constraint of a row when
2018 parameters are instanciated*/
2024 /*P->Constraint constraint matrice of polyhedron P*/
2025 for(row
=0;row
<P
->NbConstraints
;row
++) {
2026 value_assign(v
,P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
2027 for(col
=1;col
<P
->Dimension
+1;col
++) {
2028 value_multiply(tmp
,P
->Constraint
[row
][col
],list_args
[col
-1]);
2029 value_addto(v
,v
,tmp
);
2031 if (value_notzero_p(P
->Constraint
[row
][0])) {
2033 /*if v is not >=0 then this constraint is not respected */
2034 if (value_neg_p(v
)) {
2038 return P
->next
? in_domain(P
->next
, list_args
) : 0;
2043 /*if v is not = 0 then this constraint is not respected */
2044 if (value_notzero_p(v
))
2049 /*if not return before this point => all
2050 the constraints are respected */
2056 /****************************************************/
2057 /* function compute enode */
2058 /* compute the value of enode p with parameters */
2059 /* list "list_args */
2060 /* compute the polynomial or the periodic */
2061 /****************************************************/
2063 static double compute_enode(enode
*p
, Value
*list_args
) {
2075 if (p
->type
== polynomial
) {
2077 value_assign(param
,list_args
[p
->pos
-1]);
2079 /* Compute the polynomial using Horner's rule */
2080 for (i
=p
->size
-1;i
>0;i
--) {
2081 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2082 res
*=VALUE_TO_DOUBLE(param
);
2084 res
+=compute_evalue(&p
->arr
[0],list_args
);
2086 else if (p
->type
== fractional
) {
2087 double d
= compute_evalue(&p
->arr
[0], list_args
);
2088 d
-= floor(d
+1e-10);
2090 /* Compute the polynomial using Horner's rule */
2091 for (i
=p
->size
-1;i
>1;i
--) {
2092 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2095 res
+=compute_evalue(&p
->arr
[1],list_args
);
2097 else if (p
->type
== flooring
) {
2098 double d
= compute_evalue(&p
->arr
[0], list_args
);
2101 /* Compute the polynomial using Horner's rule */
2102 for (i
=p
->size
-1;i
>1;i
--) {
2103 res
+=compute_evalue(&p
->arr
[i
],list_args
);
2106 res
+=compute_evalue(&p
->arr
[1],list_args
);
2108 else if (p
->type
== periodic
) {
2109 value_assign(m
,list_args
[p
->pos
-1]);
2111 /* Choose the right element of the periodic */
2112 value_set_si(param
,p
->size
);
2113 value_pmodulus(m
,m
,param
);
2114 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
2116 else if (p
->type
== relation
) {
2117 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
2118 res
= compute_evalue(&p
->arr
[1], list_args
);
2119 else if (p
->size
> 2)
2120 res
= compute_evalue(&p
->arr
[2], list_args
);
2122 else if (p
->type
== partition
) {
2123 int dim
= EVALUE_DOMAIN(p
->arr
[0])->Dimension
;
2124 Value
*vals
= list_args
;
2127 for (i
= 0; i
< dim
; ++i
) {
2128 value_init(vals
[i
]);
2130 value_assign(vals
[i
], list_args
[i
]);
2133 for (i
= 0; i
< p
->size
/2; ++i
)
2134 if (DomainContains(EVALUE_DOMAIN(p
->arr
[2*i
]), vals
, p
->pos
, 0, 1)) {
2135 res
= compute_evalue(&p
->arr
[2*i
+1], vals
);
2139 for (i
= 0; i
< dim
; ++i
)
2140 value_clear(vals
[i
]);
2149 } /* compute_enode */
2151 /*************************************************/
2152 /* return the value of Ehrhart Polynomial */
2153 /* It returns a double, because since it is */
2154 /* a recursive function, some intermediate value */
2155 /* might not be integral */
2156 /*************************************************/
2158 double compute_evalue(evalue
*e
,Value
*list_args
) {
2162 if (value_notzero_p(e
->d
)) {
2163 if (value_notone_p(e
->d
))
2164 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
2166 res
= VALUE_TO_DOUBLE(e
->x
.n
);
2169 res
= compute_enode(e
->x
.p
,list_args
);
2171 } /* compute_evalue */
2174 /****************************************************/
2175 /* function compute_poly : */
2176 /* Check for the good validity domain */
2177 /* return the number of point in the Polyhedron */
2178 /* in allocated memory */
2179 /* Using the Ehrhart pseudo-polynomial */
2180 /****************************************************/
2181 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
2184 /* double d; int i; */
2186 tmp
= (Value
*) malloc (sizeof(Value
));
2187 assert(tmp
!= NULL
);
2189 value_set_si(*tmp
,0);
2192 return(tmp
); /* no ehrhart polynomial */
2193 if(en
->ValidityDomain
) {
2194 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
2195 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2200 return(tmp
); /* no Validity Domain */
2202 if(in_domain(en
->ValidityDomain
,list_args
)) {
2204 #ifdef EVAL_EHRHART_DEBUG
2205 Print_Domain(stdout
,en
->ValidityDomain
);
2206 print_evalue(stdout
,&en
->EP
);
2209 /* d = compute_evalue(&en->EP,list_args);
2211 printf("(double)%lf = %d\n", d, i ); */
2212 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
2218 value_set_si(*tmp
,0);
2219 return(tmp
); /* no compatible domain with the arguments */
2220 } /* compute_poly */
2222 size_t value_size(Value v
) {
2223 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
2224 * sizeof(v
[0]._mp_d
[0]);
2227 size_t domain_size(Polyhedron
*D
)
2230 size_t s
= sizeof(*D
);
2232 for (i
= 0; i
< D
->NbConstraints
; ++i
)
2233 for (j
= 0; j
< D
->Dimension
+2; ++j
)
2234 s
+= value_size(D
->Constraint
[i
][j
]);
2237 for (i = 0; i < D->NbRays; ++i)
2238 for (j = 0; j < D->Dimension+2; ++j)
2239 s += value_size(D->Ray[i][j]);
2242 return D
->next
? s
+domain_size(D
->next
) : s
;
2245 size_t enode_size(enode
*p
) {
2246 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
2249 if (p
->type
== partition
)
2250 for (i
= 0; i
< p
->size
/2; ++i
) {
2251 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
2252 s
+= evalue_size(&p
->arr
[2*i
+1]);
2255 for (i
= 0; i
< p
->size
; ++i
) {
2256 s
+= evalue_size(&p
->arr
[i
]);
2261 size_t evalue_size(evalue
*e
)
2263 size_t s
= sizeof(*e
);
2264 s
+= value_size(e
->d
);
2265 if (value_notzero_p(e
->d
))
2266 s
+= value_size(e
->x
.n
);
2268 s
+= enode_size(e
->x
.p
);
2272 static evalue
*find_second(evalue
*base
, evalue
*cst
, evalue
*e
, Value m
)
2274 evalue
*found
= NULL
;
2279 if (value_pos_p(e
->d
) || e
->x
.p
->type
!= fractional
)
2282 value_init(offset
.d
);
2283 value_init(offset
.x
.n
);
2284 poly_denom(&e
->x
.p
->arr
[0], &offset
.d
);
2285 value_lcm(m
, offset
.d
, &offset
.d
);
2286 value_set_si(offset
.x
.n
, 1);
2289 evalue_copy(©
, cst
);
2292 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2294 if (eequal(base
, &e
->x
.p
->arr
[0]))
2295 found
= &e
->x
.p
->arr
[0];
2297 value_set_si(offset
.x
.n
, -2);
2300 mpz_fdiv_r(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2302 if (eequal(base
, &e
->x
.p
->arr
[0]))
2305 free_evalue_refs(cst
);
2306 free_evalue_refs(&offset
);
2309 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2310 found
= find_second(base
, cst
, &e
->x
.p
->arr
[i
], m
);
2315 static evalue
*find_relation_pair(evalue
*e
)
2318 evalue
*found
= NULL
;
2320 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
2323 if (e
->x
.p
->type
== fractional
) {
2328 poly_denom(&e
->x
.p
->arr
[0], &m
);
2330 for (cst
= &e
->x
.p
->arr
[0]; value_zero_p(cst
->d
);
2331 cst
= &cst
->x
.p
->arr
[0])
2334 for (i
= 1; !found
&& i
< e
->x
.p
->size
; ++i
)
2335 found
= find_second(&e
->x
.p
->arr
[0], cst
, &e
->x
.p
->arr
[i
], m
);
2340 i
= e
->x
.p
->type
== relation
;
2341 for (; !found
&& i
< e
->x
.p
->size
; ++i
)
2342 found
= find_relation_pair(&e
->x
.p
->arr
[i
]);
2347 void evalue_mod2relation(evalue
*e
) {
2350 if (value_zero_p(e
->d
) && e
->x
.p
->type
== partition
) {
2353 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2354 evalue_mod2relation(&e
->x
.p
->arr
[2*i
+1]);
2355 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2356 value_clear(e
->x
.p
->arr
[2*i
].d
);
2357 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2359 if (2*i
< e
->x
.p
->size
) {
2360 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2361 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2366 if (e
->x
.p
->size
== 0) {
2368 evalue_set_si(e
, 0, 1);
2374 while ((d
= find_relation_pair(e
)) != NULL
) {
2378 value_init(split
.d
);
2379 value_set_si(split
.d
, 0);
2380 split
.x
.p
= new_enode(relation
, 3, 0);
2381 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
2382 evalue_set_si(&split
.x
.p
->arr
[2], 1, 1);
2384 ev
= &split
.x
.p
->arr
[0];
2385 value_set_si(ev
->d
, 0);
2386 ev
->x
.p
= new_enode(fractional
, 3, -1);
2387 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
2388 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
2389 evalue_copy(&ev
->x
.p
->arr
[0], d
);
2395 free_evalue_refs(&split
);
2399 static int evalue_comp(const void * a
, const void * b
)
2401 const evalue
*e1
= *(const evalue
**)a
;
2402 const evalue
*e2
= *(const evalue
**)b
;
2403 return ecmp(e1
, e2
);
2406 void evalue_combine(evalue
*e
)
2413 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2416 NALLOC(evs
, e
->x
.p
->size
/2);
2417 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2418 evs
[i
] = &e
->x
.p
->arr
[2*i
+1];
2419 qsort(evs
, e
->x
.p
->size
/2, sizeof(evs
[0]), evalue_comp
);
2420 p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
2421 for (i
= 0, k
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2422 if (k
== 0 || ecmp(&p
->arr
[2*k
-1], evs
[i
]) != 0) {
2423 value_clear(p
->arr
[2*k
].d
);
2424 value_clear(p
->arr
[2*k
+1].d
);
2425 p
->arr
[2*k
] = *(evs
[i
]-1);
2426 p
->arr
[2*k
+1] = *(evs
[i
]);
2429 Polyhedron
*D
= EVALUE_DOMAIN(*(evs
[i
]-1));
2432 value_clear((evs
[i
]-1)->d
);
2436 L
->next
= EVALUE_DOMAIN(p
->arr
[2*k
-2]);
2437 EVALUE_SET_DOMAIN(p
->arr
[2*k
-2], D
);
2438 free_evalue_refs(evs
[i
]);
2442 for (i
= 2*k
; i
< p
->size
; ++i
)
2443 value_clear(p
->arr
[i
].d
);
2450 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2452 if (value_notzero_p(e
->x
.p
->arr
[2*i
+1].d
))
2454 H
= DomainConvex(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
2457 for (k
= 0; k
< e
->x
.p
->size
/2; ++k
) {
2458 Polyhedron
*D
, *N
, **P
;
2461 P
= &EVALUE_DOMAIN(e
->x
.p
->arr
[2*k
]);
2468 if (D
->NbEq
<= H
->NbEq
) {
2474 tmp
.x
.p
= new_enode(partition
, 2, e
->x
.p
->pos
);
2475 EVALUE_SET_DOMAIN(tmp
.x
.p
->arr
[0], Polyhedron_Copy(D
));
2476 evalue_copy(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*i
+1]);
2477 reduce_evalue(&tmp
);
2478 if (value_notzero_p(tmp
.d
) ||
2479 ecmp(&tmp
.x
.p
->arr
[1], &e
->x
.p
->arr
[2*k
+1]) != 0)
2482 D
->next
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2483 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]) = D
;
2486 free_evalue_refs(&tmp
);
2492 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
2494 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
2496 value_clear(e
->x
.p
->arr
[2*i
].d
);
2497 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2499 if (2*i
< e
->x
.p
->size
) {
2500 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2501 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2508 H
= DomainConvex(D
, 0);
2509 E
= DomainDifference(H
, D
, 0);
2511 D
= DomainDifference(H
, E
, 0);
2514 EVALUE_SET_DOMAIN(p
->arr
[2*i
], D
);
2518 /* May change coefficients to become non-standard if fiddle is set
2519 * => reduce p afterwards to correct
2521 static Polyhedron
*polynomial_projection(enode
*p
, Polyhedron
*D
, Value
*d
,
2522 Matrix
**R
, int fiddle
)
2526 unsigned dim
= D
->Dimension
;
2527 Matrix
*T
= Matrix_Alloc(2, dim
+1);
2532 assert(p
->type
== fractional
);
2534 value_set_si(T
->p
[1][dim
], 1);
2536 while (value_zero_p(pp
->d
)) {
2537 assert(pp
->x
.p
->type
== polynomial
);
2538 assert(pp
->x
.p
->size
== 2);
2539 assert(value_notzero_p(pp
->x
.p
->arr
[1].d
));
2540 value_division(T
->p
[0][pp
->x
.p
->pos
-1], *d
, pp
->x
.p
->arr
[1].d
);
2541 mpz_mul_ui(twice
, pp
->x
.p
->arr
[1].x
.n
, 2);
2542 if (fiddle
&& value_gt(twice
, pp
->x
.p
->arr
[1].d
))
2543 value_subtract(pp
->x
.p
->arr
[1].x
.n
,
2544 pp
->x
.p
->arr
[1].x
.n
, pp
->x
.p
->arr
[1].d
);
2545 value_multiply(T
->p
[0][pp
->x
.p
->pos
-1],
2546 T
->p
[0][pp
->x
.p
->pos
-1], pp
->x
.p
->arr
[1].x
.n
);
2547 pp
= &pp
->x
.p
->arr
[0];
2549 value_division(T
->p
[0][dim
], *d
, pp
->d
);
2550 value_multiply(T
->p
[0][dim
], T
->p
[0][dim
], pp
->x
.n
);
2551 I
= DomainImage(D
, T
, 0);
2552 H
= DomainConvex(I
, 0);
2561 int reduce_in_domain(evalue
*e
, Polyhedron
*D
)
2571 if (value_notzero_p(e
->d
))
2576 if (p
->type
== relation
) {
2582 I
= polynomial_projection(p
->arr
[0].x
.p
, D
, &d
, &T
, 1);
2583 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2584 equal
= value_eq(min
, max
);
2585 mpz_cdiv_q(min
, min
, d
);
2586 mpz_fdiv_q(max
, max
, d
);
2588 if (bounded
&& value_gt(min
, max
)) {
2594 evalue_set_si(e
, 0, 1);
2597 free_evalue_refs(&(p
->arr
[1]));
2598 free_evalue_refs(&(p
->arr
[0]));
2604 return r
? r
: reduce_in_domain(e
, D
);
2605 } else if (bounded
&& equal
) {
2608 free_evalue_refs(&(p
->arr
[2]));
2611 free_evalue_refs(&(p
->arr
[0]));
2617 return reduce_in_domain(e
, D
);
2618 } else if (bounded
&& value_eq(min
, max
)) {
2619 /* zero for a single value */
2621 Matrix
*M
= Matrix_Alloc(1, D
->Dimension
+2);
2622 Vector_Copy(T
->p
[0], M
->p
[0]+1, D
->Dimension
+1);
2623 value_multiply(min
, min
, d
);
2624 value_subtract(M
->p
[0][D
->Dimension
+1],
2625 M
->p
[0][D
->Dimension
+1], min
);
2626 E
= DomainAddConstraints(D
, M
, 0);
2632 r
= reduce_in_domain(&p
->arr
[1], E
);
2634 r
|= reduce_in_domain(&p
->arr
[2], D
);
2636 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2644 _reduce_evalue(&p
->arr
[0].x
.p
->arr
[0], 0, 1);
2647 i
= p
->type
== relation
? 1 :
2648 p
->type
== fractional
? 1 : 0;
2649 for (; i
<p
->size
; i
++)
2650 r
|= reduce_in_domain(&p
->arr
[i
], D
);
2652 if (p
->type
!= fractional
) {
2653 if (r
&& p
->type
== polynomial
) {
2656 value_set_si(f
.d
, 0);
2657 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2658 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2659 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2660 reorder_terms(p
, &f
);
2671 I
= polynomial_projection(p
, D
, &d
, &T
, 1);
2673 bounded
= line_minmax(I
, &min
, &max
); /* frees I */
2674 mpz_fdiv_q(min
, min
, d
);
2675 mpz_fdiv_q(max
, max
, d
);
2676 value_subtract(d
, max
, min
);
2678 if (bounded
&& value_eq(min
, max
)) {
2681 value_init(inc
.x
.n
);
2682 value_set_si(inc
.d
, 1);
2683 value_oppose(inc
.x
.n
, min
);
2684 eadd(&inc
, &p
->arr
[0]);
2685 reorder_terms(p
, &p
->arr
[0]); /* frees arr[0] */
2689 free_evalue_refs(&inc
);
2691 } else if (bounded
&& value_one_p(d
) && p
->size
> 3) {
2698 value_set_si(rem
.d
, 0);
2699 rem
.x
.p
= new_enode(fractional
, 3, -1);
2700 evalue_copy(&rem
.x
.p
->arr
[0], &p
->arr
[0]);
2701 rem
.x
.p
->arr
[1] = p
->arr
[1];
2702 rem
.x
.p
->arr
[2] = p
->arr
[2];
2703 for (i
= 3; i
< p
->size
; ++i
)
2704 p
->arr
[i
-2] = p
->arr
[i
];
2708 value_init(inc
.x
.n
);
2709 value_set_si(inc
.d
, 1);
2710 value_oppose(inc
.x
.n
, min
);
2713 evalue_copy(&t
, &p
->arr
[0]);
2717 value_set_si(f
.d
, 0);
2718 f
.x
.p
= new_enode(fractional
, 3, -1);
2719 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
2720 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2721 evalue_set_si(&f
.x
.p
->arr
[2], 2, 1);
2723 value_init(factor
.d
);
2724 evalue_set_si(&factor
, -1, 1);
2730 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2731 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
2737 free_evalue_refs(&inc
);
2738 free_evalue_refs(&t
);
2739 free_evalue_refs(&f
);
2740 free_evalue_refs(&factor
);
2741 free_evalue_refs(&rem
);
2743 reduce_in_domain(e
, D
);
2747 _reduce_evalue(&p
->arr
[0], 0, 1);
2751 value_set_si(f
.d
, 0);
2752 f
.x
.p
= new_enode(fractional
, 3, -1);
2753 value_clear(f
.x
.p
->arr
[0].d
);
2754 f
.x
.p
->arr
[0] = p
->arr
[0];
2755 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2756 evalue_set_si(&f
.x
.p
->arr
[2], 1, 1);
2757 reorder_terms(p
, &f
);
2771 void evalue_range_reduction(evalue
*e
)
2774 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2777 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2778 if (reduce_in_domain(&e
->x
.p
->arr
[2*i
+1],
2779 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))) {
2780 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
2782 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
2783 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
2784 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
2785 value_clear(e
->x
.p
->arr
[2*i
].d
);
2787 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
2788 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
2796 Enumeration
* partition2enumeration(evalue
*EP
)
2799 Enumeration
*en
, *res
= NULL
;
2801 if (EVALUE_IS_ZERO(*EP
)) {
2806 for (i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
2807 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
])->Dimension
);
2808 en
= (Enumeration
*)malloc(sizeof(Enumeration
));
2811 res
->ValidityDomain
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
2812 value_clear(EP
->x
.p
->arr
[2*i
].d
);
2813 res
->EP
= EP
->x
.p
->arr
[2*i
+1];
2821 static int frac2floor_in_domain(evalue
*e
, Polyhedron
*D
)
2831 if (value_notzero_p(e
->d
))
2836 i
= p
->type
== relation
? 1 :
2837 p
->type
== fractional
? 1 : 0;
2838 for (; i
<p
->size
; i
++)
2839 r
|= frac2floor_in_domain(&p
->arr
[i
], D
);
2841 if (p
->type
!= fractional
) {
2842 if (r
&& p
->type
== polynomial
) {
2845 value_set_si(f
.d
, 0);
2846 f
.x
.p
= new_enode(polynomial
, 2, p
->pos
);
2847 evalue_set_si(&f
.x
.p
->arr
[0], 0, 1);
2848 evalue_set_si(&f
.x
.p
->arr
[1], 1, 1);
2849 reorder_terms(p
, &f
);
2858 I
= polynomial_projection(p
, D
, &d
, &T
, 0);
2861 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2864 assert(I
->NbEq
== 0); /* Should have been reduced */
2867 for (i
= 0; i
< I
->NbConstraints
; ++i
)
2868 if (value_pos_p(I
->Constraint
[i
][1]))
2871 assert(i
< I
->NbConstraints
);
2873 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
2874 mpz_cdiv_q(min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
2875 if (value_neg_p(min
)) {
2877 mpz_fdiv_q(min
, min
, d
);
2878 value_init(offset
.d
);
2879 value_set_si(offset
.d
, 1);
2880 value_init(offset
.x
.n
);
2881 value_oppose(offset
.x
.n
, min
);
2882 eadd(&offset
, &p
->arr
[0]);
2883 free_evalue_refs(&offset
);
2892 value_set_si(fl
.d
, 0);
2893 fl
.x
.p
= new_enode(flooring
, 3, -1);
2894 evalue_set_si(&fl
.x
.p
->arr
[1], 0, 1);
2895 evalue_set_si(&fl
.x
.p
->arr
[2], -1, 1);
2896 evalue_copy(&fl
.x
.p
->arr
[0], &p
->arr
[0]);
2898 eadd(&fl
, &p
->arr
[0]);
2899 reorder_terms(p
, &p
->arr
[0]);
2902 free_evalue_refs(&fl
);
2907 void evalue_frac2floor(evalue
*e
)
2910 if (value_notzero_p(e
->d
) || e
->x
.p
->type
!= partition
)
2913 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
)
2914 if (frac2floor_in_domain(&e
->x
.p
->arr
[2*i
+1],
2915 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
])))
2916 reduce_evalue(&e
->x
.p
->arr
[2*i
+1]);
2919 static Matrix
*esum_add_constraint(int nvar
, Polyhedron
*D
, Matrix
*C
,
2924 int nparam
= D
->Dimension
- nvar
;
2927 nr
= D
->NbConstraints
+ 2;
2928 nc
= D
->Dimension
+ 2 + 1;
2929 C
= Matrix_Alloc(nr
, nc
);
2930 for (i
= 0; i
< D
->NbConstraints
; ++i
) {
2931 Vector_Copy(D
->Constraint
[i
], C
->p
[i
], 1 + nvar
);
2932 Vector_Copy(D
->Constraint
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
2933 D
->Dimension
+ 1 - nvar
);
2938 nc
= C
->NbColumns
+ 1;
2939 C
= Matrix_Alloc(nr
, nc
);
2940 for (i
= 0; i
< oldC
->NbRows
; ++i
) {
2941 Vector_Copy(oldC
->p
[i
], C
->p
[i
], 1 + nvar
);
2942 Vector_Copy(oldC
->p
[i
] + 1 + nvar
, C
->p
[i
] + 1 + nvar
+ 1,
2943 oldC
->NbColumns
- 1 - nvar
);
2946 value_set_si(C
->p
[nr
-2][0], 1);
2947 value_set_si(C
->p
[nr
-2][1 + nvar
], 1);
2948 value_set_si(C
->p
[nr
-2][nc
- 1], -1);
2950 Vector_Copy(row
->p
, C
->p
[nr
-1], 1 + nvar
+ 1);
2951 Vector_Copy(row
->p
+ 1 + nvar
+ 1, C
->p
[nr
-1] + C
->NbColumns
- 1 - nparam
,
2957 static void floor2frac_r(evalue
*e
, int nvar
)
2964 if (value_notzero_p(e
->d
))
2969 assert(p
->type
== flooring
);
2970 for (i
= 1; i
< p
->size
; i
++)
2971 floor2frac_r(&p
->arr
[i
], nvar
);
2973 for (pp
= &p
->arr
[0]; value_zero_p(pp
->d
); pp
= &pp
->x
.p
->arr
[0]) {
2974 assert(pp
->x
.p
->type
== polynomial
);
2975 pp
->x
.p
->pos
-= nvar
;
2979 value_set_si(f
.d
, 0);
2980 f
.x
.p
= new_enode(fractional
, 3, -1);
2981 evalue_set_si(&f
.x
.p
->arr
[1], 0, 1);
2982 evalue_set_si(&f
.x
.p
->arr
[2], -1, 1);
2983 evalue_copy(&f
.x
.p
->arr
[0], &p
->arr
[0]);
2985 eadd(&f
, &p
->arr
[0]);
2986 reorder_terms(p
, &p
->arr
[0]);
2989 free_evalue_refs(&f
);
2992 /* Convert flooring back to fractional and shift position
2993 * of the parameters by nvar
2995 static void floor2frac(evalue
*e
, int nvar
)
2997 floor2frac_r(e
, nvar
);
3001 evalue
*esum_over_domain_cst(int nvar
, Polyhedron
*D
, Matrix
*C
)
3004 int nparam
= D
->Dimension
- nvar
;
3008 D
= Constraints2Polyhedron(C
, 0);
3012 t
= barvinok_enumerate_e(D
, 0, nparam
, 0);
3014 /* Double check that D was not unbounded. */
3015 assert(!(value_pos_p(t
->d
) && value_neg_p(t
->x
.n
)));
3023 evalue
*esum_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
3030 evalue
*factor
= NULL
;
3033 if (EVALUE_IS_ZERO(*e
))
3037 Polyhedron
*DD
= Disjoint_Domain(D
, 0, 0);
3044 res
= esum_over_domain(e
, nvar
, Q
, C
);
3047 for (Q
= DD
; Q
; Q
= DD
) {
3053 t
= esum_over_domain(e
, nvar
, Q
, C
);
3060 free_evalue_refs(t
);
3067 if (value_notzero_p(e
->d
)) {
3070 t
= esum_over_domain_cst(nvar
, D
, C
);
3072 if (!EVALUE_IS_ONE(*e
))
3078 switch (e
->x
.p
->type
) {
3080 evalue
*pp
= &e
->x
.p
->arr
[0];
3082 if (pp
->x
.p
->pos
> nvar
) {
3083 /* remainder is independent of the summated vars */
3089 floor2frac(&f
, nvar
);
3091 t
= esum_over_domain_cst(nvar
, D
, C
);
3095 free_evalue_refs(&f
);
3100 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3101 poly_denom(pp
, &row
->p
[1 + nvar
]);
3102 value_set_si(row
->p
[0], 1);
3103 for (pp
= &e
->x
.p
->arr
[0]; value_zero_p(pp
->d
);
3104 pp
= &pp
->x
.p
->arr
[0]) {
3106 assert(pp
->x
.p
->type
== polynomial
);
3108 if (pos
>= 1 + nvar
)
3110 value_assign(row
->p
[pos
], row
->p
[1+nvar
]);
3111 value_division(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].d
);
3112 value_multiply(row
->p
[pos
], row
->p
[pos
], pp
->x
.p
->arr
[1].x
.n
);
3114 value_assign(row
->p
[1 + D
->Dimension
+ 1], row
->p
[1+nvar
]);
3115 value_division(row
->p
[1 + D
->Dimension
+ 1],
3116 row
->p
[1 + D
->Dimension
+ 1],
3118 value_multiply(row
->p
[1 + D
->Dimension
+ 1],
3119 row
->p
[1 + D
->Dimension
+ 1],
3121 value_oppose(row
->p
[1 + nvar
], row
->p
[1 + nvar
]);
3125 int pos
= e
->x
.p
->pos
;
3129 value_init(factor
->d
);
3130 value_set_si(factor
->d
, 0);
3131 factor
->x
.p
= new_enode(polynomial
, 2, pos
- nvar
);
3132 evalue_set_si(&factor
->x
.p
->arr
[0], 0, 1);
3133 evalue_set_si(&factor
->x
.p
->arr
[1], 1, 1);
3137 row
= Vector_Alloc(1 + D
->Dimension
+ 1 + 1);
3138 for (i
= 0; i
< D
->NbRays
; ++i
)
3139 if (value_notzero_p(D
->Ray
[i
][pos
]))
3141 assert(i
< D
->NbRays
);
3142 if (value_neg_p(D
->Ray
[i
][pos
])) {
3144 value_init(factor
->d
);
3145 evalue_set_si(factor
, -1, 1);
3147 value_set_si(row
->p
[0], 1);
3148 value_set_si(row
->p
[pos
], 1);
3149 value_set_si(row
->p
[1 + nvar
], -1);
3156 i
= type_offset(e
->x
.p
);
3158 res
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
3163 evalue_copy(&cum
, factor
);
3167 for (; i
< e
->x
.p
->size
; ++i
) {
3171 C
= esum_add_constraint(nvar
, D
, C
, row
);
3177 Vector_Print(stderr, P_VALUE_FMT, row);
3179 Matrix_Print(stderr, P_VALUE_FMT, C);
3181 t
= esum_over_domain(&e
->x
.p
->arr
[i
], nvar
, D
, C
);
3190 free_evalue_refs(t
);
3193 if (factor
&& i
+1 < e
->x
.p
->size
)
3200 free_evalue_refs(factor
);
3201 free_evalue_refs(&cum
);
3213 evalue
*esum(evalue
*e
, int nvar
)
3221 if (nvar
== 0 || EVALUE_IS_ZERO(*e
)) {
3222 evalue_copy(res
, e
);
3226 evalue_set_si(res
, 0, 1);
3228 assert(value_zero_p(e
->d
));
3229 assert(e
->x
.p
->type
== partition
);
3231 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
3233 t
= esum_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
3234 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), 0);
3236 free_evalue_refs(t
);
3245 /* Initial silly implementation */
3246 void eor(evalue
*e1
, evalue
*res
)
3252 evalue_set_si(&mone
, -1, 1);
3254 evalue_copy(&E
, res
);
3260 free_evalue_refs(&E
);
3261 free_evalue_refs(&mone
);