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