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