4 #include "ev_operations.h"
7 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
9 void evalue_set_si(evalue
*ev
, int n
, int d
) {
10 value_set_si(ev
->d
, d
);
12 value_set_si(ev
->x
.n
, n
);
15 void evalue_set(evalue
*ev
, Value n
, Value d
) {
16 value_assign(ev
->d
, d
);
18 value_assign(ev
->x
.n
, n
);
21 void aep_evalue(evalue
*e
, int *ref
) {
26 if (value_notzero_p(e
->d
))
27 return; /* a rational number, its already reduced */
29 return; /* hum... an overflow probably occured */
31 /* First check the components of p */
32 for (i
=0;i
<p
->size
;i
++)
33 aep_evalue(&p
->arr
[i
],ref
);
40 p
->pos
= ref
[p
->pos
-1]+1;
46 void addeliminatedparams_evalue(evalue
*e
,Matrix
*CT
) {
52 if (value_notzero_p(e
->d
))
53 return; /* a rational number, its already reduced */
55 return; /* hum... an overflow probably occured */
58 ref
= (int *)malloc(sizeof(int)*(CT
->NbRows
-1));
59 for(i
=0;i
<CT
->NbRows
-1;i
++)
60 for(j
=0;j
<CT
->NbColumns
;j
++)
61 if(value_notzero_p(CT
->p
[i
][j
])) {
66 /* Transform the references in e, using ref */
70 } /* addeliminatedparams_evalue */
72 static int mod_rational_smaller(evalue
*e1
, evalue
*e2
)
78 assert(value_notzero_p(e1
->d
));
79 assert(value_notzero_p(e2
->d
));
80 value_multiply(m
, e1
->x
.n
, e2
->d
);
81 value_division(m
, m
, e1
->d
);
82 if (value_lt(m
, e2
->x
.n
))
84 else if (value_gt(m
, e2
->x
.n
))
93 static int mod_term_smaller_r(evalue
*e1
, evalue
*e2
)
95 if (value_notzero_p(e1
->d
)) {
96 if (value_zero_p(e2
->d
))
98 int r
= mod_rational_smaller(e1
, e2
);
99 return r
== -1 ? 0 : r
;
101 if (value_notzero_p(e2
->d
))
103 if (e1
->x
.p
->pos
< e2
->x
.p
->pos
)
105 else if (e1
->x
.p
->pos
> e2
->x
.p
->pos
)
108 int r
= mod_rational_smaller(&e1
->x
.p
->arr
[1], &e2
->x
.p
->arr
[1]);
110 ? mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0])
115 static int mod_term_smaller(evalue
*e1
, evalue
*e2
)
117 assert(value_zero_p(e1
->d
));
118 assert(value_zero_p(e2
->d
));
119 assert(e1
->x
.p
->type
== fractional
);
120 assert(e2
->x
.p
->type
== fractional
);
121 return mod_term_smaller_r(&e1
->x
.p
->arr
[0], &e2
->x
.p
->arr
[0]);
132 struct fixed_param
*fixed
;
137 static int relations_depth(evalue
*e
)
142 value_zero_p(e
->d
) && e
->x
.p
->type
== relation
;
143 e
= &e
->x
.p
->arr
[1], ++d
);
147 #define EVALUE_IS_ONE(ev) (value_pos_p((ev).d) && value_one_p((ev).x.n))
149 static void realloc_substitution(struct subst
*s
, int d
)
151 struct fixed_param
*n
;
154 for (i
= 0; i
< s
->n
; ++i
)
161 static int add_modulo_substitution(struct subst
*s
, evalue
*r
)
168 assert(value_zero_p(r
->d
) && r
->x
.p
->type
== relation
);
171 /* May have been reduced already */
172 if (value_notzero_p(m
->d
))
175 if (s
->n
>= s
->max
) {
176 int d
= relations_depth(r
);
177 realloc_substitution(s
, d
);
180 assert(value_zero_p(m
->d
) && m
->x
.p
->type
== fractional
);
181 assert(m
->x
.p
->size
== 3);
182 assert(EVALUE_IS_ONE(m
->x
.p
->arr
[2]));
183 assert(EVALUE_IS_ZERO(m
->x
.p
->arr
[1]));
186 assert(value_zero_p(p
->d
) && p
->x
.p
->type
== polynomial
);
187 assert(p
->x
.p
->size
== 2);
190 neg
= value_neg_p(f
->x
.n
);
192 value_init(s
->fixed
[s
->n
].m
);
193 value_assign(s
->fixed
[s
->n
].m
, f
->d
);
194 s
->fixed
[s
->n
].pos
= p
->x
.p
->pos
;
195 value_init(s
->fixed
[s
->n
].d
);
197 value_oppose(s
->fixed
[s
->n
].d
, f
->x
.n
);
199 value_assign(s
->fixed
[s
->n
].d
, f
->x
.n
);
200 value_init(s
->fixed
[s
->n
].s
.d
);
201 evalue_copy(&s
->fixed
[s
->n
].s
, &p
->x
.p
->arr
[0]);
205 evalue_set_si(&mone
, -1, 1);
206 emul(&mone
, &s
->fixed
[s
->n
].s
);
207 free_evalue_refs(&mone
);
214 void _reduce_evalue (evalue
*e
, struct subst
*s
, int fract
) {
219 if (value_notzero_p(e
->d
)) {
221 mpz_fdiv_r(e
->x
.n
, e
->x
.n
, e
->d
);
222 return; /* a rational number, its already reduced */
226 return; /* hum... an overflow probably occured */
228 /* First reduce the components of p */
229 for (i
=0; i
<p
->size
; i
++) {
230 int add
= p
->type
== relation
&& i
== 1;
232 add
= add_modulo_substitution(s
, e
);
234 if (i
== 0 && p
->type
==fractional
)
235 _reduce_evalue(&p
->arr
[i
], s
, 1);
237 _reduce_evalue(&p
->arr
[i
], s
, fract
);
241 value_clear(s
->fixed
[s
->n
].m
);
242 value_clear(s
->fixed
[s
->n
].d
);
243 free_evalue_refs(&s
->fixed
[s
->n
].s
);
247 if (p
->type
==periodic
) {
249 /* Try to reduce the period */
250 for (i
=1; i
<=(p
->size
)/2; i
++) {
251 if ((p
->size
% i
)==0) {
253 /* Can we reduce the size to i ? */
255 for (k
=j
+i
; k
<e
->x
.p
->size
; k
+=i
)
256 if (!eequal(&p
->arr
[j
], &p
->arr
[k
])) goto you_lose
;
259 for (j
=i
; j
<p
->size
; j
++) free_evalue_refs(&p
->arr
[j
]);
263 you_lose
: /* OK, lets not do it */
268 /* Try to reduce its strength */
271 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
275 else if (p
->type
==polynomial
) {
276 for (k
= 0; s
&& k
< s
->n
; ++k
) {
277 if (s
->fixed
[k
].pos
== p
->pos
) {
278 int divide
= value_notone_p(s
->fixed
[k
].d
);
284 if (value_notzero_p(s
->fixed
[k
].m
)) {
287 assert(p
->size
== 2);
288 if (!mpz_divisible_p(s
->fixed
[k
].m
, p
->arr
[1].d
))
295 value_assign(d
.d
, s
->fixed
[k
].d
);
297 if (value_notzero_p(s
->fixed
[k
].m
))
298 value_assign(d
.x
.n
, s
->fixed
[k
].m
);
300 value_set_si(d
.x
.n
, 1);
303 for (i
=p
->size
-1;i
>=1;i
--) {
304 emul(&s
->fixed
[k
].s
, &p
->arr
[i
]);
306 emul(&d
, &p
->arr
[i
]);
307 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
308 free_evalue_refs(&(p
->arr
[i
]));
311 _reduce_evalue(&p
->arr
[0], s
, fract
);
314 free_evalue_refs(&d
);
320 /* Try to reduce the degree */
321 for (i
=p
->size
-1;i
>=1;i
--) {
322 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
324 /* Zero coefficient */
325 free_evalue_refs(&(p
->arr
[i
]));
330 /* Try to reduce its strength */
333 memcpy(e
,&p
->arr
[0],sizeof(evalue
));
337 else if (p
->type
==fractional
) {
341 if (value_notzero_p(p
->arr
[0].d
)) {
343 value_assign(v
.d
, p
->arr
[0].d
);
345 mpz_fdiv_r(v
.x
.n
, p
->arr
[0].x
.n
, p
->arr
[0].d
);
349 /* reduction may have made this fractional arg smaller */
350 for (i
= 1; i
< p
->size
; ++i
)
351 if (value_zero_p(p
->arr
[i
].d
) &&
352 p
->arr
[i
].x
.p
->type
== fractional
&&
353 !mod_term_smaller(e
, &p
->arr
[i
]))
357 value_set_si(v
.d
, 0);
358 v
.x
.p
= new_enode(fractional
, 3, -1);
359 evalue_set_si(&v
.x
.p
->arr
[1], 0, 1);
360 evalue_set_si(&v
.x
.p
->arr
[2], 1, 1);
361 evalue_copy(&v
.x
.p
->arr
[0], &p
->arr
[0]);
368 for (i
=p
->size
-1;i
>=2;i
--) {
369 emul(&v
, &p
->arr
[i
]);
370 eadd(&p
->arr
[i
], &p
->arr
[i
-1]);
371 free_evalue_refs(&(p
->arr
[i
]));
374 free_evalue_refs(&v
);
375 _reduce_evalue(&p
->arr
[1], s
, fract
);
378 /* Try to reduce the degree */
379 for (i
=p
->size
-1;i
>=2;i
--) {
380 if (!(value_one_p(p
->arr
[i
].d
) && value_zero_p(p
->arr
[i
].x
.n
)))
382 /* Zero coefficient */
383 free_evalue_refs(&(p
->arr
[i
]));
388 /* Try to reduce its strength */
391 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
392 free_evalue_refs(&(p
->arr
[0]));
396 else if (p
->type
== relation
) {
397 if (p
->size
== 3 && EVALUE_IS_ZERO(p
->arr
[2])) {
398 free_evalue_refs(&(p
->arr
[2]));
401 if (p
->size
== 2 && EVALUE_IS_ZERO(p
->arr
[1])) {
402 free_evalue_refs(&(p
->arr
[1]));
403 free_evalue_refs(&(p
->arr
[0]));
404 evalue_set_si(e
, 0, 1);
406 } else if (value_notzero_p(p
->arr
[0].d
)) {
407 if (value_zero_p(p
->arr
[0].x
.n
)) {
409 memcpy(e
,&p
->arr
[1],sizeof(evalue
));
411 free_evalue_refs(&(p
->arr
[2]));
415 memcpy(e
,&p
->arr
[2],sizeof(evalue
));
417 evalue_set_si(e
, 0, 1);
418 free_evalue_refs(&(p
->arr
[1]));
420 free_evalue_refs(&(p
->arr
[0]));
425 } /* reduce_evalue */
427 static void add_substitution(struct subst
*s
, Value
*row
, unsigned dim
)
432 for (k
= 0; k
< dim
; ++k
)
433 if (value_notzero_p(row
[k
+1]))
436 Vector_Normalize_Positive(row
+1, dim
+1, k
);
437 assert(s
->n
< s
->max
);
438 value_init(s
->fixed
[s
->n
].d
);
439 value_init(s
->fixed
[s
->n
].m
);
440 value_assign(s
->fixed
[s
->n
].d
, row
[k
+1]);
441 s
->fixed
[s
->n
].pos
= k
+1;
442 value_set_si(s
->fixed
[s
->n
].m
, 0);
443 r
= &s
->fixed
[s
->n
].s
;
445 for (l
= k
+1; l
< dim
; ++l
)
446 if (value_notzero_p(row
[l
+1])) {
447 value_set_si(r
->d
, 0);
448 r
->x
.p
= new_enode(polynomial
, 2, l
+ 1);
449 value_init(r
->x
.p
->arr
[1].x
.n
);
450 value_oppose(r
->x
.p
->arr
[1].x
.n
, row
[l
+1]);
451 value_set_si(r
->x
.p
->arr
[1].d
, 1);
455 value_oppose(r
->x
.n
, row
[dim
+1]);
456 value_set_si(r
->d
, 1);
460 void reduce_evalue (evalue
*e
) {
461 if (value_notzero_p(e
->d
))
462 return; /* a rational number, its already reduced */
464 if (e
->x
.p
->type
== partition
) {
465 struct subst s
= { NULL
, 0, 0 };
468 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
470 Polyhedron
*D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
471 /* This shouldn't really happen;
472 * Empty domains should not be added.
479 D
= DomainConvex(D
, 0);
480 if (!D
->next
&& D
->NbEq
) {
484 realloc_substitution(&s
, dim
);
486 int d
= relations_depth(&e
->x
.p
->arr
[2*i
+1]);
488 NALLOC(s
.fixed
, s
.max
);
491 for (j
= 0; j
< D
->NbEq
; ++j
)
492 add_substitution(&s
, D
->Constraint
[j
], dim
);
494 if (D
!= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))
496 _reduce_evalue(&e
->x
.p
->arr
[2*i
+1], &s
, 0);
497 if (EVALUE_IS_ZERO(e
->x
.p
->arr
[2*i
+1])) {
499 free_evalue_refs(&e
->x
.p
->arr
[2*i
+1]);
500 Domain_Free(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]));
501 value_clear(e
->x
.p
->arr
[2*i
].d
);
503 e
->x
.p
->arr
[2*i
] = e
->x
.p
->arr
[e
->x
.p
->size
];
504 e
->x
.p
->arr
[2*i
+1] = e
->x
.p
->arr
[e
->x
.p
->size
+1];
509 for (j
= 0; j
< s
.n
; ++j
) {
510 value_clear(s
.fixed
[j
].d
);
511 free_evalue_refs(&s
.fixed
[j
].s
);
518 _reduce_evalue(e
, 0, 0);
521 void print_evalue(FILE *DST
,evalue
*e
,char **pname
) {
523 if(value_notzero_p(e
->d
)) {
524 if(value_notone_p(e
->d
)) {
525 value_print(DST
,VALUE_FMT
,e
->x
.n
);
527 value_print(DST
,VALUE_FMT
,e
->d
);
530 value_print(DST
,VALUE_FMT
,e
->x
.n
);
534 print_enode(DST
,e
->x
.p
,pname
);
538 void print_enode(FILE *DST
,enode
*p
,char **pname
) {
543 fprintf(DST
, "NULL");
549 for (i
=0; i
<p
->size
; i
++) {
550 print_evalue(DST
, &p
->arr
[i
], pname
);
554 fprintf(DST
, " }\n");
558 for (i
=p
->size
-1; i
>=0; i
--) {
559 print_evalue(DST
, &p
->arr
[i
], pname
);
560 if (i
==1) fprintf(DST
, " * %s + ", pname
[p
->pos
-1]);
562 fprintf(DST
, " * %s^%d + ", pname
[p
->pos
-1], i
);
564 fprintf(DST
, " )\n");
568 for (i
=0; i
<p
->size
; i
++) {
569 print_evalue(DST
, &p
->arr
[i
], pname
);
570 if (i
!=(p
->size
-1)) fprintf(DST
, ", ");
572 fprintf(DST
," ]_%s", pname
[p
->pos
-1]);
576 for (i
=p
->size
-1; i
>=1; i
--) {
577 print_evalue(DST
, &p
->arr
[i
], pname
);
581 print_evalue(DST
, &p
->arr
[0], pname
);
584 fprintf(DST
, "^%d + ", i
-1);
589 fprintf(DST
, " )\n");
593 print_evalue(DST
, &p
->arr
[0], pname
);
594 fprintf(DST
, "= 0 ] * \n");
595 print_evalue(DST
, &p
->arr
[1], pname
);
597 fprintf(DST
, " +\n [ ");
598 print_evalue(DST
, &p
->arr
[0], pname
);
599 fprintf(DST
, "!= 0 ] * \n");
600 print_evalue(DST
, &p
->arr
[2], pname
);
604 for (i
=0; i
<p
->size
/2; i
++) {
605 Print_Domain(DST
, EVALUE_DOMAIN(p
->arr
[2*i
]), pname
);
606 print_evalue(DST
, &p
->arr
[2*i
+1], pname
);
615 static void eadd_rev(evalue
*e1
, evalue
*res
)
619 evalue_copy(&ev
, e1
);
621 free_evalue_refs(res
);
625 static void eadd_rev_cst (evalue
*e1
, evalue
*res
)
629 evalue_copy(&ev
, e1
);
630 eadd(res
, &ev
.x
.p
->arr
[ev
.x
.p
->type
==fractional
]);
631 free_evalue_refs(res
);
635 struct section
{ Polyhedron
* D
; evalue E
; };
637 void eadd_partitions (evalue
*e1
,evalue
*res
)
642 s
= (struct section
*)
643 malloc((e1
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2+1) *
644 sizeof(struct section
));
648 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
649 assert(res
->x
.p
->size
>= 2);
650 fd
= DomainDifference(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
651 EVALUE_DOMAIN(res
->x
.p
->arr
[0]), 0);
653 for (i
= 1; i
< res
->x
.p
->size
/2; ++i
) {
655 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
664 value_init(s
[n
].E
.d
);
665 evalue_copy(&s
[n
].E
, &e1
->x
.p
->arr
[2*j
+1]);
669 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
670 fd
= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]);
671 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
673 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
674 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
680 fd
= DomainDifference(fd
, EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]), 0);
681 if (t
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
683 value_init(s
[n
].E
.d
);
684 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
685 eadd(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
690 s
[n
].E
= res
->x
.p
->arr
[2*i
+1];
694 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
697 if (fd
!= EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]))
698 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
699 value_clear(res
->x
.p
->arr
[2*i
].d
);
703 res
->x
.p
= new_enode(partition
, 2*n
, -1);
704 for (j
= 0; j
< n
; ++j
) {
705 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
706 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
707 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
708 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
714 static void explicit_complement(evalue
*res
)
716 enode
*rel
= new_enode(relation
, 3, 0);
718 value_clear(rel
->arr
[0].d
);
719 rel
->arr
[0] = res
->x
.p
->arr
[0];
720 value_clear(rel
->arr
[1].d
);
721 rel
->arr
[1] = res
->x
.p
->arr
[1];
722 value_set_si(rel
->arr
[2].d
, 1);
723 value_init(rel
->arr
[2].x
.n
);
724 value_set_si(rel
->arr
[2].x
.n
, 0);
729 void eadd(evalue
*e1
,evalue
*res
) {
732 if (value_notzero_p(e1
->d
) && value_notzero_p(res
->d
)) {
733 /* Add two rational numbers */
739 value_multiply(m1
,e1
->x
.n
,res
->d
);
740 value_multiply(m2
,res
->x
.n
,e1
->d
);
741 value_addto(res
->x
.n
,m1
,m2
);
742 value_multiply(res
->d
,e1
->d
,res
->d
);
743 Gcd(res
->x
.n
,res
->d
,&g
);
744 if (value_notone_p(g
)) {
745 value_division(res
->d
,res
->d
,g
);
746 value_division(res
->x
.n
,res
->x
.n
,g
);
748 value_clear(g
); value_clear(m1
); value_clear(m2
);
751 else if (value_notzero_p(e1
->d
) && value_zero_p(res
->d
)) {
752 switch (res
->x
.p
->type
) {
754 /* Add the constant to the constant term of a polynomial*/
755 eadd(e1
, &res
->x
.p
->arr
[0]);
758 /* Add the constant to all elements of a periodic number */
759 for (i
=0; i
<res
->x
.p
->size
; i
++) {
760 eadd(e1
, &res
->x
.p
->arr
[i
]);
764 fprintf(stderr
, "eadd: cannot add const with vector\n");
767 eadd(e1
, &res
->x
.p
->arr
[1]);
770 assert(EVALUE_IS_ZERO(*e1
));
771 break; /* Do nothing */
773 /* Create (zero) complement if needed */
774 if (res
->x
.p
->size
< 3 && !EVALUE_IS_ZERO(*e1
))
775 explicit_complement(res
);
776 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
777 eadd(e1
, &res
->x
.p
->arr
[i
]);
783 /* add polynomial or periodic to constant
784 * you have to exchange e1 and res, before doing addition */
786 else if (value_zero_p(e1
->d
) && value_notzero_p(res
->d
)) {
790 else { // ((e1->d==0) && (res->d==0))
791 assert(!((e1
->x
.p
->type
== partition
) ^
792 (res
->x
.p
->type
== partition
)));
793 if (e1
->x
.p
->type
== partition
) {
794 eadd_partitions(e1
, res
);
797 if (e1
->x
.p
->type
== relation
&&
798 (res
->x
.p
->type
!= relation
||
799 mod_term_smaller(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
803 if (res
->x
.p
->type
== relation
) {
804 if (e1
->x
.p
->type
== relation
&&
805 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
806 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
807 explicit_complement(res
);
808 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
809 explicit_complement(e1
);
810 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
811 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
814 if (res
->x
.p
->size
< 3)
815 explicit_complement(res
);
816 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
817 eadd(e1
, &res
->x
.p
->arr
[i
]);
820 if ((e1
->x
.p
->type
!= res
->x
.p
->type
) ) {
821 /* adding to evalues of different type. two cases are possible
822 * res is periodic and e1 is polynomial, you have to exchange
823 * e1 and res then to add e1 to the constant term of res */
824 if (e1
->x
.p
->type
== polynomial
) {
825 eadd_rev_cst(e1
, res
);
827 else if (res
->x
.p
->type
== polynomial
) {
828 /* res is polynomial and e1 is periodic,
829 add e1 to the constant term of res */
831 eadd(e1
,&res
->x
.p
->arr
[0]);
837 else if (e1
->x
.p
->pos
!= res
->x
.p
->pos
||
838 (res
->x
.p
->type
== fractional
&&
839 !eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]))) {
840 /* adding evalues of different position (i.e function of different unknowns
841 * to case are possible */
843 switch (res
->x
.p
->type
) {
845 if(mod_term_smaller(res
, e1
))
846 eadd(e1
,&res
->x
.p
->arr
[1]);
848 eadd_rev_cst(e1
, res
);
850 case polynomial
: // res and e1 are polynomials
851 // add e1 to the constant term of res
853 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
854 eadd(e1
,&res
->x
.p
->arr
[0]);
856 eadd_rev_cst(e1
, res
);
857 // value_clear(g); value_clear(m1); value_clear(m2);
859 case periodic
: // res and e1 are pointers to periodic numbers
860 //add e1 to all elements of res
862 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
863 for (i
=0;i
<res
->x
.p
->size
;i
++) {
864 eadd(e1
,&res
->x
.p
->arr
[i
]);
873 //same type , same pos and same size
874 if (e1
->x
.p
->size
== res
->x
.p
->size
) {
875 // add any element in e1 to the corresponding element in res
876 if (res
->x
.p
->type
== fractional
)
877 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
878 i
= res
->x
.p
->type
== fractional
? 1 : 0;
879 for (; i
<res
->x
.p
->size
; i
++) {
880 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
885 /* Sizes are different */
886 switch(res
->x
.p
->type
) {
889 /* VIN100: if e1-size > res-size you have to copy e1 in a */
890 /* new enode and add res to that new node. If you do not do */
891 /* that, you lose the the upper weight part of e1 ! */
893 if(e1
->x
.p
->size
> res
->x
.p
->size
)
897 if (res
->x
.p
->type
== fractional
)
898 assert(eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0]));
899 i
= res
->x
.p
->type
== fractional
? 1 : 0;
900 for (; i
<e1
->x
.p
->size
; i
++) {
901 eadd(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
908 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
911 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
912 of the sizes of e1 and res, then to copy res periodicaly in ne, after
913 to add periodicaly elements of e1 to elements of ne, and finaly to
918 value_init(ex
); value_init(ey
);value_init(ep
);
921 value_set_si(ex
,e1
->x
.p
->size
);
922 value_set_si(ey
,res
->x
.p
->size
);
923 value_assign (ep
,*Lcm(ex
,ey
));
924 p
=(int)mpz_get_si(ep
);
925 ne
= (evalue
*) malloc (sizeof(evalue
));
927 value_set_si( ne
->d
,0);
929 ne
->x
.p
=new_enode(res
->x
.p
->type
,p
, res
->x
.p
->pos
);
931 value_assign(ne
->x
.p
->arr
[i
].d
, res
->x
.p
->arr
[i
%y
].d
);
932 if (value_notzero_p(ne
->x
.p
->arr
[i
].d
)) {
933 value_init(ne
->x
.p
->arr
[i
].x
.n
);
934 value_assign(ne
->x
.p
->arr
[i
].x
.n
, res
->x
.p
->arr
[i
%y
].x
.n
);
937 ne
->x
.p
->arr
[i
].x
.p
=ecopy(res
->x
.p
->arr
[i
%y
].x
.p
);
941 eadd(&e1
->x
.p
->arr
[i
%x
], &ne
->x
.p
->arr
[i
]);
944 value_assign(res
->d
, ne
->d
);
950 fprintf(stderr
, "eadd: ?cannot add vectors of different length\n");
957 static void emul_rev (evalue
*e1
, evalue
*res
)
961 evalue_copy(&ev
, e1
);
963 free_evalue_refs(res
);
967 static void emul_poly (evalue
*e1
, evalue
*res
)
969 int i
, j
, o
= res
->x
.p
->type
== fractional
;
971 int size
=(e1
->x
.p
->size
+ res
->x
.p
->size
- o
- 1);
973 value_set_si(tmp
.d
,0);
974 tmp
.x
.p
=new_enode(res
->x
.p
->type
, size
, res
->x
.p
->pos
);
976 evalue_copy(&tmp
.x
.p
->arr
[0], &e1
->x
.p
->arr
[0]);
977 for (i
=o
; i
< e1
->x
.p
->size
; i
++) {
978 evalue_copy(&tmp
.x
.p
->arr
[i
], &e1
->x
.p
->arr
[i
]);
979 emul(&res
->x
.p
->arr
[o
], &tmp
.x
.p
->arr
[i
]);
982 evalue_set_si(&tmp
.x
.p
->arr
[i
], 0, 1);
983 for (i
=o
+1; i
<res
->x
.p
->size
; i
++)
984 for (j
=o
; j
<e1
->x
.p
->size
; j
++) {
987 evalue_copy(&ev
, &e1
->x
.p
->arr
[j
]);
988 emul(&res
->x
.p
->arr
[i
], &ev
);
989 eadd(&ev
, &tmp
.x
.p
->arr
[i
+j
-o
]);
990 free_evalue_refs(&ev
);
992 free_evalue_refs(res
);
996 void emul_partitions (evalue
*e1
,evalue
*res
)
1001 s
= (struct section
*)
1002 malloc((e1
->x
.p
->size
/2) * (res
->x
.p
->size
/2) *
1003 sizeof(struct section
));
1007 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
) {
1008 for (j
= 0; j
< e1
->x
.p
->size
/2; ++j
) {
1009 d
= DomainIntersection(EVALUE_DOMAIN(e1
->x
.p
->arr
[2*j
]),
1010 EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]), 0);
1016 /* This code is only needed because the partitions
1017 are not true partitions.
1019 for (k
= 0; k
< n
; ++k
) {
1020 if (DomainIncludes(s
[k
].D
, d
))
1022 if (DomainIncludes(d
, s
[k
].D
)) {
1023 Domain_Free(s
[k
].D
);
1024 free_evalue_refs(&s
[k
].E
);
1035 value_init(s
[n
].E
.d
);
1036 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*i
+1]);
1037 emul(&e1
->x
.p
->arr
[2*j
+1], &s
[n
].E
);
1041 Domain_Free(EVALUE_DOMAIN(res
->x
.p
->arr
[2*i
]));
1042 value_clear(res
->x
.p
->arr
[2*i
].d
);
1043 free_evalue_refs(&res
->x
.p
->arr
[2*i
+1]);
1047 res
->x
.p
= new_enode(partition
, 2*n
, -1);
1048 for (j
= 0; j
< n
; ++j
) {
1049 s
[j
].D
= DomainConstraintSimplify(s
[j
].D
, 0);
1050 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1051 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1052 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1058 static void Lcm3(Value i
, Value j
, Value
*res
)
1064 value_multiply(*res
,i
,j
);
1065 value_division(*res
, *res
, aux
);
1069 static void poly_denom(evalue
*p
, Value
*d
)
1071 value_set_si(*d
, 1);
1073 while (value_zero_p(p
->d
)) {
1074 assert(p
->x
.p
->type
== polynomial
);
1075 assert(p
->x
.p
->size
== 2);
1076 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
1077 Lcm3(*d
, p
->x
.p
->arr
[1].d
, d
);
1078 p
= &p
->x
.p
->arr
[0];
1083 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1085 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1086 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1087 * evalues is not treated here */
1089 void emul (evalue
*e1
, evalue
*res
){
1092 if((value_zero_p(e1
->d
)&&e1
->x
.p
->type
==evector
)||(value_zero_p(res
->d
)&&(res
->x
.p
->type
==evector
))) {
1093 fprintf(stderr
, "emul: do not proced on evector type !\n");
1097 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== partition
) {
1098 if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
)
1099 emul_partitions(e1
, res
);
1102 } else if (value_zero_p(res
->d
) && res
->x
.p
->type
== partition
) {
1103 for (i
= 0; i
< res
->x
.p
->size
/2; ++i
)
1104 emul(e1
, &res
->x
.p
->arr
[2*i
+1]);
1106 if (value_zero_p(res
->d
) && res
->x
.p
->type
== relation
) {
1107 if (value_zero_p(e1
->d
) && e1
->x
.p
->type
== relation
&&
1108 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1109 if (res
->x
.p
->size
< 3 && e1
->x
.p
->size
== 3)
1110 explicit_complement(res
);
1111 if (e1
->x
.p
->size
< 3 && res
->x
.p
->size
== 3)
1112 explicit_complement(e1
);
1113 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1114 emul(&e1
->x
.p
->arr
[i
], &res
->x
.p
->arr
[i
]);
1117 for (i
= 1; i
< res
->x
.p
->size
; ++i
)
1118 emul(e1
, &res
->x
.p
->arr
[i
]);
1120 if(value_zero_p(e1
->d
)&& value_zero_p(res
->d
)) {
1121 switch(e1
->x
.p
->type
) {
1123 switch(res
->x
.p
->type
) {
1125 if(e1
->x
.p
->pos
== res
->x
.p
->pos
) {
1126 /* Product of two polynomials of the same variable */
1131 /* Product of two polynomials of different variables */
1133 if(res
->x
.p
->pos
< e1
->x
.p
->pos
)
1134 for( i
=0; i
<res
->x
.p
->size
; i
++)
1135 emul(e1
, &res
->x
.p
->arr
[i
]);
1143 /* Product of a polynomial and a periodic or fractional */
1148 switch(res
->x
.p
->type
) {
1150 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
==res
->x
.p
->size
) {
1151 /* Product of two periodics of the same parameter and period */
1153 for(i
=0; i
<res
->x
.p
->size
;i
++)
1154 emul(&(e1
->x
.p
->arr
[i
]), &(res
->x
.p
->arr
[i
]));
1159 if(e1
->x
.p
->pos
==res
->x
.p
->pos
&& e1
->x
.p
->size
!=res
->x
.p
->size
) {
1160 /* Product of two periodics of the same parameter and different periods */
1164 value_init(x
); value_init(y
);value_init(z
);
1167 value_set_si(x
,e1
->x
.p
->size
);
1168 value_set_si(y
,res
->x
.p
->size
);
1169 value_assign (z
,*Lcm(x
,y
));
1170 lcm
=(int)mpz_get_si(z
);
1171 newp
= (evalue
*) malloc (sizeof(evalue
));
1172 value_init(newp
->d
);
1173 value_set_si( newp
->d
,0);
1174 newp
->x
.p
=new_enode(periodic
,lcm
, e1
->x
.p
->pos
);
1175 for(i
=0;i
<lcm
;i
++) {
1176 value_assign(newp
->x
.p
->arr
[i
].d
, res
->x
.p
->arr
[i
%iy
].d
);
1177 if (value_notzero_p(newp
->x
.p
->arr
[i
].d
)) {
1178 value_assign(newp
->x
.p
->arr
[i
].x
.n
, res
->x
.p
->arr
[i
%iy
].x
.n
);
1181 newp
->x
.p
->arr
[i
].x
.p
=ecopy(res
->x
.p
->arr
[i
%iy
].x
.p
);
1186 emul(&e1
->x
.p
->arr
[i
%ix
], &newp
->x
.p
->arr
[i
]);
1188 value_assign(res
->d
,newp
->d
);
1191 value_clear(x
); value_clear(y
);value_clear(z
);
1195 /* Product of two periodics of different parameters */
1197 for(i
=0; i
<res
->x
.p
->size
; i
++)
1198 emul(e1
, &(res
->x
.p
->arr
[i
]));
1204 /* Product of a periodic and a polynomial */
1206 for(i
=0; i
<res
->x
.p
->size
; i
++)
1207 emul(e1
, &(res
->x
.p
->arr
[i
]));
1213 switch(res
->x
.p
->type
) {
1215 for(i
=0; i
<res
->x
.p
->size
; i
++)
1216 emul(e1
, &(res
->x
.p
->arr
[i
]));
1221 if (e1
->x
.p
->pos
== res
->x
.p
->pos
&&
1222 eequal(&e1
->x
.p
->arr
[0], &res
->x
.p
->arr
[0])) {
1225 poly_denom(&e1
->x
.p
->arr
[0], &d
.d
);
1226 if (!value_two_p(d
.d
))
1230 value_set_si(d
.x
.n
, 1);
1232 /* { x }^2 == { x }/2 */
1233 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1234 assert(e1
->x
.p
->size
== 3);
1235 assert(res
->x
.p
->size
== 3);
1237 evalue_copy(&tmp
, &res
->x
.p
->arr
[2]);
1239 eadd(&res
->x
.p
->arr
[1], &tmp
);
1240 emul(&e1
->x
.p
->arr
[2], &tmp
);
1241 emul(&e1
->x
.p
->arr
[1], res
);
1242 eadd(&tmp
, &res
->x
.p
->arr
[2]);
1243 free_evalue_refs(&tmp
);
1248 if(mod_term_smaller(res
, e1
))
1249 for(i
=1; i
<res
->x
.p
->size
; i
++)
1250 emul(e1
, &(res
->x
.p
->arr
[i
]));
1265 if (value_notzero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1266 /* Product of two rational numbers */
1270 value_multiply(res
->d
,e1
->d
,res
->d
);
1271 value_multiply(res
->x
.n
,e1
->x
.n
,res
->x
.n
);
1272 Gcd(res
->x
.n
, res
->d
,&g
);
1273 if (value_notone_p(g
)) {
1274 value_division(res
->d
,res
->d
,g
);
1275 value_division(res
->x
.n
,res
->x
.n
,g
);
1281 if(value_zero_p(e1
->d
)&& value_notzero_p(res
->d
)) {
1282 /* Product of an expression (polynomial or peririodic) and a rational number */
1288 /* Product of a rationel number and an expression (polynomial or peririodic) */
1290 i
= res
->x
.p
->type
== fractional
? 1 : 0;
1291 for (; i
<res
->x
.p
->size
; i
++)
1292 emul(e1
, &res
->x
.p
->arr
[i
]);
1303 void emask(evalue
*mask
, evalue
*res
) {
1309 if (EVALUE_IS_ZERO(*res
))
1312 assert(value_zero_p(mask
->d
));
1313 assert(mask
->x
.p
->type
== partition
);
1314 assert(value_zero_p(res
->d
));
1315 assert(res
->x
.p
->type
== partition
);
1317 s
= (struct section
*)
1318 malloc((mask
->x
.p
->size
/2+1) * (res
->x
.p
->size
/2) *
1319 sizeof(struct section
));
1323 evalue_set_si(&mone
, -1, 1);
1326 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1327 assert(mask
->x
.p
->size
>= 2);
1328 fd
= DomainDifference(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1329 EVALUE_DOMAIN(mask
->x
.p
->arr
[0]), 0);
1331 for (i
= 1; i
< mask
->x
.p
->size
/2; ++i
) {
1333 fd
= DomainDifference(fd
, EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1342 value_init(s
[n
].E
.d
);
1343 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1347 for (i
= 0; i
< mask
->x
.p
->size
/2; ++i
) {
1348 if (EVALUE_IS_ONE(mask
->x
.p
->arr
[2*i
+1]))
1351 fd
= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]);
1352 eadd(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1353 emul(&mone
, &mask
->x
.p
->arr
[2*i
+1]);
1354 for (j
= 0; j
< res
->x
.p
->size
/2; ++j
) {
1356 d
= DomainIntersection(EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]),
1357 EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]), 0);
1363 fd
= DomainDifference(fd
, EVALUE_DOMAIN(res
->x
.p
->arr
[2*j
]), 0);
1364 if (t
!= EVALUE_DOMAIN(mask
->x
.p
->arr
[2*i
]))
1366 value_init(s
[n
].E
.d
);
1367 evalue_copy(&s
[n
].E
, &res
->x
.p
->arr
[2*j
+1]);
1368 emul(&mask
->x
.p
->arr
[2*i
+1], &s
[n
].E
);
1374 /* Just ignore; this may have been previously masked off */
1378 free_evalue_refs(&mone
);
1379 free_evalue_refs(mask
);
1380 free_evalue_refs(res
);
1383 evalue_set_si(res
, 0, 1);
1385 res
->x
.p
= new_enode(partition
, 2*n
, -1);
1386 for (j
= 0; j
< n
; ++j
) {
1387 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*j
], s
[j
].D
);
1388 value_clear(res
->x
.p
->arr
[2*j
+1].d
);
1389 res
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
1396 void evalue_copy(evalue
*dst
, evalue
*src
)
1398 value_assign(dst
->d
, src
->d
);
1399 if(value_notzero_p(src
->d
)) {
1400 value_init(dst
->x
.n
);
1401 value_assign(dst
->x
.n
, src
->x
.n
);
1403 dst
->x
.p
= ecopy(src
->x
.p
);
1406 enode
*new_enode(enode_type type
,int size
,int pos
) {
1412 fprintf(stderr
, "Allocating enode of size 0 !\n" );
1415 res
= (enode
*) malloc(sizeof(enode
) + (size
-1)*sizeof(evalue
));
1419 for(i
=0; i
<size
; i
++) {
1420 value_init(res
->arr
[i
].d
);
1421 value_set_si(res
->arr
[i
].d
,0);
1422 res
->arr
[i
].x
.p
= 0;
1427 enode
*ecopy(enode
*e
) {
1432 res
= new_enode(e
->type
,e
->size
,e
->pos
);
1433 for(i
=0;i
<e
->size
;++i
) {
1434 value_assign(res
->arr
[i
].d
,e
->arr
[i
].d
);
1435 if(value_zero_p(res
->arr
[i
].d
))
1436 res
->arr
[i
].x
.p
= ecopy(e
->arr
[i
].x
.p
);
1437 else if (EVALUE_IS_DOMAIN(res
->arr
[i
]))
1438 EVALUE_SET_DOMAIN(res
->arr
[i
], Domain_Copy(EVALUE_DOMAIN(e
->arr
[i
])));
1440 value_init(res
->arr
[i
].x
.n
);
1441 value_assign(res
->arr
[i
].x
.n
,e
->arr
[i
].x
.n
);
1447 int eequal(evalue
*e1
,evalue
*e2
) {
1452 if (value_ne(e1
->d
,e2
->d
))
1455 /* e1->d == e2->d */
1456 if (value_notzero_p(e1
->d
)) {
1457 if (value_ne(e1
->x
.n
,e2
->x
.n
))
1460 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1464 /* e1->d == e2->d == 0 */
1467 if (p1
->type
!= p2
->type
) return 0;
1468 if (p1
->size
!= p2
->size
) return 0;
1469 if (p1
->pos
!= p2
->pos
) return 0;
1470 for (i
=0; i
<p1
->size
; i
++)
1471 if (!eequal(&p1
->arr
[i
], &p2
->arr
[i
]) )
1476 void free_evalue_refs(evalue
*e
) {
1481 if (EVALUE_IS_DOMAIN(*e
)) {
1482 Domain_Free(EVALUE_DOMAIN(*e
));
1485 } else if (value_pos_p(e
->d
)) {
1487 /* 'e' stores a constant */
1489 value_clear(e
->x
.n
);
1492 assert(value_zero_p(e
->d
));
1495 if (!p
) return; /* null pointer */
1496 for (i
=0; i
<p
->size
; i
++) {
1497 free_evalue_refs(&(p
->arr
[i
]));
1501 } /* free_evalue_refs */
1503 static void mod2table_r(evalue
*e
, Vector
*periods
, Value m
, int p
,
1504 Vector
* val
, evalue
*res
)
1506 unsigned nparam
= periods
->Size
;
1509 double d
= compute_evalue(e
, val
->p
);
1510 d
*= VALUE_TO_DOUBLE(m
);
1515 value_assign(res
->d
, m
);
1516 value_init(res
->x
.n
);
1517 value_set_double(res
->x
.n
, d
);
1518 mpz_fdiv_r(res
->x
.n
, res
->x
.n
, m
);
1521 if (value_one_p(periods
->p
[p
]))
1522 mod2table_r(e
, periods
, m
, p
+1, val
, res
);
1527 value_assign(tmp
, periods
->p
[p
]);
1528 value_set_si(res
->d
, 0);
1529 res
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
1531 value_decrement(tmp
, tmp
);
1532 value_assign(val
->p
[p
], tmp
);
1533 mod2table_r(e
, periods
, m
, p
+1, val
,
1534 &res
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
1535 } while (value_pos_p(tmp
));
1541 static void rel2table(evalue
*e
, int zero
)
1543 if (value_pos_p(e
->d
)) {
1544 if (value_zero_p(e
->x
.n
) == zero
)
1545 value_set_si(e
->x
.n
, 1);
1547 value_set_si(e
->x
.n
, 0);
1548 value_set_si(e
->d
, 1);
1551 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
1552 rel2table(&e
->x
.p
->arr
[i
], zero
);
1556 void evalue_mod2table(evalue
*e
, int nparam
)
1561 if (EVALUE_IS_DOMAIN(*e
) || value_pos_p(e
->d
))
1564 for (i
=0; i
<p
->size
; i
++) {
1565 evalue_mod2table(&(p
->arr
[i
]), nparam
);
1567 if (p
->type
== relation
) {
1573 evalue_copy(©
, &p
->arr
[0]);
1575 rel2table(&p
->arr
[0], 1);
1576 emul(&p
->arr
[0], &p
->arr
[1]);
1578 rel2table(©
, 0);
1579 emul(©
, &p
->arr
[2]);
1580 eadd(&p
->arr
[2], &p
->arr
[1]);
1581 free_evalue_refs(&p
->arr
[2]);
1582 free_evalue_refs(©
);
1584 free_evalue_refs(&p
->arr
[0]);
1589 } else if (p
->type
== fractional
) {
1590 Vector
*periods
= Vector_Alloc(nparam
);
1591 Vector
*val
= Vector_Alloc(nparam
);
1597 value_set_si(tmp
, 1);
1598 Vector_Set(periods
->p
, 1, nparam
);
1599 Vector_Set(val
->p
, 0, nparam
);
1600 for (ev
= &p
->arr
[0]; value_zero_p(ev
->d
); ev
= &ev
->x
.p
->arr
[0]) {
1603 assert(p
->type
== polynomial
);
1604 assert(p
->size
== 2);
1605 value_assign(periods
->p
[p
->pos
-1], p
->arr
[1].d
);
1606 Lcm3(tmp
, p
->arr
[1].d
, &tmp
);
1609 mod2table_r(&p
->arr
[0], periods
, tmp
, 0, val
, &EP
);
1612 evalue_set_si(&res
, 0, 1);
1613 /* Compute the polynomial using Horner's rule */
1614 for (i
=p
->size
-1;i
>1;i
--) {
1615 eadd(&p
->arr
[i
], &res
);
1618 eadd(&p
->arr
[1], &res
);
1620 free_evalue_refs(e
);
1621 free_evalue_refs(&EP
);
1626 Vector_Free(periods
);
1628 } /* evalue_mod2table */
1630 /********************************************************/
1631 /* function in domain */
1632 /* check if the parameters in list_args */
1633 /* verifies the constraints of Domain P */
1634 /********************************************************/
1635 int in_domain(Polyhedron
*P
, Value
*list_args
) {
1638 Value v
; /* value of the constraint of a row when
1639 parameters are instanciated*/
1645 /*P->Constraint constraint matrice of polyhedron P*/
1646 for(row
=0;row
<P
->NbConstraints
;row
++) {
1647 value_assign(v
,P
->Constraint
[row
][P
->Dimension
+1]); /*constant part*/
1648 for(col
=1;col
<P
->Dimension
+1;col
++) {
1649 value_multiply(tmp
,P
->Constraint
[row
][col
],list_args
[col
-1]);
1650 value_addto(v
,v
,tmp
);
1652 if (value_notzero_p(P
->Constraint
[row
][0])) {
1654 /*if v is not >=0 then this constraint is not respected */
1655 if (value_neg_p(v
)) {
1659 return P
->next
? in_domain(P
->next
, list_args
) : 0;
1664 /*if v is not = 0 then this constraint is not respected */
1665 if (value_notzero_p(v
))
1670 /*if not return before this point => all
1671 the constraints are respected */
1677 /****************************************************/
1678 /* function compute enode */
1679 /* compute the value of enode p with parameters */
1680 /* list "list_args */
1681 /* compute the polynomial or the periodic */
1682 /****************************************************/
1684 static double compute_enode(enode
*p
, Value
*list_args
) {
1696 if (p
->type
== polynomial
) {
1698 value_assign(param
,list_args
[p
->pos
-1]);
1700 /* Compute the polynomial using Horner's rule */
1701 for (i
=p
->size
-1;i
>0;i
--) {
1702 res
+=compute_evalue(&p
->arr
[i
],list_args
);
1703 res
*=VALUE_TO_DOUBLE(param
);
1705 res
+=compute_evalue(&p
->arr
[0],list_args
);
1707 else if (p
->type
== fractional
) {
1708 double d
= compute_evalue(&p
->arr
[0], list_args
);
1709 d
-= floor(d
+1e-10);
1711 /* Compute the polynomial using Horner's rule */
1712 for (i
=p
->size
-1;i
>1;i
--) {
1713 res
+=compute_evalue(&p
->arr
[i
],list_args
);
1716 res
+=compute_evalue(&p
->arr
[1],list_args
);
1718 else if (p
->type
== periodic
) {
1719 value_assign(param
,list_args
[p
->pos
-1]);
1721 /* Choose the right element of the periodic */
1722 value_absolute(m
,param
);
1723 value_set_si(param
,p
->size
);
1724 value_modulus(m
,m
,param
);
1725 res
= compute_evalue(&p
->arr
[VALUE_TO_INT(m
)],list_args
);
1727 else if (p
->type
== relation
) {
1728 if (fabs(compute_evalue(&p
->arr
[0], list_args
)) < 1e-10)
1729 res
= compute_evalue(&p
->arr
[1], list_args
);
1730 else if (p
->size
> 2)
1731 res
= compute_evalue(&p
->arr
[2], list_args
);
1733 else if (p
->type
== partition
) {
1734 for (i
= 0; i
< p
->size
/2; ++i
)
1735 if (in_domain(EVALUE_DOMAIN(p
->arr
[2*i
]), list_args
)) {
1736 res
= compute_evalue(&p
->arr
[2*i
+1], list_args
);
1743 } /* compute_enode */
1745 /*************************************************/
1746 /* return the value of Ehrhart Polynomial */
1747 /* It returns a double, because since it is */
1748 /* a recursive function, some intermediate value */
1749 /* might not be integral */
1750 /*************************************************/
1752 double compute_evalue(evalue
*e
,Value
*list_args
) {
1756 if (value_notzero_p(e
->d
)) {
1757 if (value_notone_p(e
->d
))
1758 res
= VALUE_TO_DOUBLE(e
->x
.n
) / VALUE_TO_DOUBLE(e
->d
);
1760 res
= VALUE_TO_DOUBLE(e
->x
.n
);
1763 res
= compute_enode(e
->x
.p
,list_args
);
1765 } /* compute_evalue */
1768 /****************************************************/
1769 /* function compute_poly : */
1770 /* Check for the good validity domain */
1771 /* return the number of point in the Polyhedron */
1772 /* in allocated memory */
1773 /* Using the Ehrhart pseudo-polynomial */
1774 /****************************************************/
1775 Value
*compute_poly(Enumeration
*en
,Value
*list_args
) {
1778 /* double d; int i; */
1780 tmp
= (Value
*) malloc (sizeof(Value
));
1781 assert(tmp
!= NULL
);
1783 value_set_si(*tmp
,0);
1786 return(tmp
); /* no ehrhart polynomial */
1787 if(en
->ValidityDomain
) {
1788 if(!en
->ValidityDomain
->Dimension
) { /* no parameters */
1789 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
1794 return(tmp
); /* no Validity Domain */
1796 if(in_domain(en
->ValidityDomain
,list_args
)) {
1798 #ifdef EVAL_EHRHART_DEBUG
1799 Print_Domain(stdout
,en
->ValidityDomain
);
1800 print_evalue(stdout
,&en
->EP
);
1803 /* d = compute_evalue(&en->EP,list_args);
1805 printf("(double)%lf = %d\n", d, i ); */
1806 value_set_double(*tmp
,compute_evalue(&en
->EP
,list_args
)+.25);
1812 value_set_si(*tmp
,0);
1813 return(tmp
); /* no compatible domain with the arguments */
1814 } /* compute_poly */
1816 size_t value_size(Value v
) {
1817 return (v
[0]._mp_size
> 0 ? v
[0]._mp_size
: -v
[0]._mp_size
)
1818 * sizeof(v
[0]._mp_d
[0]);
1821 size_t domain_size(Polyhedron
*D
)
1824 size_t s
= sizeof(*D
);
1826 for (i
= 0; i
< D
->NbConstraints
; ++i
)
1827 for (j
= 0; j
< D
->Dimension
+2; ++j
)
1828 s
+= value_size(D
->Constraint
[i
][j
]);
1830 for (i
= 0; i
< D
->NbRays
; ++i
)
1831 for (j
= 0; j
< D
->Dimension
+2; ++j
)
1832 s
+= value_size(D
->Ray
[i
][j
]);
1837 size_t enode_size(enode
*p
) {
1838 size_t s
= sizeof(*p
) - sizeof(p
->arr
[0]);
1841 if (p
->type
== partition
)
1842 for (i
= 0; i
< p
->size
/2; ++i
) {
1843 s
+= domain_size(EVALUE_DOMAIN(p
->arr
[2*i
]));
1844 s
+= evalue_size(&p
->arr
[2*i
+1]);
1847 for (i
= 0; i
< p
->size
; ++i
) {
1848 s
+= evalue_size(&p
->arr
[i
]);
1853 size_t evalue_size(evalue
*e
)
1855 size_t s
= sizeof(*e
);
1856 s
+= value_size(e
->d
);
1857 if (value_notzero_p(e
->d
))
1858 s
+= value_size(e
->x
.n
);
1860 s
+= enode_size(e
->x
.p
);