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