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