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