?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
small memory leak
[barvinok.git] / ev_operations.c
blob8afea7f3b56c3c28427ee72aa648e2755d3f314f
1 #include <assert.h>
2 #include <math.h>
3 #include <stdlib.h>
4
5 #include "ev_operations.h"
6 #include "barvinok.h"
7 #include "util.h"
8
9 #ifndef value_pmodulus
10 #define value_pmodulus(ref,val1,val2) (mpz_fdiv_r((ref),(val1),(val2)))
11 #endif
12
13 #define ALLOC(p) p = (typeof(p))malloc(sizeof(*p))
14 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
15
16 void evalue_set_si(evalue *ev, int n, int d) {
17 value_set_si(ev->d, d);
18 value_init(ev->x.n);
19 value_set_si(ev->x.n, n);
20 }
21
22 void evalue_set(evalue *ev, Value n, Value d) {
23 value_assign(ev->d, d);
24 value_init(ev->x.n);
25 value_assign(ev->x.n, n);
26 }
27
28 void aep_evalue(evalue *e, int *ref) {
29
30 enode *p;
31 int i;
32
33 if (value_notzero_p(e->d))
34 return; /* a rational number, its already reduced */
35 if(!(p = e->x.p))
36 return; /* hum... an overflow probably occured */
37
38 /* First check the components of p */
39 for (i=0;i<p->size;i++)
40 aep_evalue(&p->arr[i],ref);
41
42 /* Then p itself */
43 switch (p->type) {
44 case polynomial:
45 case periodic:
46 case evector:
47 p->pos = ref[p->pos-1]+1;
48 }
49 return;
50 } /* aep_evalue */
51
52 /** Comments */
53 void addeliminatedparams_evalue(evalue *e,Matrix *CT) {
54
55 enode *p;
56 int i, j;
57 int *ref;
58
59 if (value_notzero_p(e->d))
60 return; /* a rational number, its already reduced */
61 if(!(p = e->x.p))
62 return; /* hum... an overflow probably occured */
63
64 /* Compute ref */
65 ref = (int *)malloc(sizeof(int)*(CT->NbRows-1));
66 for(i=0;i<CT->NbRows-1;i++)
67 for(j=0;j<CT->NbColumns;j++)
68 if(value_notzero_p(CT->p[i][j])) {
69 ref[i] = j;
70 break;
71 }
72
73 /* Transform the references in e, using ref */
74 aep_evalue(e,ref);
75 free( ref );
76 return;
77 } /* addeliminatedparams_evalue */
78
79 void addeliminatedparams_enum(evalue *e, Matrix *CT, Polyhedron *CEq,
80 unsigned MaxRays, unsigned nparam)
81 {
82 enode *p;
83 int i;
84
85 if (CT->NbRows == CT->NbColumns)
86 return;
87
88 if (EVALUE_IS_ZERO(*e))
89 return;
90
91 if (value_notzero_p(e->d)) {
92 evalue res;
93 value_init(res.d);
94 value_set_si(res.d, 0);
95 res.x.p = new_enode(partition, 2, nparam);
96 EVALUE_SET_DOMAIN(res.x.p->arr[0],
97 DomainConstraintSimplify(Polyhedron_Copy(CEq), MaxRays));
98 value_clear(res.x.p->arr[1].d);
99 res.x.p->arr[1] = *e;
100 *e = res;
101 return;
102 }
103
104 p = e->x.p;
105 assert(p);
106 assert(p->type == partition);
107 p->pos = nparam;
108
109 for (i=0; i<p->size/2; i++) {
110 Polyhedron *D = EVALUE_DOMAIN(p->arr[2*i]);
111 Polyhedron *T = DomainPreimage(D, CT, MaxRays);
112 Domain_Free(D);
113 D = T;
114 T = DomainIntersection(D, CEq, MaxRays);
115 Domain_Free(D);
116 EVALUE_SET_DOMAIN(p->arr[2*i], T);
117 addeliminatedparams_evalue(&p->arr[2*i+1], CT);
118 }
119 }
120
121 static int mod_rational_smaller(evalue *e1, evalue *e2)
122 {
123 int r;
124 Value m;
125 Value m2;
126 value_init(m);
127 value_init(m2);
128
129 assert(value_notzero_p(e1->d));
130 assert(value_notzero_p(e2->d));
131 value_multiply(m, e1->x.n, e2->d);
132 value_multiply(m2, e2->x.n, e1->d);
133 if (value_lt(m, m2))
134 r = 1;
135 else if (value_gt(m, m2))
136 r = 0;
137 else
138 r = -1;
139 value_clear(m);
140 value_clear(m2);
141
142 return r;
143 }
144
145 static int mod_term_smaller_r(evalue *e1, evalue *e2)
146 {
147 if (value_notzero_p(e1->d)) {
148 int r;
149 if (value_zero_p(e2->d))
150 return 1;
151 r = mod_rational_smaller(e1, e2);
152 return r == -1 ? 0 : r;
153 }
154 if (value_notzero_p(e2->d))
155 return 0;
156 if (e1->x.p->pos < e2->x.p->pos)
157 return 1;
158 else if (e1->x.p->pos > e2->x.p->pos)
159 return 0;
160 else {
161 int r = mod_rational_smaller(&e1->x.p->arr[1], &e2->x.p->arr[1]);
162 return r == -1
163 ? mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0])
164 : r;
165 }
166 }
167
168 static int mod_term_smaller(evalue *e1, evalue *e2)
169 {
170 assert(value_zero_p(e1->d));
171 assert(value_zero_p(e2->d));
172 assert(e1->x.p->type == fractional || e1->x.p->type == flooring);
173 assert(e2->x.p->type == fractional || e2->x.p->type == flooring);
174 return mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0]);
175 }
176
177 /* Negative pos means inequality */
178 /* s is negative of substitution if m is not zero */
179 struct fixed_param {
180 int pos;
181 evalue s;
182 Value d;
183 Value m;
184 };
185
186 struct subst {
187 struct fixed_param *fixed;
188 int n;
189 int max;
190 };
191
192 static int relations_depth(evalue *e)
193 {
194 int d;
195
196 for (d = 0;
197 value_zero_p(e->d) && e->x.p->type == relation;
198 e = &e->x.p->arr[1], ++d);
199 return d;
200 }
201
202 static void poly_denom(evalue *p, Value *d)
203 {
204 value_set_si(*d, 1);
205
206 while (value_zero_p(p->d)) {
207 assert(p->x.p->type == polynomial);
208 assert(p->x.p->size == 2);
209 assert(value_notzero_p(p->x.p->arr[1].d));
210 value_lcm(*d, p->x.p->arr[1].d, d);
211 p = &p->x.p->arr[0];
212 }
213 value_lcm(*d, p->d, d);
214 }
215
216 #define EVALUE_IS_ONE(ev) (value_one_p((ev).d) && value_one_p((ev).x.n))
217
218 static void realloc_substitution(struct subst *s, int d)
219 {
220 struct fixed_param *n;
221 int i;
222 NALLOC(n, s->max+d);
223 for (i = 0; i < s->n; ++i)
224 n[i] = s->fixed[i];
225 free(s->fixed);
226 s->fixed = n;
227 s->max += d;
228 }
229
230 static int add_modulo_substitution(struct subst *s, evalue *r)
231 {
232 evalue *p;
233 evalue *f;
234 evalue *m;
235
236 assert(value_zero_p(r->d) && r->x.p->type == relation);
237 m = &r->x.p->arr[0];
238
239 /* May have been reduced already */
240 if (value_notzero_p(m->d))
241 return 0;
242
243 assert(value_zero_p(m->d) && m->x.p->type == fractional);
244 assert(m->x.p->size == 3);
245
246 /* fractional was inverted during reduction
247 * invert it back and move constant in
248 */
249 if (!EVALUE_IS_ONE(m->x.p->arr[2])) {
250 assert(value_pos_p(m->x.p->arr[2].d));
251 assert(value_mone_p(m->x.p->arr[2].x.n));
252 value_set_si(m->x.p->arr[2].x.n, 1);
253 value_increment(m->x.p->arr[1].x.n, m->x.p->arr[1].x.n);
254 assert(value_eq(m->x.p->arr[1].x.n, m->x.p->arr[1].d));
255 value_set_si(m->x.p->arr[1].x.n, 1);
256 eadd(&m->x.p->arr[1], &m->x.p->arr[0]);
257 value_set_si(m->x.p->arr[1].x.n, 0);
258 value_set_si(m->x.p->arr[1].d, 1);
259 }
260
261 /* Oops. Nested identical relations. */
262 if (!EVALUE_IS_ZERO(m->x.p->arr[1]))
263 return 0;
264
265 if (s->n >= s->max) {
266 int d = relations_depth(r);
267 realloc_substitution(s, d);
268 }
269
270 p = &m->x.p->arr[0];
271 assert(value_zero_p(p->d) && p->x.p->type == polynomial);
272 assert(p->x.p->size == 2);
273 f = &p->x.p->arr[1];
274
275 assert(value_pos_p(f->x.n));
276
277 value_init(s->fixed[s->n].m);
278 value_assign(s->fixed[s->n].m, f->d);
279 s->fixed[s->n].pos = p->x.p->pos;
280 value_init(s->fixed[s->n].d);
281 value_assign(s->fixed[s->n].d, f->x.n);
282 value_init(s->fixed[s->n].s.d);
283 evalue_copy(&s->fixed[s->n].s, &p->x.p->arr[0]);
284 ++s->n;
285
286 return 1;
287 }
288
289 static int type_offset(enode *p)
290 {
291 return p->type == fractional ? 1 :
292 p->type == flooring ? 1 : 0;
293 }
294
295 static void reorder_terms(enode *p, evalue *v)
296 {
297 int i;
298 int offset = type_offset(p);
299
300 for (i = p->size-1; i >= offset+1; i--) {
301 emul(v, &p->arr[i]);
302 eadd(&p->arr[i], &p->arr[i-1]);
303 free_evalue_refs(&(p->arr[i]));
304 }
305 p->size = offset+1;
306 free_evalue_refs(v);
307 }
308
309 void _reduce_evalue (evalue *e, struct subst *s, int fract) {
310
311 enode *p;
312 int i, j, k;
313 int add;
314
315 if (value_notzero_p(e->d)) {
316 if (fract)
317 mpz_fdiv_r(e->x.n, e->x.n, e->d);
318 return; /* a rational number, its already reduced */
319 }
320
321 if(!(p = e->x.p))
322 return; /* hum... an overflow probably occured */
323
324 /* First reduce the components of p */
325 add = p->type == relation;
326 for (i=0; i<p->size; i++) {
327 if (add && i == 1)
328 add = add_modulo_substitution(s, e);
329
330 if (i == 0 && p->type==fractional)
331 _reduce_evalue(&p->arr[i], s, 1);
332 else
333 _reduce_evalue(&p->arr[i], s, fract);
334
335 if (add && i == p->size-1) {
336 --s->n;
337 value_clear(s->fixed[s->n].m);
338 value_clear(s->fixed[s->n].d);
339 free_evalue_refs(&s->fixed[s->n].s);
340 } else if (add && i == 1)
341 s->fixed[s->n-1].pos *= -1;
342 }
343
344 if (p->type==periodic) {
345
346 /* Try to reduce the period */
347 for (i=1; i<=(p->size)/2; i++) {
348 if ((p->size % i)==0) {
349
350 /* Can we reduce the size to i ? */
351 for (j=0; j<i; j++)
352 for (k=j+i; k<e->x.p->size; k+=i)
353 if (!eequal(&p->arr[j], &p->arr[k])) goto you_lose;
354
355 /* OK, lets do it */
356 for (j=i; j<p->size; j++) free_evalue_refs(&p->arr[j]);
357 p->size = i;
358 break;
359
360 you_lose: /* OK, lets not do it */
361 continue;
362 }
363 }
364
365 /* Try to reduce its strength */
366 if (p->size == 1) {
367 value_clear(e->d);
368 memcpy(e,&p->arr[0],sizeof(evalue));
369 free(p);
370 }
371 }
372 else if (p->type==polynomial) {
373 for (k = 0; s && k < s->n; ++k) {
374 if (s->fixed[k].pos == p->pos) {
375 int divide = value_notone_p(s->fixed[k].d);
376 evalue d;
377
378 if (value_notzero_p(s->fixed[k].m)) {
379 if (!fract)
380 continue;
381 assert(p->size == 2);
382 if (divide && value_ne(s->fixed[k].d, p->arr[1].x.n))
383 continue;
384 if (!mpz_divisible_p(s->fixed[k].m, p->arr[1].d))
385 continue;
386 divide = 1;
387 }
388
389 if (divide) {
390 value_init(d.d);
391 value_assign(d.d, s->fixed[k].d);
392 value_init(d.x.n);
393 if (value_notzero_p(s->fixed[k].m))
394 value_oppose(d.x.n, s->fixed[k].m);
395 else
396 value_set_si(d.x.n, 1);
397 }
398
399 for (i=p->size-1;i>=1;i--) {
400 emul(&s->fixed[k].s, &p->arr[i]);
401 if (divide)
402 emul(&d, &p->arr[i]);
403 eadd(&p->arr[i], &p->arr[i-1]);
404 free_evalue_refs(&(p->arr[i]));
405 }
406 p->size = 1;
407 _reduce_evalue(&p->arr[0], s, fract);
408
409 if (divide)
410 free_evalue_refs(&d);
411
412 break;
413 }
414 }
415
416 /* Try to reduce the degree */
417 for (i=p->size-1;i>=1;i--) {
418 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
419 break;
420 /* Zero coefficient */
421 free_evalue_refs(&(p->arr[i]));
422 }
423 if (i+1<p->size)
424 p->size = i+1;
425
426 /* Try to reduce its strength */
427 if (p->size == 1) {
428 value_clear(e->d);
429 memcpy(e,&p->arr[0],sizeof(evalue));
430 free(p);
431 }
432 }
433 else if (p->type==fractional) {
434 int reorder = 0;
435 evalue v;
436
437 if (value_notzero_p(p->arr[0].d)) {
438 value_init(v.d);
439 value_assign(v.d, p->arr[0].d);
440 value_init(v.x.n);
441 mpz_fdiv_r(v.x.n, p->arr[0].x.n, p->arr[0].d);
442
443 reorder = 1;
444 } else {
445 evalue *f, *base;
446 evalue *pp = &p->arr[0];
447 assert(value_zero_p(pp->d) && pp->x.p->type == polynomial);
448 assert(pp->x.p->size == 2);
449
450 /* search for exact duplicate among the modulo inequalities */
451 do {
452 f = &pp->x.p->arr[1];
453 for (k = 0; s && k < s->n; ++k) {
454 if (-s->fixed[k].pos == pp->x.p->pos &&
455 value_eq(s->fixed[k].d, f->x.n) &&
456 value_eq(s->fixed[k].m, f->d) &&
457 eequal(&s->fixed[k].s, &pp->x.p->arr[0]))
458 break;
459 }
460 if (k < s->n) {
461 Value g;
462 value_init(g);
463
464 /* replace { E/m } by { (E-1)/m } + 1/m */
465 poly_denom(pp, &g);
466 if (reorder) {
467 evalue extra;
468 value_init(extra.d);
469 evalue_set_si(&extra, 1, 1);
470 value_assign(extra.d, g);
471 eadd(&extra, &v.x.p->arr[1]);
472 free_evalue_refs(&extra);
473
474 /* We've been going in circles; stop now */
475 if (value_ge(v.x.p->arr[1].x.n, v.x.p->arr[1].d)) {
476 free_evalue_refs(&v);
477 value_init(v.d);
478 evalue_set_si(&v, 0, 1);
479 break;
480 }
481 } else {
482 value_init(v.d);
483 value_set_si(v.d, 0);
484 v.x.p = new_enode(fractional, 3, -1);
485 evalue_set_si(&v.x.p->arr[1], 1, 1);
486 value_assign(v.x.p->arr[1].d, g);
487 evalue_set_si(&v.x.p->arr[2], 1, 1);
488 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
489 }
490
491 for (f = &v.x.p->arr[0]; value_zero_p(f->d);
492 f = &f->x.p->arr[0])
493 ;
494 value_division(f->d, g, f->d);
495 value_multiply(f->x.n, f->x.n, f->d);
496 value_assign(f->d, g);
497 value_decrement(f->x.n, f->x.n);
498 mpz_fdiv_r(f->x.n, f->x.n, f->d);
499
500 Gcd(f->d, f->x.n, &g);
501 value_division(f->d, f->d, g);
502 value_division(f->x.n, f->x.n, g);
503
504 value_clear(g);
505 pp = &v.x.p->arr[0];
506
507 reorder = 1;
508 }
509 } while (k < s->n);
510
511 /* reduction may have made this fractional arg smaller */
512 i = reorder ? p->size : 1;
513 for ( ; i < p->size; ++i)
514 if (value_zero_p(p->arr[i].d) &&
515 p->arr[i].x.p->type == fractional &&
516 !mod_term_smaller(e, &p->arr[i]))
517 break;
518 if (i < p->size) {
519 value_init(v.d);
520 value_set_si(v.d, 0);
521 v.x.p = new_enode(fractional, 3, -1);
522 evalue_set_si(&v.x.p->arr[1], 0, 1);
523 evalue_set_si(&v.x.p->arr[2], 1, 1);
524 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
525
526 reorder = 1;
527 }
528
529 if (!reorder) {
530 int invert = 0;
531 Value twice;
532 value_init(twice);
533
534 for (pp = &p->arr[0]; value_zero_p(pp->d);
535 pp = &pp->x.p->arr[0]) {
536 f = &pp->x.p->arr[1];
537 assert(value_pos_p(f->d));
538 mpz_mul_ui(twice, f->x.n, 2);
539 if (value_lt(twice, f->d))
540 break;
541 if (value_eq(twice, f->d))
542 continue;
543 invert = 1;
544 break;
545 }
546
547 if (invert) {
548 value_init(v.d);
549 value_set_si(v.d, 0);
550 v.x.p = new_enode(fractional, 3, -1);
551 evalue_set_si(&v.x.p->arr[1], 0, 1);
552 poly_denom(&p->arr[0], &twice);
553 value_assign(v.x.p->arr[1].d, twice);
554 value_decrement(v.x.p->arr[1].x.n, twice);
555 evalue_set_si(&v.x.p->arr[2], -1, 1);
556 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
557
558 for (pp = &v.x.p->arr[0]; value_zero_p(pp->d);
559 pp = &pp->x.p->arr[0]) {
560 f = &pp->x.p->arr[1];
561 value_oppose(f->x.n, f->x.n);
562 mpz_fdiv_r(f->x.n, f->x.n, f->d);
563 }
564 value_division(pp->d, twice, pp->d);
565 value_multiply(pp->x.n, pp->x.n, pp->d);
566 value_assign(pp->d, twice);
567 value_oppose(pp->x.n, pp->x.n);
568 value_decrement(pp->x.n, pp->x.n);
569 mpz_fdiv_r(pp->x.n, pp->x.n, pp->d);
570
571 reorder = 1;
572 }
573
574 value_clear(twice);
575 }
576 }
577
578 if (reorder) {
579 reorder_terms(p, &v);
580 _reduce_evalue(&p->arr[1], s, fract);
581 }
582
583 /* Try to reduce the degree */
584 for (i=p->size-1;i>=2;i--) {
585 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
586 break;
587 /* Zero coefficient */
588 free_evalue_refs(&(p->arr[i]));
589 }
590 if (i+1<p->size)
591 p->size = i+1;
592
593 /* Try to reduce its strength */
594 if (p->size == 2) {
595 value_clear(e->d);
596 memcpy(e,&p->arr[1],sizeof(evalue));
597 free_evalue_refs(&(p->arr[0]));
598 free(p);
599 }
600 }
601 else if (p->type == flooring) {
602 /* Try to reduce the degree */
603 for (i=p->size-1;i>=2;i--) {
604 if (!EVALUE_IS_ZERO(p->arr[i]))
605 break;
606 /* Zero coefficient */
607 free_evalue_refs(&(p->arr[i]));
608 }
609 if (i+1<p->size)
610 p->size = i+1;
611
612 /* Try to reduce its strength */
613 if (p->size == 2) {
614 value_clear(e->d);
615 memcpy(e,&p->arr[1],sizeof(evalue));
616 free_evalue_refs(&(p->arr[0]));
617 free(p);
618 }
619 }
620 else if (p->type == relation) {
621 if (p->size == 3 && eequal(&p->arr[1], &p->arr[2])) {
622 free_evalue_refs(&(p->arr[2]));
623 free_evalue_refs(&(p->arr[0]));
624 p->size = 2;
625 value_clear(e->d);
626 *e = p->arr[1];
627 free(p);
628 return;
629 }
630 if (p->size == 3 && EVALUE_IS_ZERO(p->arr[2])) {
631 free_evalue_refs(&(p->arr[2]));
632 p->size = 2;
633 }
634 if (p->size == 2 && EVALUE_IS_ZERO(p->arr[1])) {
635 free_evalue_refs(&(p->arr[1]));
636 free_evalue_refs(&(p->arr[0]));
637 evalue_set_si(e, 0, 1);
638 free(p);
639 } else {
640 int reduced = 0;
641 evalue *m;
642 m = &p->arr[0];
643
644 /* Relation was reduced by means of an identical
645 * inequality => remove
646 */
647 if (value_zero_p(m->d) && !EVALUE_IS_ZERO(m->x.p->arr[1]))
648 reduced = 1;
649
650 if (reduced || value_notzero_p(p->arr[0].d)) {
651 if (!reduced && value_zero_p(p->arr[0].x.n)) {
652 value_clear(e->d);
653 memcpy(e,&p->arr[1],sizeof(evalue));
654 if (p->size == 3)
655 free_evalue_refs(&(p->arr[2]));
656 } else {
657 if (p->size == 3) {
658 value_clear(e->d);
659 memcpy(e,&p->arr[2],sizeof(evalue));
660 } else
661 evalue_set_si(e, 0, 1);
662 free_evalue_refs(&(p->arr[1]));
663 }
664 free_evalue_refs(&(p->arr[0]));
665 free(p);
666 }
667 }
668 }
669 return;
670 } /* reduce_evalue */
671
672 static void add_substitution(struct subst *s, Value *row, unsigned dim)
673 {
674 int k, l;
675 evalue *r;
676
677 for (k = 0; k < dim; ++k)
678 if (value_notzero_p(row[k+1]))
679 break;
680
681 Vector_Normalize_Positive(row+1, dim+1, k);
682 assert(s->n < s->max);
683 value_init(s->fixed[s->n].d);
684 value_init(s->fixed[s->n].m);
685 value_assign(s->fixed[s->n].d, row[k+1]);
686 s->fixed[s->n].pos = k+1;
687 value_set_si(s->fixed[s->n].m, 0);
688 r = &s->fixed[s->n].s;
689 value_init(r->d);
690 for (l = k+1; l < dim; ++l)
691 if (value_notzero_p(row[l+1])) {
692 value_set_si(r->d, 0);
693 r->x.p = new_enode(polynomial, 2, l + 1);
694 value_init(r->x.p->arr[1].x.n);
695 value_oppose(r->x.p->arr[1].x.n, row[l+1]);
696 value_set_si(r->x.p->arr[1].d, 1);
697 r = &r->x.p->arr[0];
698 }
699 value_init(r->x.n);
700 value_oppose(r->x.n, row[dim+1]);
701 value_set_si(r->d, 1);
702 ++s->n;
703 }
704
705 void reduce_evalue (evalue *e) {
706 struct subst s = { NULL, 0, 0 };
707
708 if (value_notzero_p(e->d))
709 return; /* a rational number, its already reduced */
710
711 if (e->x.p->type == partition) {
712 int i;
713 unsigned dim = -1;
714 for (i = 0; i < e->x.p->size/2; ++i) {
715 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
716 s.n = 0;
717 /* This shouldn't really happen;
718 * Empty domains should not be added.
719 */
720 if (emptyQ(D))
721 goto discard;
722
723 dim = D->Dimension;
724 if (D->next)
725 D = DomainConvex(D, 0);
726 if (!D->next && D->NbEq) {
727 int j, k;
728 if (s.max < dim) {
729 if (s.max != 0)
730 realloc_substitution(&s, dim);
731 else {
732 int d = relations_depth(&e->x.p->arr[2*i+1]);
733 s.max = dim+d;
734 NALLOC(s.fixed, s.max);
735 }
736 }
737 for (j = 0; j < D->NbEq; ++j)
738 add_substitution(&s, D->Constraint[j], dim);
739 }
740 if (D != EVALUE_DOMAIN(e->x.p->arr[2*i]))
741 Domain_Free(D);
742 _reduce_evalue(&e->x.p->arr[2*i+1], &s, 0);
743 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
744 discard:
745 free_evalue_refs(&e->x.p->arr[2*i+1]);
746 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
747 value_clear(e->x.p->arr[2*i].d);
748 e->x.p->size -= 2;
749 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
750 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
751 --i;
752 }
753 if (s.n != 0) {
754 int j;
755 for (j = 0; j < s.n; ++j) {
756 value_clear(s.fixed[j].d);
757 value_clear(s.fixed[j].m);
758 free_evalue_refs(&s.fixed[j].s);
759 }
760 }
761 }
762 if (e->x.p->size == 0) {
763 free(e->x.p);
764 evalue_set_si(e, 0, 1);
765 }
766 } else
767 _reduce_evalue(e, &s, 0);
768 if (s.max != 0)
769 free(s.fixed);
770 }
771
772 void print_evalue(FILE *DST,evalue *e,char **pname) {
773
774 if(value_notzero_p(e->d)) {
775 if(value_notone_p(e->d)) {
776 value_print(DST,VALUE_FMT,e->x.n);
777 fprintf(DST,"/");
778 value_print(DST,VALUE_FMT,e->d);
779 }
780 else {
781 value_print(DST,VALUE_FMT,e->x.n);
782 }
783 }
784 else
785 print_enode(DST,e->x.p,pname);
786 return;
787 } /* print_evalue */
788
789 void print_enode(FILE *DST,enode *p,char **pname) {
790
791 int i;
792
793 if (!p) {
794 fprintf(DST, "NULL");
795 return;
796 }
797 switch (p->type) {
798 case evector:
799 fprintf(DST, "{ ");
800 for (i=0; i<p->size; i++) {
801 print_evalue(DST, &p->arr[i], pname);
802 if (i!=(p->size-1))
803 fprintf(DST, ", ");
804 }
805 fprintf(DST, " }\n");
806 break;
807 case polynomial:
808 fprintf(DST, "( ");
809 for (i=p->size-1; i>=0; i--) {
810 print_evalue(DST, &p->arr[i], pname);
811 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
812 else if (i>1)
813 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
814 }
815 fprintf(DST, " )\n");
816 break;
817 case periodic:
818 fprintf(DST, "[ ");
819 for (i=0; i<p->size; i++) {
820 print_evalue(DST, &p->arr[i], pname);
821 if (i!=(p->size-1)) fprintf(DST, ", ");
822 }
823 fprintf(DST," ]_%s", pname[p->pos-1]);
824 break;
825 case flooring:
826 case fractional:
827 fprintf(DST, "( ");
828 for (i=p->size-1; i>=1; i--) {
829 print_evalue(DST, &p->arr[i], pname);
830 if (i >= 2) {
831 fprintf(DST, " * ");
832 fprintf(DST, p->type == flooring ? "[" : "{");
833 print_evalue(DST, &p->arr[0], pname);
834 fprintf(DST, p->type == flooring ? "]" : "}");
835 if (i>2)
836 fprintf(DST, "^%d + ", i-1);
837 else
838 fprintf(DST, " + ");
839 }
840 }
841 fprintf(DST, " )\n");
842 break;
843 case relation:
844 fprintf(DST, "[ ");
845 print_evalue(DST, &p->arr[0], pname);
846 fprintf(DST, "= 0 ] * \n");
847 print_evalue(DST, &p->arr[1], pname);
848 if (p->size > 2) {
849 fprintf(DST, " +\n [ ");
850 print_evalue(DST, &p->arr[0], pname);
851 fprintf(DST, "!= 0 ] * \n");
852 print_evalue(DST, &p->arr[2], pname);
853 }
854 break;
855 case partition: {
856 char **names = pname;
857 int maxdim = EVALUE_DOMAIN(p->arr[0])->Dimension;
858 if (p->pos < maxdim) {
859 NALLOC(names, maxdim);
860 for (i = 0; i < p->pos; ++i)
861 names[i] = pname[i];
862 for ( ; i < maxdim; ++i) {
863 NALLOC(names[i], 10);
864 snprintf(names[i], 10, "_p%d", i);
865 }
866 }
867
868 for (i=0; i<p->size/2; i++) {
869 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), names);
870 print_evalue(DST, &p->arr[2*i+1], names);
871 }
872
873 if (p->pos < maxdim) {
874 for (i = p->pos ; i < maxdim; ++i)
875 free(names[i]);
876 free(names);
877 }
878
879 break;
880 }
881 default:
882 assert(0);
883 }
884 return;
885 } /* print_enode */
886
887 static void eadd_rev(evalue *e1, evalue *res)
888 {
889 evalue ev;
890 value_init(ev.d);
891 evalue_copy(&ev, e1);
892 eadd(res, &ev);
893 free_evalue_refs(res);
894 *res = ev;
895 }
896
897 static void eadd_rev_cst (evalue *e1, evalue *res)
898 {
899 evalue ev;
900 value_init(ev.d);
901 evalue_copy(&ev, e1);
902 eadd(res, &ev.x.p->arr[type_offset(ev.x.p)]);
903 free_evalue_refs(res);
904 *res = ev;
905 }
906
907 struct section { Polyhedron * D; evalue E; };
908
909 void eadd_partitions (evalue *e1,evalue *res)
910 {
911 int n, i, j;
912 Polyhedron *d, *fd;
913 struct section *s;
914 s = (struct section *)
915 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
916 sizeof(struct section));
917 assert(s);
918 assert(e1->x.p->pos == res->x.p->pos);
919 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
920 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
921
922 n = 0;
923 for (j = 0; j < e1->x.p->size/2; ++j) {
924 assert(res->x.p->size >= 2);
925 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
926 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
927 if (!emptyQ(fd))
928 for (i = 1; i < res->x.p->size/2; ++i) {
929 Polyhedron *t = fd;
930 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
931 Domain_Free(t);
932 if (emptyQ(fd))
933 break;
934 }
935 if (emptyQ(fd)) {
936 Domain_Free(fd);
937 continue;
938 }
939 value_init(s[n].E.d);
940 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
941 s[n].D = fd;
942 ++n;
943 }
944 for (i = 0; i < res->x.p->size/2; ++i) {
945 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
946 for (j = 0; j < e1->x.p->size/2; ++j) {
947 Polyhedron *t;
948 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
949 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
950 if (emptyQ(d)) {
951 Domain_Free(d);
952 continue;
953 }
954 t = fd;
955 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
956 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
957 Domain_Free(t);
958 value_init(s[n].E.d);
959 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
960 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
961 s[n].D = d;
962 ++n;
963 }
964 if (!emptyQ(fd)) {
965 s[n].E = res->x.p->arr[2*i+1];
966 s[n].D = fd;
967 ++n;
968 } else {
969 free_evalue_refs(&res->x.p->arr[2*i+1]);
970 Domain_Free(fd);
971 }
972 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
973 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
974 value_clear(res->x.p->arr[2*i].d);
975 }
976
977 free(res->x.p);
978 assert(n > 0);
979 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
980 for (j = 0; j < n; ++j) {
981 s[j].D = DomainConstraintSimplify(s[j].D, 0);
982 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
983 value_clear(res->x.p->arr[2*j+1].d);
984 res->x.p->arr[2*j+1] = s[j].E;
985 }
986
987 free(s);
988 }
989
990 static void explicit_complement(evalue *res)
991 {
992 enode *rel = new_enode(relation, 3, 0);
993 assert(rel);
994 value_clear(rel->arr[0].d);
995 rel->arr[0] = res->x.p->arr[0];
996 value_clear(rel->arr[1].d);
997 rel->arr[1] = res->x.p->arr[1];
998 value_set_si(rel->arr[2].d, 1);
999 value_init(rel->arr[2].x.n);
1000 value_set_si(rel->arr[2].x.n, 0);
1001 free(res->x.p);
1002 res->x.p = rel;
1003 }
1004
1005 void eadd(evalue *e1,evalue *res) {
1006
1007 int i;
1008 if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
1009 /* Add two rational numbers */
1010 Value g,m1,m2;
1011 value_init(g);
1012 value_init(m1);
1013 value_init(m2);
1014
1015 value_multiply(m1,e1->x.n,res->d);
1016 value_multiply(m2,res->x.n,e1->d);
1017 value_addto(res->x.n,m1,m2);
1018 value_multiply(res->d,e1->d,res->d);
1019 Gcd(res->x.n,res->d,&g);
1020 if (value_notone_p(g)) {
1021 value_division(res->d,res->d,g);
1022 value_division(res->x.n,res->x.n,g);
1023 }
1024 value_clear(g); value_clear(m1); value_clear(m2);
1025 return ;
1026 }
1027 else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
1028 switch (res->x.p->type) {
1029 case polynomial:
1030 /* Add the constant to the constant term of a polynomial*/
1031 eadd(e1, &res->x.p->arr[0]);
1032 return ;
1033 case periodic:
1034 /* Add the constant to all elements of a periodic number */
1035 for (i=0; i<res->x.p->size; i++) {
1036 eadd(e1, &res->x.p->arr[i]);
1037 }
1038 return ;
1039 case evector:
1040 fprintf(stderr, "eadd: cannot add const with vector\n");
1041 return;
1042 case flooring:
1043 case fractional:
1044 eadd(e1, &res->x.p->arr[1]);
1045 return ;
1046 case partition:
1047 assert(EVALUE_IS_ZERO(*e1));
1048 break; /* Do nothing */
1049 case relation:
1050 /* Create (zero) complement if needed */
1051 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
1052 explicit_complement(res);
1053 for (i = 1; i < res->x.p->size; ++i)
1054 eadd(e1, &res->x.p->arr[i]);
1055 break;
1056 default:
1057 assert(0);
1058 }
1059 }
1060 /* add polynomial or periodic to constant
1061 * you have to exchange e1 and res, before doing addition */
1062
1063 else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
1064 eadd_rev(e1, res);
1065 return;
1066 }
1067 else { // ((e1->d==0) && (res->d==0))
1068 assert(!((e1->x.p->type == partition) ^
1069 (res->x.p->type == partition)));
1070 if (e1->x.p->type == partition) {
1071 eadd_partitions(e1, res);
1072 return;
1073 }
1074 if (e1->x.p->type == relation &&
1075 (res->x.p->type != relation ||
1076 mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1077 eadd_rev(e1, res);
1078 return;
1079 }
1080 if (res->x.p->type == relation) {
1081 if (e1->x.p->type == relation &&
1082 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1083 if (res->x.p->size < 3 && e1->x.p->size == 3)
1084 explicit_complement(res);
1085 if (e1->x.p->size < 3 && res->x.p->size == 3)
1086 explicit_complement(e1);
1087 for (i = 1; i < res->x.p->size; ++i)
1088 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1089 return;
1090 }
1091 if (res->x.p->size < 3)
1092 explicit_complement(res);
1093 for (i = 1; i < res->x.p->size; ++i)
1094 eadd(e1, &res->x.p->arr[i]);
1095 return;
1096 }
1097 if ((e1->x.p->type != res->x.p->type) ) {
1098 /* adding to evalues of different type. two cases are possible
1099 * res is periodic and e1 is polynomial, you have to exchange
1100 * e1 and res then to add e1 to the constant term of res */
1101 if (e1->x.p->type == polynomial) {
1102 eadd_rev_cst(e1, res);
1103 }
1104 else if (res->x.p->type == polynomial) {
1105 /* res is polynomial and e1 is periodic,
1106 add e1 to the constant term of res */
1107
1108 eadd(e1,&res->x.p->arr[0]);
1109 } else
1110 assert(0);
1111
1112 return;
1113 }
1114 else if (e1->x.p->pos != res->x.p->pos ||
1115 ((res->x.p->type == fractional ||
1116 res->x.p->type == flooring) &&
1117 !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1118 /* adding evalues of different position (i.e function of different unknowns
1119 * to case are possible */
1120
1121 switch (res->x.p->type) {
1122 case flooring:
1123 case fractional:
1124 if(mod_term_smaller(res, e1))
1125 eadd(e1,&res->x.p->arr[1]);
1126 else
1127 eadd_rev_cst(e1, res);
1128 return;
1129 case polynomial: // res and e1 are polynomials
1130 // add e1 to the constant term of res
1131
1132 if(res->x.p->pos < e1->x.p->pos)
1133 eadd(e1,&res->x.p->arr[0]);
1134 else
1135 eadd_rev_cst(e1, res);
1136 // value_clear(g); value_clear(m1); value_clear(m2);
1137 return;
1138 case periodic: // res and e1 are pointers to periodic numbers
1139 //add e1 to all elements of res
1140
1141 if(res->x.p->pos < e1->x.p->pos)
1142 for (i=0;i<res->x.p->size;i++) {
1143 eadd(e1,&res->x.p->arr[i]);
1144 }
1145 else
1146 eadd_rev(e1, res);
1147 return;
1148 default:
1149 assert(0);
1150 }
1151 }
1152
1153
1154 //same type , same pos and same size
1155 if (e1->x.p->size == res->x.p->size) {
1156 // add any element in e1 to the corresponding element in res
1157 i = type_offset(res->x.p);
1158 if (i == 1)
1159 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1160 for (; i<res->x.p->size; i++) {
1161 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1162 }
1163 return ;
1164 }
1165
1166 /* Sizes are different */
1167 switch(res->x.p->type) {
1168 case polynomial:
1169 case flooring:
1170 case fractional:
1171 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1172 /* new enode and add res to that new node. If you do not do */
1173 /* that, you lose the the upper weight part of e1 ! */
1174
1175 if(e1->x.p->size > res->x.p->size)
1176 eadd_rev(e1, res);
1177 else {
1178 i = type_offset(res->x.p);
1179 if (i == 1)
1180 assert(eequal(&e1->x.p->arr[0],
1181 &res->x.p->arr[0]));
1182 for (; i<e1->x.p->size ; i++) {
1183 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1184 }
1185
1186 return ;
1187 }
1188 break;
1189
1190 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1191 case periodic:
1192 {
1193 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1194 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1195 to add periodicaly elements of e1 to elements of ne, and finaly to
1196 return ne. */
1197 int x,y,p;
1198 Value ex, ey ,ep;
1199 evalue *ne;
1200 value_init(ex); value_init(ey);value_init(ep);
1201 x=e1->x.p->size;
1202 y= res->x.p->size;
1203 value_set_si(ex,e1->x.p->size);
1204 value_set_si(ey,res->x.p->size);
1205 value_assign (ep,*Lcm(ex,ey));
1206 p=(int)mpz_get_si(ep);
1207 ne= (evalue *) malloc (sizeof(evalue));
1208 value_init(ne->d);
1209 value_set_si( ne->d,0);
1210
1211 ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
1212 for(i=0;i<p;i++) {
1213 evalue_copy(&ne->x.p->arr[i], &res->x.p->arr[i%y]);
1214 }
1215 for(i=0;i<p;i++) {
1216 eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
1217 }
1218
1219 value_assign(res->d, ne->d);
1220 res->x.p=ne->x.p;
1221
1222 return ;
1223 }
1224 case evector:
1225 fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
1226 return ;
1227 default:
1228 assert(0);
1229 }
1230 }
1231 return ;
1232 }/* eadd */
1233
1234 static void emul_rev (evalue *e1, evalue *res)
1235 {
1236 evalue ev;
1237 value_init(ev.d);
1238 evalue_copy(&ev, e1);
1239 emul(res, &ev);
1240 free_evalue_refs(res);
1241 *res = ev;
1242 }
1243
1244 static void emul_poly (evalue *e1, evalue *res)
1245 {
1246 int i, j, o = type_offset(res->x.p);
1247 evalue tmp;
1248 int size=(e1->x.p->size + res->x.p->size - o - 1);
1249 value_init(tmp.d);
1250 value_set_si(tmp.d,0);
1251 tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
1252 if (o)
1253 evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
1254 for (i=o; i < e1->x.p->size; i++) {
1255 evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
1256 emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
1257 }
1258 for (; i<size; i++)
1259 evalue_set_si(&tmp.x.p->arr[i], 0, 1);
1260 for (i=o+1; i<res->x.p->size; i++)
1261 for (j=o; j<e1->x.p->size; j++) {
1262 evalue ev;
1263 value_init(ev.d);
1264 evalue_copy(&ev, &e1->x.p->arr[j]);
1265 emul(&res->x.p->arr[i], &ev);
1266 eadd(&ev, &tmp.x.p->arr[i+j-o]);
1267 free_evalue_refs(&ev);
1268 }
1269 free_evalue_refs(res);
1270 *res = tmp;
1271 }
1272
1273 void emul_partitions (evalue *e1,evalue *res)
1274 {
1275 int n, i, j, k;
1276 Polyhedron *d;
1277 struct section *s;
1278 s = (struct section *)
1279 malloc((e1->x.p->size/2) * (res->x.p->size/2) *
1280 sizeof(struct section));
1281 assert(s);
1282 assert(e1->x.p->pos == res->x.p->pos);
1283 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
1284 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1285
1286 n = 0;
1287 for (i = 0; i < res->x.p->size/2; ++i) {
1288 for (j = 0; j < e1->x.p->size/2; ++j) {
1289 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1290 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1291 if (emptyQ(d)) {
1292 Domain_Free(d);
1293 continue;
1294 }
1295
1296 /* This code is only needed because the partitions
1297 are not true partitions.
1298 */
1299 for (k = 0; k < n; ++k) {
1300 if (DomainIncludes(s[k].D, d))
1301 break;
1302 if (DomainIncludes(d, s[k].D)) {
1303 Domain_Free(s[k].D);
1304 free_evalue_refs(&s[k].E);
1305 if (n > k)
1306 s[k] = s[--n];
1307 --k;
1308 }
1309 }
1310 if (k < n) {
1311 Domain_Free(d);
1312 continue;
1313 }
1314
1315 value_init(s[n].E.d);
1316 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1317 emul(&e1->x.p->arr[2*j+1], &s[n].E);
1318 s[n].D = d;
1319 ++n;
1320 }
1321 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1322 value_clear(res->x.p->arr[2*i].d);
1323 free_evalue_refs(&res->x.p->arr[2*i+1]);
1324 }
1325
1326 free(res->x.p);
1327 if (n == 0)
1328 evalue_set_si(res, 0, 1);
1329 else {
1330 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1331 for (j = 0; j < n; ++j) {
1332 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1333 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1334 value_clear(res->x.p->arr[2*j+1].d);
1335 res->x.p->arr[2*j+1] = s[j].E;
1336 }
1337 }
1338
1339 free(s);
1340 }
1341
1342 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1343
1344 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1345 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1346 * evalues is not treated here */
1347
1348 void emul (evalue *e1, evalue *res ){
1349 int i,j;
1350
1351 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
1352 fprintf(stderr, "emul: do not proced on evector type !\n");
1353 return;
1354 }
1355
1356 if (EVALUE_IS_ZERO(*res))
1357 return;
1358
1359 if (value_zero_p(e1->d) && e1->x.p->type == partition) {
1360 if (value_zero_p(res->d) && res->x.p->type == partition)
1361 emul_partitions(e1, res);
1362 else
1363 emul_rev(e1, res);
1364 } else if (value_zero_p(res->d) && res->x.p->type == partition) {
1365 for (i = 0; i < res->x.p->size/2; ++i)
1366 emul(e1, &res->x.p->arr[2*i+1]);
1367 } else
1368 if (value_zero_p(res->d) && res->x.p->type == relation) {
1369 if (value_zero_p(e1->d) && e1->x.p->type == relation &&
1370 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1371 if (res->x.p->size < 3 && e1->x.p->size == 3)
1372 explicit_complement(res);
1373 if (e1->x.p->size < 3 && res->x.p->size == 3)
1374 explicit_complement(e1);
1375 for (i = 1; i < res->x.p->size; ++i)
1376 emul(&e1->x.p->arr[i], &res->x.p->arr[i]);
1377 return;
1378 }
1379 for (i = 1; i < res->x.p->size; ++i)
1380 emul(e1, &res->x.p->arr[i]);
1381 } else
1382 if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
1383 switch(e1->x.p->type) {
1384 case polynomial:
1385 switch(res->x.p->type) {
1386 case polynomial:
1387 if(e1->x.p->pos == res->x.p->pos) {
1388 /* Product of two polynomials of the same variable */
1389 emul_poly(e1, res);
1390 return;
1391 }
1392 else {
1393 /* Product of two polynomials of different variables */
1394
1395 if(res->x.p->pos < e1->x.p->pos)
1396 for( i=0; i<res->x.p->size ; i++)
1397 emul(e1, &res->x.p->arr[i]);
1398 else
1399 emul_rev(e1, res);
1400
1401 return ;
1402 }
1403 case periodic:
1404 case flooring:
1405 case fractional:
1406 /* Product of a polynomial and a periodic or fractional */
1407 emul_rev(e1, res);
1408 return;
1409 default:
1410 assert(0);
1411 }
1412 case periodic:
1413 switch(res->x.p->type) {
1414 case periodic:
1415 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
1416 /* Product of two periodics of the same parameter and period */
1417
1418 for(i=0; i<res->x.p->size;i++)
1419 emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
1420
1421 return;
1422 }
1423 else{
1424 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
1425 /* Product of two periodics of the same parameter and different periods */
1426 evalue *newp;
1427 Value x,y,z;
1428 int ix,iy,lcm;
1429 value_init(x); value_init(y);value_init(z);
1430 ix=e1->x.p->size;
1431 iy=res->x.p->size;
1432 value_set_si(x,e1->x.p->size);
1433 value_set_si(y,res->x.p->size);
1434 value_assign (z,*Lcm(x,y));
1435 lcm=(int)mpz_get_si(z);
1436 newp= (evalue *) malloc (sizeof(evalue));
1437 value_init(newp->d);
1438 value_set_si( newp->d,0);
1439 newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
1440 for(i=0;i<lcm;i++) {
1441 evalue_copy(&newp->x.p->arr[i],
1442 &res->x.p->arr[i%iy]);
1443 }
1444 for(i=0;i<lcm;i++)
1445 emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
1446
1447 value_assign(res->d,newp->d);
1448 res->x.p=newp->x.p;
1449
1450 value_clear(x); value_clear(y);value_clear(z);
1451 return ;
1452 }
1453 else {
1454 /* Product of two periodics of different parameters */
1455
1456 if(res->x.p->pos < e1->x.p->pos)
1457 for(i=0; i<res->x.p->size; i++)
1458 emul(e1, &(res->x.p->arr[i]));
1459 else
1460 emul_rev(e1, res);
1461
1462 return;
1463 }
1464 }
1465 case polynomial:
1466 /* Product of a periodic and a polynomial */
1467
1468 for(i=0; i<res->x.p->size ; i++)
1469 emul(e1, &(res->x.p->arr[i]));
1470
1471 return;
1472
1473 }
1474 case flooring:
1475 case fractional:
1476 switch(res->x.p->type) {
1477 case polynomial:
1478 for(i=0; i<res->x.p->size ; i++)
1479 emul(e1, &(res->x.p->arr[i]));
1480 return;
1481 default:
1482 case periodic:
1483 assert(0);
1484 case flooring:
1485 case fractional:
1486 assert(e1->x.p->type == res->x.p->type);
1487 if (e1->x.p->pos == res->x.p->pos &&
1488 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1489 evalue d;
1490 value_init(d.d);
1491 poly_denom(&e1->x.p->arr[0], &d.d);
1492 if (e1->x.p->type != fractional || !value_two_p(d.d))
1493 emul_poly(e1, res);
1494 else {
1495 evalue tmp;
1496 value_init(d.x.n);
1497 value_set_si(d.x.n, 1);
1498 /* { x }^2 == { x }/2 */
1499 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1500 assert(e1->x.p->size == 3);
1501 assert(res->x.p->size == 3);
1502 value_init(tmp.d);
1503 evalue_copy(&tmp, &res->x.p->arr[2]);
1504 emul(&d, &tmp);
1505 eadd(&res->x.p->arr[1], &tmp);
1506 emul(&e1->x.p->arr[2], &tmp);
1507 emul(&e1->x.p->arr[1], res);
1508 eadd(&tmp, &res->x.p->arr[2]);
1509 free_evalue_refs(&tmp);
1510 value_clear(d.x.n);
1511 }
1512 value_clear(d.d);
1513 } else {
1514 if(mod_term_smaller(res, e1))
1515 for(i=1; i<res->x.p->size ; i++)
1516 emul(e1, &(res->x.p->arr[i]));
1517 else
1518 emul_rev(e1, res);
1519 return;
1520 }
1521 }
1522 break;
1523 case relation:
1524 emul_rev(e1, res);
1525 break;
1526 default:
1527 assert(0);
1528 }
1529 }
1530 else {
1531 if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
1532 /* Product of two rational numbers */
1533
1534 Value g;
1535 value_init(g);
1536 value_multiply(res->d,e1->d,res->d);
1537 value_multiply(res->x.n,e1->x.n,res->x.n );
1538 Gcd(res->x.n, res->d,&g);
1539 if (value_notone_p(g)) {
1540 value_division(res->d,res->d,g);
1541 value_division(res->x.n,res->x.n,g);
1542 }
1543 value_clear(g);
1544 return ;
1545 }
1546 else {
1547 if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
1548 /* Product of an expression (polynomial or peririodic) and a rational number */
1549
1550 emul_rev(e1, res);
1551 return ;
1552 }
1553 else {
1554 /* Product of a rationel number and an expression (polynomial or peririodic) */
1555
1556 i = type_offset(res->x.p);
1557 for (; i<res->x.p->size; i++)
1558 emul(e1, &res->x.p->arr[i]);
1559
1560 return ;
1561 }
1562 }
1563 }
1564
1565 return ;
1566 }
1567
1568 /* Frees mask content ! */
1569 void emask(evalue *mask, evalue *res) {
1570 int n, i, j;
1571 Polyhedron *d, *fd;
1572 struct section *s;
1573 evalue mone;
1574 int pos;
1575
1576 if (EVALUE_IS_ZERO(*res)) {
1577 free_evalue_refs(mask);
1578 return;
1579 }
1580
1581 assert(value_zero_p(mask->d));
1582 assert(mask->x.p->type == partition);
1583 assert(value_zero_p(res->d));
1584 assert(res->x.p->type == partition);
1585 assert(mask->x.p->pos == res->x.p->pos);
1586 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1587 assert(mask->x.p->pos == EVALUE_DOMAIN(mask->x.p->arr[0])->Dimension);
1588 pos = res->x.p->pos;
1589
1590 s = (struct section *)
1591 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1592 sizeof(struct section));
1593 assert(s);
1594
1595 value_init(mone.d);
1596 evalue_set_si(&mone, -1, 1);
1597
1598 n = 0;
1599 for (j = 0; j < res->x.p->size/2; ++j) {
1600 assert(mask->x.p->size >= 2);
1601 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1602 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1603 if (!emptyQ(fd))
1604 for (i = 1; i < mask->x.p->size/2; ++i) {
1605 Polyhedron *t = fd;
1606 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1607 Domain_Free(t);
1608 if (emptyQ(fd))
1609 break;
1610 }
1611 if (emptyQ(fd)) {
1612 Domain_Free(fd);
1613 continue;
1614 }
1615 value_init(s[n].E.d);
1616 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1617 s[n].D = fd;
1618 ++n;
1619 }
1620 for (i = 0; i < mask->x.p->size/2; ++i) {
1621 if (EVALUE_IS_ONE(mask->x.p->arr[2*i+1]))
1622 continue;
1623
1624 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1625 eadd(&mone, &mask->x.p->arr[2*i+1]);
1626 emul(&mone, &mask->x.p->arr[2*i+1]);
1627 for (j = 0; j < res->x.p->size/2; ++j) {
1628 Polyhedron *t;
1629 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1630 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1631 if (emptyQ(d)) {
1632 Domain_Free(d);
1633 continue;
1634 }
1635 t = fd;
1636 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1637 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1638 Domain_Free(t);
1639 value_init(s[n].E.d);
1640 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1641 emul(&mask->x.p->arr[2*i+1], &s[n].E);
1642 s[n].D = d;
1643 ++n;
1644 }
1645
1646 if (!emptyQ(fd)) {
1647 /* Just ignore; this may have been previously masked off */
1648 }
1649 if (fd != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1650 Domain_Free(fd);
1651 }
1652
1653 free_evalue_refs(&mone);
1654 free_evalue_refs(mask);
1655 free_evalue_refs(res);
1656 value_init(res->d);
1657 if (n == 0)
1658 evalue_set_si(res, 0, 1);
1659 else {
1660 res->x.p = new_enode(partition, 2*n, pos);
1661 for (j = 0; j < n; ++j) {
1662 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1663 value_clear(res->x.p->arr[2*j+1].d);
1664 res->x.p->arr[2*j+1] = s[j].E;
1665 }
1666 }
1667
1668 free(s);
1669 }
1670
1671 void evalue_copy(evalue *dst, evalue *src)
1672 {
1673 value_assign(dst->d, src->d);
1674 if(value_notzero_p(src->d)) {
1675 value_init(dst->x.n);
1676 value_assign(dst->x.n, src->x.n);
1677 } else
1678 dst->x.p = ecopy(src->x.p);
1679 }
1680
1681 enode *new_enode(enode_type type,int size,int pos) {
1682
1683 enode *res;
1684 int i;
1685
1686 if(size == 0) {
1687 fprintf(stderr, "Allocating enode of size 0 !\n" );
1688 return NULL;
1689 }
1690 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1691 res->type = type;
1692 res->size = size;
1693 res->pos = pos;
1694 for(i=0; i<size; i++) {
1695 value_init(res->arr[i].d);
1696 value_set_si(res->arr[i].d,0);
1697 res->arr[i].x.p = 0;
1698 }
1699 return res;
1700 } /* new_enode */
1701
1702 enode *ecopy(enode *e) {
1703
1704 enode *res;
1705 int i;
1706
1707 res = new_enode(e->type,e->size,e->pos);
1708 for(i=0;i<e->size;++i) {
1709 value_assign(res->arr[i].d,e->arr[i].d);
1710 if(value_zero_p(res->arr[i].d))
1711 res->arr[i].x.p = ecopy(e->arr[i].x.p);
1712 else if (EVALUE_IS_DOMAIN(res->arr[i]))
1713 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1714 else {
1715 value_init(res->arr[i].x.n);
1716 value_assign(res->arr[i].x.n,e->arr[i].x.n);
1717 }
1718 }
1719 return(res);
1720 } /* ecopy */
1721
1722 int ecmp(const evalue *e1, const evalue *e2)
1723 {
1724 enode *p1, *p2;
1725 int i;
1726 int r;
1727
1728 if (value_notzero_p(e1->d) && value_notzero_p(e2->d)) {
1729 Value m, m2;
1730 value_init(m);
1731 value_init(m2);
1732 value_multiply(m, e1->x.n, e2->d);
1733 value_multiply(m2, e2->x.n, e1->d);
1734
1735 if (value_lt(m, m2))
1736 r = -1;
1737 else if (value_gt(m, m2))
1738 r = 1;
1739 else
1740 r = 0;
1741
1742 value_clear(m);
1743 value_clear(m2);
1744
1745 return r;
1746 }
1747 if (value_notzero_p(e1->d))
1748 return -1;
1749 if (value_notzero_p(e2->d))
1750 return 1;
1751
1752 p1 = e1->x.p;
1753 p2 = e2->x.p;
1754
1755 if (p1->type != p2->type)
1756 return p1->type - p2->type;
1757 if (p1->pos != p2->pos)
1758 return p1->pos - p2->pos;
1759 if (p1->size != p2->size)
1760 return p1->size - p2->size;
1761
1762 for (i = p1->size-1; i >= 0; --i)
1763 if ((r = ecmp(&p1->arr[i], &p2->arr[i])) != 0)
1764 return r;
1765
1766 return 0;
1767 }
1768
1769 int eequal(evalue *e1,evalue *e2) {
1770
1771 int i;
1772 enode *p1, *p2;
1773
1774 if (value_ne(e1->d,e2->d))
1775 return 0;
1776
1777 /* e1->d == e2->d */
1778 if (value_notzero_p(e1->d)) {
1779 if (value_ne(e1->x.n,e2->x.n))
1780 return 0;
1781
1782 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1783 return 1;
1784 }
1785
1786 /* e1->d == e2->d == 0 */
1787 p1 = e1->x.p;
1788 p2 = e2->x.p;
1789 if (p1->type != p2->type) return 0;
1790 if (p1->size != p2->size) return 0;
1791 if (p1->pos != p2->pos) return 0;
1792 for (i=0; i<p1->size; i++)
1793 if (!eequal(&p1->arr[i], &p2->arr[i]) )
1794 return 0;
1795 return 1;
1796 } /* eequal */
1797
1798 void free_evalue_refs(evalue *e) {
1799
1800 enode *p;
1801 int i;
1802
1803 if (EVALUE_IS_DOMAIN(*e)) {
1804 Domain_Free(EVALUE_DOMAIN(*e));
1805 value_clear(e->d);
1806 return;
1807 } else if (value_pos_p(e->d)) {
1808
1809 /* 'e' stores a constant */
1810 value_clear(e->d);
1811 value_clear(e->x.n);
1812 return;
1813 }
1814 assert(value_zero_p(e->d));
1815 value_clear(e->d);
1816 p = e->x.p;
1817 if (!p) return; /* null pointer */
1818 for (i=0; i<p->size; i++) {
1819 free_evalue_refs(&(p->arr[i]));
1820 }
1821 free(p);
1822 return;
1823 } /* free_evalue_refs */
1824
1825 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
1826 Vector * val, evalue *res)
1827 {
1828 unsigned nparam = periods->Size;
1829
1830 if (p == nparam) {
1831 double d = compute_evalue(e, val->p);
1832 d *= VALUE_TO_DOUBLE(m);
1833 if (d > 0)
1834 d += .25;
1835 else
1836 d -= .25;
1837 value_assign(res->d, m);
1838 value_init(res->x.n);
1839 value_set_double(res->x.n, d);
1840 mpz_fdiv_r(res->x.n, res->x.n, m);
1841 return;
1842 }
1843 if (value_one_p(periods->p[p]))
1844 mod2table_r(e, periods, m, p+1, val, res);
1845 else {
1846 Value tmp;
1847 value_init(tmp);
1848
1849 value_assign(tmp, periods->p[p]);
1850 value_set_si(res->d, 0);
1851 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
1852 do {
1853 value_decrement(tmp, tmp);
1854 value_assign(val->p[p], tmp);
1855 mod2table_r(e, periods, m, p+1, val,
1856 &res->x.p->arr[VALUE_TO_INT(tmp)]);
1857 } while (value_pos_p(tmp));
1858
1859 value_clear(tmp);
1860 }
1861 }
1862
1863 static void rel2table(evalue *e, int zero)
1864 {
1865 if (value_pos_p(e->d)) {
1866 if (value_zero_p(e->x.n) == zero)
1867 value_set_si(e->x.n, 1);
1868 else
1869 value_set_si(e->x.n, 0);
1870 value_set_si(e->d, 1);
1871 } else {
1872 int i;
1873 for (i = 0; i < e->x.p->size; ++i)
1874 rel2table(&e->x.p->arr[i], zero);
1875 }
1876 }
1877
1878 void evalue_mod2table(evalue *e, int nparam)
1879 {
1880 enode *p;
1881 int i;
1882
1883 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
1884 return;
1885 p = e->x.p;
1886 for (i=0; i<p->size; i++) {
1887 evalue_mod2table(&(p->arr[i]), nparam);
1888 }
1889 if (p->type == relation) {
1890 evalue copy;
1891
1892 if (p->size > 2) {
1893 value_init(copy.d);
1894 evalue_copy(&copy, &p->arr[0]);
1895 }
1896 rel2table(&p->arr[0], 1);
1897 emul(&p->arr[0], &p->arr[1]);
1898 if (p->size > 2) {
1899 rel2table(&copy, 0);
1900 emul(&copy, &p->arr[2]);
1901 eadd(&p->arr[2], &p->arr[1]);
1902 free_evalue_refs(&p->arr[2]);
1903 free_evalue_refs(&copy);
1904 }
1905 free_evalue_refs(&p->arr[0]);
1906 value_clear(e->d);
1907 *e = p->arr[1];
1908 free(p);
1909 } else if (p->type == fractional) {
1910 Vector *periods = Vector_Alloc(nparam);
1911 Vector *val = Vector_Alloc(nparam);
1912 Value tmp;
1913 evalue *ev;
1914 evalue EP, res;
1915
1916 value_init(tmp);
1917 value_set_si(tmp, 1);
1918 Vector_Set(periods->p, 1, nparam);
1919 Vector_Set(val->p, 0, nparam);
1920 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
1921 enode *p = ev->x.p;
1922
1923 assert(p->type == polynomial);
1924 assert(p->size == 2);
1925 value_assign(periods->p[p->pos-1], p->arr[1].d);
1926 value_lcm(tmp, p->arr[1].d, &tmp);
1927 }
1928 value_lcm(tmp, ev->d, &tmp);
1929 value_init(EP.d);
1930 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
1931
1932 value_init(res.d);
1933 evalue_set_si(&res, 0, 1);
1934 /* Compute the polynomial using Horner's rule */
1935 for (i=p->size-1;i>1;i--) {
1936 eadd(&p->arr[i], &res);
1937 emul(&EP, &res);
1938 }
1939 eadd(&p->arr[1], &res);
1940
1941 free_evalue_refs(e);
1942 free_evalue_refs(&EP);
1943 *e = res;
1944
1945 value_clear(tmp);
1946 Vector_Free(val);
1947 Vector_Free(periods);
1948 }
1949 } /* evalue_mod2table */
1950
1951 /********************************************************/
1952 /* function in domain */
1953 /* check if the parameters in list_args */
1954 /* verifies the constraints of Domain P */
1955 /********************************************************/
1956 int in_domain(Polyhedron *P, Value *list_args) {
1957
1958 int col,row;
1959 Value v; /* value of the constraint of a row when
1960 parameters are instanciated*/
1961 Value tmp;
1962
1963 value_init(v);
1964 value_init(tmp);
1965
1966 /*P->Constraint constraint matrice of polyhedron P*/
1967 for(row=0;row<P->NbConstraints;row++) {
1968 value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
1969 for(col=1;col<P->Dimension+1;col++) {
1970 value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
1971 value_addto(v,v,tmp);
1972 }
1973 if (value_notzero_p(P->Constraint[row][0])) {
1974
1975 /*if v is not >=0 then this constraint is not respected */
1976 if (value_neg_p(v)) {
1977 next:
1978 value_clear(v);
1979 value_clear(tmp);
1980 return P->next ? in_domain(P->next, list_args) : 0;
1981 }
1982 }
1983 else {
1984
1985 /*if v is not = 0 then this constraint is not respected */
1986 if (value_notzero_p(v))
1987 goto next;
1988 }
1989 }
1990
1991 /*if not return before this point => all
1992 the constraints are respected */
1993 value_clear(v);
1994 value_clear(tmp);
1995 return 1;
1996 } /* in_domain */
1997
1998 /****************************************************/
1999 /* function compute enode */
2000 /* compute the value of enode p with parameters */
2001 /* list "list_args */
2002 /* compute the polynomial or the periodic */
2003 /****************************************************/
2004
2005 static double compute_enode(enode *p, Value *list_args) {
2006
2007 int i;
2008 Value m, param;
2009 double res=0.0;
2010
2011 if (!p)
2012 return(0.);
2013
2014 value_init(m);
2015 value_init(param);
2016
2017 if (p->type == polynomial) {
2018 if (p->size > 1)
2019 value_assign(param,list_args[p->pos-1]);
2020
2021 /* Compute the polynomial using Horner's rule */
2022 for (i=p->size-1;i>0;i--) {
2023 res +=compute_evalue(&p->arr[i],list_args);
2024 res *=VALUE_TO_DOUBLE(param);
2025 }
2026 res +=compute_evalue(&p->arr[0],list_args);
2027 }
2028 else if (p->type == fractional) {
2029 double d = compute_evalue(&p->arr[0], list_args);
2030 d -= floor(d+1e-10);
2031
2032 /* Compute the polynomial using Horner's rule */
2033 for (i=p->size-1;i>1;i--) {
2034 res +=compute_evalue(&p->arr[i],list_args);
2035 res *=d;
2036 }
2037 res +=compute_evalue(&p->arr[1],list_args);
2038 }
2039 else if (p->type == flooring) {
2040 double d = compute_evalue(&p->arr[0], list_args);
2041 d = floor(d+1e-10);
2042
2043 /* Compute the polynomial using Horner's rule */
2044 for (i=p->size-1;i>1;i--) {
2045 res +=compute_evalue(&p->arr[i],list_args);
2046 res *=d;
2047 }
2048 res +=compute_evalue(&p->arr[1],list_args);
2049 }
2050 else if (p->type == periodic) {
2051 value_assign(m,list_args[p->pos-1]);
2052
2053 /* Choose the right element of the periodic */
2054 value_set_si(param,p->size);
2055 value_pmodulus(m,m,param);
2056 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
2057 }
2058 else if (p->type == relation) {
2059 if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
2060 res = compute_evalue(&p->arr[1], list_args);
2061 else if (p->size > 2)
2062 res = compute_evalue(&p->arr[2], list_args);
2063 }
2064 else if (p->type == partition) {
2065 int dim = EVALUE_DOMAIN(p->arr[0])->Dimension;
2066 Value *vals = list_args;
2067 if (p->pos < dim) {
2068 NALLOC(vals, dim);
2069 for (i = 0; i < dim; ++i) {
2070 value_init(vals[i]);
2071 if (i < p->pos)
2072 value_assign(vals[i], list_args[i]);
2073 }
2074 }
2075 for (i = 0; i < p->size/2; ++i)
2076 if (DomainContains(EVALUE_DOMAIN(p->arr[2*i]), vals, p->pos, 0, 1)) {
2077 res = compute_evalue(&p->arr[2*i+1], vals);
2078 break;
2079 }
2080 if (p->pos < dim) {
2081 for (i = 0; i < dim; ++i)
2082 value_clear(vals[i]);
2083 free(vals);
2084 }
2085 }
2086 else
2087 assert(0);
2088 value_clear(m);
2089 value_clear(param);
2090 return res;
2091 } /* compute_enode */
2092
2093 /*************************************************/
2094 /* return the value of Ehrhart Polynomial */
2095 /* It returns a double, because since it is */
2096 /* a recursive function, some intermediate value */
2097 /* might not be integral */
2098 /*************************************************/
2099
2100 double compute_evalue(evalue *e,Value *list_args) {
2101
2102 double res;
2103
2104 if (value_notzero_p(e->d)) {
2105 if (value_notone_p(e->d))
2106 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
2107 else
2108 res = VALUE_TO_DOUBLE(e->x.n);
2109 }
2110 else
2111 res = compute_enode(e->x.p,list_args);
2112 return res;
2113 } /* compute_evalue */
2114
2115
2116 /****************************************************/
2117 /* function compute_poly : */
2118 /* Check for the good validity domain */
2119 /* return the number of point in the Polyhedron */
2120 /* in allocated memory */
2121 /* Using the Ehrhart pseudo-polynomial */
2122 /****************************************************/
2123 Value *compute_poly(Enumeration *en,Value *list_args) {
2124
2125 Value *tmp;
2126 /* double d; int i; */
2127
2128 tmp = (Value *) malloc (sizeof(Value));
2129 assert(tmp != NULL);
2130 value_init(*tmp);
2131 value_set_si(*tmp,0);
2132
2133 if(!en)
2134 return(tmp); /* no ehrhart polynomial */
2135 if(en->ValidityDomain) {
2136 if(!en->ValidityDomain->Dimension) { /* no parameters */
2137 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2138 return(tmp);
2139 }
2140 }
2141 else
2142 return(tmp); /* no Validity Domain */
2143 while(en) {
2144 if(in_domain(en->ValidityDomain,list_args)) {
2145
2146 #ifdef EVAL_EHRHART_DEBUG
2147 Print_Domain(stdout,en->ValidityDomain);
2148 print_evalue(stdout,&en->EP);
2149 #endif
2150
2151 /* d = compute_evalue(&en->EP,list_args);
2152 i = d;
2153 printf("(double)%lf = %d\n", d, i ); */
2154 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2155 return(tmp);
2156 }
2157 else
2158 en=en->next;
2159 }
2160 value_set_si(*tmp,0);
2161 return(tmp); /* no compatible domain with the arguments */
2162 } /* compute_poly */
2163
2164 size_t value_size(Value v) {
2165 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2166 * sizeof(v[0]._mp_d[0]);
2167 }
2168
2169 size_t domain_size(Polyhedron *D)
2170 {
2171 int i, j;
2172 size_t s = sizeof(*D);
2173
2174 for (i = 0; i < D->NbConstraints; ++i)
2175 for (j = 0; j < D->Dimension+2; ++j)
2176 s += value_size(D->Constraint[i][j]);
2177
2178 /*
2179 for (i = 0; i < D->NbRays; ++i)
2180 for (j = 0; j < D->Dimension+2; ++j)
2181 s += value_size(D->Ray[i][j]);
2182 */
2183
2184 return D->next ? s+domain_size(D->next) : s;
2185 }
2186
2187 size_t enode_size(enode *p) {
2188 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2189 int i;
2190
2191 if (p->type == partition)
2192 for (i = 0; i < p->size/2; ++i) {
2193 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2194 s += evalue_size(&p->arr[2*i+1]);
2195 }
2196 else
2197 for (i = 0; i < p->size; ++i) {
2198 s += evalue_size(&p->arr[i]);
2199 }
2200 return s;
2201 }
2202
2203 size_t evalue_size(evalue *e)
2204 {
2205 size_t s = sizeof(*e);
2206 s += value_size(e->d);
2207 if (value_notzero_p(e->d))
2208 s += value_size(e->x.n);
2209 else
2210 s += enode_size(e->x.p);
2211 return s;
2212 }
2213
2214 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2215 {
2216 evalue *found = NULL;
2217 evalue offset;
2218 evalue copy;
2219 int i;
2220
2221 if (value_pos_p(e->d) || e->x.p->type != fractional)
2222 return NULL;
2223
2224 value_init(offset.d);
2225 value_init(offset.x.n);
2226 poly_denom(&e->x.p->arr[0], &offset.d);
2227 value_lcm(m, offset.d, &offset.d);
2228 value_set_si(offset.x.n, 1);
2229
2230 value_init(copy.d);
2231 evalue_copy(&copy, cst);
2232
2233 eadd(&offset, cst);
2234 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2235
2236 if (eequal(base, &e->x.p->arr[0]))
2237 found = &e->x.p->arr[0];
2238 else {
2239 value_set_si(offset.x.n, -2);
2240
2241 eadd(&offset, cst);
2242 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2243
2244 if (eequal(base, &e->x.p->arr[0]))
2245 found = base;
2246 }
2247 free_evalue_refs(cst);
2248 free_evalue_refs(&offset);
2249 *cst = copy;
2250
2251 for (i = 1; !found && i < e->x.p->size; ++i)
2252 found = find_second(base, cst, &e->x.p->arr[i], m);
2253
2254 return found;
2255 }
2256
2257 static evalue *find_relation_pair(evalue *e)
2258 {
2259 int i;
2260 evalue *found = NULL;
2261
2262 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2263 return NULL;
2264
2265 if (e->x.p->type == fractional) {
2266 Value m;
2267 evalue *cst;
2268
2269 value_init(m);
2270 poly_denom(&e->x.p->arr[0], &m);
2271
2272 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2273 cst = &cst->x.p->arr[0])
2274 ;
2275
2276 for (i = 1; !found && i < e->x.p->size; ++i)
2277 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2278
2279 value_clear(m);
2280 }
2281
2282 i = e->x.p->type == relation;
2283 for (; !found && i < e->x.p->size; ++i)
2284 found = find_relation_pair(&e->x.p->arr[i]);
2285
2286 return found;
2287 }
2288
2289 void evalue_mod2relation(evalue *e) {
2290 evalue *d;
2291
2292 if (value_zero_p(e->d) && e->x.p->type == partition) {
2293 int i;
2294
2295 for (i = 0; i < e->x.p->size/2; ++i) {
2296 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2297 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2298 value_clear(e->x.p->arr[2*i].d);
2299 free_evalue_refs(&e->x.p->arr[2*i+1]);
2300 e->x.p->size -= 2;
2301 if (2*i < e->x.p->size) {
2302 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2303 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2304 }
2305 --i;
2306 }
2307 }
2308 if (e->x.p->size == 0) {
2309 free(e->x.p);
2310 evalue_set_si(e, 0, 1);
2311 }
2312
2313 return;
2314 }
2315
2316 while ((d = find_relation_pair(e)) != NULL) {
2317 evalue split;
2318 evalue *ev;
2319
2320 value_init(split.d);
2321 value_set_si(split.d, 0);
2322 split.x.p = new_enode(relation, 3, 0);
2323 evalue_set_si(&split.x.p->arr[1], 1, 1);
2324 evalue_set_si(&split.x.p->arr[2], 1, 1);
2325
2326 ev = &split.x.p->arr[0];
2327 value_set_si(ev->d, 0);
2328 ev->x.p = new_enode(fractional, 3, -1);
2329 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2330 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2331 evalue_copy(&ev->x.p->arr[0], d);
2332
2333 emul(&split, e);
2334
2335 reduce_evalue(e);
2336
2337 free_evalue_refs(&split);
2338 }
2339 }
2340
2341 static int evalue_comp(const void * a, const void * b)
2342 {
2343 const evalue *e1 = *(const evalue **)a;
2344 const evalue *e2 = *(const evalue **)b;
2345 return ecmp(e1, e2);
2346 }
2347
2348 void evalue_combine(evalue *e)
2349 {
2350 evalue **evs;
2351 int i, k;
2352 enode *p;
2353 evalue tmp;
2354
2355 if (value_notzero_p(e->d) || e->x.p->type != partition)
2356 return;
2357
2358 NALLOC(evs, e->x.p->size/2);
2359 for (i = 0; i < e->x.p->size/2; ++i)
2360 evs[i] = &e->x.p->arr[2*i+1];
2361 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2362 p = new_enode(partition, e->x.p->size, e->x.p->pos);
2363 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2364 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2365 value_clear(p->arr[2*k].d);
2366 value_clear(p->arr[2*k+1].d);
2367 p->arr[2*k] = *(evs[i]-1);
2368 p->arr[2*k+1] = *(evs[i]);
2369 ++k;
2370 } else {
2371 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2372 Polyhedron *L = D;
2373
2374 value_clear((evs[i]-1)->d);
2375
2376 while (L->next)
2377 L = L->next;
2378 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2379 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2380 free_evalue_refs(evs[i]);
2381 }
2382 }
2383
2384 for (i = 2*k ; i < p->size; ++i)
2385 value_clear(p->arr[i].d);
2386
2387 free(evs);
2388 free(e->x.p);
2389 p->size = 2*k;
2390 e->x.p = p;
2391
2392 for (i = 0; i < e->x.p->size/2; ++i) {
2393 Polyhedron *H;
2394 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2395 continue;
2396 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2397 if (H == NULL)
2398 continue;
2399 for (k = 0; k < e->x.p->size/2; ++k) {
2400 Polyhedron *D, *N, **P;
2401 if (i == k)
2402 continue;
2403 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2404 D = *P;
2405 if (D == NULL)
2406 continue;
2407 for (; D; D = N) {
2408 *P = D;
2409 N = D->next;
2410 if (D->NbEq <= H->NbEq) {
2411 P = &D->next;
2412 continue;
2413 }
2414
2415 value_init(tmp.d);
2416 tmp.x.p = new_enode(partition, 2, e->x.p->pos);
2417 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2418 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2419 reduce_evalue(&tmp);
2420 if (value_notzero_p(tmp.d) ||
2421 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2422 P = &D->next;
2423 else {
2424 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2425 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2426 *P = NULL;
2427 }
2428 free_evalue_refs(&tmp);
2429 }
2430 }
2431 Polyhedron_Free(H);
2432 }
2433
2434 for (i = 0; i < e->x.p->size/2; ++i) {
2435 Polyhedron *H, *E;
2436 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2437 if (!D) {
2438 value_clear(e->x.p->arr[2*i].d);
2439 free_evalue_refs(&e->x.p->arr[2*i+1]);
2440 e->x.p->size -= 2;
2441 if (2*i < e->x.p->size) {
2442 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2443 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2444 }
2445 --i;
2446 continue;
2447 }
2448 if (!D->next)
2449 continue;
2450 H = DomainConvex(D, 0);
2451 E = DomainDifference(H, D, 0);
2452 Domain_Free(D);
2453 D = DomainDifference(H, E, 0);
2454 Domain_Free(H);
2455 Domain_Free(E);
2456 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2457 }
2458 }
2459
2460 /* May change coefficients to become non-standard if fiddle is set
2461 * => reduce p afterwards to correct
2462 */
2463 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2464 Matrix **R, int fiddle)
2465 {
2466 Polyhedron *I, *H;
2467 evalue *pp;
2468 unsigned dim = D->Dimension;
2469 Matrix *T = Matrix_Alloc(2, dim+1);
2470 Value twice;
2471 value_init(twice);
2472 assert(T);
2473
2474 assert(p->type == fractional);
2475 pp = &p->arr[0];
2476 value_set_si(T->p[1][dim], 1);
2477 poly_denom(pp, d);
2478 while (value_zero_p(pp->d)) {
2479 assert(pp->x.p->type == polynomial);
2480 assert(pp->x.p->size == 2);
2481 assert(value_notzero_p(pp->x.p->arr[1].d));
2482 value_division(T->p[0][pp->x.p->pos-1], *d, pp->x.p->arr[1].d);
2483 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2484 if (fiddle && value_gt(twice, pp->x.p->arr[1].d))
2485 value_substract(pp->x.p->arr[1].x.n,
2486 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2487 value_multiply(T->p[0][pp->x.p->pos-1],
2488 T->p[0][pp->x.p->pos-1], pp->x.p->arr[1].x.n);
2489 pp = &pp->x.p->arr[0];
2490 }
2491 value_division(T->p[0][dim], *d, pp->d);
2492 value_multiply(T->p[0][dim], T->p[0][dim], pp->x.n);
2493 I = DomainImage(D, T, 0);
2494 H = DomainConvex(I, 0);
2495 Domain_Free(I);
2496 *R = T;
2497
2498 value_clear(twice);
2499
2500 return H;
2501 }
2502
2503 static int reduce_in_domain(evalue *e, Polyhedron *D)
2504 {
2505 int i;
2506 enode *p;
2507 Value d, min, max;
2508 int r = 0;
2509 Polyhedron *I;
2510 Matrix *T;
2511 int bounded;
2512
2513 if (value_notzero_p(e->d))
2514 return r;
2515
2516 p = e->x.p;
2517
2518 if (p->type == relation) {
2519 int equal;
2520 value_init(d);
2521 value_init(min);
2522 value_init(max);
2523
2524 I = polynomial_projection(p->arr[0].x.p, D, &d, &T, 1);
2525 bounded = line_minmax(I, &min, &max); /* frees I */
2526 equal = value_eq(min, max);
2527 mpz_cdiv_q(min, min, d);
2528 mpz_fdiv_q(max, max, d);
2529
2530 if (bounded && value_gt(min, max)) {
2531 /* Never zero */
2532 if (p->size == 3) {
2533 value_clear(e->d);
2534 *e = p->arr[2];
2535 } else {
2536 evalue_set_si(e, 0, 1);
2537 r = 1;
2538 }
2539 free_evalue_refs(&(p->arr[1]));
2540 free_evalue_refs(&(p->arr[0]));
2541 free(p);
2542 value_clear(d);
2543 value_clear(min);
2544 value_clear(max);
2545 Matrix_Free(T);
2546 return r ? r : reduce_in_domain(e, D);
2547 } else if (bounded && equal) {
2548 /* Always zero */
2549 if (p->size == 3)
2550 free_evalue_refs(&(p->arr[2]));
2551 value_clear(e->d);
2552 *e = p->arr[1];
2553 free_evalue_refs(&(p->arr[0]));
2554 free(p);
2555 value_clear(d);
2556 value_clear(min);
2557 value_clear(max);
2558 Matrix_Free(T);
2559 return reduce_in_domain(e, D);
2560 } else if (bounded && value_eq(min, max)) {
2561 /* zero for a single value */
2562 Polyhedron *E;
2563 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2564 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2565 value_multiply(min, min, d);
2566 value_substract(M->p[0][D->Dimension+1],
2567 M->p[0][D->Dimension+1], min);
2568 E = DomainAddConstraints(D, M, 0);
2569 value_clear(d);
2570 value_clear(min);
2571 value_clear(max);
2572 Matrix_Free(T);
2573 Matrix_Free(M);
2574 r = reduce_in_domain(&p->arr[1], E);
2575 if (p->size == 3)
2576 r |= reduce_in_domain(&p->arr[2], D);
2577 Domain_Free(E);
2578 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2579 return r;
2580 }
2581
2582 value_clear(d);
2583 value_clear(min);
2584 value_clear(max);
2585 Matrix_Free(T);
2586 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2587 }
2588
2589 i = p->type == relation ? 1 :
2590 p->type == fractional ? 1 : 0;
2591 for (; i<p->size; i++)
2592 r |= reduce_in_domain(&p->arr[i], D);
2593
2594 if (p->type != fractional) {
2595 if (r && p->type == polynomial) {
2596 evalue f;
2597 value_init(f.d);
2598 value_set_si(f.d, 0);
2599 f.x.p = new_enode(polynomial, 2, p->pos);
2600 evalue_set_si(&f.x.p->arr[0], 0, 1);
2601 evalue_set_si(&f.x.p->arr[1], 1, 1);
2602 reorder_terms(p, &f);
2603 value_clear(e->d);
2604 *e = p->arr[0];
2605 free(p);
2606 }
2607 return r;
2608 }
2609
2610 value_init(d);
2611 value_init(min);
2612 value_init(max);
2613 I = polynomial_projection(p, D, &d, &T, 1);
2614 Matrix_Free(T);
2615 bounded = line_minmax(I, &min, &max); /* frees I */
2616 mpz_fdiv_q(min, min, d);
2617 mpz_fdiv_q(max, max, d);
2618 value_substract(d, max, min);
2619
2620 if (bounded && value_eq(min, max)) {
2621 evalue inc;
2622 value_init(inc.d);
2623 value_init(inc.x.n);
2624 value_set_si(inc.d, 1);
2625 value_oppose(inc.x.n, min);
2626 eadd(&inc, &p->arr[0]);
2627 reorder_terms(p, &p->arr[0]); /* frees arr[0] */
2628 value_clear(e->d);
2629 *e = p->arr[1];
2630 free(p);
2631 free_evalue_refs(&inc);
2632 r = 1;
2633 } else if (bounded && value_one_p(d) && p->size > 3) {
2634 evalue rem;
2635 evalue inc;
2636 evalue t;
2637 evalue f;
2638 evalue factor;
2639 value_init(rem.d);
2640 value_set_si(rem.d, 0);
2641 rem.x.p = new_enode(fractional, 3, -1);
2642 evalue_copy(&rem.x.p->arr[0], &p->arr[0]);
2643 rem.x.p->arr[1] = p->arr[1];
2644 rem.x.p->arr[2] = p->arr[2];
2645 for (i = 3; i < p->size; ++i)
2646 p->arr[i-2] = p->arr[i];
2647 p->size -= 2;
2648
2649 value_init(inc.d);
2650 value_init(inc.x.n);
2651 value_set_si(inc.d, 1);
2652 value_oppose(inc.x.n, min);
2653
2654 value_init(t.d);
2655 evalue_copy(&t, &p->arr[0]);
2656 eadd(&inc, &t);
2657
2658 value_init(f.d);
2659 value_set_si(f.d, 0);
2660 f.x.p = new_enode(fractional, 3, -1);
2661 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2662 evalue_set_si(&f.x.p->arr[1], 1, 1);
2663 evalue_set_si(&f.x.p->arr[2], 2, 1);
2664
2665 value_init(factor.d);
2666 evalue_set_si(&factor, -1, 1);
2667 emul(&t, &factor);
2668
2669 eadd(&f, &factor);
2670 emul(&t, &factor);
2671
2672 evalue_set_si(&f.x.p->arr[1], 0, 1);
2673 evalue_set_si(&f.x.p->arr[2], -1, 1);
2674 eadd(&f, &factor);
2675
2676 emul(&factor, e);
2677 eadd(&rem, e);
2678
2679 free_evalue_refs(&inc);
2680 free_evalue_refs(&t);
2681 free_evalue_refs(&f);
2682 free_evalue_refs(&factor);
2683 free_evalue_refs(&rem);
2684
2685 reduce_in_domain(e, D);
2686
2687 r = 1;
2688 } else {
2689 _reduce_evalue(&p->arr[0], 0, 1);
2690 if (r == 1) {
2691 evalue f;
2692 value_init(f.d);
2693 value_set_si(f.d, 0);
2694 f.x.p = new_enode(fractional, 3, -1);
2695 value_clear(f.x.p->arr[0].d);
2696 f.x.p->arr[0] = p->arr[0];
2697 evalue_set_si(&f.x.p->arr[1], 0, 1);
2698 evalue_set_si(&f.x.p->arr[2], 1, 1);
2699 reorder_terms(p, &f);
2700 value_clear(e->d);
2701 *e = p->arr[1];
2702 free(p);
2703 }
2704 }
2705
2706 value_clear(d);
2707 value_clear(min);
2708 value_clear(max);
2709
2710 return r;
2711 }
2712
2713 void evalue_range_reduction(evalue *e)
2714 {
2715 int i;
2716 if (value_notzero_p(e->d) || e->x.p->type != partition)
2717 return;
2718
2719 for (i = 0; i < e->x.p->size/2; ++i)
2720 if (reduce_in_domain(&e->x.p->arr[2*i+1],
2721 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
2722 reduce_evalue(&e->x.p->arr[2*i+1]);
2723
2724 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2725 free_evalue_refs(&e->x.p->arr[2*i+1]);
2726 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
2727 value_clear(e->x.p->arr[2*i].d);
2728 e->x.p->size -= 2;
2729 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2730 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2731 --i;
2732 }
2733 }
2734 }
2735
2736 /* Frees EP
2737 */
2738 Enumeration* partition2enumeration(evalue *EP)
2739 {
2740 int i;
2741 Enumeration *en, *res = NULL;
2742
2743 if (EVALUE_IS_ZERO(*EP)) {
2744 free(EP);
2745 return res;
2746 }
2747
2748 for (i = 0; i < EP->x.p->size/2; ++i) {
2749 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[2*i])->Dimension);
2750 en = (Enumeration *)malloc(sizeof(Enumeration));
2751 en->next = res;
2752 res = en;
2753 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2754 value_clear(EP->x.p->arr[2*i].d);
2755 res->EP = EP->x.p->arr[2*i+1];
2756 }
2757 free(EP->x.p);
2758 value_clear(EP->d);
2759 free(EP);
2760 return res;
2761 }
2762
2763 static int frac2floor_in_domain(evalue *e, Polyhedron *D)
2764 {
2765 enode *p;
2766 int r = 0;
2767 int i;
2768 Polyhedron *I;
2769 Matrix *T;
2770 Value d, min;
2771 evalue fl;
2772
2773 if (value_notzero_p(e->d))
2774 return r;
2775
2776 p = e->x.p;
2777
2778 i = p->type == relation ? 1 :
2779 p->type == fractional ? 1 : 0;
2780 for (; i<p->size; i++)
2781 r |= frac2floor_in_domain(&p->arr[i], D);
2782
2783 if (p->type != fractional) {
2784 if (r && p->type == polynomial) {
2785 evalue f;
2786 value_init(f.d);
2787 value_set_si(f.d, 0);
2788 f.x.p = new_enode(polynomial, 2, p->pos);
2789 evalue_set_si(&f.x.p->arr[0], 0, 1);
2790 evalue_set_si(&f.x.p->arr[1], 1, 1);
2791 reorder_terms(p, &f);
2792 value_clear(e->d);
2793 *e = p->arr[0];
2794 free(p);
2795 }
2796 return r;
2797 }
2798
2799 value_init(d);
2800 I = polynomial_projection(p, D, &d, &T, 0);
2801
2802 /*
2803 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2804 */
2805
2806 assert(I->NbEq == 0); /* Should have been reduced */
2807
2808 /* Find minimum */
2809 for (i = 0; i < I->NbConstraints; ++i)
2810 if (value_pos_p(I->Constraint[i][1]))
2811 break;
2812
2813 assert(i < I->NbConstraints);
2814 value_init(min);
2815 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
2816 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
2817 if (value_neg_p(min)) {
2818 evalue offset;
2819 mpz_fdiv_q(min, min, d);
2820 value_init(offset.d);
2821 value_set_si(offset.d, 1);
2822 value_init(offset.x.n);
2823 value_oppose(offset.x.n, min);
2824 eadd(&offset, &p->arr[0]);
2825 free_evalue_refs(&offset);
2826 }
2827
2828 Polyhedron_Free(I);
2829 Matrix_Free(T);
2830 value_clear(min);
2831 value_clear(d);
2832
2833 value_init(fl.d);
2834 value_set_si(fl.d, 0);
2835 fl.x.p = new_enode(flooring, 3, -1);
2836 evalue_set_si(&fl.x.p->arr[1], 0, 1);
2837 evalue_set_si(&fl.x.p->arr[2], -1, 1);
2838 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
2839
2840 eadd(&fl, &p->arr[0]);
2841 reorder_terms(p, &p->arr[0]);
2842 *e = p->arr[1];
2843 free(p);
2844 free_evalue_refs(&fl);
2845
2846 return 1;
2847 }
2848
2849 void evalue_frac2floor(evalue *e)
2850 {
2851 int i;
2852 if (value_notzero_p(e->d) || e->x.p->type != partition)
2853 return;
2854
2855 for (i = 0; i < e->x.p->size/2; ++i)
2856 if (frac2floor_in_domain(&e->x.p->arr[2*i+1],
2857 EVALUE_DOMAIN(e->x.p->arr[2*i])))
2858 reduce_evalue(&e->x.p->arr[2*i+1]);
2859 }
2860
2861 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
2862 Vector *row)
2863 {
2864 int nr, nc;
2865 int i;
2866 int nparam = D->Dimension - nvar;
2867
2868 if (C == 0) {
2869 nr = D->NbConstraints + 2;
2870 nc = D->Dimension + 2 + 1;
2871 C = Matrix_Alloc(nr, nc);
2872 for (i = 0; i < D->NbConstraints; ++i) {
2873 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
2874 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2875 D->Dimension + 1 - nvar);
2876 }
2877 } else {
2878 Matrix *oldC = C;
2879 nr = C->NbRows + 2;
2880 nc = C->NbColumns + 1;
2881 C = Matrix_Alloc(nr, nc);
2882 for (i = 0; i < oldC->NbRows; ++i) {
2883 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
2884 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2885 oldC->NbColumns - 1 - nvar);
2886 }
2887 }
2888 value_set_si(C->p[nr-2][0], 1);
2889 value_set_si(C->p[nr-2][1 + nvar], 1);
2890 value_set_si(C->p[nr-2][nc - 1], -1);
2891
2892 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
2893 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
2894 1 + nparam);
2895
2896 return C;
2897 }
2898
2899 static void floor2frac_r(evalue *e, int nvar)
2900 {
2901 enode *p;
2902 int i;
2903 evalue f;
2904 evalue *pp;
2905
2906 if (value_notzero_p(e->d))
2907 return;
2908
2909 p = e->x.p;
2910
2911 assert(p->type == flooring);
2912 for (i = 1; i < p->size; i++)
2913 floor2frac_r(&p->arr[i], nvar);
2914
2915 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
2916 assert(pp->x.p->type == polynomial);
2917 pp->x.p->pos -= nvar;
2918 }
2919
2920 value_init(f.d);
2921 value_set_si(f.d, 0);
2922 f.x.p = new_enode(fractional, 3, -1);
2923 evalue_set_si(&f.x.p->arr[1], 0, 1);
2924 evalue_set_si(&f.x.p->arr[2], -1, 1);
2925 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2926
2927 eadd(&f, &p->arr[0]);
2928 reorder_terms(p, &p->arr[0]);
2929 *e = p->arr[1];
2930 free(p);
2931 free_evalue_refs(&f);
2932 }
2933
2934 /* Convert flooring back to fractional and shift position
2935 * of the parameters by nvar
2936 */
2937 static void floor2frac(evalue *e, int nvar)
2938 {
2939 floor2frac_r(e, nvar);
2940 reduce_evalue(e);
2941 }
2942
2943 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
2944 {
2945 evalue *t;
2946 int nparam = D->Dimension - nvar;
2947
2948 if (C != 0) {
2949 C = Matrix_Copy(C);
2950 D = Constraints2Polyhedron(C, 0);
2951 Matrix_Free(C);
2952 }
2953
2954 t = barvinok_enumerate_e(D, 0, nparam, 0);
2955
2956 /* Double check that D was not unbounded. */
2957 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
2958
2959 if (C != 0)
2960 Polyhedron_Free(D);
2961
2962 return t;
2963 }
2964
2965 evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
2966 Matrix *C)
2967 {
2968 Vector *row = NULL;
2969 int i;
2970 evalue *res;
2971 Matrix *origC;
2972 evalue *factor = NULL;
2973 evalue cum;
2974
2975 if (EVALUE_IS_ZERO(*e))
2976 return 0;
2977
2978 if (D->next) {
2979 Polyhedron *DD = Disjoint_Domain(D, 0, 0);
2980 Polyhedron *Q;
2981
2982 Q = DD;
2983 DD = Q->next;
2984 Q->next = 0;
2985
2986 res = esum_over_domain(e, nvar, Q, C);
2987 Polyhedron_Free(Q);
2988
2989 for (Q = DD; Q; Q = DD) {
2990 evalue *t;
2991
2992 DD = Q->next;
2993 Q->next = 0;
2994
2995 t = esum_over_domain(e, nvar, Q, C);
2996 Polyhedron_Free(Q);
2997
2998 if (!res)
2999 res = t;
3000 else if (t) {
3001 eadd(t, res);
3002 free_evalue_refs(t);
3003 free(t);
3004 }
3005 }
3006 return res;
3007 }
3008
3009 if (value_notzero_p(e->d)) {
3010 evalue *t;
3011
3012 t = esum_over_domain_cst(nvar, D, C);
3013
3014 if (!EVALUE_IS_ONE(*e))
3015 emul(e, t);
3016
3017 return t;
3018 }
3019
3020 switch (e->x.p->type) {
3021 case flooring: {
3022 evalue *pp = &e->x.p->arr[0];
3023
3024 if (pp->x.p->pos > nvar) {
3025 /* remainder is independent of the summated vars */
3026 evalue f;
3027 evalue *t;
3028
3029 value_init(f.d);
3030 evalue_copy(&f, e);
3031 floor2frac(&f, nvar);
3032
3033 t = esum_over_domain_cst(nvar, D, C);
3034
3035 emul(&f, t);
3036
3037 free_evalue_refs(&f);
3038
3039 return t;
3040 }
3041
3042 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3043 poly_denom(pp, &row->p[1 + nvar]);
3044 value_set_si(row->p[0], 1);
3045 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3046 pp = &pp->x.p->arr[0]) {
3047 int pos;
3048 assert(pp->x.p->type == polynomial);
3049 pos = pp->x.p->pos;
3050 if (pos >= 1 + nvar)
3051 ++pos;
3052 value_assign(row->p[pos], row->p[1+nvar]);
3053 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3054 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3055 }
3056 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3057 value_division(row->p[1 + D->Dimension + 1],
3058 row->p[1 + D->Dimension + 1],
3059 pp->d);
3060 value_multiply(row->p[1 + D->Dimension + 1],
3061 row->p[1 + D->Dimension + 1],
3062 pp->x.n);
3063 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3064 break;
3065 }
3066 case polynomial: {
3067 int pos = e->x.p->pos;
3068
3069 if (pos > nvar) {
3070 ALLOC(factor);
3071 value_init(factor->d);
3072 value_set_si(factor->d, 0);
3073 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3074 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3075 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3076 break;
3077 }
3078
3079 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3080 for (i = 0; i < D->NbRays; ++i)
3081 if (value_notzero_p(D->Ray[i][pos]))
3082 break;
3083 assert(i < D->NbRays);
3084 if (value_neg_p(D->Ray[i][pos])) {
3085 ALLOC(factor);
3086 value_init(factor->d);
3087 evalue_set_si(factor, -1, 1);
3088 }
3089 value_set_si(row->p[0], 1);
3090 value_set_si(row->p[pos], 1);
3091 value_set_si(row->p[1 + nvar], -1);
3092 break;
3093 }
3094 default:
3095 assert(0);
3096 }
3097
3098 i = type_offset(e->x.p);
3099
3100 res = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3101 ++i;
3102
3103 if (factor) {
3104 value_init(cum.d);
3105 evalue_copy(&cum, factor);
3106 }
3107
3108 origC = C;
3109 for (; i < e->x.p->size; ++i) {
3110 evalue *t;
3111 if (row) {
3112 Matrix *prevC = C;
3113 C = esum_add_constraint(nvar, D, C, row);
3114 if (prevC != origC)
3115 Matrix_Free(prevC);
3116 }
3117 /*
3118 if (row)
3119 Vector_Print(stderr, P_VALUE_FMT, row);
3120 if (C)
3121 Matrix_Print(stderr, P_VALUE_FMT, C);
3122 */
3123 t = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3124
3125 if (t && factor)
3126 emul(&cum, t);
3127
3128 if (!res)
3129 res = t;
3130 else if (t) {
3131 eadd(t, res);
3132 free_evalue_refs(t);
3133 free(t);
3134 }
3135 if (factor && i+1 < e->x.p->size)
3136 emul(factor, &cum);
3137 }
3138 if (C != origC)
3139 Matrix_Free(C);
3140
3141 if (factor) {
3142 free_evalue_refs(factor);
3143 free_evalue_refs(&cum);
3144 free(factor);
3145 }
3146
3147 if (row)
3148 Vector_Free(row);
3149
3150 reduce_evalue(res);
3151
3152 return res;
3153 }
3154
3155 evalue *esum(evalue *e, int nvar)
3156 {
3157 int i;
3158 evalue *res;
3159 ALLOC(res);
3160 value_init(res->d);
3161
3162 assert(nvar >= 0);
3163 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3164 evalue_copy(res, e);
3165 return res;
3166 }
3167
3168 evalue_set_si(res, 0, 1);
3169
3170 assert(value_zero_p(e->d));
3171 assert(e->x.p->type == partition);
3172
3173 for (i = 0; i < e->x.p->size/2; ++i) {
3174 evalue *t;
3175 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3176 EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
3177 eadd(t, res);
3178 free_evalue_refs(t);
3179 free(t);
3180 }
3181
3182 reduce_evalue(res);
3183
3184 return res;
3185 }
3186
3187 /* Initial silly implementation */
3188 void eor(evalue *e1, evalue *res)
3189 {
3190 evalue E;
3191 evalue mone;
3192 value_init(E.d);
3193 value_init(mone.d);
3194 evalue_set_si(&mone, -1, 1);
3195
3196 evalue_copy(&E, res);
3197 eadd(e1, &E);
3198 emul(e1, res);
3199 emul(&mone, res);
3200 eadd(&E, res);
3201
3202 free_evalue_refs(&E);
3203 free_evalue_refs(&mone);
3204 }
lattice point enumeration library