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