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