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