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