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