?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
remove debugging output in unfringe
[barvinok.git] / ev_operations.c
blobb2acdb9774d0150ccaeacf8ed37d7cba7fc94040
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
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 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 (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 (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 (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 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 (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 (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 char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
3132 evalue *t;
3133 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3134 EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
3135 Polyhedron_Print(stderr, P_VALUE_FMT, EVALUE_DOMAIN(e->x.p->arr[2*i]));
3136 print_evalue(stderr, t, test);
3137 eadd(t, res);
3138 free_evalue_refs(t);
3139 free(t);
3140 }
3141
3142 reduce_evalue(res);
3143
3144 return res;
3145 }
lattice point enumeration library