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