?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
remove polyhedron_range
[barvinok.git] / evalue.c
blobd96d312d3a04d9a62835027476c02d5ca415488e
1 /***********************************************************************/
2 /* copyright 1997, Doran Wilde */
3 /* copyright 1997-2000, Vincent Loechner */
4 /* copyright 2003-2006, Sven Verdoolaege */
5 /* Permission is granted to copy, use, and distribute */
6 /* for any commercial or noncommercial purpose under the terms */
7 /* of the GNU General Public license, version 2, June 1991 */
8 /* (see file : LICENSE). */
9 /***********************************************************************/
10
11 #include <assert.h>
12 #include <math.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <barvinok/evalue.h>
16 #include <barvinok/barvinok.h>
17 #include <barvinok/util.h>
18 #include "summate.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 void evalue_set_reduce(evalue *ev, Value n, Value d) {
46 value_init(ev->x.n);
47 value_gcd(ev->x.n, n, d);
48 value_divexact(ev->d, d, ev->x.n);
49 value_divexact(ev->x.n, n, ev->x.n);
50 }
51
52 evalue* evalue_zero()
53 {
54 evalue *EP = ALLOC(evalue);
55 value_init(EP->d);
56 evalue_set_si(EP, 0, 1);
57 return EP;
58 }
59
60 evalue *evalue_nan()
61 {
62 evalue *EP = ALLOC(evalue);
63 value_init(EP->d);
64 value_set_si(EP->d, -2);
65 EP->x.p = NULL;
66 return EP;
67 }
68
69 /* returns an evalue that corresponds to
70 *
71 * x_var
72 */
73 evalue *evalue_var(int var)
74 {
75 evalue *EP = ALLOC(evalue);
76 value_init(EP->d);
77 value_set_si(EP->d,0);
78 EP->x.p = new_enode(polynomial, 2, var + 1);
79 evalue_set_si(&EP->x.p->arr[0], 0, 1);
80 evalue_set_si(&EP->x.p->arr[1], 1, 1);
81 return EP;
82 }
83
84 void aep_evalue(evalue *e, int *ref) {
85
86 enode *p;
87 int i;
88
89 if (value_notzero_p(e->d))
90 return; /* a rational number, its already reduced */
91 if(!(p = e->x.p))
92 return; /* hum... an overflow probably occured */
93
94 /* First check the components of p */
95 for (i=0;i<p->size;i++)
96 aep_evalue(&p->arr[i],ref);
97
98 /* Then p itself */
99 switch (p->type) {
100 case polynomial:
101 case periodic:
102 case evector:
103 p->pos = ref[p->pos-1]+1;
104 }
105 return;
106 } /* aep_evalue */
107
108 /** Comments */
109 void addeliminatedparams_evalue(evalue *e,Matrix *CT) {
110
111 enode *p;
112 int i, j;
113 int *ref;
114
115 if (value_notzero_p(e->d))
116 return; /* a rational number, its already reduced */
117 if(!(p = e->x.p))
118 return; /* hum... an overflow probably occured */
119
120 /* Compute ref */
121 ref = (int *)malloc(sizeof(int)*(CT->NbRows-1));
122 for(i=0;i<CT->NbRows-1;i++)
123 for(j=0;j<CT->NbColumns;j++)
124 if(value_notzero_p(CT->p[i][j])) {
125 ref[i] = j;
126 break;
127 }
128
129 /* Transform the references in e, using ref */
130 aep_evalue(e,ref);
131 free( ref );
132 return;
133 } /* addeliminatedparams_evalue */
134
135 static void addeliminatedparams_partition(enode *p, Matrix *CT, Polyhedron *CEq,
136 unsigned nparam, unsigned MaxRays)
137 {
138 int i;
139 assert(p->type == partition);
140 p->pos = nparam;
141
142 for (i = 0; i < p->size/2; i++) {
143 Polyhedron *D = EVALUE_DOMAIN(p->arr[2*i]);
144 Polyhedron *T = DomainPreimage(D, CT, MaxRays);
145 Domain_Free(D);
146 if (CEq) {
147 D = T;
148 T = DomainIntersection(D, CEq, MaxRays);
149 Domain_Free(D);
150 }
151 EVALUE_SET_DOMAIN(p->arr[2*i], T);
152 }
153 }
154
155 void addeliminatedparams_enum(evalue *e, Matrix *CT, Polyhedron *CEq,
156 unsigned MaxRays, unsigned nparam)
157 {
158 enode *p;
159 int i;
160
161 if (CT->NbRows == CT->NbColumns)
162 return;
163
164 if (EVALUE_IS_ZERO(*e))
165 return;
166
167 if (value_notzero_p(e->d)) {
168 evalue res;
169 value_init(res.d);
170 value_set_si(res.d, 0);
171 res.x.p = new_enode(partition, 2, nparam);
172 EVALUE_SET_DOMAIN(res.x.p->arr[0],
173 DomainConstraintSimplify(Polyhedron_Copy(CEq), MaxRays));
174 value_clear(res.x.p->arr[1].d);
175 res.x.p->arr[1] = *e;
176 *e = res;
177 return;
178 }
179
180 p = e->x.p;
181 assert(p);
182
183 addeliminatedparams_partition(p, CT, CEq, nparam, MaxRays);
184 for (i = 0; i < p->size/2; i++)
185 addeliminatedparams_evalue(&p->arr[2*i+1], CT);
186 }
187
188 static int mod_rational_cmp(evalue *e1, evalue *e2)
189 {
190 int r;
191 Value m;
192 Value m2;
193 value_init(m);
194 value_init(m2);
195
196 assert(value_notzero_p(e1->d));
197 assert(value_notzero_p(e2->d));
198 value_multiply(m, e1->x.n, e2->d);
199 value_multiply(m2, e2->x.n, e1->d);
200 if (value_lt(m, m2))
201 r = -1;
202 else if (value_gt(m, m2))
203 r = 1;
204 else
205 r = 0;
206 value_clear(m);
207 value_clear(m2);
208
209 return r;
210 }
211
212 static int mod_term_cmp_r(evalue *e1, evalue *e2)
213 {
214 if (value_notzero_p(e1->d)) {
215 int r;
216 if (value_zero_p(e2->d))
217 return -1;
218 return mod_rational_cmp(e1, e2);
219 }
220 if (value_notzero_p(e2->d))
221 return 1;
222 if (e1->x.p->pos < e2->x.p->pos)
223 return -1;
224 else if (e1->x.p->pos > e2->x.p->pos)
225 return 1;
226 else {
227 int r = mod_rational_cmp(&e1->x.p->arr[1], &e2->x.p->arr[1]);
228 return r == 0
229 ? mod_term_cmp_r(&e1->x.p->arr[0], &e2->x.p->arr[0])
230 : r;
231 }
232 }
233
234 static int mod_term_cmp(const evalue *e1, const evalue *e2)
235 {
236 assert(value_zero_p(e1->d));
237 assert(value_zero_p(e2->d));
238 assert(e1->x.p->type == fractional || e1->x.p->type == flooring);
239 assert(e2->x.p->type == fractional || e2->x.p->type == flooring);
240 return mod_term_cmp_r(&e1->x.p->arr[0], &e2->x.p->arr[0]);
241 }
242
243 static void check_order(const evalue *e)
244 {
245 int i;
246 evalue *a;
247
248 if (value_notzero_p(e->d))
249 return;
250
251 switch (e->x.p->type) {
252 case partition:
253 for (i = 0; i < e->x.p->size/2; ++i)
254 check_order(&e->x.p->arr[2*i+1]);
255 break;
256 case relation:
257 for (i = 1; 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 relation:
263 assert(mod_term_cmp(&e->x.p->arr[0], &a->x.p->arr[0]) < 0);
264 break;
265 case partition:
266 assert(0);
267 }
268 check_order(a);
269 }
270 break;
271 case polynomial:
272 for (i = 0; i < e->x.p->size; ++i) {
273 a = &e->x.p->arr[i];
274 if (value_notzero_p(a->d))
275 continue;
276 switch (a->x.p->type) {
277 case polynomial:
278 assert(e->x.p->pos < a->x.p->pos);
279 break;
280 case relation:
281 case partition:
282 assert(0);
283 }
284 check_order(a);
285 }
286 break;
287 case fractional:
288 case flooring:
289 for (i = 1; i < e->x.p->size; ++i) {
290 a = &e->x.p->arr[i];
291 if (value_notzero_p(a->d))
292 continue;
293 switch (a->x.p->type) {
294 case polynomial:
295 case relation:
296 case partition:
297 assert(0);
298 }
299 }
300 break;
301 }
302 }
303
304 /* Negative pos means inequality */
305 /* s is negative of substitution if m is not zero */
306 struct fixed_param {
307 int pos;
308 evalue s;
309 Value d;
310 Value m;
311 };
312
313 struct subst {
314 struct fixed_param *fixed;
315 int n;
316 int max;
317 };
318
319 static int relations_depth(evalue *e)
320 {
321 int d;
322
323 for (d = 0;
324 value_zero_p(e->d) && e->x.p->type == relation;
325 e = &e->x.p->arr[1], ++d);
326 return d;
327 }
328
329 static void poly_denom_not_constant(evalue **pp, Value *d)
330 {
331 evalue *p = *pp;
332 value_set_si(*d, 1);
333
334 while (value_zero_p(p->d)) {
335 assert(p->x.p->type == polynomial);
336 assert(p->x.p->size == 2);
337 assert(value_notzero_p(p->x.p->arr[1].d));
338 value_lcm(*d, *d, p->x.p->arr[1].d);
339 p = &p->x.p->arr[0];
340 }
341 *pp = p;
342 }
343
344 static void poly_denom(evalue *p, Value *d)
345 {
346 poly_denom_not_constant(&p, d);
347 value_lcm(*d, *d, p->d);
348 }
349
350 static void realloc_substitution(struct subst *s, int d)
351 {
352 struct fixed_param *n;
353 int i;
354 NALLOC(n, s->max+d);
355 for (i = 0; i < s->n; ++i)
356 n[i] = s->fixed[i];
357 free(s->fixed);
358 s->fixed = n;
359 s->max += d;
360 }
361
362 static int add_modulo_substitution(struct subst *s, evalue *r)
363 {
364 evalue *p;
365 evalue *f;
366 evalue *m;
367
368 assert(value_zero_p(r->d) && r->x.p->type == relation);
369 m = &r->x.p->arr[0];
370
371 /* May have been reduced already */
372 if (value_notzero_p(m->d))
373 return 0;
374
375 assert(value_zero_p(m->d) && m->x.p->type == fractional);
376 assert(m->x.p->size == 3);
377
378 /* fractional was inverted during reduction
379 * invert it back and move constant in
380 */
381 if (!EVALUE_IS_ONE(m->x.p->arr[2])) {
382 assert(value_pos_p(m->x.p->arr[2].d));
383 assert(value_mone_p(m->x.p->arr[2].x.n));
384 value_set_si(m->x.p->arr[2].x.n, 1);
385 value_increment(m->x.p->arr[1].x.n, m->x.p->arr[1].x.n);
386 assert(value_eq(m->x.p->arr[1].x.n, m->x.p->arr[1].d));
387 value_set_si(m->x.p->arr[1].x.n, 1);
388 eadd(&m->x.p->arr[1], &m->x.p->arr[0]);
389 value_set_si(m->x.p->arr[1].x.n, 0);
390 value_set_si(m->x.p->arr[1].d, 1);
391 }
392
393 /* Oops. Nested identical relations. */
394 if (!EVALUE_IS_ZERO(m->x.p->arr[1]))
395 return 0;
396
397 if (s->n >= s->max) {
398 int d = relations_depth(r);
399 realloc_substitution(s, d);
400 }
401
402 p = &m->x.p->arr[0];
403 assert(value_zero_p(p->d) && p->x.p->type == polynomial);
404 assert(p->x.p->size == 2);
405 f = &p->x.p->arr[1];
406
407 assert(value_pos_p(f->x.n));
408
409 value_init(s->fixed[s->n].m);
410 value_assign(s->fixed[s->n].m, f->d);
411 s->fixed[s->n].pos = p->x.p->pos;
412 value_init(s->fixed[s->n].d);
413 value_assign(s->fixed[s->n].d, f->x.n);
414 value_init(s->fixed[s->n].s.d);
415 evalue_copy(&s->fixed[s->n].s, &p->x.p->arr[0]);
416 ++s->n;
417
418 return 1;
419 }
420
421 static int type_offset(enode *p)
422 {
423 return p->type == fractional ? 1 :
424 p->type == flooring ? 1 :
425 p->type == relation ? 1 : 0;
426 }
427
428 static void reorder_terms_about(enode *p, evalue *v)
429 {
430 int i;
431 int offset = type_offset(p);
432
433 for (i = p->size-1; i >= offset+1; i--) {
434 emul(v, &p->arr[i]);
435 eadd(&p->arr[i], &p->arr[i-1]);
436 free_evalue_refs(&(p->arr[i]));
437 }
438 p->size = offset+1;
439 free_evalue_refs(v);
440 }
441
442 void evalue_reorder_terms(evalue *e)
443 {
444 enode *p;
445 evalue f;
446 int offset;
447
448 assert(value_zero_p(e->d));
449 p = e->x.p;
450 assert(p->type == fractional ||
451 p->type == flooring ||
452 p->type == polynomial); /* for now */
453
454 offset = type_offset(p);
455 value_init(f.d);
456 value_set_si(f.d, 0);
457 f.x.p = new_enode(p->type, offset+2, p->pos);
458 if (offset == 1) {
459 value_clear(f.x.p->arr[0].d);
460 f.x.p->arr[0] = p->arr[0];
461 }
462 evalue_set_si(&f.x.p->arr[offset], 0, 1);
463 evalue_set_si(&f.x.p->arr[offset+1], 1, 1);
464 reorder_terms_about(p, &f);
465 value_clear(e->d);
466 *e = p->arr[offset];
467 free(p);
468 }
469
470 static void evalue_reduce_size(evalue *e)
471 {
472 int i, offset;
473 enode *p;
474 assert(value_zero_p(e->d));
475
476 p = e->x.p;
477 offset = type_offset(p);
478
479 /* Try to reduce the degree */
480 for (i = p->size-1; i >= offset+1; i--) {
481 if (!EVALUE_IS_ZERO(p->arr[i]))
482 break;
483 free_evalue_refs(&p->arr[i]);
484 }
485 if (i+1 < p->size)
486 p->size = i+1;
487
488 /* Try to reduce its strength */
489 if (p->type == relation) {
490 if (p->size == 1) {
491 free_evalue_refs(&p->arr[0]);
492 evalue_set_si(e, 0, 1);
493 free(p);
494 }
495 } else if (p->size == offset+1) {
496 value_clear(e->d);
497 memcpy(e, &p->arr[offset], sizeof(evalue));
498 if (offset == 1)
499 free_evalue_refs(&p->arr[0]);
500 free(p);
501 }
502 }
503
504 void _reduce_evalue (evalue *e, struct subst *s, int fract) {
505
506 enode *p;
507 int i, j, k;
508 int add;
509
510 if (value_notzero_p(e->d)) {
511 if (fract)
512 mpz_fdiv_r(e->x.n, e->x.n, e->d);
513 return; /* a rational number, its already reduced */
514 }
515
516 if(!(p = e->x.p))
517 return; /* hum... an overflow probably occured */
518
519 /* First reduce the components of p */
520 add = p->type == relation;
521 for (i=0; i<p->size; i++) {
522 if (add && i == 1)
523 add = add_modulo_substitution(s, e);
524
525 if (i == 0 && p->type==fractional)
526 _reduce_evalue(&p->arr[i], s, 1);
527 else
528 _reduce_evalue(&p->arr[i], s, fract);
529
530 if (add && i == p->size-1) {
531 --s->n;
532 value_clear(s->fixed[s->n].m);
533 value_clear(s->fixed[s->n].d);
534 free_evalue_refs(&s->fixed[s->n].s);
535 } else if (add && i == 1)
536 s->fixed[s->n-1].pos *= -1;
537 }
538
539 if (p->type==periodic) {
540
541 /* Try to reduce the period */
542 for (i=1; i<=(p->size)/2; i++) {
543 if ((p->size % i)==0) {
544
545 /* Can we reduce the size to i ? */
546 for (j=0; j<i; j++)
547 for (k=j+i; k<e->x.p->size; k+=i)
548 if (!eequal(&p->arr[j], &p->arr[k])) goto you_lose;
549
550 /* OK, lets do it */
551 for (j=i; j<p->size; j++) free_evalue_refs(&p->arr[j]);
552 p->size = i;
553 break;
554
555 you_lose: /* OK, lets not do it */
556 continue;
557 }
558 }
559
560 /* Try to reduce its strength */
561 if (p->size == 1) {
562 value_clear(e->d);
563 memcpy(e,&p->arr[0],sizeof(evalue));
564 free(p);
565 }
566 }
567 else if (p->type==polynomial) {
568 for (k = 0; s && k < s->n; ++k) {
569 if (s->fixed[k].pos == p->pos) {
570 int divide = value_notone_p(s->fixed[k].d);
571 evalue d;
572
573 if (value_notzero_p(s->fixed[k].m)) {
574 if (!fract)
575 continue;
576 assert(p->size == 2);
577 if (divide && value_ne(s->fixed[k].d, p->arr[1].x.n))
578 continue;
579 if (!mpz_divisible_p(s->fixed[k].m, p->arr[1].d))
580 continue;
581 divide = 1;
582 }
583
584 if (divide) {
585 value_init(d.d);
586 value_assign(d.d, s->fixed[k].d);
587 value_init(d.x.n);
588 if (value_notzero_p(s->fixed[k].m))
589 value_oppose(d.x.n, s->fixed[k].m);
590 else
591 value_set_si(d.x.n, 1);
592 }
593
594 for (i=p->size-1;i>=1;i--) {
595 emul(&s->fixed[k].s, &p->arr[i]);
596 if (divide)
597 emul(&d, &p->arr[i]);
598 eadd(&p->arr[i], &p->arr[i-1]);
599 free_evalue_refs(&(p->arr[i]));
600 }
601 p->size = 1;
602 _reduce_evalue(&p->arr[0], s, fract);
603
604 if (divide)
605 free_evalue_refs(&d);
606
607 break;
608 }
609 }
610
611 evalue_reduce_size(e);
612 }
613 else if (p->type==fractional) {
614 int reorder = 0;
615 evalue v;
616
617 if (value_notzero_p(p->arr[0].d)) {
618 value_init(v.d);
619 value_assign(v.d, p->arr[0].d);
620 value_init(v.x.n);
621 mpz_fdiv_r(v.x.n, p->arr[0].x.n, p->arr[0].d);
622
623 reorder = 1;
624 } else {
625 evalue *f, *base;
626 evalue *pp = &p->arr[0];
627 assert(value_zero_p(pp->d) && pp->x.p->type == polynomial);
628 assert(pp->x.p->size == 2);
629
630 /* search for exact duplicate among the modulo inequalities */
631 do {
632 f = &pp->x.p->arr[1];
633 for (k = 0; s && k < s->n; ++k) {
634 if (-s->fixed[k].pos == pp->x.p->pos &&
635 value_eq(s->fixed[k].d, f->x.n) &&
636 value_eq(s->fixed[k].m, f->d) &&
637 eequal(&s->fixed[k].s, &pp->x.p->arr[0]))
638 break;
639 }
640 if (k < s->n) {
641 Value g;
642 value_init(g);
643
644 /* replace { E/m } by { (E-1)/m } + 1/m */
645 poly_denom(pp, &g);
646 if (reorder) {
647 evalue extra;
648 value_init(extra.d);
649 evalue_set_si(&extra, 1, 1);
650 value_assign(extra.d, g);
651 eadd(&extra, &v.x.p->arr[1]);
652 free_evalue_refs(&extra);
653
654 /* We've been going in circles; stop now */
655 if (value_ge(v.x.p->arr[1].x.n, v.x.p->arr[1].d)) {
656 free_evalue_refs(&v);
657 value_init(v.d);
658 evalue_set_si(&v, 0, 1);
659 break;
660 }
661 } else {
662 value_init(v.d);
663 value_set_si(v.d, 0);
664 v.x.p = new_enode(fractional, 3, -1);
665 evalue_set_si(&v.x.p->arr[1], 1, 1);
666 value_assign(v.x.p->arr[1].d, g);
667 evalue_set_si(&v.x.p->arr[2], 1, 1);
668 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
669 }
670
671 for (f = &v.x.p->arr[0]; value_zero_p(f->d);
672 f = &f->x.p->arr[0])
673 ;
674 value_division(f->d, g, f->d);
675 value_multiply(f->x.n, f->x.n, f->d);
676 value_assign(f->d, g);
677 value_decrement(f->x.n, f->x.n);
678 mpz_fdiv_r(f->x.n, f->x.n, f->d);
679
680 value_gcd(g, f->d, f->x.n);
681 value_division(f->d, f->d, g);
682 value_division(f->x.n, f->x.n, g);
683
684 value_clear(g);
685 pp = &v.x.p->arr[0];
686
687 reorder = 1;
688 }
689 } while (k < s->n);
690
691 /* reduction may have made this fractional arg smaller */
692 i = reorder ? p->size : 1;
693 for ( ; i < p->size; ++i)
694 if (value_zero_p(p->arr[i].d) &&
695 p->arr[i].x.p->type == fractional &&
696 mod_term_cmp(e, &p->arr[i]) >= 0)
697 break;
698 if (i < p->size) {
699 value_init(v.d);
700 value_set_si(v.d, 0);
701 v.x.p = new_enode(fractional, 3, -1);
702 evalue_set_si(&v.x.p->arr[1], 0, 1);
703 evalue_set_si(&v.x.p->arr[2], 1, 1);
704 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
705
706 reorder = 1;
707 }
708
709 if (!reorder) {
710 Value m;
711 Value r;
712 evalue *pp = &p->arr[0];
713 value_init(m);
714 value_init(r);
715 poly_denom_not_constant(&pp, &m);
716 mpz_fdiv_r(r, m, pp->d);
717 if (value_notzero_p(r)) {
718 value_init(v.d);
719 value_set_si(v.d, 0);
720 v.x.p = new_enode(fractional, 3, -1);
721
722 value_multiply(r, m, pp->x.n);
723 value_multiply(v.x.p->arr[1].d, m, pp->d);
724 value_init(v.x.p->arr[1].x.n);
725 mpz_fdiv_r(v.x.p->arr[1].x.n, r, pp->d);
726 mpz_fdiv_q(r, r, pp->d);
727
728 evalue_set_si(&v.x.p->arr[2], 1, 1);
729 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
730 pp = &v.x.p->arr[0];
731 while (value_zero_p(pp->d))
732 pp = &pp->x.p->arr[0];
733
734 value_assign(pp->d, m);
735 value_assign(pp->x.n, r);
736
737 value_gcd(r, pp->d, pp->x.n);
738 value_division(pp->d, pp->d, r);
739 value_division(pp->x.n, pp->x.n, r);
740
741 reorder = 1;
742 }
743 value_clear(m);
744 value_clear(r);
745 }
746
747 if (!reorder) {
748 int invert = 0;
749 Value twice;
750 value_init(twice);
751
752 for (pp = &p->arr[0]; value_zero_p(pp->d);
753 pp = &pp->x.p->arr[0]) {
754 f = &pp->x.p->arr[1];
755 assert(value_pos_p(f->d));
756 mpz_mul_ui(twice, f->x.n, 2);
757 if (value_lt(twice, f->d))
758 break;
759 if (value_eq(twice, f->d))
760 continue;
761 invert = 1;
762 break;
763 }
764
765 if (invert) {
766 value_init(v.d);
767 value_set_si(v.d, 0);
768 v.x.p = new_enode(fractional, 3, -1);
769 evalue_set_si(&v.x.p->arr[1], 0, 1);
770 poly_denom(&p->arr[0], &twice);
771 value_assign(v.x.p->arr[1].d, twice);
772 value_decrement(v.x.p->arr[1].x.n, twice);
773 evalue_set_si(&v.x.p->arr[2], -1, 1);
774 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
775
776 for (pp = &v.x.p->arr[0]; value_zero_p(pp->d);
777 pp = &pp->x.p->arr[0]) {
778 f = &pp->x.p->arr[1];
779 value_oppose(f->x.n, f->x.n);
780 mpz_fdiv_r(f->x.n, f->x.n, f->d);
781 }
782 value_division(pp->d, twice, pp->d);
783 value_multiply(pp->x.n, pp->x.n, pp->d);
784 value_assign(pp->d, twice);
785 value_oppose(pp->x.n, pp->x.n);
786 value_decrement(pp->x.n, pp->x.n);
787 mpz_fdiv_r(pp->x.n, pp->x.n, pp->d);
788
789 /* Maybe we should do this during reduction of
790 * the constant.
791 */
792 value_gcd(twice, pp->d, pp->x.n);
793 value_division(pp->d, pp->d, twice);
794 value_division(pp->x.n, pp->x.n, twice);
795
796 reorder = 1;
797 }
798
799 value_clear(twice);
800 }
801 }
802
803 if (reorder) {
804 reorder_terms_about(p, &v);
805 _reduce_evalue(&p->arr[1], s, fract);
806 }
807
808 evalue_reduce_size(e);
809 }
810 else if (p->type == flooring) {
811 /* Replace floor(constant) by its value */
812 if (value_notzero_p(p->arr[0].d)) {
813 evalue v;
814 value_init(v.d);
815 value_set_si(v.d, 1);
816 value_init(v.x.n);
817 mpz_fdiv_q(v.x.n, p->arr[0].x.n, p->arr[0].d);
818 reorder_terms_about(p, &v);
819 _reduce_evalue(&p->arr[1], s, fract);
820 }
821 evalue_reduce_size(e);
822 }
823 else if (p->type == relation) {
824 if (p->size == 3 && eequal(&p->arr[1], &p->arr[2])) {
825 free_evalue_refs(&(p->arr[2]));
826 free_evalue_refs(&(p->arr[0]));
827 p->size = 2;
828 value_clear(e->d);
829 *e = p->arr[1];
830 free(p);
831 return;
832 }
833 evalue_reduce_size(e);
834 if (value_notzero_p(e->d) || p != e->x.p)
835 return;
836 else {
837 int reduced = 0;
838 evalue *m;
839 m = &p->arr[0];
840
841 /* Relation was reduced by means of an identical
842 * inequality => remove
843 */
844 if (value_zero_p(m->d) && !EVALUE_IS_ZERO(m->x.p->arr[1]))
845 reduced = 1;
846
847 if (reduced || value_notzero_p(p->arr[0].d)) {
848 if (!reduced && value_zero_p(p->arr[0].x.n)) {
849 value_clear(e->d);
850 memcpy(e,&p->arr[1],sizeof(evalue));
851 if (p->size == 3)
852 free_evalue_refs(&(p->arr[2]));
853 } else {
854 if (p->size == 3) {
855 value_clear(e->d);
856 memcpy(e,&p->arr[2],sizeof(evalue));
857 } else
858 evalue_set_si(e, 0, 1);
859 free_evalue_refs(&(p->arr[1]));
860 }
861 free_evalue_refs(&(p->arr[0]));
862 free(p);
863 }
864 }
865 }
866 return;
867 } /* reduce_evalue */
868
869 static void add_substitution(struct subst *s, Value *row, unsigned dim)
870 {
871 int k, l;
872 evalue *r;
873
874 for (k = 0; k < dim; ++k)
875 if (value_notzero_p(row[k+1]))
876 break;
877
878 Vector_Normalize_Positive(row+1, dim+1, k);
879 assert(s->n < s->max);
880 value_init(s->fixed[s->n].d);
881 value_init(s->fixed[s->n].m);
882 value_assign(s->fixed[s->n].d, row[k+1]);
883 s->fixed[s->n].pos = k+1;
884 value_set_si(s->fixed[s->n].m, 0);
885 r = &s->fixed[s->n].s;
886 value_init(r->d);
887 for (l = k+1; l < dim; ++l)
888 if (value_notzero_p(row[l+1])) {
889 value_set_si(r->d, 0);
890 r->x.p = new_enode(polynomial, 2, l + 1);
891 value_init(r->x.p->arr[1].x.n);
892 value_oppose(r->x.p->arr[1].x.n, row[l+1]);
893 value_set_si(r->x.p->arr[1].d, 1);
894 r = &r->x.p->arr[0];
895 }
896 value_init(r->x.n);
897 value_oppose(r->x.n, row[dim+1]);
898 value_set_si(r->d, 1);
899 ++s->n;
900 }
901
902 static void _reduce_evalue_in_domain(evalue *e, Polyhedron *D, struct subst *s)
903 {
904 unsigned dim;
905 Polyhedron *orig = D;
906
907 s->n = 0;
908 dim = D->Dimension;
909 if (D->next)
910 D = DomainConvex(D, 0);
911 /* We don't perform any substitutions if the domain is a union.
912 * We may therefore miss out on some possible simplifications,
913 * e.g., if a variable is always even in the whole union,
914 * while there is a relation in the evalue that evaluates
915 * to zero for even values of the variable.
916 */
917 if (!D->next && D->NbEq) {
918 int j, k;
919 if (s->max < dim) {
920 if (s->max != 0)
921 realloc_substitution(s, dim);
922 else {
923 int d = relations_depth(e);
924 s->max = dim+d;
925 NALLOC(s->fixed, s->max);
926 }
927 }
928 for (j = 0; j < D->NbEq; ++j)
929 add_substitution(s, D->Constraint[j], dim);
930 }
931 if (D != orig)
932 Domain_Free(D);
933 _reduce_evalue(e, s, 0);
934 if (s->n != 0) {
935 int j;
936 for (j = 0; j < s->n; ++j) {
937 value_clear(s->fixed[j].d);
938 value_clear(s->fixed[j].m);
939 free_evalue_refs(&s->fixed[j].s);
940 }
941 }
942 }
943
944 void reduce_evalue_in_domain(evalue *e, Polyhedron *D)
945 {
946 struct subst s = { NULL, 0, 0 };
947 POL_ENSURE_VERTICES(D);
948 if (emptyQ(D)) {
949 if (EVALUE_IS_ZERO(*e))
950 return;
951 free_evalue_refs(e);
952 value_init(e->d);
953 evalue_set_si(e, 0, 1);
954 return;
955 }
956 _reduce_evalue_in_domain(e, D, &s);
957 if (s.max != 0)
958 free(s.fixed);
959 }
960
961 void reduce_evalue (evalue *e) {
962 struct subst s = { NULL, 0, 0 };
963
964 if (value_notzero_p(e->d))
965 return; /* a rational number, its already reduced */
966
967 if (e->x.p->type == partition) {
968 int i;
969 for (i = 0; i < e->x.p->size/2; ++i) {
970 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
971
972 /* This shouldn't really happen;
973 * Empty domains should not be added.
974 */
975 POL_ENSURE_VERTICES(D);
976 if (!emptyQ(D))
977 _reduce_evalue_in_domain(&e->x.p->arr[2*i+1], D, &s);
978
979 if (emptyQ(D) || EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
980 free_evalue_refs(&e->x.p->arr[2*i+1]);
981 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
982 value_clear(e->x.p->arr[2*i].d);
983 e->x.p->size -= 2;
984 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
985 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
986 --i;
987 }
988 }
989 if (e->x.p->size == 0) {
990 free(e->x.p);
991 evalue_set_si(e, 0, 1);
992 }
993 } else
994 _reduce_evalue(e, &s, 0);
995 if (s.max != 0)
996 free(s.fixed);
997 }
998
999 static void print_evalue_r(FILE *DST, const evalue *e, const char **pname)
1000 {
1001 if (EVALUE_IS_NAN(*e)) {
1002 fprintf(DST, "NaN");
1003 return;
1004 }
1005
1006 if(value_notzero_p(e->d)) {
1007 if(value_notone_p(e->d)) {
1008 value_print(DST,VALUE_FMT,e->x.n);
1009 fprintf(DST,"/");
1010 value_print(DST,VALUE_FMT,e->d);
1011 }
1012 else {
1013 value_print(DST,VALUE_FMT,e->x.n);
1014 }
1015 }
1016 else
1017 print_enode(DST,e->x.p,pname);
1018 return;
1019 } /* print_evalue */
1020
1021 void print_evalue(FILE *DST, const evalue *e, const char **pname)
1022 {
1023 print_evalue_r(DST, e, pname);
1024 if (value_notzero_p(e->d))
1025 fprintf(DST, "\n");
1026 }
1027
1028 void print_enode(FILE *DST, enode *p, const char **pname)
1029 {
1030 int i;
1031
1032 if (!p) {
1033 fprintf(DST, "NULL");
1034 return;
1035 }
1036 switch (p->type) {
1037 case evector:
1038 fprintf(DST, "{ ");
1039 for (i=0; i<p->size; i++) {
1040 print_evalue_r(DST, &p->arr[i], pname);
1041 if (i!=(p->size-1))
1042 fprintf(DST, ", ");
1043 }
1044 fprintf(DST, " }\n");
1045 break;
1046 case polynomial:
1047 fprintf(DST, "( ");
1048 for (i=p->size-1; i>=0; i--) {
1049 print_evalue_r(DST, &p->arr[i], pname);
1050 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
1051 else if (i>1)
1052 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
1053 }
1054 fprintf(DST, " )\n");
1055 break;
1056 case periodic:
1057 fprintf(DST, "[ ");
1058 for (i=0; i<p->size; i++) {
1059 print_evalue_r(DST, &p->arr[i], pname);
1060 if (i!=(p->size-1)) fprintf(DST, ", ");
1061 }
1062 fprintf(DST," ]_%s", pname[p->pos-1]);
1063 break;
1064 case flooring:
1065 case fractional:
1066 fprintf(DST, "( ");
1067 for (i=p->size-1; i>=1; i--) {
1068 print_evalue_r(DST, &p->arr[i], pname);
1069 if (i >= 2) {
1070 fprintf(DST, " * ");
1071 fprintf(DST, p->type == flooring ? "[" : "{");
1072 print_evalue_r(DST, &p->arr[0], pname);
1073 fprintf(DST, p->type == flooring ? "]" : "}");
1074 if (i>2)
1075 fprintf(DST, "^%d + ", i-1);
1076 else
1077 fprintf(DST, " + ");
1078 }
1079 }
1080 fprintf(DST, " )\n");
1081 break;
1082 case relation:
1083 fprintf(DST, "[ ");
1084 print_evalue_r(DST, &p->arr[0], pname);
1085 fprintf(DST, "= 0 ] * \n");
1086 print_evalue_r(DST, &p->arr[1], pname);
1087 if (p->size > 2) {
1088 fprintf(DST, " +\n [ ");
1089 print_evalue_r(DST, &p->arr[0], pname);
1090 fprintf(DST, "!= 0 ] * \n");
1091 print_evalue_r(DST, &p->arr[2], pname);
1092 }
1093 break;
1094 case partition: {
1095 char **new_names = NULL;
1096 const char **names = pname;
1097 int maxdim = EVALUE_DOMAIN(p->arr[0])->Dimension;
1098 if (!pname || p->pos < maxdim) {
1099 new_names = ALLOCN(char *, maxdim);
1100 for (i = 0; i < p->pos; ++i) {
1101 if (pname)
1102 new_names[i] = (char *)pname[i];
1103 else {
1104 new_names[i] = ALLOCN(char, 10);
1105 snprintf(new_names[i], 10, "%c", 'P'+i);
1106 }
1107 }
1108 for ( ; i < maxdim; ++i) {
1109 new_names[i] = ALLOCN(char, 10);
1110 snprintf(new_names[i], 10, "_p%d", i);
1111 }
1112 names = (const char**)new_names;
1113 }
1114
1115 for (i=0; i<p->size/2; i++) {
1116 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), names);
1117 print_evalue_r(DST, &p->arr[2*i+1], names);
1118 if (value_notzero_p(p->arr[2*i+1].d))
1119 fprintf(DST, "\n");
1120 }
1121
1122 if (!pname || p->pos < maxdim) {
1123 for (i = pname ? p->pos : 0; i < maxdim; ++i)
1124 free(new_names[i]);
1125 free(new_names);
1126 }
1127
1128 break;
1129 }
1130 default:
1131 assert(0);
1132 }
1133 return;
1134 } /* print_enode */
1135
1136 /* Returns
1137 * 0 if toplevels of e1 and e2 are at the same level
1138 * <0 if toplevel of e1 should be outside of toplevel of e2
1139 * >0 if toplevel of e2 should be outside of toplevel of e1
1140 */
1141 static int evalue_level_cmp(const evalue *e1, const evalue *e2)
1142 {
1143 if (value_notzero_p(e1->d) && value_notzero_p(e2->d))
1144 return 0;
1145 if (value_notzero_p(e1->d))
1146 return 1;
1147 if (value_notzero_p(e2->d))
1148 return -1;
1149 if (e1->x.p->type == partition && e2->x.p->type == partition)
1150 return 0;
1151 if (e1->x.p->type == partition)
1152 return -1;
1153 if (e2->x.p->type == partition)
1154 return 1;
1155 if (e1->x.p->type == relation && e2->x.p->type == relation) {
1156 if (eequal(&e1->x.p->arr[0], &e2->x.p->arr[0]))
1157 return 0;
1158 return mod_term_cmp(&e1->x.p->arr[0], &e2->x.p->arr[0]);
1159 }
1160 if (e1->x.p->type == relation)
1161 return -1;
1162 if (e2->x.p->type == relation)
1163 return 1;
1164 if (e1->x.p->type == polynomial && e2->x.p->type == polynomial)
1165 return e1->x.p->pos - e2->x.p->pos;
1166 if (e1->x.p->type == polynomial)
1167 return -1;
1168 if (e2->x.p->type == polynomial)
1169 return 1;
1170 if (e1->x.p->type == periodic && e2->x.p->type == periodic)
1171 return e1->x.p->pos - e2->x.p->pos;
1172 assert(e1->x.p->type != periodic);
1173 assert(e2->x.p->type != periodic);
1174 assert(e1->x.p->type == e2->x.p->type);
1175 if (eequal(&e1->x.p->arr[0], &e2->x.p->arr[0]))
1176 return 0;
1177 return mod_term_cmp(e1, e2);
1178 }
1179
1180 static void eadd_rev(const evalue *e1, evalue *res)
1181 {
1182 evalue ev;
1183 value_init(ev.d);
1184 evalue_copy(&ev, e1);
1185 eadd(res, &ev);
1186 free_evalue_refs(res);
1187 *res = ev;
1188 }
1189
1190 static void eadd_rev_cst(const evalue *e1, evalue *res)
1191 {
1192 evalue ev;
1193 value_init(ev.d);
1194 evalue_copy(&ev, e1);
1195 eadd(res, &ev.x.p->arr[type_offset(ev.x.p)]);
1196 free_evalue_refs(res);
1197 *res = ev;
1198 }
1199
1200 struct section { Polyhedron * D; evalue E; };
1201
1202 void eadd_partitions(const evalue *e1, evalue *res)
1203 {
1204 int n, i, j;
1205 Polyhedron *d, *fd;
1206 struct section *s;
1207 s = (struct section *)
1208 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
1209 sizeof(struct section));
1210 assert(s);
1211 assert(e1->x.p->pos == res->x.p->pos);
1212 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
1213 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1214
1215 n = 0;
1216 for (j = 0; j < e1->x.p->size/2; ++j) {
1217 assert(res->x.p->size >= 2);
1218 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1219 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
1220 if (!emptyQ(fd))
1221 for (i = 1; i < res->x.p->size/2; ++i) {
1222 Polyhedron *t = fd;
1223 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1224 Domain_Free(t);
1225 if (emptyQ(fd))
1226 break;
1227 }
1228 fd = DomainConstraintSimplify(fd, 0);
1229 if (emptyQ(fd)) {
1230 Domain_Free(fd);
1231 continue;
1232 }
1233 value_init(s[n].E.d);
1234 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
1235 s[n].D = fd;
1236 ++n;
1237 }
1238 for (i = 0; i < res->x.p->size/2; ++i) {
1239 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
1240 for (j = 0; j < e1->x.p->size/2; ++j) {
1241 Polyhedron *t;
1242 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1243 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1244 d = DomainConstraintSimplify(d, 0);
1245 if (emptyQ(d)) {
1246 Domain_Free(d);
1247 continue;
1248 }
1249 t = fd;
1250 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
1251 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
1252 Domain_Free(t);
1253 value_init(s[n].E.d);
1254 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1255 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
1256 s[n].D = d;
1257 ++n;
1258 }
1259 if (!emptyQ(fd)) {
1260 s[n].E = res->x.p->arr[2*i+1];
1261 s[n].D = fd;
1262 ++n;
1263 } else {
1264 free_evalue_refs(&res->x.p->arr[2*i+1]);
1265 Domain_Free(fd);
1266 }
1267 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
1268 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1269 value_clear(res->x.p->arr[2*i].d);
1270 }
1271
1272 free(res->x.p);
1273 assert(n > 0);
1274 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1275 for (j = 0; j < n; ++j) {
1276 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1277 value_clear(res->x.p->arr[2*j+1].d);
1278 res->x.p->arr[2*j+1] = s[j].E;
1279 }
1280
1281 free(s);
1282 }
1283
1284 static void explicit_complement(evalue *res)
1285 {
1286 enode *rel = new_enode(relation, 3, 0);
1287 assert(rel);
1288 value_clear(rel->arr[0].d);
1289 rel->arr[0] = res->x.p->arr[0];
1290 value_clear(rel->arr[1].d);
1291 rel->arr[1] = res->x.p->arr[1];
1292 value_set_si(rel->arr[2].d, 1);
1293 value_init(rel->arr[2].x.n);
1294 value_set_si(rel->arr[2].x.n, 0);
1295 free(res->x.p);
1296 res->x.p = rel;
1297 }
1298
1299 static void reduce_constant(evalue *e)
1300 {
1301 Value g;
1302 value_init(g);
1303
1304 value_gcd(g, e->x.n, e->d);
1305 if (value_notone_p(g)) {
1306 value_division(e->d, e->d,g);
1307 value_division(e->x.n, e->x.n,g);
1308 }
1309 value_clear(g);
1310 }
1311
1312 /* Add two rational numbers */
1313 static void eadd_rationals(const evalue *e1, evalue *res)
1314 {
1315 if (value_eq(e1->d, res->d))
1316 value_addto(res->x.n, res->x.n, e1->x.n);
1317 else {
1318 value_multiply(res->x.n, res->x.n, e1->d);
1319 value_addmul(res->x.n, e1->x.n, res->d);
1320 value_multiply(res->d,e1->d,res->d);
1321 }
1322 reduce_constant(res);
1323 }
1324
1325 static void eadd_relations(const evalue *e1, evalue *res)
1326 {
1327 int i;
1328
1329 if (res->x.p->size < 3 && e1->x.p->size == 3)
1330 explicit_complement(res);
1331 for (i = 1; i < e1->x.p->size; ++i)
1332 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1333 }
1334
1335 static void eadd_arrays(const evalue *e1, evalue *res, int n)
1336 {
1337 int i;
1338
1339 // add any element in e1 to the corresponding element in res
1340 i = type_offset(res->x.p);
1341 if (i == 1)
1342 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1343 for (; i < n; i++)
1344 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1345 }
1346
1347 static void eadd_poly(const evalue *e1, evalue *res)
1348 {
1349 if (e1->x.p->size > res->x.p->size)
1350 eadd_rev(e1, res);
1351 else
1352 eadd_arrays(e1, res, e1->x.p->size);
1353 }
1354
1355 /*
1356 * Product or sum of two periodics of the same parameter
1357 * and different periods
1358 */
1359 static void combine_periodics(const evalue *e1, evalue *res,
1360 void (*op)(const evalue *, evalue*))
1361 {
1362 Value es, rs;
1363 int i, size;
1364 enode *p;
1365
1366 value_init(es);
1367 value_init(rs);
1368 value_set_si(es, e1->x.p->size);
1369 value_set_si(rs, res->x.p->size);
1370 value_lcm(rs, es, rs);
1371 size = (int)mpz_get_si(rs);
1372 value_clear(es);
1373 value_clear(rs);
1374 p = new_enode(periodic, size, e1->x.p->pos);
1375 for (i = 0; i < res->x.p->size; i++) {
1376 value_clear(p->arr[i].d);
1377 p->arr[i] = res->x.p->arr[i];
1378 }
1379 for (i = res->x.p->size; i < size; i++)
1380 evalue_copy(&p->arr[i], &res->x.p->arr[i % res->x.p->size]);
1381 for (i = 0; i < size; i++)
1382 op(&e1->x.p->arr[i % e1->x.p->size], &p->arr[i]);
1383 free(res->x.p);
1384 res->x.p = p;
1385 }
1386
1387 static void eadd_periodics(const evalue *e1, evalue *res)
1388 {
1389 int i;
1390 int x, y, p;
1391 evalue *ne;
1392
1393 if (e1->x.p->size == res->x.p->size) {
1394 eadd_arrays(e1, res, e1->x.p->size);
1395 return;
1396 }
1397
1398 combine_periodics(e1, res, eadd);
1399 }
1400
1401 void evalue_assign(evalue *dst, const evalue *src)
1402 {
1403 if (value_pos_p(dst->d) && value_pos_p(src->d)) {
1404 value_assign(dst->d, src->d);
1405 value_assign(dst->x.n, src->x.n);
1406 return;
1407 }
1408 free_evalue_refs(dst);
1409 value_init(dst->d);
1410 evalue_copy(dst, src);
1411 }
1412
1413 void eadd(const evalue *e1, evalue *res)
1414 {
1415 int cmp;
1416
1417 if (EVALUE_IS_ZERO(*e1))
1418 return;
1419
1420 if (EVALUE_IS_NAN(*res))
1421 return;
1422
1423 if (EVALUE_IS_NAN(*e1)) {
1424 evalue_assign(res, e1);
1425 return;
1426 }
1427
1428 if (EVALUE_IS_ZERO(*res)) {
1429 evalue_assign(res, e1);
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 combine_periodics(e1, res, emul);
1671 }
1672
1673 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1674
1675 static void emul_fractionals(const evalue *e1, evalue *res)
1676 {
1677 evalue d;
1678 value_init(d.d);
1679 poly_denom(&e1->x.p->arr[0], &d.d);
1680 if (!value_two_p(d.d))
1681 emul_poly(e1, res);
1682 else {
1683 evalue tmp;
1684 value_init(d.x.n);
1685 value_set_si(d.x.n, 1);
1686 /* { x }^2 == { x }/2 */
1687 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1688 assert(e1->x.p->size == 3);
1689 assert(res->x.p->size == 3);
1690 value_init(tmp.d);
1691 evalue_copy(&tmp, &res->x.p->arr[2]);
1692 emul(&d, &tmp);
1693 eadd(&res->x.p->arr[1], &tmp);
1694 emul(&e1->x.p->arr[2], &tmp);
1695 emul(&e1->x.p->arr[1], &res->x.p->arr[1]);
1696 emul(&e1->x.p->arr[1], &res->x.p->arr[2]);
1697 eadd(&tmp, &res->x.p->arr[2]);
1698 free_evalue_refs(&tmp);
1699 value_clear(d.x.n);
1700 }
1701 value_clear(d.d);
1702 }
1703
1704 /* Computes the product of two evalues "e1" and "res" and puts
1705 * the result in "res". You need to make a copy of "res"
1706 * before calling this function if you still need it afterward.
1707 * The vector type of evalues is not treated here
1708 */
1709 void emul(const evalue *e1, evalue *res)
1710 {
1711 int cmp;
1712
1713 assert(!(value_zero_p(e1->d) && e1->x.p->type == evector));
1714 assert(!(value_zero_p(res->d) && res->x.p->type == evector));
1715
1716 if (EVALUE_IS_ZERO(*res))
1717 return;
1718
1719 if (EVALUE_IS_ONE(*e1))
1720 return;
1721
1722 if (EVALUE_IS_ZERO(*e1)) {
1723 evalue_assign(res, e1);
1724 return;
1725 }
1726
1727 if (EVALUE_IS_NAN(*res))
1728 return;
1729
1730 if (EVALUE_IS_NAN(*e1)) {
1731 evalue_assign(res, e1);
1732 return;
1733 }
1734
1735 cmp = evalue_level_cmp(res, e1);
1736 if (cmp > 0) {
1737 emul_rev(e1, res);
1738 } else if (cmp == 0) {
1739 if (value_notzero_p(e1->d)) {
1740 emul_rationals(e1, res);
1741 } else {
1742 switch (e1->x.p->type) {
1743 case partition:
1744 emul_partitions(e1, res);
1745 break;
1746 case relation:
1747 emul_relations(e1, res);
1748 break;
1749 case polynomial:
1750 case flooring:
1751 emul_poly(e1, res);
1752 break;
1753 case periodic:
1754 emul_periodics(e1, res);
1755 break;
1756 case fractional:
1757 emul_fractionals(e1, res);
1758 break;
1759 }
1760 }
1761 } else {
1762 int i;
1763 switch (res->x.p->type) {
1764 case partition:
1765 for (i = 0; i < res->x.p->size/2; ++i)
1766 emul(e1, &res->x.p->arr[2*i+1]);
1767 break;
1768 case relation:
1769 case polynomial:
1770 case periodic:
1771 case flooring:
1772 case fractional:
1773 for (i = type_offset(res->x.p); i < res->x.p->size; ++i)
1774 emul(e1, &res->x.p->arr[i]);
1775 break;
1776 }
1777 }
1778 }
1779
1780 /* Frees mask content ! */
1781 void emask(evalue *mask, evalue *res) {
1782 int n, i, j;
1783 Polyhedron *d, *fd;
1784 struct section *s;
1785 evalue mone;
1786 int pos;
1787
1788 if (EVALUE_IS_ZERO(*res)) {
1789 free_evalue_refs(mask);
1790 return;
1791 }
1792
1793 assert(value_zero_p(mask->d));
1794 assert(mask->x.p->type == partition);
1795 assert(value_zero_p(res->d));
1796 assert(res->x.p->type == partition);
1797 assert(mask->x.p->pos == res->x.p->pos);
1798 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1799 assert(mask->x.p->pos == EVALUE_DOMAIN(mask->x.p->arr[0])->Dimension);
1800 pos = res->x.p->pos;
1801
1802 s = (struct section *)
1803 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1804 sizeof(struct section));
1805 assert(s);
1806
1807 value_init(mone.d);
1808 evalue_set_si(&mone, -1, 1);
1809
1810 n = 0;
1811 for (j = 0; j < res->x.p->size/2; ++j) {
1812 assert(mask->x.p->size >= 2);
1813 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1814 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1815 if (!emptyQ(fd))
1816 for (i = 1; i < mask->x.p->size/2; ++i) {
1817 Polyhedron *t = fd;
1818 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1819 Domain_Free(t);
1820 if (emptyQ(fd))
1821 break;
1822 }
1823 if (emptyQ(fd)) {
1824 Domain_Free(fd);
1825 continue;
1826 }
1827 value_init(s[n].E.d);
1828 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1829 s[n].D = fd;
1830 ++n;
1831 }
1832 for (i = 0; i < mask->x.p->size/2; ++i) {
1833 if (EVALUE_IS_ONE(mask->x.p->arr[2*i+1]))
1834 continue;
1835
1836 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1837 eadd(&mone, &mask->x.p->arr[2*i+1]);
1838 emul(&mone, &mask->x.p->arr[2*i+1]);
1839 for (j = 0; j < res->x.p->size/2; ++j) {
1840 Polyhedron *t;
1841 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1842 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1843 if (emptyQ(d)) {
1844 Domain_Free(d);
1845 continue;
1846 }
1847 t = fd;
1848 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1849 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1850 Domain_Free(t);
1851 value_init(s[n].E.d);
1852 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1853 emul(&mask->x.p->arr[2*i+1], &s[n].E);
1854 s[n].D = d;
1855 ++n;
1856 }
1857
1858 if (!emptyQ(fd)) {
1859 /* Just ignore; this may have been previously masked off */
1860 }
1861 if (fd != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1862 Domain_Free(fd);
1863 }
1864
1865 free_evalue_refs(&mone);
1866 free_evalue_refs(mask);
1867 free_evalue_refs(res);
1868 value_init(res->d);
1869 if (n == 0)
1870 evalue_set_si(res, 0, 1);
1871 else {
1872 res->x.p = new_enode(partition, 2*n, pos);
1873 for (j = 0; j < n; ++j) {
1874 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1875 value_clear(res->x.p->arr[2*j+1].d);
1876 res->x.p->arr[2*j+1] = s[j].E;
1877 }
1878 }
1879
1880 free(s);
1881 }
1882
1883 void evalue_copy(evalue *dst, const evalue *src)
1884 {
1885 value_assign(dst->d, src->d);
1886 if (EVALUE_IS_NAN(*dst)) {
1887 dst->x.p = NULL;
1888 return;
1889 }
1890 if (value_pos_p(src->d)) {
1891 value_init(dst->x.n);
1892 value_assign(dst->x.n, src->x.n);
1893 } else
1894 dst->x.p = ecopy(src->x.p);
1895 }
1896
1897 evalue *evalue_dup(const evalue *e)
1898 {
1899 evalue *res = ALLOC(evalue);
1900 value_init(res->d);
1901 evalue_copy(res, e);
1902 return res;
1903 }
1904
1905 enode *new_enode(enode_type type,int size,int pos) {
1906
1907 enode *res;
1908 int i;
1909
1910 if(size == 0) {
1911 fprintf(stderr, "Allocating enode of size 0 !\n" );
1912 return NULL;
1913 }
1914 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1915 res->type = type;
1916 res->size = size;
1917 res->pos = pos;
1918 for(i=0; i<size; i++) {
1919 value_init(res->arr[i].d);
1920 value_set_si(res->arr[i].d,0);
1921 res->arr[i].x.p = 0;
1922 }
1923 return res;
1924 } /* new_enode */
1925
1926 enode *ecopy(enode *e) {
1927
1928 enode *res;
1929 int i;
1930
1931 res = new_enode(e->type,e->size,e->pos);
1932 for(i=0;i<e->size;++i) {
1933 value_assign(res->arr[i].d,e->arr[i].d);
1934 if(value_zero_p(res->arr[i].d))
1935 res->arr[i].x.p = ecopy(e->arr[i].x.p);
1936 else if (EVALUE_IS_DOMAIN(res->arr[i]))
1937 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1938 else {
1939 value_init(res->arr[i].x.n);
1940 value_assign(res->arr[i].x.n,e->arr[i].x.n);
1941 }
1942 }
1943 return(res);
1944 } /* ecopy */
1945
1946 int ecmp(const evalue *e1, const evalue *e2)
1947 {
1948 enode *p1, *p2;
1949 int i;
1950 int r;
1951
1952 if (value_notzero_p(e1->d) && value_notzero_p(e2->d)) {
1953 Value m, m2;
1954 value_init(m);
1955 value_init(m2);
1956 value_multiply(m, e1->x.n, e2->d);
1957 value_multiply(m2, e2->x.n, e1->d);
1958
1959 if (value_lt(m, m2))
1960 r = -1;
1961 else if (value_gt(m, m2))
1962 r = 1;
1963 else
1964 r = 0;
1965
1966 value_clear(m);
1967 value_clear(m2);
1968
1969 return r;
1970 }
1971 if (value_notzero_p(e1->d))
1972 return -1;
1973 if (value_notzero_p(e2->d))
1974 return 1;
1975
1976 p1 = e1->x.p;
1977 p2 = e2->x.p;
1978
1979 if (p1->type != p2->type)
1980 return p1->type - p2->type;
1981 if (p1->pos != p2->pos)
1982 return p1->pos - p2->pos;
1983 if (p1->size != p2->size)
1984 return p1->size - p2->size;
1985
1986 for (i = p1->size-1; i >= 0; --i)
1987 if ((r = ecmp(&p1->arr[i], &p2->arr[i])) != 0)
1988 return r;
1989
1990 return 0;
1991 }
1992
1993 int eequal(const evalue *e1, const evalue *e2)
1994 {
1995 int i;
1996 enode *p1, *p2;
1997
1998 if (value_ne(e1->d,e2->d))
1999 return 0;
2000
2001 if (EVALUE_IS_DOMAIN(*e1))
2002 return PolyhedronIncludes(EVALUE_DOMAIN(*e2), EVALUE_DOMAIN(*e1)) &&
2003 PolyhedronIncludes(EVALUE_DOMAIN(*e1), EVALUE_DOMAIN(*e2));
2004
2005 if (EVALUE_IS_NAN(*e1))
2006 return 1;
2007
2008 assert(value_posz_p(e1->d));
2009
2010 /* e1->d == e2->d */
2011 if (value_notzero_p(e1->d)) {
2012 if (value_ne(e1->x.n,e2->x.n))
2013 return 0;
2014
2015 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2016 return 1;
2017 }
2018
2019 /* e1->d == e2->d == 0 */
2020 p1 = e1->x.p;
2021 p2 = e2->x.p;
2022 if (p1->type != p2->type) return 0;
2023 if (p1->size != p2->size) return 0;
2024 if (p1->pos != p2->pos) return 0;
2025 for (i=0; i<p1->size; i++)
2026 if (!eequal(&p1->arr[i], &p2->arr[i]) )
2027 return 0;
2028 return 1;
2029 } /* eequal */
2030
2031 void free_evalue_refs(evalue *e) {
2032
2033 enode *p;
2034 int i;
2035
2036 if (EVALUE_IS_NAN(*e)) {
2037 value_clear(e->d);
2038 return;
2039 }
2040
2041 if (EVALUE_IS_DOMAIN(*e)) {
2042 Domain_Free(EVALUE_DOMAIN(*e));
2043 value_clear(e->d);
2044 return;
2045 } else if (value_pos_p(e->d)) {
2046
2047 /* 'e' stores a constant */
2048 value_clear(e->d);
2049 value_clear(e->x.n);
2050 return;
2051 }
2052 assert(value_zero_p(e->d));
2053 value_clear(e->d);
2054 p = e->x.p;
2055 if (!p) return; /* null pointer */
2056 for (i=0; i<p->size; i++) {
2057 free_evalue_refs(&(p->arr[i]));
2058 }
2059 free(p);
2060 return;
2061 } /* free_evalue_refs */
2062
2063 void evalue_free(evalue *e)
2064 {
2065 free_evalue_refs(e);
2066 free(e);
2067 }
2068
2069 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
2070 Vector * val, evalue *res)
2071 {
2072 unsigned nparam = periods->Size;
2073
2074 if (p == nparam) {
2075 double d = compute_evalue(e, val->p);
2076 d *= VALUE_TO_DOUBLE(m);
2077 if (d > 0)
2078 d += .25;
2079 else
2080 d -= .25;
2081 value_assign(res->d, m);
2082 value_init(res->x.n);
2083 value_set_double(res->x.n, d);
2084 mpz_fdiv_r(res->x.n, res->x.n, m);
2085 return;
2086 }
2087 if (value_one_p(periods->p[p]))
2088 mod2table_r(e, periods, m, p+1, val, res);
2089 else {
2090 Value tmp;
2091 value_init(tmp);
2092
2093 value_assign(tmp, periods->p[p]);
2094 value_set_si(res->d, 0);
2095 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
2096 do {
2097 value_decrement(tmp, tmp);
2098 value_assign(val->p[p], tmp);
2099 mod2table_r(e, periods, m, p+1, val,
2100 &res->x.p->arr[VALUE_TO_INT(tmp)]);
2101 } while (value_pos_p(tmp));
2102
2103 value_clear(tmp);
2104 }
2105 }
2106
2107 static void rel2table(evalue *e, int zero)
2108 {
2109 if (value_pos_p(e->d)) {
2110 if (value_zero_p(e->x.n) == zero)
2111 value_set_si(e->x.n, 1);
2112 else
2113 value_set_si(e->x.n, 0);
2114 value_set_si(e->d, 1);
2115 } else {
2116 int i;
2117 for (i = 0; i < e->x.p->size; ++i)
2118 rel2table(&e->x.p->arr[i], zero);
2119 }
2120 }
2121
2122 void evalue_mod2table(evalue *e, int nparam)
2123 {
2124 enode *p;
2125 int i;
2126
2127 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2128 return;
2129 p = e->x.p;
2130 for (i=0; i<p->size; i++) {
2131 evalue_mod2table(&(p->arr[i]), nparam);
2132 }
2133 if (p->type == relation) {
2134 evalue copy;
2135
2136 if (p->size > 2) {
2137 value_init(copy.d);
2138 evalue_copy(&copy, &p->arr[0]);
2139 }
2140 rel2table(&p->arr[0], 1);
2141 emul(&p->arr[0], &p->arr[1]);
2142 if (p->size > 2) {
2143 rel2table(&copy, 0);
2144 emul(&copy, &p->arr[2]);
2145 eadd(&p->arr[2], &p->arr[1]);
2146 free_evalue_refs(&p->arr[2]);
2147 free_evalue_refs(&copy);
2148 }
2149 free_evalue_refs(&p->arr[0]);
2150 value_clear(e->d);
2151 *e = p->arr[1];
2152 free(p);
2153 } else if (p->type == fractional) {
2154 Vector *periods = Vector_Alloc(nparam);
2155 Vector *val = Vector_Alloc(nparam);
2156 Value tmp;
2157 evalue *ev;
2158 evalue EP, res;
2159
2160 value_init(tmp);
2161 value_set_si(tmp, 1);
2162 Vector_Set(periods->p, 1, nparam);
2163 Vector_Set(val->p, 0, nparam);
2164 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
2165 enode *p = ev->x.p;
2166
2167 assert(p->type == polynomial);
2168 assert(p->size == 2);
2169 value_assign(periods->p[p->pos-1], p->arr[1].d);
2170 value_lcm(tmp, tmp, p->arr[1].d);
2171 }
2172 value_lcm(tmp, tmp, ev->d);
2173 value_init(EP.d);
2174 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
2175
2176 value_init(res.d);
2177 evalue_set_si(&res, 0, 1);
2178 /* Compute the polynomial using Horner's rule */
2179 for (i=p->size-1;i>1;i--) {
2180 eadd(&p->arr[i], &res);
2181 emul(&EP, &res);
2182 }
2183 eadd(&p->arr[1], &res);
2184
2185 free_evalue_refs(e);
2186 free_evalue_refs(&EP);
2187 *e = res;
2188
2189 value_clear(tmp);
2190 Vector_Free(val);
2191 Vector_Free(periods);
2192 }
2193 } /* evalue_mod2table */
2194
2195 /********************************************************/
2196 /* function in domain */
2197 /* check if the parameters in list_args */
2198 /* verifies the constraints of Domain P */
2199 /********************************************************/
2200 int in_domain(Polyhedron *P, Value *list_args)
2201 {
2202 int row, in = 1;
2203 Value v; /* value of the constraint of a row when
2204 parameters are instantiated*/
2205
2206 if (P->Dimension == 0)
2207 return !emptyQ(P);
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 evalue_reorder_terms(&rem);
3012 evalue_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 evalue_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 assert(value_zero_p(EP->d));
3076 assert(EP->x.p->type == partition);
3077 assert(EP->x.p->size >= 2);
3078
3079 for (i = 0; i < EP->x.p->size/2; ++i) {
3080 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[2*i])->Dimension);
3081 en = (Enumeration *)malloc(sizeof(Enumeration));
3082 en->next = res;
3083 res = en;
3084 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
3085 value_clear(EP->x.p->arr[2*i].d);
3086 res->EP = EP->x.p->arr[2*i+1];
3087 }
3088 free(EP->x.p);
3089 value_clear(EP->d);
3090 free(EP);
3091 return res;
3092 }
3093
3094 int evalue_frac2floor_in_domain3(evalue *e, Polyhedron *D, int shift)
3095 {
3096 enode *p;
3097 int r = 0;
3098 int i;
3099 Polyhedron *I;
3100 Value d, min;
3101 evalue fl;
3102
3103 if (value_notzero_p(e->d))
3104 return r;
3105
3106 p = e->x.p;
3107
3108 i = p->type == relation ? 1 :
3109 p->type == fractional ? 1 : 0;
3110 for (; i<p->size; i++)
3111 r |= evalue_frac2floor_in_domain3(&p->arr[i], D, shift);
3112
3113 if (p->type != fractional) {
3114 if (r && p->type == polynomial) {
3115 evalue f;
3116 value_init(f.d);
3117 value_set_si(f.d, 0);
3118 f.x.p = new_enode(polynomial, 2, p->pos);
3119 evalue_set_si(&f.x.p->arr[0], 0, 1);
3120 evalue_set_si(&f.x.p->arr[1], 1, 1);
3121 reorder_terms_about(p, &f);
3122 value_clear(e->d);
3123 *e = p->arr[0];
3124 free(p);
3125 }
3126 return r;
3127 }
3128
3129 if (shift) {
3130 value_init(d);
3131 I = polynomial_projection(p, D, &d, NULL);
3132
3133 /*
3134 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3135 */
3136
3137 assert(I->NbEq == 0); /* Should have been reduced */
3138
3139 /* Find minimum */
3140 for (i = 0; i < I->NbConstraints; ++i)
3141 if (value_pos_p(I->Constraint[i][1]))
3142 break;
3143
3144 if (i < I->NbConstraints) {
3145 value_init(min);
3146 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
3147 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
3148 if (value_neg_p(min)) {
3149 evalue offset;
3150 mpz_fdiv_q(min, min, d);
3151 value_init(offset.d);
3152 value_set_si(offset.d, 1);
3153 value_init(offset.x.n);
3154 value_oppose(offset.x.n, min);
3155 eadd(&offset, &p->arr[0]);
3156 free_evalue_refs(&offset);
3157 }
3158 value_clear(min);
3159 }
3160
3161 Polyhedron_Free(I);
3162 value_clear(d);
3163 }
3164
3165 value_init(fl.d);
3166 value_set_si(fl.d, 0);
3167 fl.x.p = new_enode(flooring, 3, -1);
3168 evalue_set_si(&fl.x.p->arr[1], 0, 1);
3169 evalue_set_si(&fl.x.p->arr[2], -1, 1);
3170 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
3171
3172 eadd(&fl, &p->arr[0]);
3173 reorder_terms_about(p, &p->arr[0]);
3174 value_clear(e->d);
3175 *e = p->arr[1];
3176 free(p);
3177 free_evalue_refs(&fl);
3178
3179 return 1;
3180 }
3181
3182 int evalue_frac2floor_in_domain(evalue *e, Polyhedron *D)
3183 {
3184 return evalue_frac2floor_in_domain3(e, D, 1);
3185 }
3186
3187 void evalue_frac2floor2(evalue *e, int shift)
3188 {
3189 int i;
3190 if (value_notzero_p(e->d) || e->x.p->type != partition) {
3191 if (!shift) {
3192 if (evalue_frac2floor_in_domain3(e, NULL, 0))
3193 reduce_evalue(e);
3194 }
3195 return;
3196 }
3197
3198 for (i = 0; i < e->x.p->size/2; ++i)
3199 if (evalue_frac2floor_in_domain3(&e->x.p->arr[2*i+1],
3200 EVALUE_DOMAIN(e->x.p->arr[2*i]), shift))
3201 reduce_evalue(&e->x.p->arr[2*i+1]);
3202 }
3203
3204 void evalue_frac2floor(evalue *e)
3205 {
3206 evalue_frac2floor2(e, 1);
3207 }
3208
3209 /* Add a new variable with lower bound 1 and upper bound specified
3210 * by row. If negative is true, then the new variable has upper
3211 * bound -1 and lower bound specified by row.
3212 */
3213 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
3214 Vector *row, int negative)
3215 {
3216 int nr, nc;
3217 int i;
3218 int nparam = D->Dimension - nvar;
3219
3220 if (C == 0) {
3221 nr = D->NbConstraints + 2;
3222 nc = D->Dimension + 2 + 1;
3223 C = Matrix_Alloc(nr, nc);
3224 for (i = 0; i < D->NbConstraints; ++i) {
3225 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
3226 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3227 D->Dimension + 1 - nvar);
3228 }
3229 } else {
3230 Matrix *oldC = C;
3231 nr = C->NbRows + 2;
3232 nc = C->NbColumns + 1;
3233 C = Matrix_Alloc(nr, nc);
3234 for (i = 0; i < oldC->NbRows; ++i) {
3235 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
3236 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3237 oldC->NbColumns - 1 - nvar);
3238 }
3239 }
3240 value_set_si(C->p[nr-2][0], 1);
3241 if (negative)
3242 value_set_si(C->p[nr-2][1 + nvar], -1);
3243 else
3244 value_set_si(C->p[nr-2][1 + nvar], 1);
3245 value_set_si(C->p[nr-2][nc - 1], -1);
3246
3247 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
3248 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
3249 1 + nparam);
3250
3251 return C;
3252 }
3253
3254 static void floor2frac_r(evalue *e, int nvar)
3255 {
3256 enode *p;
3257 int i;
3258 evalue f;
3259 evalue *pp;
3260
3261 if (value_notzero_p(e->d))
3262 return;
3263
3264 p = e->x.p;
3265
3266 for (i = type_offset(p); i < p->size; i++)
3267 floor2frac_r(&p->arr[i], nvar);
3268
3269 if (p->type != flooring)
3270 return;
3271
3272 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
3273 assert(pp->x.p->type == polynomial);
3274 pp->x.p->pos -= nvar;
3275 }
3276
3277 value_init(f.d);
3278 value_set_si(f.d, 0);
3279 f.x.p = new_enode(fractional, 3, -1);
3280 evalue_set_si(&f.x.p->arr[1], 0, 1);
3281 evalue_set_si(&f.x.p->arr[2], -1, 1);
3282 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3283
3284 eadd(&f, &p->arr[0]);
3285 reorder_terms_about(p, &p->arr[0]);
3286 value_clear(e->d);
3287 *e = p->arr[1];
3288 free(p);
3289 free_evalue_refs(&f);
3290 }
3291
3292 /* Convert flooring back to fractional and shift position
3293 * of the parameters by nvar
3294 */
3295 static void floor2frac(evalue *e, int nvar)
3296 {
3297 floor2frac_r(e, nvar);
3298 reduce_evalue(e);
3299 }
3300
3301 int evalue_floor2frac(evalue *e)
3302 {
3303 int i;
3304 int r = 0;
3305
3306 if (value_notzero_p(e->d))
3307 return 0;
3308
3309 if (e->x.p->type == partition) {
3310 for (i = 0; i < e->x.p->size/2; ++i)
3311 if (evalue_floor2frac(&e->x.p->arr[2*i+1]))
3312 reduce_evalue(&e->x.p->arr[2*i+1]);
3313 return 0;
3314 }
3315
3316 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
3317 r |= evalue_floor2frac(&e->x.p->arr[i]);
3318
3319 if (e->x.p->type == flooring) {
3320 floor2frac(e, 0);
3321 return 1;
3322 }
3323
3324 if (r)
3325 evalue_reorder_terms(e);
3326
3327 return r;
3328 }
3329
3330 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
3331 {
3332 evalue *t;
3333 int nparam = D->Dimension - nvar;
3334
3335 if (C != 0) {
3336 C = Matrix_Copy(C);
3337 D = Constraints2Polyhedron(C, 0);
3338 Matrix_Free(C);
3339 }
3340
3341 t = barvinok_enumerate_e(D, 0, nparam, 0);
3342
3343 /* Double check that D was not unbounded. */
3344 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3345
3346 if (C != 0)
3347 Polyhedron_Free(D);
3348
3349 return t;
3350 }
3351
3352 static void domain_signs(Polyhedron *D, int *signs)
3353 {
3354 int j, k;
3355
3356 POL_ENSURE_VERTICES(D);
3357 for (j = 0; j < D->Dimension; ++j) {
3358 signs[j] = 0;
3359 for (k = 0; k < D->NbRays; ++k) {
3360 signs[j] = value_sign(D->Ray[k][1+j]);
3361 if (signs[j])
3362 break;
3363 }
3364 }
3365 }
3366
3367 static evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3368 int *signs, Matrix *C, unsigned MaxRays)
3369 {
3370 Vector *row = NULL;
3371 int i, offset;
3372 evalue *res;
3373 Matrix *origC;
3374 evalue *factor = NULL;
3375 evalue cum;
3376 int negative = 0;
3377 int allocated = 0;
3378
3379 if (EVALUE_IS_ZERO(*e))
3380 return 0;
3381
3382 if (D->next) {
3383 Polyhedron *DD = Disjoint_Domain(D, 0, MaxRays);
3384 Polyhedron *Q;
3385
3386 Q = DD;
3387 DD = Q->next;
3388 Q->next = 0;
3389
3390 res = esum_over_domain(e, nvar, Q, signs, C, MaxRays);
3391 Polyhedron_Free(Q);
3392
3393 for (Q = DD; Q; Q = DD) {
3394 evalue *t;
3395
3396 DD = Q->next;
3397 Q->next = 0;
3398
3399 t = esum_over_domain(e, nvar, Q, signs, C, MaxRays);
3400 Polyhedron_Free(Q);
3401
3402 if (!res)
3403 res = t;
3404 else if (t) {
3405 eadd(t, res);
3406 evalue_free(t);
3407 }
3408 }
3409 return res;
3410 }
3411
3412 if (value_notzero_p(e->d)) {
3413 evalue *t;
3414
3415 t = esum_over_domain_cst(nvar, D, C);
3416
3417 if (!EVALUE_IS_ONE(*e))
3418 emul(e, t);
3419
3420 return t;
3421 }
3422
3423 if (!signs) {
3424 allocated = 1;
3425 signs = ALLOCN(int, D->Dimension);
3426 domain_signs(D, signs);
3427 }
3428
3429 switch (e->x.p->type) {
3430 case flooring: {
3431 evalue *pp = &e->x.p->arr[0];
3432
3433 if (pp->x.p->pos > nvar) {
3434 /* remainder is independent of the summated vars */
3435 evalue f;
3436 evalue *t;
3437
3438 value_init(f.d);
3439 evalue_copy(&f, e);
3440 floor2frac(&f, nvar);
3441
3442 t = esum_over_domain_cst(nvar, D, C);
3443
3444 emul(&f, t);
3445
3446 free_evalue_refs(&f);
3447
3448 if (allocated)
3449 free(signs);
3450
3451 return t;
3452 }
3453
3454 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3455 poly_denom(pp, &row->p[1 + nvar]);
3456 value_set_si(row->p[0], 1);
3457 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3458 pp = &pp->x.p->arr[0]) {
3459 int pos;
3460 assert(pp->x.p->type == polynomial);
3461 pos = pp->x.p->pos;
3462 if (pos >= 1 + nvar)
3463 ++pos;
3464 value_assign(row->p[pos], row->p[1+nvar]);
3465 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3466 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3467 }
3468 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3469 value_division(row->p[1 + D->Dimension + 1],
3470 row->p[1 + D->Dimension + 1],
3471 pp->d);
3472 value_multiply(row->p[1 + D->Dimension + 1],
3473 row->p[1 + D->Dimension + 1],
3474 pp->x.n);
3475 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3476 break;
3477 }
3478 case polynomial: {
3479 int pos = e->x.p->pos;
3480
3481 if (pos > nvar) {
3482 factor = ALLOC(evalue);
3483 value_init(factor->d);
3484 value_set_si(factor->d, 0);
3485 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3486 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3487 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3488 break;
3489 }
3490
3491 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3492 negative = signs[pos-1] < 0;
3493 value_set_si(row->p[0], 1);
3494 if (negative) {
3495 value_set_si(row->p[pos], -1);
3496 value_set_si(row->p[1 + nvar], 1);
3497 } else {
3498 value_set_si(row->p[pos], 1);
3499 value_set_si(row->p[1 + nvar], -1);
3500 }
3501 break;
3502 }
3503 default:
3504 assert(0);
3505 }
3506
3507 offset = type_offset(e->x.p);
3508
3509 res = esum_over_domain(&e->x.p->arr[offset], nvar, D, signs, C, MaxRays);
3510
3511 if (factor) {
3512 value_init(cum.d);
3513 evalue_copy(&cum, factor);
3514 }
3515
3516 origC = C;
3517 for (i = 1; offset+i < e->x.p->size; ++i) {
3518 evalue *t;
3519 if (row) {
3520 Matrix *prevC = C;
3521 C = esum_add_constraint(nvar, D, C, row, negative);
3522 if (prevC != origC)
3523 Matrix_Free(prevC);
3524 }
3525 /*
3526 if (row)
3527 Vector_Print(stderr, P_VALUE_FMT, row);
3528 if (C)
3529 Matrix_Print(stderr, P_VALUE_FMT, C);
3530 */
3531 t = esum_over_domain(&e->x.p->arr[offset+i], nvar, D, signs, C, MaxRays);
3532
3533 if (t) {
3534 if (factor)
3535 emul(&cum, t);
3536 if (negative && (i % 2))
3537 evalue_negate(t);
3538 }
3539
3540 if (!res)
3541 res = t;
3542 else if (t) {
3543 eadd(t, res);
3544 evalue_free(t);
3545 }
3546 if (factor && offset+i+1 < e->x.p->size)
3547 emul(factor, &cum);
3548 }
3549 if (C != origC)
3550 Matrix_Free(C);
3551
3552 if (factor) {
3553 free_evalue_refs(&cum);
3554 evalue_free(factor);
3555 }
3556
3557 if (row)
3558 Vector_Free(row);
3559
3560 reduce_evalue(res);
3561
3562 if (allocated)
3563 free(signs);
3564
3565 return res;
3566 }
3567
3568 static void Polyhedron_Insert(Polyhedron ***next, Polyhedron *Q)
3569 {
3570 if (emptyQ(Q))
3571 Polyhedron_Free(Q);
3572 else {
3573 **next = Q;
3574 *next = &(Q->next);
3575 }
3576 }
3577
3578 static Polyhedron *Polyhedron_Split_Into_Orthants(Polyhedron *P,
3579 unsigned MaxRays)
3580 {
3581 int i = 0;
3582 Polyhedron *D = P;
3583 Vector *c = Vector_Alloc(1 + P->Dimension + 1);
3584 value_set_si(c->p[0], 1);
3585
3586 if (P->Dimension == 0)
3587 return Polyhedron_Copy(P);
3588
3589 for (i = 0; i < P->Dimension; ++i) {
3590 Polyhedron *L = NULL;
3591 Polyhedron **next = &L;
3592 Polyhedron *I;
3593
3594 for (I = D; I; I = I->next) {
3595 Polyhedron *Q;
3596 value_set_si(c->p[1+i], 1);
3597 value_set_si(c->p[1+P->Dimension], 0);
3598 Q = AddConstraints(c->p, 1, I, MaxRays);
3599 Polyhedron_Insert(&next, Q);
3600 value_set_si(c->p[1+i], -1);
3601 value_set_si(c->p[1+P->Dimension], -1);
3602 Q = AddConstraints(c->p, 1, I, MaxRays);
3603 Polyhedron_Insert(&next, Q);
3604 value_set_si(c->p[1+i], 0);
3605 }
3606 if (D != P)
3607 Domain_Free(D);
3608 D = L;
3609 }
3610 Vector_Free(c);
3611 return D;
3612 }
3613
3614 /* Make arguments of all floors non-negative */
3615 static void shift_floor_in_domain(evalue *e, Polyhedron *D)
3616 {
3617 Value d, m;
3618 Polyhedron *I;
3619 int i;
3620 enode *p;
3621
3622 if (value_notzero_p(e->d))
3623 return;
3624
3625 p = e->x.p;
3626
3627 for (i = type_offset(p); i < p->size; ++i)
3628 shift_floor_in_domain(&p->arr[i], D);
3629
3630 if (p->type != flooring)
3631 return;
3632
3633 value_init(d);
3634 value_init(m);
3635
3636 I = polynomial_projection(p, D, &d, NULL);
3637 assert(I->NbEq == 0); /* Should have been reduced */
3638
3639 for (i = 0; i < I->NbConstraints; ++i)
3640 if (value_pos_p(I->Constraint[i][1]))
3641 break;
3642 assert(i < I->NbConstraints);
3643 if (i < I->NbConstraints) {
3644 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
3645 mpz_fdiv_q(m, I->Constraint[i][2], I->Constraint[i][1]);
3646 if (value_neg_p(m)) {
3647 /* replace [e] by [e-m]+m such that e-m >= 0 */
3648 evalue f;
3649
3650 value_init(f.d);
3651 value_init(f.x.n);
3652 value_set_si(f.d, 1);
3653 value_oppose(f.x.n, m);
3654 eadd(&f, &p->arr[0]);
3655 value_clear(f.x.n);
3656
3657 value_set_si(f.d, 0);
3658 f.x.p = new_enode(flooring, 3, -1);
3659 value_clear(f.x.p->arr[0].d);
3660 f.x.p->arr[0] = p->arr[0];
3661 evalue_set_si(&f.x.p->arr[2], 1, 1);
3662 value_set_si(f.x.p->arr[1].d, 1);
3663 value_init(f.x.p->arr[1].x.n);
3664 value_assign(f.x.p->arr[1].x.n, m);
3665 reorder_terms_about(p, &f);
3666 value_clear(e->d);
3667 *e = p->arr[1];
3668 free(p);
3669 }
3670 }
3671 Polyhedron_Free(I);
3672 value_clear(d);
3673 value_clear(m);
3674 }
3675
3676 evalue *box_summate(Polyhedron *P, evalue *E, unsigned nvar, unsigned MaxRays)
3677 {
3678 evalue *sum = evalue_zero();
3679 Polyhedron *D = Polyhedron_Split_Into_Orthants(P, MaxRays);
3680 for (P = D; P; P = P->next) {
3681 evalue *t;
3682 evalue *fe = evalue_dup(E);
3683 Polyhedron *next = P->next;
3684 P->next = NULL;
3685 reduce_evalue_in_domain(fe, P);
3686 evalue_frac2floor2(fe, 0);
3687 shift_floor_in_domain(fe, P);
3688 t = esum_over_domain(fe, nvar, P, NULL, NULL, MaxRays);
3689 if (t) {
3690 eadd(t, sum);
3691 evalue_free(t);
3692 }
3693 evalue_free(fe);
3694 P->next = next;
3695 }
3696 Domain_Free(D);
3697 return sum;
3698 }
3699
3700 /* Initial silly implementation */
3701 void eor(evalue *e1, evalue *res)
3702 {
3703 evalue E;
3704 evalue mone;
3705 value_init(E.d);
3706 value_init(mone.d);
3707 evalue_set_si(&mone, -1, 1);
3708
3709 evalue_copy(&E, res);
3710 eadd(e1, &E);
3711 emul(e1, res);
3712 emul(&mone, res);
3713 eadd(&E, res);
3714
3715 free_evalue_refs(&E);
3716 free_evalue_refs(&mone);
3717 }
3718
3719 /* computes denominator of polynomial evalue
3720 * d should point to a value initialized to 1
3721 */
3722 void evalue_denom(const evalue *e, Value *d)
3723 {
3724 int i, offset;
3725
3726 if (value_notzero_p(e->d)) {
3727 value_lcm(*d, *d, e->d);
3728 return;
3729 }
3730 assert(e->x.p->type == polynomial ||
3731 e->x.p->type == fractional ||
3732 e->x.p->type == flooring);
3733 offset = type_offset(e->x.p);
3734 for (i = e->x.p->size-1; i >= offset; --i)
3735 evalue_denom(&e->x.p->arr[i], d);
3736 }
3737
3738 /* Divides the evalue e by the integer n */
3739 void evalue_div(evalue *e, Value n)
3740 {
3741 int i, offset;
3742
3743 if (value_one_p(n) || 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->d, e->d, n);
3750 value_gcd(gc, e->x.n, e->d);
3751 if (value_notone_p(gc)) {
3752 value_division(e->d, e->d, gc);
3753 value_division(e->x.n, e->x.n, gc);
3754 }
3755 value_clear(gc);
3756 return;
3757 }
3758 if (e->x.p->type == partition) {
3759 for (i = 0; i < e->x.p->size/2; ++i)
3760 evalue_div(&e->x.p->arr[2*i+1], n);
3761 return;
3762 }
3763 offset = type_offset(e->x.p);
3764 for (i = e->x.p->size-1; i >= offset; --i)
3765 evalue_div(&e->x.p->arr[i], n);
3766 }
3767
3768 /* Multiplies the evalue e by the integer n */
3769 void evalue_mul(evalue *e, Value n)
3770 {
3771 int i, offset;
3772
3773 if (value_one_p(n) || EVALUE_IS_ZERO(*e))
3774 return;
3775
3776 if (value_notzero_p(e->d)) {
3777 Value gc;
3778 value_init(gc);
3779 value_multiply(e->x.n, e->x.n, n);
3780 value_gcd(gc, e->x.n, e->d);
3781 if (value_notone_p(gc)) {
3782 value_division(e->d, e->d, gc);
3783 value_division(e->x.n, e->x.n, gc);
3784 }
3785 value_clear(gc);
3786 return;
3787 }
3788 if (e->x.p->type == partition) {
3789 for (i = 0; i < e->x.p->size/2; ++i)
3790 evalue_mul(&e->x.p->arr[2*i+1], n);
3791 return;
3792 }
3793 offset = type_offset(e->x.p);
3794 for (i = e->x.p->size-1; i >= offset; --i)
3795 evalue_mul(&e->x.p->arr[i], n);
3796 }
3797
3798 /* Multiplies the evalue e by the n/d */
3799 void evalue_mul_div(evalue *e, Value n, Value d)
3800 {
3801 int i, offset;
3802
3803 if ((value_one_p(n) && value_one_p(d)) || EVALUE_IS_ZERO(*e))
3804 return;
3805
3806 if (value_notzero_p(e->d)) {
3807 Value gc;
3808 value_init(gc);
3809 value_multiply(e->x.n, e->x.n, n);
3810 value_multiply(e->d, e->d, d);
3811 value_gcd(gc, e->x.n, e->d);
3812 if (value_notone_p(gc)) {
3813 value_division(e->d, e->d, gc);
3814 value_division(e->x.n, e->x.n, gc);
3815 }
3816 value_clear(gc);
3817 return;
3818 }
3819 if (e->x.p->type == partition) {
3820 for (i = 0; i < e->x.p->size/2; ++i)
3821 evalue_mul_div(&e->x.p->arr[2*i+1], n, d);
3822 return;
3823 }
3824 offset = type_offset(e->x.p);
3825 for (i = e->x.p->size-1; i >= offset; --i)
3826 evalue_mul_div(&e->x.p->arr[i], n, d);
3827 }
3828
3829 void evalue_negate(evalue *e)
3830 {
3831 int i, offset;
3832
3833 if (value_notzero_p(e->d)) {
3834 value_oppose(e->x.n, e->x.n);
3835 return;
3836 }
3837 if (e->x.p->type == partition) {
3838 for (i = 0; i < e->x.p->size/2; ++i)
3839 evalue_negate(&e->x.p->arr[2*i+1]);
3840 return;
3841 }
3842 offset = type_offset(e->x.p);
3843 for (i = e->x.p->size-1; i >= offset; --i)
3844 evalue_negate(&e->x.p->arr[i]);
3845 }
3846
3847 void evalue_add_constant(evalue *e, const Value cst)
3848 {
3849 int i;
3850
3851 if (value_zero_p(e->d)) {
3852 if (e->x.p->type == partition) {
3853 for (i = 0; i < e->x.p->size/2; ++i)
3854 evalue_add_constant(&e->x.p->arr[2*i+1], cst);
3855 return;
3856 }
3857 if (e->x.p->type == relation) {
3858 for (i = 1; i < e->x.p->size; ++i)
3859 evalue_add_constant(&e->x.p->arr[i], cst);
3860 return;
3861 }
3862 do {
3863 e = &e->x.p->arr[type_offset(e->x.p)];
3864 } while (value_zero_p(e->d));
3865 }
3866 value_addmul(e->x.n, cst, e->d);
3867 }
3868
3869 static void evalue_frac2polynomial_r(evalue *e, int *signs, int sign, int in_frac)
3870 {
3871 int i, offset;
3872 Value d;
3873 enode *p;
3874 int sign_odd = sign;
3875
3876 if (value_notzero_p(e->d)) {
3877 if (in_frac && sign * value_sign(e->x.n) < 0) {
3878 value_set_si(e->x.n, 0);
3879 value_set_si(e->d, 1);
3880 }
3881 return;
3882 }
3883
3884 if (e->x.p->type == relation) {
3885 for (i = e->x.p->size-1; i >= 1; --i)
3886 evalue_frac2polynomial_r(&e->x.p->arr[i], signs, sign, in_frac);
3887 return;
3888 }
3889
3890 if (e->x.p->type == polynomial)
3891 sign_odd *= signs[e->x.p->pos-1];
3892 offset = type_offset(e->x.p);
3893 evalue_frac2polynomial_r(&e->x.p->arr[offset], signs, sign, in_frac);
3894 in_frac |= e->x.p->type == fractional;
3895 for (i = e->x.p->size-1; i > offset; --i)
3896 evalue_frac2polynomial_r(&e->x.p->arr[i], signs,
3897 (i - offset) % 2 ? sign_odd : sign, in_frac);
3898
3899 if (e->x.p->type != fractional)
3900 return;
3901
3902 /* replace { a/m } by (m-1)/m if sign != 0
3903 * and by (m-1)/(2m) if sign == 0
3904 */
3905 value_init(d);
3906 value_set_si(d, 1);
3907 evalue_denom(&e->x.p->arr[0], &d);
3908 free_evalue_refs(&e->x.p->arr[0]);
3909 value_init(e->x.p->arr[0].d);
3910 value_init(e->x.p->arr[0].x.n);
3911 if (sign == 0)
3912 value_addto(e->x.p->arr[0].d, d, d);
3913 else
3914 value_assign(e->x.p->arr[0].d, d);
3915 value_decrement(e->x.p->arr[0].x.n, d);
3916 value_clear(d);
3917
3918 p = e->x.p;
3919 reorder_terms_about(p, &p->arr[0]);
3920 value_clear(e->d);
3921 *e = p->arr[1];
3922 free(p);
3923 }
3924
3925 /* Approximate the evalue in fractional representation by a polynomial.
3926 * If sign > 0, the result is an upper bound;
3927 * if sign < 0, the result is a lower bound;
3928 * if sign = 0, the result is an intermediate approximation.
3929 */
3930 void evalue_frac2polynomial(evalue *e, int sign, unsigned MaxRays)
3931 {
3932 int i, dim;
3933 int *signs;
3934
3935 if (value_notzero_p(e->d))
3936 return;
3937 assert(e->x.p->type == partition);
3938 /* make sure all variables in the domains have a fixed sign */
3939 if (sign) {
3940 evalue_split_domains_into_orthants(e, MaxRays);
3941 if (EVALUE_IS_ZERO(*e))
3942 return;
3943 }
3944
3945 assert(e->x.p->size >= 2);
3946 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3947
3948 signs = ALLOCN(int, dim);
3949
3950 if (!sign)
3951 for (i = 0; i < dim; ++i)
3952 signs[i] = 0;
3953 for (i = 0; i < e->x.p->size/2; ++i) {
3954 if (sign)
3955 domain_signs(EVALUE_DOMAIN(e->x.p->arr[2*i]), signs);
3956 evalue_frac2polynomial_r(&e->x.p->arr[2*i+1], signs, sign, 0);
3957 }
3958
3959 free(signs);
3960 }
3961
3962 /* Split the domains of e (which is assumed to be a partition)
3963 * such that each resulting domain lies entirely in one orthant.
3964 */
3965 void evalue_split_domains_into_orthants(evalue *e, unsigned MaxRays)
3966 {
3967 int i, dim;
3968 assert(value_zero_p(e->d));
3969 assert(e->x.p->type == partition);
3970 assert(e->x.p->size >= 2);
3971 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3972
3973 for (i = 0; i < dim; ++i) {
3974 evalue split;
3975 Matrix *C, *C2;
3976 C = Matrix_Alloc(1, 1 + dim + 1);
3977 value_set_si(C->p[0][0], 1);
3978 value_init(split.d);
3979 value_set_si(split.d, 0);
3980 split.x.p = new_enode(partition, 4, dim);
3981 value_set_si(C->p[0][1+i], 1);
3982 C2 = Matrix_Copy(C);
3983 EVALUE_SET_DOMAIN(split.x.p->arr[0], Constraints2Polyhedron(C2, MaxRays));
3984 Matrix_Free(C2);
3985 evalue_set_si(&split.x.p->arr[1], 1, 1);
3986 value_set_si(C->p[0][1+i], -1);
3987 value_set_si(C->p[0][1+dim], -1);
3988 EVALUE_SET_DOMAIN(split.x.p->arr[2], Constraints2Polyhedron(C, MaxRays));
3989 evalue_set_si(&split.x.p->arr[3], 1, 1);
3990 emul(&split, e);
3991 free_evalue_refs(&split);
3992 Matrix_Free(C);
3993 }
3994 }
3995
3996 static evalue *find_fractional_with_max_periods(evalue *e, Polyhedron *D,
3997 int max_periods,
3998 Matrix **TT,
3999 Value *min, Value *max)
4000 {
4001 Matrix *T;
4002 evalue *f = NULL;
4003 Value d;
4004 int i;
4005
4006 if (value_notzero_p(e->d))
4007 return NULL;
4008
4009 if (e->x.p->type == fractional) {
4010 Polyhedron *I;
4011 int bounded;
4012
4013 value_init(d);
4014 I = polynomial_projection(e->x.p, D, &d, &T);
4015 bounded = line_minmax(I, min, max); /* frees I */
4016 if (bounded) {
4017 Value mp;
4018 value_init(mp);
4019 value_set_si(mp, max_periods);
4020 mpz_fdiv_q(*min, *min, d);
4021 mpz_fdiv_q(*max, *max, d);
4022 value_assign(T->p[1][D->Dimension], d);
4023 value_subtract(d, *max, *min);
4024 if (value_ge(d, mp))
4025 Matrix_Free(T);
4026 else
4027 f = evalue_dup(&e->x.p->arr[0]);
4028 value_clear(mp);
4029 } else
4030 Matrix_Free(T);
4031 value_clear(d);
4032 if (f) {
4033 *TT = T;
4034 return f;
4035 }
4036 }
4037
4038 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
4039 if ((f = find_fractional_with_max_periods(&e->x.p->arr[i], D, max_periods,
4040 TT, min, max)))
4041 return f;
4042
4043 return NULL;
4044 }
4045
4046 static void replace_fract_by_affine(evalue *e, evalue *f, Value val)
4047 {
4048 int i, offset;
4049
4050 if (value_notzero_p(e->d))
4051 return;
4052
4053 offset = type_offset(e->x.p);
4054 for (i = e->x.p->size-1; i >= offset; --i)
4055 replace_fract_by_affine(&e->x.p->arr[i], f, val);
4056
4057 if (e->x.p->type != fractional)
4058 return;
4059
4060 if (!eequal(&e->x.p->arr[0], f))
4061 return;
4062
4063 replace_by_affine(e, val);
4064 }
4065
4066 /* Look for fractional parts that can be removed by splitting the corresponding
4067 * domain into at most max_periods parts.
4068 * We use a very simply strategy that looks for the first fractional part
4069 * that satisfies the condition, performs the split and then continues
4070 * looking for other fractional parts in the split domains until no
4071 * such fractional part can be found anymore.
4072 */
4073 void evalue_split_periods(evalue *e, int max_periods, unsigned int MaxRays)
4074 {
4075 int i, j, n;
4076 Value min;
4077 Value max;
4078 Value d;
4079
4080 if (EVALUE_IS_ZERO(*e))
4081 return;
4082 if (value_notzero_p(e->d) || e->x.p->type != partition) {
4083 fprintf(stderr,
4084 "WARNING: evalue_split_periods called on incorrect evalue type\n");
4085 return;
4086 }
4087
4088 value_init(min);
4089 value_init(max);
4090 value_init(d);
4091
4092 for (i = 0; i < e->x.p->size/2; ++i) {
4093 enode *p;
4094 evalue *f;
4095 Matrix *T;
4096 Matrix *M;
4097 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
4098 Polyhedron *E;
4099 f = find_fractional_with_max_periods(&e->x.p->arr[2*i+1], D, max_periods,
4100 &T, &min, &max);
4101 if (!f)
4102 continue;
4103
4104 M = Matrix_Alloc(2, 2+D->Dimension);
4105
4106 value_subtract(d, max, min);
4107 n = VALUE_TO_INT(d)+1;
4108
4109 value_set_si(M->p[0][0], 1);
4110 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
4111 value_multiply(d, max, T->p[1][D->Dimension]);
4112 value_subtract(M->p[0][1+D->Dimension], M->p[0][1+D->Dimension], d);
4113 value_set_si(d, -1);
4114 value_set_si(M->p[1][0], 1);
4115 Vector_Scale(T->p[0], M->p[1]+1, d, D->Dimension+1);
4116 value_addmul(M->p[1][1+D->Dimension], max, T->p[1][D->Dimension]);
4117 value_addto(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension],
4118 T->p[1][D->Dimension]);
4119 value_decrement(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension]);
4120
4121 p = new_enode(partition, e->x.p->size + (n-1)*2, e->x.p->pos);
4122 for (j = 0; j < 2*i; ++j) {
4123 value_clear(p->arr[j].d);
4124 p->arr[j] = e->x.p->arr[j];
4125 }
4126 for (j = 2*i+2; j < e->x.p->size; ++j) {
4127 value_clear(p->arr[j+2*(n-1)].d);
4128 p->arr[j+2*(n-1)] = e->x.p->arr[j];
4129 }
4130 for (j = n-1; j >= 0; --j) {
4131 if (j == 0) {
4132 value_clear(p->arr[2*i+1].d);
4133 p->arr[2*i+1] = e->x.p->arr[2*i+1];
4134 } else
4135 evalue_copy(&p->arr[2*(i+j)+1], &e->x.p->arr[2*i+1]);
4136 if (j != n-1) {
4137 value_subtract(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension],
4138 T->p[1][D->Dimension]);
4139 value_addto(M->p[0][1+D->Dimension], M->p[0][1+D->Dimension],
4140 T->p[1][D->Dimension]);
4141 }
4142 replace_fract_by_affine(&p->arr[2*(i+j)+1], f, max);
4143 E = DomainAddConstraints(D, M, MaxRays);
4144 EVALUE_SET_DOMAIN(p->arr[2*(i+j)], E);
4145 if (evalue_range_reduction_in_domain(&p->arr[2*(i+j)+1], E))
4146 reduce_evalue(&p->arr[2*(i+j)+1]);
4147 value_decrement(max, max);
4148 }
4149 value_clear(e->x.p->arr[2*i].d);
4150 Domain_Free(D);
4151 Matrix_Free(M);
4152 Matrix_Free(T);
4153 evalue_free(f);
4154 free(e->x.p);
4155 e->x.p = p;
4156 --i;
4157 }
4158
4159 value_clear(d);
4160 value_clear(min);
4161 value_clear(max);
4162 }
4163
4164 void evalue_extract_affine(const evalue *e, Value *coeff, Value *cst, Value *d)
4165 {
4166 value_set_si(*d, 1);
4167 evalue_denom(e, d);
4168 for ( ; value_zero_p(e->d); e = &e->x.p->arr[0]) {
4169 evalue *c;
4170 assert(e->x.p->type == polynomial);
4171 assert(e->x.p->size == 2);
4172 c = &e->x.p->arr[1];
4173 value_multiply(coeff[e->x.p->pos-1], *d, c->x.n);
4174 value_division(coeff[e->x.p->pos-1], coeff[e->x.p->pos-1], c->d);
4175 }
4176 value_multiply(*cst, *d, e->x.n);
4177 value_division(*cst, *cst, e->d);
4178 }
4179
4180 /* returns an evalue that corresponds to
4181 *
4182 * c/den x_param
4183 */
4184 static evalue *term(int param, Value c, Value den)
4185 {
4186 evalue *EP = ALLOC(evalue);
4187 value_init(EP->d);
4188 value_set_si(EP->d,0);
4189 EP->x.p = new_enode(polynomial, 2, param + 1);
4190 evalue_set_si(&EP->x.p->arr[0], 0, 1);
4191 evalue_set_reduce(&EP->x.p->arr[1], c, den);
4192 return EP;
4193 }
4194
4195 evalue *affine2evalue(Value *coeff, Value denom, int nvar)
4196 {
4197 int i;
4198 evalue *E = ALLOC(evalue);
4199 value_init(E->d);
4200 evalue_set_reduce(E, coeff[nvar], denom);
4201 for (i = 0; i < nvar; ++i) {
4202 evalue *t;
4203 if (value_zero_p(coeff[i]))
4204 continue;
4205 t = term(i, coeff[i], denom);
4206 eadd(t, E);
4207 evalue_free(t);
4208 }
4209 return E;
4210 }
4211
4212 void evalue_substitute(evalue *e, evalue **subs)
4213 {
4214 int i, offset;
4215 evalue *v;
4216 enode *p;
4217
4218 if (value_notzero_p(e->d))
4219 return;
4220
4221 p = e->x.p;
4222 assert(p->type != partition);
4223
4224 for (i = 0; i < p->size; ++i)
4225 evalue_substitute(&p->arr[i], subs);
4226
4227 if (p->type == relation) {
4228 /* For relation a ? b : c, compute (a' ? 1) * b' + (a' ? 0 : 1) * c' */
4229 if (p->size == 3) {
4230 v = ALLOC(evalue);
4231 value_init(v->d);
4232 value_set_si(v->d, 0);
4233 v->x.p = new_enode(relation, 3, 0);
4234 evalue_copy(&v->x.p->arr[0], &p->arr[0]);
4235 evalue_set_si(&v->x.p->arr[1], 0, 1);
4236 evalue_set_si(&v->x.p->arr[2], 1, 1);
4237 emul(v, &p->arr[2]);
4238 evalue_free(v);
4239 }
4240 v = ALLOC(evalue);
4241 value_init(v->d);
4242 value_set_si(v->d, 0);
4243 v->x.p = new_enode(relation, 2, 0);
4244 value_clear(v->x.p->arr[0].d);
4245 v->x.p->arr[0] = p->arr[0];
4246 evalue_set_si(&v->x.p->arr[1], 1, 1);
4247 emul(v, &p->arr[1]);
4248 evalue_free(v);
4249 if (p->size == 3) {
4250 eadd(&p->arr[2], &p->arr[1]);
4251 free_evalue_refs(&p->arr[2]);
4252 }
4253 value_clear(e->d);
4254 *e = p->arr[1];
4255 free(p);
4256 return;
4257 }
4258
4259 if (p->type == polynomial)
4260 v = subs[p->pos-1];
4261 else {
4262 v = ALLOC(evalue);
4263 value_init(v->d);
4264 value_set_si(v->d, 0);
4265 v->x.p = new_enode(p->type, 3, -1);
4266 value_clear(v->x.p->arr[0].d);
4267 v->x.p->arr[0] = p->arr[0];
4268 evalue_set_si(&v->x.p->arr[1], 0, 1);
4269 evalue_set_si(&v->x.p->arr[2], 1, 1);
4270 }
4271
4272 offset = type_offset(p);
4273
4274 for (i = p->size-1; i >= offset+1; i--) {
4275 emul(v, &p->arr[i]);
4276 eadd(&p->arr[i], &p->arr[i-1]);
4277 free_evalue_refs(&(p->arr[i]));
4278 }
4279
4280 if (p->type != polynomial)
4281 evalue_free(v);
4282
4283 value_clear(e->d);
4284 *e = p->arr[offset];
4285 free(p);
4286 }
4287
4288 /* evalue e is given in terms of "new" parameter; CP maps the new
4289 * parameters back to the old parameters.
4290 * Transforms e such that it refers back to the old parameters and
4291 * adds appropriate constraints to the domain.
4292 * In particular, if CP maps the new parameters onto an affine
4293 * subspace of the old parameters, then the corresponding equalities
4294 * are added to the domain.
4295 * Also, if any of the new parameters was a rational combination
4296 * of the old parameters $p' = (<a, p> + c)/m$, then modulo
4297 * constraints ${<a, p> + c)/m} = 0$ are added to ensure
4298 * the new evalue remains non-zero only for integer parameters
4299 * of the new parameters (which have been removed by the substitution).
4300 */
4301 void evalue_backsubstitute(evalue *e, Matrix *CP, unsigned MaxRays)
4302 {
4303 Matrix *eq;
4304 Matrix *inv;
4305 evalue **subs;
4306 enode *p;
4307 int i, j;
4308 unsigned nparam = CP->NbColumns-1;
4309 Polyhedron *CEq;
4310 Value gcd;
4311
4312 if (EVALUE_IS_ZERO(*e))
4313 return;
4314
4315 assert(value_zero_p(e->d));
4316 p = e->x.p;
4317 assert(p->type == partition);
4318
4319 inv = left_inverse(CP, &eq);
4320 subs = ALLOCN(evalue *, nparam);
4321 for (i = 0; i < nparam; ++i)
4322 subs[i] = affine2evalue(inv->p[i], inv->p[nparam][inv->NbColumns-1],
4323 inv->NbColumns-1);
4324
4325 CEq = Constraints2Polyhedron(eq, MaxRays);
4326 addeliminatedparams_partition(p, inv, CEq, inv->NbColumns-1, MaxRays);
4327 Polyhedron_Free(CEq);
4328
4329 for (i = 0; i < p->size/2; ++i)
4330 evalue_substitute(&p->arr[2*i+1], subs);
4331
4332 for (i = 0; i < nparam; ++i)
4333 evalue_free(subs[i]);
4334 free(subs);
4335
4336 value_init(gcd);
4337 for (i = 0; i < inv->NbRows-1; ++i) {
4338 Vector_Gcd(inv->p[i], inv->NbColumns, &gcd);
4339 value_gcd(gcd, gcd, inv->p[inv->NbRows-1][inv->NbColumns-1]);
4340 if (value_eq(gcd, inv->p[inv->NbRows-1][inv->NbColumns-1]))
4341 continue;
4342 Vector_AntiScale(inv->p[i], inv->p[i], gcd, inv->NbColumns);
4343 value_divexact(gcd, inv->p[inv->NbRows-1][inv->NbColumns-1], gcd);
4344
4345 for (j = 0; j < p->size/2; ++j) {
4346 evalue *arg = affine2evalue(inv->p[i], gcd, inv->NbColumns-1);
4347 evalue *ev;
4348 evalue rel;
4349
4350 value_init(rel.d);
4351 value_set_si(rel.d, 0);
4352 rel.x.p = new_enode(relation, 2, 0);
4353 value_clear(rel.x.p->arr[1].d);
4354 rel.x.p->arr[1] = p->arr[2*j+1];
4355 ev = &rel.x.p->arr[0];
4356 value_set_si(ev->d, 0);
4357 ev->x.p = new_enode(fractional, 3, -1);
4358 evalue_set_si(&ev->x.p->arr[1], 0, 1);
4359 evalue_set_si(&ev->x.p->arr[2], 1, 1);
4360 value_clear(ev->x.p->arr[0].d);
4361 ev->x.p->arr[0] = *arg;
4362 free(arg);
4363
4364 p->arr[2*j+1] = rel;
4365 }
4366 }
4367 value_clear(gcd);
4368
4369 Matrix_Free(eq);
4370 Matrix_Free(inv);
4371 }
4372
4373 /* Computes
4374 *
4375 * \sum_{i=0}^n c_i/d X^i
4376 *
4377 * where d is the last element in the vector c.
4378 */
4379 evalue *evalue_polynomial(Vector *c, const evalue* X)
4380 {
4381 unsigned dim = c->Size-2;
4382 evalue EC;
4383 evalue *EP = ALLOC(evalue);
4384 int i;
4385
4386 value_init(EP->d);
4387
4388 if (EVALUE_IS_ZERO(*X) || dim == 0) {
4389 evalue_set(EP, c->p[0], c->p[dim+1]);
4390 reduce_constant(EP);
4391 return EP;
4392 }
4393
4394 evalue_set(EP, c->p[dim], c->p[dim+1]);
4395
4396 value_init(EC.d);
4397 evalue_set(&EC, c->p[dim], c->p[dim+1]);
4398
4399 for (i = dim-1; i >= 0; --i) {
4400 emul(X, EP);
4401 value_assign(EC.x.n, c->p[i]);
4402 eadd(&EC, EP);
4403 }
4404 free_evalue_refs(&EC);
4405 return EP;
4406 }
4407
4408 /* Create an evalue from an array of pairs of domains and evalues. */
4409 evalue *evalue_from_section_array(struct evalue_section *s, int n)
4410 {
4411 int i;
4412 evalue *res;
4413
4414 res = ALLOC(evalue);
4415 value_init(res->d);
4416
4417 if (n == 0)
4418 evalue_set_si(res, 0, 1);
4419 else {
4420 value_set_si(res->d, 0);
4421 res->x.p = new_enode(partition, 2*n, s[0].D->Dimension);
4422 for (i = 0; i < n; ++i) {
4423 EVALUE_SET_DOMAIN(res->x.p->arr[2*i], s[i].D);
4424 value_clear(res->x.p->arr[2*i+1].d);
4425 res->x.p->arr[2*i+1] = *s[i].E;
4426 free(s[i].E);
4427 }
4428 }
4429 return res;
4430 }
4431
4432 /* shift variables (>= first, 0-based) in polynomial n up (may be negative) */
4433 void evalue_shift_variables(evalue *e, int first, int n)
4434 {
4435 int i;
4436 if (value_notzero_p(e->d))
4437 return;
4438 assert(e->x.p->type == polynomial ||
4439 e->x.p->type == flooring ||
4440 e->x.p->type == fractional);
4441 if (e->x.p->type == polynomial && e->x.p->pos >= first+1) {
4442 assert(e->x.p->pos + n >= 1);
4443 e->x.p->pos += n;
4444 }
4445 for (i = 0; i < e->x.p->size; ++i)
4446 evalue_shift_variables(&e->x.p->arr[i], first, n);
4447 }
4448
4449 static const evalue *outer_floor(evalue *e, const evalue *outer)
4450 {
4451 int i;
4452
4453 if (value_notzero_p(e->d))
4454 return outer;
4455 switch (e->x.p->type) {
4456 case flooring:
4457 if (!outer || evalue_level_cmp(outer, &e->x.p->arr[0]) > 0)
4458 return &e->x.p->arr[0];
4459 else
4460 return outer;
4461 case polynomial:
4462 case fractional:
4463 case relation:
4464 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
4465 outer = outer_floor(&e->x.p->arr[i], outer);
4466 return outer;
4467 case partition:
4468 for (i = 0; i < e->x.p->size/2; ++i)
4469 outer = outer_floor(&e->x.p->arr[2*i+1], outer);
4470 return outer;
4471 default:
4472 assert(0);
4473 }
4474 }
4475
4476 /* Find and return outermost floor argument or NULL if e has no floors */
4477 evalue *evalue_outer_floor(evalue *e)
4478 {
4479 const evalue *floor = outer_floor(e, NULL);
4480 return floor ? evalue_dup(floor): NULL;
4481 }
4482
4483 static void evalue_set_to_zero(evalue *e)
4484 {
4485 if (EVALUE_IS_ZERO(*e))
4486 return;
4487 if (value_zero_p(e->d)) {
4488 free_evalue_refs(e);
4489 value_init(e->d);
4490 value_init(e->x.n);
4491 }
4492 value_set_si(e->d, 1);
4493 value_set_si(e->x.n, 0);
4494 }
4495
4496 /* Replace (outer) floor with argument "floor" by variable "var" (0-based)
4497 * and drop terms not containing the floor.
4498 * Returns true if e contains the floor.
4499 */
4500 int evalue_replace_floor(evalue *e, const evalue *floor, int var)
4501 {
4502 int i;
4503 int contains = 0;
4504 int reorder = 0;
4505
4506 if (value_notzero_p(e->d))
4507 return 0;
4508 switch (e->x.p->type) {
4509 case flooring:
4510 if (!eequal(floor, &e->x.p->arr[0]))
4511 return 0;
4512 e->x.p->type = polynomial;
4513 e->x.p->pos = 1 + var;
4514 e->x.p->size--;
4515 free_evalue_refs(&e->x.p->arr[0]);
4516 for (i = 0; i < e->x.p->size; ++i)
4517 e->x.p->arr[i] = e->x.p->arr[i+1];
4518 evalue_set_to_zero(&e->x.p->arr[0]);
4519 return 1;
4520 case polynomial:
4521 case fractional:
4522 case relation:
4523 for (i = type_offset(e->x.p); i < e->x.p->size; ++i) {
4524 int c = evalue_replace_floor(&e->x.p->arr[i], floor, var);
4525 contains |= c;
4526 if (!c)
4527 evalue_set_to_zero(&e->x.p->arr[i]);
4528 if (c && !reorder && evalue_level_cmp(&e->x.p->arr[i], e) < 0)
4529 reorder = 1;
4530 }
4531 evalue_reduce_size(e);
4532 if (reorder)
4533 evalue_reorder_terms(e);
4534 return contains;
4535 case partition:
4536 default:
4537 assert(0);
4538 }
4539 }
4540
4541 /* Replace (outer) floor with argument "floor" by variable zero */
4542 void evalue_drop_floor(evalue *e, const evalue *floor)
4543 {
4544 int i;
4545 enode *p;
4546
4547 if (value_notzero_p(e->d))
4548 return;
4549 switch (e->x.p->type) {
4550 case flooring:
4551 if (!eequal(floor, &e->x.p->arr[0]))
4552 return;
4553 p = e->x.p;
4554 free_evalue_refs(&p->arr[0]);
4555 for (i = 2; i < p->size; ++i)
4556 free_evalue_refs(&p->arr[i]);
4557 value_clear(e->d);
4558 *e = p->arr[1];
4559 free(p);
4560 break;
4561 case polynomial:
4562 case fractional:
4563 case relation:
4564 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
4565 evalue_drop_floor(&e->x.p->arr[i], floor);
4566 evalue_reduce_size(e);
4567 break;
4568 case partition:
4569 default:
4570 assert(0);
4571 }
4572 }
lattice point enumeration library