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