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