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