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