?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
store "real" dimension in pos field of partition
[barvinok.git] / ev_operations.c
blobf2d4067875a03938b4bf526ede407760d9924b88
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 for (i=0; i<p->size/2; i++) {
852 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), pname);
853 print_evalue(DST, &p->arr[2*i+1], pname);
854 }
855 break;
856 default:
857 assert(0);
858 }
859 return;
860 } /* print_enode */
861
862 static void eadd_rev(evalue *e1, evalue *res)
863 {
864 evalue ev;
865 value_init(ev.d);
866 evalue_copy(&ev, e1);
867 eadd(res, &ev);
868 free_evalue_refs(res);
869 *res = ev;
870 }
871
872 static void eadd_rev_cst (evalue *e1, evalue *res)
873 {
874 evalue ev;
875 value_init(ev.d);
876 evalue_copy(&ev, e1);
877 eadd(res, &ev.x.p->arr[type_offset(ev.x.p)]);
878 free_evalue_refs(res);
879 *res = ev;
880 }
881
882 struct section { Polyhedron * D; evalue E; };
883
884 void eadd_partitions (evalue *e1,evalue *res)
885 {
886 int n, i, j;
887 Polyhedron *d, *fd;
888 struct section *s;
889 s = (struct section *)
890 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
891 sizeof(struct section));
892 assert(s);
893 assert(e1->x.p->pos == res->x.p->pos);
894 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
895
896 n = 0;
897 for (j = 0; j < e1->x.p->size/2; ++j) {
898 assert(res->x.p->size >= 2);
899 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
900 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
901 if (!emptyQ(fd))
902 for (i = 1; i < res->x.p->size/2; ++i) {
903 Polyhedron *t = fd;
904 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
905 Domain_Free(t);
906 if (emptyQ(fd))
907 break;
908 }
909 if (emptyQ(fd)) {
910 Domain_Free(fd);
911 continue;
912 }
913 value_init(s[n].E.d);
914 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
915 s[n].D = fd;
916 ++n;
917 }
918 for (i = 0; i < res->x.p->size/2; ++i) {
919 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
920 for (j = 0; j < e1->x.p->size/2; ++j) {
921 Polyhedron *t;
922 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
923 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
924 if (emptyQ(d)) {
925 Domain_Free(d);
926 continue;
927 }
928 t = fd;
929 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
930 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
931 Domain_Free(t);
932 value_init(s[n].E.d);
933 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
934 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
935 s[n].D = d;
936 ++n;
937 }
938 if (!emptyQ(fd)) {
939 s[n].E = res->x.p->arr[2*i+1];
940 s[n].D = fd;
941 ++n;
942 } else {
943 free_evalue_refs(&res->x.p->arr[2*i+1]);
944 Domain_Free(fd);
945 }
946 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
947 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
948 value_clear(res->x.p->arr[2*i].d);
949 }
950
951 free(res->x.p);
952 assert(n > 0);
953 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
954 for (j = 0; j < n; ++j) {
955 s[j].D = DomainConstraintSimplify(s[j].D, 0);
956 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
957 value_clear(res->x.p->arr[2*j+1].d);
958 res->x.p->arr[2*j+1] = s[j].E;
959 }
960
961 free(s);
962 }
963
964 static void explicit_complement(evalue *res)
965 {
966 enode *rel = new_enode(relation, 3, 0);
967 assert(rel);
968 value_clear(rel->arr[0].d);
969 rel->arr[0] = res->x.p->arr[0];
970 value_clear(rel->arr[1].d);
971 rel->arr[1] = res->x.p->arr[1];
972 value_set_si(rel->arr[2].d, 1);
973 value_init(rel->arr[2].x.n);
974 value_set_si(rel->arr[2].x.n, 0);
975 free(res->x.p);
976 res->x.p = rel;
977 }
978
979 void eadd(evalue *e1,evalue *res) {
980
981 int i;
982 if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
983 /* Add two rational numbers */
984 Value g,m1,m2;
985 value_init(g);
986 value_init(m1);
987 value_init(m2);
988
989 value_multiply(m1,e1->x.n,res->d);
990 value_multiply(m2,res->x.n,e1->d);
991 value_addto(res->x.n,m1,m2);
992 value_multiply(res->d,e1->d,res->d);
993 Gcd(res->x.n,res->d,&g);
994 if (value_notone_p(g)) {
995 value_division(res->d,res->d,g);
996 value_division(res->x.n,res->x.n,g);
997 }
998 value_clear(g); value_clear(m1); value_clear(m2);
999 return ;
1000 }
1001 else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
1002 switch (res->x.p->type) {
1003 case polynomial:
1004 /* Add the constant to the constant term of a polynomial*/
1005 eadd(e1, &res->x.p->arr[0]);
1006 return ;
1007 case periodic:
1008 /* Add the constant to all elements of a periodic number */
1009 for (i=0; i<res->x.p->size; i++) {
1010 eadd(e1, &res->x.p->arr[i]);
1011 }
1012 return ;
1013 case evector:
1014 fprintf(stderr, "eadd: cannot add const with vector\n");
1015 return;
1016 case flooring:
1017 case fractional:
1018 eadd(e1, &res->x.p->arr[1]);
1019 return ;
1020 case partition:
1021 assert(EVALUE_IS_ZERO(*e1));
1022 break; /* Do nothing */
1023 case relation:
1024 /* Create (zero) complement if needed */
1025 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
1026 explicit_complement(res);
1027 for (i = 1; i < res->x.p->size; ++i)
1028 eadd(e1, &res->x.p->arr[i]);
1029 break;
1030 default:
1031 assert(0);
1032 }
1033 }
1034 /* add polynomial or periodic to constant
1035 * you have to exchange e1 and res, before doing addition */
1036
1037 else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
1038 eadd_rev(e1, res);
1039 return;
1040 }
1041 else { // ((e1->d==0) && (res->d==0))
1042 assert(!((e1->x.p->type == partition) ^
1043 (res->x.p->type == partition)));
1044 if (e1->x.p->type == partition) {
1045 eadd_partitions(e1, res);
1046 return;
1047 }
1048 if (e1->x.p->type == relation &&
1049 (res->x.p->type != relation ||
1050 mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1051 eadd_rev(e1, res);
1052 return;
1053 }
1054 if (res->x.p->type == relation) {
1055 if (e1->x.p->type == relation &&
1056 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1057 if (res->x.p->size < 3 && e1->x.p->size == 3)
1058 explicit_complement(res);
1059 if (e1->x.p->size < 3 && res->x.p->size == 3)
1060 explicit_complement(e1);
1061 for (i = 1; i < res->x.p->size; ++i)
1062 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1063 return;
1064 }
1065 if (res->x.p->size < 3)
1066 explicit_complement(res);
1067 for (i = 1; i < res->x.p->size; ++i)
1068 eadd(e1, &res->x.p->arr[i]);
1069 return;
1070 }
1071 if ((e1->x.p->type != res->x.p->type) ) {
1072 /* adding to evalues of different type. two cases are possible
1073 * res is periodic and e1 is polynomial, you have to exchange
1074 * e1 and res then to add e1 to the constant term of res */
1075 if (e1->x.p->type == polynomial) {
1076 eadd_rev_cst(e1, res);
1077 }
1078 else if (res->x.p->type == polynomial) {
1079 /* res is polynomial and e1 is periodic,
1080 add e1 to the constant term of res */
1081
1082 eadd(e1,&res->x.p->arr[0]);
1083 } else
1084 assert(0);
1085
1086 return;
1087 }
1088 else if (e1->x.p->pos != res->x.p->pos ||
1089 ((res->x.p->type == fractional ||
1090 res->x.p->type == flooring) &&
1091 !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1092 /* adding evalues of different position (i.e function of different unknowns
1093 * to case are possible */
1094
1095 switch (res->x.p->type) {
1096 case flooring:
1097 case fractional:
1098 if(mod_term_smaller(res, e1))
1099 eadd(e1,&res->x.p->arr[1]);
1100 else
1101 eadd_rev_cst(e1, res);
1102 return;
1103 case polynomial: // res and e1 are polynomials
1104 // add e1 to the constant term of res
1105
1106 if(res->x.p->pos < e1->x.p->pos)
1107 eadd(e1,&res->x.p->arr[0]);
1108 else
1109 eadd_rev_cst(e1, res);
1110 // value_clear(g); value_clear(m1); value_clear(m2);
1111 return;
1112 case periodic: // res and e1 are pointers to periodic numbers
1113 //add e1 to all elements of res
1114
1115 if(res->x.p->pos < e1->x.p->pos)
1116 for (i=0;i<res->x.p->size;i++) {
1117 eadd(e1,&res->x.p->arr[i]);
1118 }
1119 else
1120 eadd_rev(e1, res);
1121 return;
1122 default:
1123 assert(0);
1124 }
1125 }
1126
1127
1128 //same type , same pos and same size
1129 if (e1->x.p->size == res->x.p->size) {
1130 // add any element in e1 to the corresponding element in res
1131 i = type_offset(res->x.p);
1132 if (i == 1)
1133 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1134 for (; i<res->x.p->size; i++) {
1135 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1136 }
1137 return ;
1138 }
1139
1140 /* Sizes are different */
1141 switch(res->x.p->type) {
1142 case polynomial:
1143 case flooring:
1144 case fractional:
1145 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1146 /* new enode and add res to that new node. If you do not do */
1147 /* that, you lose the the upper weight part of e1 ! */
1148
1149 if(e1->x.p->size > res->x.p->size)
1150 eadd_rev(e1, res);
1151 else {
1152 i = type_offset(res->x.p);
1153 if (i == 1)
1154 assert(eequal(&e1->x.p->arr[0],
1155 &res->x.p->arr[0]));
1156 for (; i<e1->x.p->size ; i++) {
1157 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1158 }
1159
1160 return ;
1161 }
1162 break;
1163
1164 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1165 case periodic:
1166 {
1167 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1168 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1169 to add periodicaly elements of e1 to elements of ne, and finaly to
1170 return ne. */
1171 int x,y,p;
1172 Value ex, ey ,ep;
1173 evalue *ne;
1174 value_init(ex); value_init(ey);value_init(ep);
1175 x=e1->x.p->size;
1176 y= res->x.p->size;
1177 value_set_si(ex,e1->x.p->size);
1178 value_set_si(ey,res->x.p->size);
1179 value_assign (ep,*Lcm(ex,ey));
1180 p=(int)mpz_get_si(ep);
1181 ne= (evalue *) malloc (sizeof(evalue));
1182 value_init(ne->d);
1183 value_set_si( ne->d,0);
1184
1185 ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
1186 for(i=0;i<p;i++) {
1187 evalue_copy(&ne->x.p->arr[i], &res->x.p->arr[i%y]);
1188 }
1189 for(i=0;i<p;i++) {
1190 eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
1191 }
1192
1193 value_assign(res->d, ne->d);
1194 res->x.p=ne->x.p;
1195
1196 return ;
1197 }
1198 case evector:
1199 fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
1200 return ;
1201 default:
1202 assert(0);
1203 }
1204 }
1205 return ;
1206 }/* eadd */
1207
1208 static void emul_rev (evalue *e1, evalue *res)
1209 {
1210 evalue ev;
1211 value_init(ev.d);
1212 evalue_copy(&ev, e1);
1213 emul(res, &ev);
1214 free_evalue_refs(res);
1215 *res = ev;
1216 }
1217
1218 static void emul_poly (evalue *e1, evalue *res)
1219 {
1220 int i, j, o = type_offset(res->x.p);
1221 evalue tmp;
1222 int size=(e1->x.p->size + res->x.p->size - o - 1);
1223 value_init(tmp.d);
1224 value_set_si(tmp.d,0);
1225 tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
1226 if (o)
1227 evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
1228 for (i=o; i < e1->x.p->size; i++) {
1229 evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
1230 emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
1231 }
1232 for (; i<size; i++)
1233 evalue_set_si(&tmp.x.p->arr[i], 0, 1);
1234 for (i=o+1; i<res->x.p->size; i++)
1235 for (j=o; j<e1->x.p->size; j++) {
1236 evalue ev;
1237 value_init(ev.d);
1238 evalue_copy(&ev, &e1->x.p->arr[j]);
1239 emul(&res->x.p->arr[i], &ev);
1240 eadd(&ev, &tmp.x.p->arr[i+j-o]);
1241 free_evalue_refs(&ev);
1242 }
1243 free_evalue_refs(res);
1244 *res = tmp;
1245 }
1246
1247 void emul_partitions (evalue *e1,evalue *res)
1248 {
1249 int n, i, j, k;
1250 Polyhedron *d;
1251 struct section *s;
1252 s = (struct section *)
1253 malloc((e1->x.p->size/2) * (res->x.p->size/2) *
1254 sizeof(struct section));
1255 assert(s);
1256 fprintf(stderr, "%d %d\n", e1->x.p->pos, res->x.p->pos);
1257 assert(e1->x.p->pos == res->x.p->pos);
1258 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
1259
1260 n = 0;
1261 for (i = 0; i < res->x.p->size/2; ++i) {
1262 for (j = 0; j < e1->x.p->size/2; ++j) {
1263 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1264 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1265 if (emptyQ(d)) {
1266 Domain_Free(d);
1267 continue;
1268 }
1269
1270 /* This code is only needed because the partitions
1271 are not true partitions.
1272 */
1273 for (k = 0; k < n; ++k) {
1274 if (DomainIncludes(s[k].D, d))
1275 break;
1276 if (DomainIncludes(d, s[k].D)) {
1277 Domain_Free(s[k].D);
1278 free_evalue_refs(&s[k].E);
1279 if (n > k)
1280 s[k] = s[--n];
1281 --k;
1282 }
1283 }
1284 if (k < n) {
1285 Domain_Free(d);
1286 continue;
1287 }
1288
1289 value_init(s[n].E.d);
1290 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1291 emul(&e1->x.p->arr[2*j+1], &s[n].E);
1292 s[n].D = d;
1293 ++n;
1294 }
1295 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1296 value_clear(res->x.p->arr[2*i].d);
1297 free_evalue_refs(&res->x.p->arr[2*i+1]);
1298 }
1299
1300 free(res->x.p);
1301 if (n == 0)
1302 evalue_set_si(res, 0, 1);
1303 else {
1304 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1305 for (j = 0; j < n; ++j) {
1306 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1307 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1308 value_clear(res->x.p->arr[2*j+1].d);
1309 res->x.p->arr[2*j+1] = s[j].E;
1310 }
1311 }
1312
1313 free(s);
1314 }
1315
1316 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1317
1318 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1319 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1320 * evalues is not treated here */
1321
1322 void emul (evalue *e1, evalue *res ){
1323 int i,j;
1324
1325 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
1326 fprintf(stderr, "emul: do not proced on evector type !\n");
1327 return;
1328 }
1329
1330 if (EVALUE_IS_ZERO(*res))
1331 return;
1332
1333 if (value_zero_p(e1->d) && e1->x.p->type == partition) {
1334 if (value_zero_p(res->d) && res->x.p->type == partition)
1335 emul_partitions(e1, res);
1336 else
1337 emul_rev(e1, res);
1338 } else if (value_zero_p(res->d) && res->x.p->type == partition) {
1339 for (i = 0; i < res->x.p->size/2; ++i)
1340 emul(e1, &res->x.p->arr[2*i+1]);
1341 } else
1342 if (value_zero_p(res->d) && res->x.p->type == relation) {
1343 if (value_zero_p(e1->d) && e1->x.p->type == relation &&
1344 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1345 if (res->x.p->size < 3 && e1->x.p->size == 3)
1346 explicit_complement(res);
1347 if (e1->x.p->size < 3 && res->x.p->size == 3)
1348 explicit_complement(e1);
1349 for (i = 1; i < res->x.p->size; ++i)
1350 emul(&e1->x.p->arr[i], &res->x.p->arr[i]);
1351 return;
1352 }
1353 for (i = 1; i < res->x.p->size; ++i)
1354 emul(e1, &res->x.p->arr[i]);
1355 } else
1356 if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
1357 switch(e1->x.p->type) {
1358 case polynomial:
1359 switch(res->x.p->type) {
1360 case polynomial:
1361 if(e1->x.p->pos == res->x.p->pos) {
1362 /* Product of two polynomials of the same variable */
1363 emul_poly(e1, res);
1364 return;
1365 }
1366 else {
1367 /* Product of two polynomials of different variables */
1368
1369 if(res->x.p->pos < e1->x.p->pos)
1370 for( i=0; i<res->x.p->size ; i++)
1371 emul(e1, &res->x.p->arr[i]);
1372 else
1373 emul_rev(e1, res);
1374
1375 return ;
1376 }
1377 case periodic:
1378 case flooring:
1379 case fractional:
1380 /* Product of a polynomial and a periodic or fractional */
1381 emul_rev(e1, res);
1382 return;
1383 default:
1384 assert(0);
1385 }
1386 case periodic:
1387 switch(res->x.p->type) {
1388 case periodic:
1389 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
1390 /* Product of two periodics of the same parameter and period */
1391
1392 for(i=0; i<res->x.p->size;i++)
1393 emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
1394
1395 return;
1396 }
1397 else{
1398 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
1399 /* Product of two periodics of the same parameter and different periods */
1400 evalue *newp;
1401 Value x,y,z;
1402 int ix,iy,lcm;
1403 value_init(x); value_init(y);value_init(z);
1404 ix=e1->x.p->size;
1405 iy=res->x.p->size;
1406 value_set_si(x,e1->x.p->size);
1407 value_set_si(y,res->x.p->size);
1408 value_assign (z,*Lcm(x,y));
1409 lcm=(int)mpz_get_si(z);
1410 newp= (evalue *) malloc (sizeof(evalue));
1411 value_init(newp->d);
1412 value_set_si( newp->d,0);
1413 newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
1414 for(i=0;i<lcm;i++) {
1415 evalue_copy(&newp->x.p->arr[i],
1416 &res->x.p->arr[i%iy]);
1417 }
1418 for(i=0;i<lcm;i++)
1419 emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
1420
1421 value_assign(res->d,newp->d);
1422 res->x.p=newp->x.p;
1423
1424 value_clear(x); value_clear(y);value_clear(z);
1425 return ;
1426 }
1427 else {
1428 /* Product of two periodics of different parameters */
1429
1430 if(res->x.p->pos < e1->x.p->pos)
1431 for(i=0; i<res->x.p->size; i++)
1432 emul(e1, &(res->x.p->arr[i]));
1433 else
1434 emul_rev(e1, res);
1435
1436 return;
1437 }
1438 }
1439 case polynomial:
1440 /* Product of a periodic and a polynomial */
1441
1442 for(i=0; i<res->x.p->size ; i++)
1443 emul(e1, &(res->x.p->arr[i]));
1444
1445 return;
1446
1447 }
1448 case flooring:
1449 case fractional:
1450 switch(res->x.p->type) {
1451 case polynomial:
1452 for(i=0; i<res->x.p->size ; i++)
1453 emul(e1, &(res->x.p->arr[i]));
1454 return;
1455 default:
1456 case periodic:
1457 assert(0);
1458 case flooring:
1459 case fractional:
1460 assert(e1->x.p->type == res->x.p->type);
1461 if (e1->x.p->pos == res->x.p->pos &&
1462 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1463 evalue d;
1464 value_init(d.d);
1465 poly_denom(&e1->x.p->arr[0], &d.d);
1466 if (e1->x.p->type != fractional || !value_two_p(d.d))
1467 emul_poly(e1, res);
1468 else {
1469 value_init(d.x.n);
1470 value_set_si(d.x.n, 1);
1471 evalue tmp;
1472 /* { x }^2 == { x }/2 */
1473 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1474 assert(e1->x.p->size == 3);
1475 assert(res->x.p->size == 3);
1476 value_init(tmp.d);
1477 evalue_copy(&tmp, &res->x.p->arr[2]);
1478 emul(&d, &tmp);
1479 eadd(&res->x.p->arr[1], &tmp);
1480 emul(&e1->x.p->arr[2], &tmp);
1481 emul(&e1->x.p->arr[1], res);
1482 eadd(&tmp, &res->x.p->arr[2]);
1483 free_evalue_refs(&tmp);
1484 value_clear(d.x.n);
1485 }
1486 value_clear(d.d);
1487 } else {
1488 if(mod_term_smaller(res, e1))
1489 for(i=1; i<res->x.p->size ; i++)
1490 emul(e1, &(res->x.p->arr[i]));
1491 else
1492 emul_rev(e1, res);
1493 return;
1494 }
1495 }
1496 break;
1497 case relation:
1498 emul_rev(e1, res);
1499 break;
1500 default:
1501 assert(0);
1502 }
1503 }
1504 else {
1505 if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
1506 /* Product of two rational numbers */
1507
1508 Value g;
1509 value_init(g);
1510 value_multiply(res->d,e1->d,res->d);
1511 value_multiply(res->x.n,e1->x.n,res->x.n );
1512 Gcd(res->x.n, res->d,&g);
1513 if (value_notone_p(g)) {
1514 value_division(res->d,res->d,g);
1515 value_division(res->x.n,res->x.n,g);
1516 }
1517 value_clear(g);
1518 return ;
1519 }
1520 else {
1521 if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
1522 /* Product of an expression (polynomial or peririodic) and a rational number */
1523
1524 emul_rev(e1, res);
1525 return ;
1526 }
1527 else {
1528 /* Product of a rationel number and an expression (polynomial or peririodic) */
1529
1530 i = type_offset(res->x.p);
1531 for (; i<res->x.p->size; i++)
1532 emul(e1, &res->x.p->arr[i]);
1533
1534 return ;
1535 }
1536 }
1537 }
1538
1539 return ;
1540 }
1541
1542 /* Frees mask content ! */
1543 void emask(evalue *mask, evalue *res) {
1544 int n, i, j;
1545 Polyhedron *d, *fd;
1546 struct section *s;
1547 evalue mone;
1548 int pos;
1549
1550 if (EVALUE_IS_ZERO(*res)) {
1551 free_evalue_refs(mask);
1552 return;
1553 }
1554
1555 assert(value_zero_p(mask->d));
1556 assert(mask->x.p->type == partition);
1557 assert(value_zero_p(res->d));
1558 assert(res->x.p->type == partition);
1559 assert(mask->x.p->pos == res->x.p->pos);
1560 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1561 pos = res->x.p->pos;
1562
1563 s = (struct section *)
1564 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1565 sizeof(struct section));
1566 assert(s);
1567
1568 value_init(mone.d);
1569 evalue_set_si(&mone, -1, 1);
1570
1571 n = 0;
1572 for (j = 0; j < res->x.p->size/2; ++j) {
1573 assert(mask->x.p->size >= 2);
1574 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1575 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1576 if (!emptyQ(fd))
1577 for (i = 1; i < mask->x.p->size/2; ++i) {
1578 Polyhedron *t = fd;
1579 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1580 Domain_Free(t);
1581 if (emptyQ(fd))
1582 break;
1583 }
1584 if (emptyQ(fd)) {
1585 Domain_Free(fd);
1586 continue;
1587 }
1588 value_init(s[n].E.d);
1589 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1590 s[n].D = fd;
1591 ++n;
1592 }
1593 for (i = 0; i < mask->x.p->size/2; ++i) {
1594 if (EVALUE_IS_ONE(mask->x.p->arr[2*i+1]))
1595 continue;
1596
1597 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1598 eadd(&mone, &mask->x.p->arr[2*i+1]);
1599 emul(&mone, &mask->x.p->arr[2*i+1]);
1600 for (j = 0; j < res->x.p->size/2; ++j) {
1601 Polyhedron *t;
1602 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1603 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1604 if (emptyQ(d)) {
1605 Domain_Free(d);
1606 continue;
1607 }
1608 t = fd;
1609 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1610 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1611 Domain_Free(t);
1612 value_init(s[n].E.d);
1613 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1614 emul(&mask->x.p->arr[2*i+1], &s[n].E);
1615 s[n].D = d;
1616 ++n;
1617 }
1618
1619 if (!emptyQ(fd)) {
1620 /* Just ignore; this may have been previously masked off */
1621 }
1622 if (fd != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1623 Domain_Free(fd);
1624 }
1625
1626 free_evalue_refs(&mone);
1627 free_evalue_refs(mask);
1628 free_evalue_refs(res);
1629 value_init(res->d);
1630 if (n == 0)
1631 evalue_set_si(res, 0, 1);
1632 else {
1633 res->x.p = new_enode(partition, 2*n, pos);
1634 for (j = 0; j < n; ++j) {
1635 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1636 value_clear(res->x.p->arr[2*j+1].d);
1637 res->x.p->arr[2*j+1] = s[j].E;
1638 }
1639 }
1640
1641 free(s);
1642 }
1643
1644 void evalue_copy(evalue *dst, evalue *src)
1645 {
1646 value_assign(dst->d, src->d);
1647 if(value_notzero_p(src->d)) {
1648 value_init(dst->x.n);
1649 value_assign(dst->x.n, src->x.n);
1650 } else
1651 dst->x.p = ecopy(src->x.p);
1652 }
1653
1654 enode *new_enode(enode_type type,int size,int pos) {
1655
1656 enode *res;
1657 int i;
1658
1659 if(size == 0) {
1660 fprintf(stderr, "Allocating enode of size 0 !\n" );
1661 return NULL;
1662 }
1663 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1664 res->type = type;
1665 res->size = size;
1666 res->pos = pos;
1667 for(i=0; i<size; i++) {
1668 value_init(res->arr[i].d);
1669 value_set_si(res->arr[i].d,0);
1670 res->arr[i].x.p = 0;
1671 }
1672 return res;
1673 } /* new_enode */
1674
1675 enode *ecopy(enode *e) {
1676
1677 enode *res;
1678 int i;
1679
1680 res = new_enode(e->type,e->size,e->pos);
1681 for(i=0;i<e->size;++i) {
1682 value_assign(res->arr[i].d,e->arr[i].d);
1683 if(value_zero_p(res->arr[i].d))
1684 res->arr[i].x.p = ecopy(e->arr[i].x.p);
1685 else if (EVALUE_IS_DOMAIN(res->arr[i]))
1686 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1687 else {
1688 value_init(res->arr[i].x.n);
1689 value_assign(res->arr[i].x.n,e->arr[i].x.n);
1690 }
1691 }
1692 return(res);
1693 } /* ecopy */
1694
1695 int ecmp(const evalue *e1, const evalue *e2)
1696 {
1697 enode *p1, *p2;
1698 int i;
1699 int r;
1700
1701 if (value_notzero_p(e1->d) && value_notzero_p(e2->d)) {
1702 Value m, m2;
1703 value_init(m);
1704 value_init(m2);
1705 value_multiply(m, e1->x.n, e2->d);
1706 value_multiply(m2, e2->x.n, e1->d);
1707
1708 if (value_lt(m, m2))
1709 r = -1;
1710 else if (value_gt(m, m2))
1711 r = 1;
1712 else
1713 r = 0;
1714
1715 value_clear(m);
1716 value_clear(m2);
1717
1718 return r;
1719 }
1720 if (value_notzero_p(e1->d))
1721 return -1;
1722 if (value_notzero_p(e2->d))
1723 return 1;
1724
1725 p1 = e1->x.p;
1726 p2 = e2->x.p;
1727
1728 if (p1->type != p2->type)
1729 return p1->type - p2->type;
1730 if (p1->pos != p2->pos)
1731 return p1->pos - p2->pos;
1732 if (p1->size != p2->size)
1733 return p1->size - p2->size;
1734
1735 for (i = p1->size-1; i >= 0; --i)
1736 if ((r = ecmp(&p1->arr[i], &p2->arr[i])) != 0)
1737 return r;
1738
1739 return 0;
1740 }
1741
1742 int eequal(evalue *e1,evalue *e2) {
1743
1744 int i;
1745 enode *p1, *p2;
1746
1747 if (value_ne(e1->d,e2->d))
1748 return 0;
1749
1750 /* e1->d == e2->d */
1751 if (value_notzero_p(e1->d)) {
1752 if (value_ne(e1->x.n,e2->x.n))
1753 return 0;
1754
1755 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1756 return 1;
1757 }
1758
1759 /* e1->d == e2->d == 0 */
1760 p1 = e1->x.p;
1761 p2 = e2->x.p;
1762 if (p1->type != p2->type) return 0;
1763 if (p1->size != p2->size) return 0;
1764 if (p1->pos != p2->pos) return 0;
1765 for (i=0; i<p1->size; i++)
1766 if (!eequal(&p1->arr[i], &p2->arr[i]) )
1767 return 0;
1768 return 1;
1769 } /* eequal */
1770
1771 void free_evalue_refs(evalue *e) {
1772
1773 enode *p;
1774 int i;
1775
1776 if (EVALUE_IS_DOMAIN(*e)) {
1777 Domain_Free(EVALUE_DOMAIN(*e));
1778 value_clear(e->d);
1779 return;
1780 } else if (value_pos_p(e->d)) {
1781
1782 /* 'e' stores a constant */
1783 value_clear(e->d);
1784 value_clear(e->x.n);
1785 return;
1786 }
1787 assert(value_zero_p(e->d));
1788 value_clear(e->d);
1789 p = e->x.p;
1790 if (!p) return; /* null pointer */
1791 for (i=0; i<p->size; i++) {
1792 free_evalue_refs(&(p->arr[i]));
1793 }
1794 free(p);
1795 return;
1796 } /* free_evalue_refs */
1797
1798 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
1799 Vector * val, evalue *res)
1800 {
1801 unsigned nparam = periods->Size;
1802
1803 if (p == nparam) {
1804 double d = compute_evalue(e, val->p);
1805 d *= VALUE_TO_DOUBLE(m);
1806 if (d > 0)
1807 d += .25;
1808 else
1809 d -= .25;
1810 value_assign(res->d, m);
1811 value_init(res->x.n);
1812 value_set_double(res->x.n, d);
1813 mpz_fdiv_r(res->x.n, res->x.n, m);
1814 return;
1815 }
1816 if (value_one_p(periods->p[p]))
1817 mod2table_r(e, periods, m, p+1, val, res);
1818 else {
1819 Value tmp;
1820 value_init(tmp);
1821
1822 value_assign(tmp, periods->p[p]);
1823 value_set_si(res->d, 0);
1824 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
1825 do {
1826 value_decrement(tmp, tmp);
1827 value_assign(val->p[p], tmp);
1828 mod2table_r(e, periods, m, p+1, val,
1829 &res->x.p->arr[VALUE_TO_INT(tmp)]);
1830 } while (value_pos_p(tmp));
1831
1832 value_clear(tmp);
1833 }
1834 }
1835
1836 static void rel2table(evalue *e, int zero)
1837 {
1838 if (value_pos_p(e->d)) {
1839 if (value_zero_p(e->x.n) == zero)
1840 value_set_si(e->x.n, 1);
1841 else
1842 value_set_si(e->x.n, 0);
1843 value_set_si(e->d, 1);
1844 } else {
1845 int i;
1846 for (i = 0; i < e->x.p->size; ++i)
1847 rel2table(&e->x.p->arr[i], zero);
1848 }
1849 }
1850
1851 void evalue_mod2table(evalue *e, int nparam)
1852 {
1853 enode *p;
1854 int i;
1855
1856 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
1857 return;
1858 p = e->x.p;
1859 for (i=0; i<p->size; i++) {
1860 evalue_mod2table(&(p->arr[i]), nparam);
1861 }
1862 if (p->type == relation) {
1863 evalue copy;
1864
1865 if (p->size > 2) {
1866 value_init(copy.d);
1867 evalue_copy(&copy, &p->arr[0]);
1868 }
1869 rel2table(&p->arr[0], 1);
1870 emul(&p->arr[0], &p->arr[1]);
1871 if (p->size > 2) {
1872 rel2table(&copy, 0);
1873 emul(&copy, &p->arr[2]);
1874 eadd(&p->arr[2], &p->arr[1]);
1875 free_evalue_refs(&p->arr[2]);
1876 free_evalue_refs(&copy);
1877 }
1878 free_evalue_refs(&p->arr[0]);
1879 value_clear(e->d);
1880 *e = p->arr[1];
1881 free(p);
1882 } else if (p->type == fractional) {
1883 Vector *periods = Vector_Alloc(nparam);
1884 Vector *val = Vector_Alloc(nparam);
1885 Value tmp;
1886 evalue *ev;
1887 evalue EP, res;
1888
1889 value_init(tmp);
1890 value_set_si(tmp, 1);
1891 Vector_Set(periods->p, 1, nparam);
1892 Vector_Set(val->p, 0, nparam);
1893 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
1894 enode *p = ev->x.p;
1895
1896 assert(p->type == polynomial);
1897 assert(p->size == 2);
1898 value_assign(periods->p[p->pos-1], p->arr[1].d);
1899 value_lcm(tmp, p->arr[1].d, &tmp);
1900 }
1901 value_init(EP.d);
1902 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
1903
1904 value_init(res.d);
1905 evalue_set_si(&res, 0, 1);
1906 /* Compute the polynomial using Horner's rule */
1907 for (i=p->size-1;i>1;i--) {
1908 eadd(&p->arr[i], &res);
1909 emul(&EP, &res);
1910 }
1911 eadd(&p->arr[1], &res);
1912
1913 free_evalue_refs(e);
1914 free_evalue_refs(&EP);
1915 *e = res;
1916
1917 value_clear(tmp);
1918 Vector_Free(val);
1919 Vector_Free(periods);
1920 }
1921 } /* evalue_mod2table */
1922
1923 /********************************************************/
1924 /* function in domain */
1925 /* check if the parameters in list_args */
1926 /* verifies the constraints of Domain P */
1927 /********************************************************/
1928 int in_domain(Polyhedron *P, Value *list_args) {
1929
1930 int col,row;
1931 Value v; /* value of the constraint of a row when
1932 parameters are instanciated*/
1933 Value tmp;
1934
1935 value_init(v);
1936 value_init(tmp);
1937
1938 /*P->Constraint constraint matrice of polyhedron P*/
1939 for(row=0;row<P->NbConstraints;row++) {
1940 value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
1941 for(col=1;col<P->Dimension+1;col++) {
1942 value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
1943 value_addto(v,v,tmp);
1944 }
1945 if (value_notzero_p(P->Constraint[row][0])) {
1946
1947 /*if v is not >=0 then this constraint is not respected */
1948 if (value_neg_p(v)) {
1949 next:
1950 value_clear(v);
1951 value_clear(tmp);
1952 return P->next ? in_domain(P->next, list_args) : 0;
1953 }
1954 }
1955 else {
1956
1957 /*if v is not = 0 then this constraint is not respected */
1958 if (value_notzero_p(v))
1959 goto next;
1960 }
1961 }
1962
1963 /*if not return before this point => all
1964 the constraints are respected */
1965 value_clear(v);
1966 value_clear(tmp);
1967 return 1;
1968 } /* in_domain */
1969
1970 /****************************************************/
1971 /* function compute enode */
1972 /* compute the value of enode p with parameters */
1973 /* list "list_args */
1974 /* compute the polynomial or the periodic */
1975 /****************************************************/
1976
1977 static double compute_enode(enode *p, Value *list_args) {
1978
1979 int i;
1980 Value m, param;
1981 double res=0.0;
1982
1983 if (!p)
1984 return(0.);
1985
1986 value_init(m);
1987 value_init(param);
1988
1989 if (p->type == polynomial) {
1990 if (p->size > 1)
1991 value_assign(param,list_args[p->pos-1]);
1992
1993 /* Compute the polynomial using Horner's rule */
1994 for (i=p->size-1;i>0;i--) {
1995 res +=compute_evalue(&p->arr[i],list_args);
1996 res *=VALUE_TO_DOUBLE(param);
1997 }
1998 res +=compute_evalue(&p->arr[0],list_args);
1999 }
2000 else if (p->type == fractional) {
2001 double d = compute_evalue(&p->arr[0], list_args);
2002 d -= floor(d+1e-10);
2003
2004 /* Compute the polynomial using Horner's rule */
2005 for (i=p->size-1;i>1;i--) {
2006 res +=compute_evalue(&p->arr[i],list_args);
2007 res *=d;
2008 }
2009 res +=compute_evalue(&p->arr[1],list_args);
2010 }
2011 else if (p->type == flooring) {
2012 double d = compute_evalue(&p->arr[0], list_args);
2013 d = floor(d+1e-10);
2014
2015 /* Compute the polynomial using Horner's rule */
2016 for (i=p->size-1;i>1;i--) {
2017 res +=compute_evalue(&p->arr[i],list_args);
2018 res *=d;
2019 }
2020 res +=compute_evalue(&p->arr[1],list_args);
2021 }
2022 else if (p->type == periodic) {
2023 value_assign(param,list_args[p->pos-1]);
2024
2025 /* Choose the right element of the periodic */
2026 value_absolute(m,param);
2027 value_set_si(param,p->size);
2028 value_modulus(m,m,param);
2029 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
2030 }
2031 else if (p->type == relation) {
2032 if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
2033 res = compute_evalue(&p->arr[1], list_args);
2034 else if (p->size > 2)
2035 res = compute_evalue(&p->arr[2], list_args);
2036 }
2037 else if (p->type == partition) {
2038 for (i = 0; i < p->size/2; ++i)
2039 if (in_domain(EVALUE_DOMAIN(p->arr[2*i]), list_args)) {
2040 res = compute_evalue(&p->arr[2*i+1], list_args);
2041 break;
2042 }
2043 }
2044 else
2045 assert(0);
2046 value_clear(m);
2047 value_clear(param);
2048 return res;
2049 } /* compute_enode */
2050
2051 /*************************************************/
2052 /* return the value of Ehrhart Polynomial */
2053 /* It returns a double, because since it is */
2054 /* a recursive function, some intermediate value */
2055 /* might not be integral */
2056 /*************************************************/
2057
2058 double compute_evalue(evalue *e,Value *list_args) {
2059
2060 double res;
2061
2062 if (value_notzero_p(e->d)) {
2063 if (value_notone_p(e->d))
2064 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
2065 else
2066 res = VALUE_TO_DOUBLE(e->x.n);
2067 }
2068 else
2069 res = compute_enode(e->x.p,list_args);
2070 return res;
2071 } /* compute_evalue */
2072
2073
2074 /****************************************************/
2075 /* function compute_poly : */
2076 /* Check for the good validity domain */
2077 /* return the number of point in the Polyhedron */
2078 /* in allocated memory */
2079 /* Using the Ehrhart pseudo-polynomial */
2080 /****************************************************/
2081 Value *compute_poly(Enumeration *en,Value *list_args) {
2082
2083 Value *tmp;
2084 /* double d; int i; */
2085
2086 tmp = (Value *) malloc (sizeof(Value));
2087 assert(tmp != NULL);
2088 value_init(*tmp);
2089 value_set_si(*tmp,0);
2090
2091 if(!en)
2092 return(tmp); /* no ehrhart polynomial */
2093 if(en->ValidityDomain) {
2094 if(!en->ValidityDomain->Dimension) { /* no parameters */
2095 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2096 return(tmp);
2097 }
2098 }
2099 else
2100 return(tmp); /* no Validity Domain */
2101 while(en) {
2102 if(in_domain(en->ValidityDomain,list_args)) {
2103
2104 #ifdef EVAL_EHRHART_DEBUG
2105 Print_Domain(stdout,en->ValidityDomain);
2106 print_evalue(stdout,&en->EP);
2107 #endif
2108
2109 /* d = compute_evalue(&en->EP,list_args);
2110 i = d;
2111 printf("(double)%lf = %d\n", d, i ); */
2112 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2113 return(tmp);
2114 }
2115 else
2116 en=en->next;
2117 }
2118 value_set_si(*tmp,0);
2119 return(tmp); /* no compatible domain with the arguments */
2120 } /* compute_poly */
2121
2122 size_t value_size(Value v) {
2123 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2124 * sizeof(v[0]._mp_d[0]);
2125 }
2126
2127 size_t domain_size(Polyhedron *D)
2128 {
2129 int i, j;
2130 size_t s = sizeof(*D);
2131
2132 for (i = 0; i < D->NbConstraints; ++i)
2133 for (j = 0; j < D->Dimension+2; ++j)
2134 s += value_size(D->Constraint[i][j]);
2135
2136 /*
2137 for (i = 0; i < D->NbRays; ++i)
2138 for (j = 0; j < D->Dimension+2; ++j)
2139 s += value_size(D->Ray[i][j]);
2140 */
2141
2142 return D->next ? s+domain_size(D->next) : s;
2143 }
2144
2145 size_t enode_size(enode *p) {
2146 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2147 int i;
2148
2149 if (p->type == partition)
2150 for (i = 0; i < p->size/2; ++i) {
2151 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2152 s += evalue_size(&p->arr[2*i+1]);
2153 }
2154 else
2155 for (i = 0; i < p->size; ++i) {
2156 s += evalue_size(&p->arr[i]);
2157 }
2158 return s;
2159 }
2160
2161 size_t evalue_size(evalue *e)
2162 {
2163 size_t s = sizeof(*e);
2164 s += value_size(e->d);
2165 if (value_notzero_p(e->d))
2166 s += value_size(e->x.n);
2167 else
2168 s += enode_size(e->x.p);
2169 return s;
2170 }
2171
2172 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2173 {
2174 evalue *found = NULL;
2175 evalue offset;
2176 evalue copy;
2177 int i;
2178
2179 if (value_pos_p(e->d) || e->x.p->type != fractional)
2180 return NULL;
2181
2182 value_init(offset.d);
2183 value_init(offset.x.n);
2184 poly_denom(&e->x.p->arr[0], &offset.d);
2185 value_lcm(m, offset.d, &offset.d);
2186 value_set_si(offset.x.n, 1);
2187
2188 value_init(copy.d);
2189 evalue_copy(&copy, cst);
2190
2191 eadd(&offset, cst);
2192 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2193
2194 if (eequal(base, &e->x.p->arr[0]))
2195 found = &e->x.p->arr[0];
2196 else {
2197 value_set_si(offset.x.n, -2);
2198
2199 eadd(&offset, cst);
2200 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2201
2202 if (eequal(base, &e->x.p->arr[0]))
2203 found = base;
2204 }
2205 free_evalue_refs(cst);
2206 free_evalue_refs(&offset);
2207 *cst = copy;
2208
2209 for (i = 1; !found && i < e->x.p->size; ++i)
2210 found = find_second(base, cst, &e->x.p->arr[i], m);
2211
2212 return found;
2213 }
2214
2215 static evalue *find_relation_pair(evalue *e)
2216 {
2217 int i;
2218 evalue *found = NULL;
2219
2220 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2221 return NULL;
2222
2223 if (e->x.p->type == fractional) {
2224 Value m;
2225 evalue *cst;
2226
2227 value_init(m);
2228 poly_denom(&e->x.p->arr[0], &m);
2229
2230 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2231 cst = &cst->x.p->arr[0])
2232 ;
2233
2234 for (i = 1; !found && i < e->x.p->size; ++i)
2235 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2236
2237 value_clear(m);
2238 }
2239
2240 i = e->x.p->type == relation;
2241 for (; !found && i < e->x.p->size; ++i)
2242 found = find_relation_pair(&e->x.p->arr[i]);
2243
2244 return found;
2245 }
2246
2247 void evalue_mod2relation(evalue *e) {
2248 evalue *d;
2249
2250 if (value_zero_p(e->d) && e->x.p->type == partition) {
2251 int i;
2252
2253 for (i = 0; i < e->x.p->size/2; ++i) {
2254 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2255 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2256 value_clear(e->x.p->arr[2*i].d);
2257 free_evalue_refs(&e->x.p->arr[2*i+1]);
2258 e->x.p->size -= 2;
2259 if (2*i < e->x.p->size) {
2260 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2261 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2262 }
2263 --i;
2264 }
2265 }
2266 if (e->x.p->size == 0) {
2267 free(e->x.p);
2268 evalue_set_si(e, 0, 1);
2269 }
2270
2271 return;
2272 }
2273
2274 while ((d = find_relation_pair(e)) != NULL) {
2275 evalue split;
2276 evalue *ev;
2277
2278 value_init(split.d);
2279 value_set_si(split.d, 0);
2280 split.x.p = new_enode(relation, 3, 0);
2281 evalue_set_si(&split.x.p->arr[1], 1, 1);
2282 evalue_set_si(&split.x.p->arr[2], 1, 1);
2283
2284 ev = &split.x.p->arr[0];
2285 value_set_si(ev->d, 0);
2286 ev->x.p = new_enode(fractional, 3, -1);
2287 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2288 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2289 evalue_copy(&ev->x.p->arr[0], d);
2290
2291 emul(&split, e);
2292
2293 reduce_evalue(e);
2294
2295 free_evalue_refs(&split);
2296 }
2297 }
2298
2299 static int evalue_comp(const void * a, const void * b)
2300 {
2301 const evalue *e1 = *(const evalue **)a;
2302 const evalue *e2 = *(const evalue **)b;
2303 return ecmp(e1, e2);
2304 }
2305
2306 void evalue_combine(evalue *e)
2307 {
2308 evalue **evs;
2309 int i, k;
2310 enode *p;
2311 evalue tmp;
2312
2313 if (value_notzero_p(e->d) || e->x.p->type != partition)
2314 return;
2315
2316 NALLOC(evs, e->x.p->size/2);
2317 for (i = 0; i < e->x.p->size/2; ++i)
2318 evs[i] = &e->x.p->arr[2*i+1];
2319 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2320 p = new_enode(partition, e->x.p->size, e->x.p->pos);
2321 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2322 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2323 value_clear(p->arr[2*k].d);
2324 value_clear(p->arr[2*k+1].d);
2325 p->arr[2*k] = *(evs[i]-1);
2326 p->arr[2*k+1] = *(evs[i]);
2327 ++k;
2328 } else {
2329 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2330 Polyhedron *L = D;
2331
2332 value_clear((evs[i]-1)->d);
2333
2334 while (L->next)
2335 L = L->next;
2336 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2337 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2338 free_evalue_refs(evs[i]);
2339 }
2340 }
2341
2342 for (i = 2*k ; i < p->size; ++i)
2343 value_clear(p->arr[i].d);
2344
2345 free(evs);
2346 free(e->x.p);
2347 p->size = 2*k;
2348 e->x.p = p;
2349
2350 for (i = 0; i < e->x.p->size/2; ++i) {
2351 Polyhedron *H;
2352 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2353 continue;
2354 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2355 if (H == NULL)
2356 continue;
2357 for (k = 0; k < e->x.p->size/2; ++k) {
2358 Polyhedron *D, *N, **P;
2359 if (i == k)
2360 continue;
2361 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2362 D = *P;
2363 if (D == NULL)
2364 continue;
2365 for (; D; D = N) {
2366 *P = D;
2367 N = D->next;
2368 if (D->NbEq <= H->NbEq) {
2369 P = &D->next;
2370 continue;
2371 }
2372
2373 value_init(tmp.d);
2374 tmp.x.p = new_enode(partition, 2, e->x.p->pos);
2375 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2376 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2377 reduce_evalue(&tmp);
2378 if (value_notzero_p(tmp.d) ||
2379 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2380 P = &D->next;
2381 else {
2382 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2383 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2384 *P = NULL;
2385 }
2386 free_evalue_refs(&tmp);
2387 }
2388 }
2389 Polyhedron_Free(H);
2390 }
2391
2392 for (i = 0; i < e->x.p->size/2; ++i) {
2393 Polyhedron *H, *E;
2394 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2395 if (!D) {
2396 value_clear(e->x.p->arr[2*i].d);
2397 free_evalue_refs(&e->x.p->arr[2*i+1]);
2398 e->x.p->size -= 2;
2399 if (2*i < e->x.p->size) {
2400 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2401 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2402 }
2403 --i;
2404 continue;
2405 }
2406 if (!D->next)
2407 continue;
2408 H = DomainConvex(D, 0);
2409 E = DomainDifference(H, D, 0);
2410 Domain_Free(D);
2411 D = DomainDifference(H, E, 0);
2412 Domain_Free(H);
2413 Domain_Free(E);
2414 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2415 }
2416 }
2417
2418 /* May change coefficients to become non-standard if fiddle is set
2419 * => reduce p afterwards to correct
2420 */
2421 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2422 Matrix **R, int fiddle)
2423 {
2424 Polyhedron *I, *H;
2425 evalue *pp;
2426 unsigned dim = D->Dimension;
2427 Matrix *T = Matrix_Alloc(2, dim+1);
2428 Value twice;
2429 value_init(twice);
2430 assert(T);
2431
2432 assert(p->type == fractional);
2433 pp = &p->arr[0];
2434 value_set_si(T->p[1][dim], 1);
2435 poly_denom(pp, d);
2436 while (value_zero_p(pp->d)) {
2437 assert(pp->x.p->type == polynomial);
2438 assert(pp->x.p->size == 2);
2439 assert(value_notzero_p(pp->x.p->arr[1].d));
2440 value_division(T->p[0][pp->x.p->pos-1], *d, pp->x.p->arr[1].d);
2441 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2442 if (fiddle && value_gt(twice, pp->x.p->arr[1].d))
2443 value_substract(pp->x.p->arr[1].x.n,
2444 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2445 value_multiply(T->p[0][pp->x.p->pos-1],
2446 T->p[0][pp->x.p->pos-1], pp->x.p->arr[1].x.n);
2447 pp = &pp->x.p->arr[0];
2448 }
2449 value_division(T->p[0][dim], *d, pp->d);
2450 value_multiply(T->p[0][dim], T->p[0][dim], pp->x.n);
2451 I = DomainImage(D, T, 0);
2452 H = DomainConvex(I, 0);
2453 Domain_Free(I);
2454 *R = T;
2455
2456 value_clear(twice);
2457
2458 return H;
2459 }
2460
2461 static int reduce_in_domain(evalue *e, Polyhedron *D)
2462 {
2463 int i;
2464 enode *p;
2465 Value d, min, max;
2466 int r = 0;
2467 Polyhedron *I;
2468 Matrix *T;
2469 int bounded;
2470
2471 if (value_notzero_p(e->d))
2472 return r;
2473
2474 p = e->x.p;
2475
2476 if (p->type == relation) {
2477 int equal;
2478 value_init(d);
2479 value_init(min);
2480 value_init(max);
2481
2482 I = polynomial_projection(p->arr[0].x.p, D, &d, &T, 1);
2483 bounded = line_minmax(I, &min, &max); /* frees I */
2484 equal = value_eq(min, max);
2485 mpz_cdiv_q(min, min, d);
2486 mpz_fdiv_q(max, max, d);
2487
2488 if (bounded && value_gt(min, max)) {
2489 /* Never zero */
2490 if (p->size == 3) {
2491 value_clear(e->d);
2492 *e = p->arr[2];
2493 } else {
2494 evalue_set_si(e, 0, 1);
2495 r = 1;
2496 }
2497 free_evalue_refs(&(p->arr[1]));
2498 free_evalue_refs(&(p->arr[0]));
2499 free(p);
2500 value_clear(d);
2501 value_clear(min);
2502 value_clear(max);
2503 Matrix_Free(T);
2504 return r ? r : reduce_in_domain(e, D);
2505 } else if (bounded && equal) {
2506 /* Always zero */
2507 if (p->size == 3)
2508 free_evalue_refs(&(p->arr[2]));
2509 value_clear(e->d);
2510 *e = p->arr[1];
2511 free_evalue_refs(&(p->arr[0]));
2512 free(p);
2513 value_clear(d);
2514 value_clear(min);
2515 value_clear(max);
2516 Matrix_Free(T);
2517 return reduce_in_domain(e, D);
2518 } else if (bounded && value_eq(min, max)) {
2519 /* zero for a single value */
2520 Polyhedron *E;
2521 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2522 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2523 value_multiply(min, min, d);
2524 value_substract(M->p[0][D->Dimension+1],
2525 M->p[0][D->Dimension+1], min);
2526 E = DomainAddConstraints(D, M, 0);
2527 value_clear(d);
2528 value_clear(min);
2529 value_clear(max);
2530 Matrix_Free(T);
2531 Matrix_Free(M);
2532 r = reduce_in_domain(&p->arr[1], E);
2533 if (p->size == 3)
2534 r |= reduce_in_domain(&p->arr[2], D);
2535 Domain_Free(E);
2536 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2537 return r;
2538 }
2539
2540 value_clear(d);
2541 value_clear(min);
2542 value_clear(max);
2543 Matrix_Free(T);
2544 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2545 }
2546
2547 i = p->type == relation ? 1 :
2548 p->type == fractional ? 1 : 0;
2549 for (; i<p->size; i++)
2550 r |= reduce_in_domain(&p->arr[i], D);
2551
2552 if (p->type != fractional) {
2553 if (r && p->type == polynomial) {
2554 evalue f;
2555 value_init(f.d);
2556 value_set_si(f.d, 0);
2557 f.x.p = new_enode(polynomial, 2, p->pos);
2558 evalue_set_si(&f.x.p->arr[0], 0, 1);
2559 evalue_set_si(&f.x.p->arr[1], 1, 1);
2560 reorder_terms(p, &f);
2561 value_clear(e->d);
2562 *e = p->arr[0];
2563 free(p);
2564 }
2565 return r;
2566 }
2567
2568 value_init(d);
2569 value_init(min);
2570 value_init(max);
2571 I = polynomial_projection(p, D, &d, &T, 1);
2572 Matrix_Free(T);
2573 bounded = line_minmax(I, &min, &max); /* frees I */
2574 mpz_fdiv_q(min, min, d);
2575 mpz_fdiv_q(max, max, d);
2576 value_substract(d, max, min);
2577
2578 if (bounded && value_eq(min, max)) {
2579 evalue inc;
2580 value_init(inc.d);
2581 value_init(inc.x.n);
2582 value_set_si(inc.d, 1);
2583 value_oppose(inc.x.n, min);
2584 eadd(&inc, &p->arr[0]);
2585 reorder_terms(p, &p->arr[0]); /* frees arr[0] */
2586 value_clear(e->d);
2587 *e = p->arr[1];
2588 free(p);
2589 free_evalue_refs(&inc);
2590 r = 1;
2591 } else if (bounded && value_one_p(d) && p->size > 3) {
2592 evalue rem;
2593 value_init(rem.d);
2594 value_set_si(rem.d, 0);
2595 rem.x.p = new_enode(fractional, 3, -1);
2596 evalue_copy(&rem.x.p->arr[0], &p->arr[0]);
2597 rem.x.p->arr[1] = p->arr[1];
2598 rem.x.p->arr[2] = p->arr[2];
2599 for (i = 3; i < p->size; ++i)
2600 p->arr[i-2] = p->arr[i];
2601 p->size -= 2;
2602
2603 evalue inc;
2604 value_init(inc.d);
2605 value_init(inc.x.n);
2606 value_set_si(inc.d, 1);
2607 value_oppose(inc.x.n, min);
2608
2609 evalue t;
2610 value_init(t.d);
2611 evalue_copy(&t, &p->arr[0]);
2612 eadd(&inc, &t);
2613
2614 evalue f;
2615 value_init(f.d);
2616 value_set_si(f.d, 0);
2617 f.x.p = new_enode(fractional, 3, -1);
2618 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2619 evalue_set_si(&f.x.p->arr[1], 1, 1);
2620 evalue_set_si(&f.x.p->arr[2], 2, 1);
2621
2622 evalue factor;
2623 value_init(factor.d);
2624 evalue_set_si(&factor, -1, 1);
2625 emul(&t, &factor);
2626
2627 eadd(&f, &factor);
2628 emul(&t, &factor);
2629
2630 evalue_set_si(&f.x.p->arr[1], 0, 1);
2631 evalue_set_si(&f.x.p->arr[2], -1, 1);
2632 eadd(&f, &factor);
2633
2634 emul(&factor, e);
2635 eadd(&rem, e);
2636
2637 free_evalue_refs(&inc);
2638 free_evalue_refs(&t);
2639 free_evalue_refs(&f);
2640 free_evalue_refs(&factor);
2641 free_evalue_refs(&rem);
2642
2643 reduce_in_domain(e, D);
2644
2645 r = 1;
2646 } else {
2647 _reduce_evalue(&p->arr[0], 0, 1);
2648 if (r == 1) {
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 value_clear(f.x.p->arr[0].d);
2654 f.x.p->arr[0] = p->arr[0];
2655 evalue_set_si(&f.x.p->arr[1], 0, 1);
2656 evalue_set_si(&f.x.p->arr[2], 1, 1);
2657 reorder_terms(p, &f);
2658 value_clear(e->d);
2659 *e = p->arr[1];
2660 free(p);
2661 }
2662 }
2663
2664 value_clear(d);
2665 value_clear(min);
2666 value_clear(max);
2667
2668 return r;
2669 }
2670
2671 void evalue_range_reduction(evalue *e)
2672 {
2673 int i;
2674 if (value_notzero_p(e->d) || e->x.p->type != partition)
2675 return;
2676
2677 for (i = 0; i < e->x.p->size/2; ++i)
2678 if (reduce_in_domain(&e->x.p->arr[2*i+1],
2679 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
2680 reduce_evalue(&e->x.p->arr[2*i+1]);
2681
2682 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2683 free_evalue_refs(&e->x.p->arr[2*i+1]);
2684 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
2685 value_clear(e->x.p->arr[2*i].d);
2686 e->x.p->size -= 2;
2687 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2688 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2689 --i;
2690 }
2691 }
2692 }
2693
2694 /* Frees EP
2695 */
2696 Enumeration* partition2enumeration(evalue *EP)
2697 {
2698 int i;
2699 Enumeration *en, *res = NULL;
2700
2701 if (EVALUE_IS_ZERO(*EP)) {
2702 free(EP);
2703 return res;
2704 }
2705
2706 for (i = 0; i < EP->x.p->size/2; ++i) {
2707 en = (Enumeration *)malloc(sizeof(Enumeration));
2708 en->next = res;
2709 res = en;
2710 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2711 value_clear(EP->x.p->arr[2*i].d);
2712 res->EP = EP->x.p->arr[2*i+1];
2713 }
2714 free(EP->x.p);
2715 value_clear(EP->d);
2716 free(EP);
2717 return res;
2718 }
2719
2720 static int frac2floor_in_domain(evalue *e, Polyhedron *D)
2721 {
2722 enode *p;
2723 int r = 0;
2724 int i;
2725 Polyhedron *I;
2726 Matrix *T;
2727 Value d, min;
2728 evalue fl;
2729
2730 if (value_notzero_p(e->d))
2731 return r;
2732
2733 p = e->x.p;
2734
2735 i = p->type == relation ? 1 :
2736 p->type == fractional ? 1 : 0;
2737 for (; i<p->size; i++)
2738 r |= frac2floor_in_domain(&p->arr[i], D);
2739
2740 if (p->type != fractional) {
2741 if (r && p->type == polynomial) {
2742 evalue f;
2743 value_init(f.d);
2744 value_set_si(f.d, 0);
2745 f.x.p = new_enode(polynomial, 2, p->pos);
2746 evalue_set_si(&f.x.p->arr[0], 0, 1);
2747 evalue_set_si(&f.x.p->arr[1], 1, 1);
2748 reorder_terms(p, &f);
2749 value_clear(e->d);
2750 *e = p->arr[0];
2751 free(p);
2752 }
2753 return r;
2754 }
2755
2756 value_init(d);
2757 I = polynomial_projection(p, D, &d, &T, 0);
2758
2759 /*
2760 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2761 */
2762
2763 assert(I->NbEq == 0); /* Should have been reduced */
2764
2765 /* Find minimum */
2766 for (i = 0; i < I->NbConstraints; ++i)
2767 if (value_pos_p(I->Constraint[i][1]))
2768 break;
2769
2770 assert(i < I->NbConstraints);
2771 value_init(min);
2772 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
2773 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
2774 if (value_neg_p(min)) {
2775 evalue offset;
2776 mpz_fdiv_q(min, min, d);
2777 value_init(offset.d);
2778 value_set_si(offset.d, 1);
2779 value_init(offset.x.n);
2780 value_oppose(offset.x.n, min);
2781 eadd(&offset, &p->arr[0]);
2782 free_evalue_refs(&offset);
2783 }
2784
2785 Polyhedron_Free(I);
2786 Matrix_Free(T);
2787 value_clear(min);
2788 value_clear(d);
2789
2790 value_init(fl.d);
2791 value_set_si(fl.d, 0);
2792 fl.x.p = new_enode(flooring, 3, -1);
2793 evalue_set_si(&fl.x.p->arr[1], 0, 1);
2794 evalue_set_si(&fl.x.p->arr[2], -1, 1);
2795 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
2796
2797 eadd(&fl, &p->arr[0]);
2798 reorder_terms(p, &p->arr[0]);
2799 *e = p->arr[1];
2800 free(p);
2801 free_evalue_refs(&fl);
2802
2803 return 1;
2804 }
2805
2806 void evalue_frac2floor(evalue *e)
2807 {
2808 int i;
2809 if (value_notzero_p(e->d) || e->x.p->type != partition)
2810 return;
2811
2812 for (i = 0; i < e->x.p->size/2; ++i)
2813 if (frac2floor_in_domain(&e->x.p->arr[2*i+1],
2814 EVALUE_DOMAIN(e->x.p->arr[2*i])))
2815 reduce_evalue(&e->x.p->arr[2*i+1]);
2816 }
2817
2818 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
2819 Vector *row)
2820 {
2821 int nr, nc;
2822 int i;
2823 int nparam = D->Dimension - nvar;
2824
2825 if (C == 0) {
2826 nr = D->NbConstraints + 2;
2827 nc = D->Dimension + 2 + 1;
2828 C = Matrix_Alloc(nr, nc);
2829 for (i = 0; i < D->NbConstraints; ++i) {
2830 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
2831 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2832 D->Dimension + 1 - nvar);
2833 }
2834 } else {
2835 Matrix *oldC = C;
2836 nr = C->NbRows + 2;
2837 nc = C->NbColumns + 1;
2838 C = Matrix_Alloc(nr, nc);
2839 for (i = 0; i < oldC->NbRows; ++i) {
2840 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
2841 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2842 oldC->NbColumns - 1 - nvar);
2843 }
2844 }
2845 value_set_si(C->p[nr-2][0], 1);
2846 value_set_si(C->p[nr-2][1 + nvar], 1);
2847 value_set_si(C->p[nr-2][nc - 1], -1);
2848
2849 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
2850 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
2851 1 + nparam);
2852
2853 return C;
2854 }
2855
2856 static void floor2frac_r(evalue *e, int nvar)
2857 {
2858 enode *p;
2859 int i;
2860 evalue f;
2861 evalue *pp;
2862
2863 if (value_notzero_p(e->d))
2864 return;
2865
2866 p = e->x.p;
2867
2868 assert(p->type == flooring);
2869 for (i = 1; i < p->size; i++)
2870 floor2frac_r(&p->arr[i], nvar);
2871
2872 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
2873 assert(pp->x.p->type == polynomial);
2874 pp->x.p->pos -= nvar;
2875 }
2876
2877 value_init(f.d);
2878 value_set_si(f.d, 0);
2879 f.x.p = new_enode(fractional, 3, -1);
2880 evalue_set_si(&f.x.p->arr[1], 0, 1);
2881 evalue_set_si(&f.x.p->arr[2], -1, 1);
2882 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2883
2884 eadd(&f, &p->arr[0]);
2885 reorder_terms(p, &p->arr[0]);
2886 *e = p->arr[1];
2887 free(p);
2888 free_evalue_refs(&f);
2889 }
2890
2891 /* Convert flooring back to fractional and shift position
2892 * of the parameters by nvar
2893 */
2894 static void floor2frac(evalue *e, int nvar)
2895 {
2896 floor2frac_r(e, nvar);
2897 reduce_evalue(e);
2898 }
2899
2900 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
2901 {
2902 evalue *t;
2903 int nparam = D->Dimension - nvar;
2904
2905 if (C != 0) {
2906 C = Matrix_Copy(C);
2907 D = Constraints2Polyhedron(C, 0);
2908 Matrix_Free(C);
2909 }
2910
2911 t = barvinok_enumerate_e(D, 0, nparam, 0);
2912
2913 /* Double check that D was not unbounded. */
2914 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
2915
2916 if (C != 0)
2917 Polyhedron_Free(D);
2918
2919 return t;
2920 }
2921
2922 evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
2923 Matrix *C)
2924 {
2925 Vector *row = NULL;
2926 int i;
2927 evalue *res;
2928 Matrix *origC;
2929 evalue *factor = NULL;
2930 evalue cum;
2931
2932 if (EVALUE_IS_ZERO(*e))
2933 return 0;
2934
2935 if (D->next) {
2936 Polyhedron *DD = Disjoint_Domain(D, 0, 0);
2937 Polyhedron *Q;
2938
2939 Q = DD;
2940 DD = Q->next;
2941 Q->next = 0;
2942
2943 res = esum_over_domain(e, nvar, Q, C);
2944 Polyhedron_Free(Q);
2945
2946 for (Q = DD; Q; Q = DD) {
2947 evalue *t;
2948
2949 DD = Q->next;
2950 Q->next = 0;
2951
2952 t = esum_over_domain(e, nvar, Q, C);
2953 Polyhedron_Free(Q);
2954
2955 if (!res)
2956 res = t;
2957 else if (t) {
2958 eadd(t, res);
2959 free_evalue_refs(t);
2960 free(t);
2961 }
2962 }
2963 return res;
2964 }
2965
2966 if (value_notzero_p(e->d)) {
2967 evalue *t;
2968
2969 t = esum_over_domain_cst(nvar, D, C);
2970
2971 if (!EVALUE_IS_ONE(*e))
2972 emul(e, t);
2973
2974 return t;
2975 }
2976
2977 switch (e->x.p->type) {
2978 case flooring: {
2979 evalue *pp = &e->x.p->arr[0];
2980
2981 if (pp->x.p->pos > nvar) {
2982 /* remainder is independent of the summated vars */
2983 evalue f;
2984 evalue *t;
2985
2986 value_init(f.d);
2987 evalue_copy(&f, e);
2988 floor2frac(&f, nvar);
2989
2990 t = esum_over_domain_cst(nvar, D, C);
2991
2992 emul(&f, t);
2993
2994 free_evalue_refs(&f);
2995
2996 return t;
2997 }
2998
2999 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3000 poly_denom(pp, &row->p[1 + nvar]);
3001 value_set_si(row->p[0], 1);
3002 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3003 pp = &pp->x.p->arr[0]) {
3004 int pos;
3005 assert(pp->x.p->type == polynomial);
3006 pos = pp->x.p->pos;
3007 if (pos >= 1 + nvar)
3008 ++pos;
3009 value_assign(row->p[pos], row->p[1+nvar]);
3010 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3011 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3012 }
3013 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3014 value_division(row->p[1 + D->Dimension + 1],
3015 row->p[1 + D->Dimension + 1],
3016 pp->d);
3017 value_multiply(row->p[1 + D->Dimension + 1],
3018 row->p[1 + D->Dimension + 1],
3019 pp->x.n);
3020 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3021 break;
3022 }
3023 case polynomial: {
3024 int pos = e->x.p->pos;
3025
3026 if (pos > nvar) {
3027 ALLOC(factor);
3028 value_init(factor->d);
3029 value_set_si(factor->d, 0);
3030 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3031 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3032 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3033 break;
3034 }
3035
3036 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3037 for (i = 0; i < D->NbRays; ++i)
3038 if (value_notzero_p(D->Ray[i][pos]))
3039 break;
3040 assert(i < D->NbRays);
3041 if (value_neg_p(D->Ray[i][pos])) {
3042 ALLOC(factor);
3043 value_init(factor->d);
3044 evalue_set_si(factor, -1, 1);
3045 }
3046 value_set_si(row->p[0], 1);
3047 value_set_si(row->p[pos], 1);
3048 value_set_si(row->p[1 + nvar], -1);
3049 break;
3050 }
3051 default:
3052 assert(0);
3053 }
3054
3055 i = type_offset(e->x.p);
3056
3057 res = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3058 ++i;
3059
3060 if (factor) {
3061 value_init(cum.d);
3062 evalue_copy(&cum, factor);
3063 }
3064
3065 origC = C;
3066 for (; i < e->x.p->size; ++i) {
3067 evalue *t;
3068 if (row) {
3069 Matrix *prevC = C;
3070 C = esum_add_constraint(nvar, D, C, row);
3071 if (prevC != origC)
3072 Matrix_Free(prevC);
3073 }
3074 /*
3075 if (row)
3076 Vector_Print(stderr, P_VALUE_FMT, row);
3077 if (C)
3078 Matrix_Print(stderr, P_VALUE_FMT, C);
3079 */
3080 t = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3081
3082 if (t && factor)
3083 emul(&cum, t);
3084
3085 if (!res)
3086 res = t;
3087 else if (t) {
3088 eadd(t, res);
3089 free_evalue_refs(t);
3090 free(t);
3091 }
3092 if (factor && i+1 < e->x.p->size)
3093 emul(factor, &cum);
3094 }
3095 if (C != origC)
3096 Matrix_Free(C);
3097
3098 if (factor) {
3099 free_evalue_refs(factor);
3100 free_evalue_refs(&cum);
3101 free(factor);
3102 }
3103
3104 if (row)
3105 Vector_Free(row);
3106
3107 reduce_evalue(res);
3108
3109 return res;
3110 }
3111
3112 evalue *esum(evalue *e, int nvar)
3113 {
3114 int i;
3115 evalue *res;
3116 ALLOC(res);
3117 value_init(res->d);
3118
3119 assert(nvar >= 0);
3120 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3121 evalue_copy(res, e);
3122 return res;
3123 }
3124
3125 evalue_set_si(res, 0, 1);
3126
3127 assert(value_zero_p(e->d));
3128 assert(e->x.p->type == partition);
3129
3130 for (i = 0; i < e->x.p->size/2; ++i) {
3131 evalue *t;
3132 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3133 EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
3134 eadd(t, res);
3135 free_evalue_refs(t);
3136 free(t);
3137 }
3138
3139 reduce_evalue(res);
3140
3141 return res;
3142 }
lattice point enumeration library