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