?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
bernstein: bernsteinExpansion: accept list of polynomials
[barvinok.git] / evalue.c
blob0c27ef3c4ff583547af27131dfb18bda7037be77
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,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(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 size_t value_size(Value v) {
2356 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2357 * sizeof(v[0]._mp_d[0]);
2358 }
2359
2360 size_t domain_size(Polyhedron *D)
2361 {
2362 int i, j;
2363 size_t s = sizeof(*D);
2364
2365 for (i = 0; i < D->NbConstraints; ++i)
2366 for (j = 0; j < D->Dimension+2; ++j)
2367 s += value_size(D->Constraint[i][j]);
2368
2369 /*
2370 for (i = 0; i < D->NbRays; ++i)
2371 for (j = 0; j < D->Dimension+2; ++j)
2372 s += value_size(D->Ray[i][j]);
2373 */
2374
2375 return D->next ? s+domain_size(D->next) : s;
2376 }
2377
2378 size_t enode_size(enode *p) {
2379 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2380 int i;
2381
2382 if (p->type == partition)
2383 for (i = 0; i < p->size/2; ++i) {
2384 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2385 s += evalue_size(&p->arr[2*i+1]);
2386 }
2387 else
2388 for (i = 0; i < p->size; ++i) {
2389 s += evalue_size(&p->arr[i]);
2390 }
2391 return s;
2392 }
2393
2394 size_t evalue_size(evalue *e)
2395 {
2396 size_t s = sizeof(*e);
2397 s += value_size(e->d);
2398 if (value_notzero_p(e->d))
2399 s += value_size(e->x.n);
2400 else
2401 s += enode_size(e->x.p);
2402 return s;
2403 }
2404
2405 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2406 {
2407 evalue *found = NULL;
2408 evalue offset;
2409 evalue copy;
2410 int i;
2411
2412 if (value_pos_p(e->d) || e->x.p->type != fractional)
2413 return NULL;
2414
2415 value_init(offset.d);
2416 value_init(offset.x.n);
2417 poly_denom(&e->x.p->arr[0], &offset.d);
2418 value_lcm(m, offset.d, &offset.d);
2419 value_set_si(offset.x.n, 1);
2420
2421 value_init(copy.d);
2422 evalue_copy(&copy, cst);
2423
2424 eadd(&offset, cst);
2425 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2426
2427 if (eequal(base, &e->x.p->arr[0]))
2428 found = &e->x.p->arr[0];
2429 else {
2430 value_set_si(offset.x.n, -2);
2431
2432 eadd(&offset, cst);
2433 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2434
2435 if (eequal(base, &e->x.p->arr[0]))
2436 found = base;
2437 }
2438 free_evalue_refs(cst);
2439 free_evalue_refs(&offset);
2440 *cst = copy;
2441
2442 for (i = 1; !found && i < e->x.p->size; ++i)
2443 found = find_second(base, cst, &e->x.p->arr[i], m);
2444
2445 return found;
2446 }
2447
2448 static evalue *find_relation_pair(evalue *e)
2449 {
2450 int i;
2451 evalue *found = NULL;
2452
2453 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2454 return NULL;
2455
2456 if (e->x.p->type == fractional) {
2457 Value m;
2458 evalue *cst;
2459
2460 value_init(m);
2461 poly_denom(&e->x.p->arr[0], &m);
2462
2463 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2464 cst = &cst->x.p->arr[0])
2465 ;
2466
2467 for (i = 1; !found && i < e->x.p->size; ++i)
2468 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2469
2470 value_clear(m);
2471 }
2472
2473 i = e->x.p->type == relation;
2474 for (; !found && i < e->x.p->size; ++i)
2475 found = find_relation_pair(&e->x.p->arr[i]);
2476
2477 return found;
2478 }
2479
2480 void evalue_mod2relation(evalue *e) {
2481 evalue *d;
2482
2483 if (value_zero_p(e->d) && e->x.p->type == partition) {
2484 int i;
2485
2486 for (i = 0; i < e->x.p->size/2; ++i) {
2487 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2488 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2489 value_clear(e->x.p->arr[2*i].d);
2490 free_evalue_refs(&e->x.p->arr[2*i+1]);
2491 e->x.p->size -= 2;
2492 if (2*i < e->x.p->size) {
2493 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2494 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2495 }
2496 --i;
2497 }
2498 }
2499 if (e->x.p->size == 0) {
2500 free(e->x.p);
2501 evalue_set_si(e, 0, 1);
2502 }
2503
2504 return;
2505 }
2506
2507 while ((d = find_relation_pair(e)) != NULL) {
2508 evalue split;
2509 evalue *ev;
2510
2511 value_init(split.d);
2512 value_set_si(split.d, 0);
2513 split.x.p = new_enode(relation, 3, 0);
2514 evalue_set_si(&split.x.p->arr[1], 1, 1);
2515 evalue_set_si(&split.x.p->arr[2], 1, 1);
2516
2517 ev = &split.x.p->arr[0];
2518 value_set_si(ev->d, 0);
2519 ev->x.p = new_enode(fractional, 3, -1);
2520 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2521 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2522 evalue_copy(&ev->x.p->arr[0], d);
2523
2524 emul(&split, e);
2525
2526 reduce_evalue(e);
2527
2528 free_evalue_refs(&split);
2529 }
2530 }
2531
2532 static int evalue_comp(const void * a, const void * b)
2533 {
2534 const evalue *e1 = *(const evalue **)a;
2535 const evalue *e2 = *(const evalue **)b;
2536 return ecmp(e1, e2);
2537 }
2538
2539 void evalue_combine(evalue *e)
2540 {
2541 evalue **evs;
2542 int i, k;
2543 enode *p;
2544 evalue tmp;
2545
2546 if (value_notzero_p(e->d) || e->x.p->type != partition)
2547 return;
2548
2549 NALLOC(evs, e->x.p->size/2);
2550 for (i = 0; i < e->x.p->size/2; ++i)
2551 evs[i] = &e->x.p->arr[2*i+1];
2552 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2553 p = new_enode(partition, e->x.p->size, e->x.p->pos);
2554 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2555 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2556 value_clear(p->arr[2*k].d);
2557 value_clear(p->arr[2*k+1].d);
2558 p->arr[2*k] = *(evs[i]-1);
2559 p->arr[2*k+1] = *(evs[i]);
2560 ++k;
2561 } else {
2562 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2563 Polyhedron *L = D;
2564
2565 value_clear((evs[i]-1)->d);
2566
2567 while (L->next)
2568 L = L->next;
2569 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2570 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2571 free_evalue_refs(evs[i]);
2572 }
2573 }
2574
2575 for (i = 2*k ; i < p->size; ++i)
2576 value_clear(p->arr[i].d);
2577
2578 free(evs);
2579 free(e->x.p);
2580 p->size = 2*k;
2581 e->x.p = p;
2582
2583 for (i = 0; i < e->x.p->size/2; ++i) {
2584 Polyhedron *H;
2585 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2586 continue;
2587 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2588 if (H == NULL)
2589 continue;
2590 for (k = 0; k < e->x.p->size/2; ++k) {
2591 Polyhedron *D, *N, **P;
2592 if (i == k)
2593 continue;
2594 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2595 D = *P;
2596 if (D == NULL)
2597 continue;
2598 for (; D; D = N) {
2599 *P = D;
2600 N = D->next;
2601 if (D->NbEq <= H->NbEq) {
2602 P = &D->next;
2603 continue;
2604 }
2605
2606 value_init(tmp.d);
2607 tmp.x.p = new_enode(partition, 2, e->x.p->pos);
2608 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2609 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2610 reduce_evalue(&tmp);
2611 if (value_notzero_p(tmp.d) ||
2612 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2613 P = &D->next;
2614 else {
2615 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2616 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2617 *P = NULL;
2618 }
2619 free_evalue_refs(&tmp);
2620 }
2621 }
2622 Polyhedron_Free(H);
2623 }
2624
2625 for (i = 0; i < e->x.p->size/2; ++i) {
2626 Polyhedron *H, *E;
2627 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2628 if (!D) {
2629 value_clear(e->x.p->arr[2*i].d);
2630 free_evalue_refs(&e->x.p->arr[2*i+1]);
2631 e->x.p->size -= 2;
2632 if (2*i < e->x.p->size) {
2633 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2634 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2635 }
2636 --i;
2637 continue;
2638 }
2639 if (!D->next)
2640 continue;
2641 H = DomainConvex(D, 0);
2642 E = DomainDifference(H, D, 0);
2643 Domain_Free(D);
2644 D = DomainDifference(H, E, 0);
2645 Domain_Free(H);
2646 Domain_Free(E);
2647 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2648 }
2649 }
2650
2651 /* May change coefficients to become non-standard if fiddle is set
2652 * => reduce p afterwards to correct
2653 */
2654 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2655 Matrix **R, int fiddle)
2656 {
2657 Polyhedron *I, *H;
2658 evalue *pp;
2659 unsigned dim = D->Dimension;
2660 Matrix *T = Matrix_Alloc(2, dim+1);
2661 Value twice;
2662 value_init(twice);
2663 assert(T);
2664
2665 assert(p->type == fractional);
2666 pp = &p->arr[0];
2667 value_set_si(T->p[1][dim], 1);
2668 poly_denom(pp, d);
2669 while (value_zero_p(pp->d)) {
2670 assert(pp->x.p->type == polynomial);
2671 assert(pp->x.p->size == 2);
2672 assert(value_notzero_p(pp->x.p->arr[1].d));
2673 value_division(T->p[0][pp->x.p->pos-1], *d, pp->x.p->arr[1].d);
2674 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2675 if (fiddle && value_gt(twice, pp->x.p->arr[1].d))
2676 value_subtract(pp->x.p->arr[1].x.n,
2677 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2678 value_multiply(T->p[0][pp->x.p->pos-1],
2679 T->p[0][pp->x.p->pos-1], pp->x.p->arr[1].x.n);
2680 pp = &pp->x.p->arr[0];
2681 }
2682 value_division(T->p[0][dim], *d, pp->d);
2683 value_multiply(T->p[0][dim], T->p[0][dim], pp->x.n);
2684 I = DomainImage(D, T, 0);
2685 H = DomainConvex(I, 0);
2686 Domain_Free(I);
2687 *R = T;
2688
2689 value_clear(twice);
2690
2691 return H;
2692 }
2693
2694 int evalue_range_reduction_in_domain(evalue *e, Polyhedron *D)
2695 {
2696 int i;
2697 enode *p;
2698 Value d, min, max;
2699 int r = 0;
2700 Polyhedron *I;
2701 Matrix *T;
2702 int bounded;
2703
2704 if (value_notzero_p(e->d))
2705 return r;
2706
2707 p = e->x.p;
2708
2709 if (p->type == relation) {
2710 int equal;
2711 value_init(d);
2712 value_init(min);
2713 value_init(max);
2714
2715 I = polynomial_projection(p->arr[0].x.p, D, &d, &T, 1);
2716 bounded = line_minmax(I, &min, &max); /* frees I */
2717 equal = value_eq(min, max);
2718 mpz_cdiv_q(min, min, d);
2719 mpz_fdiv_q(max, max, d);
2720
2721 if (bounded && value_gt(min, max)) {
2722 /* Never zero */
2723 if (p->size == 3) {
2724 value_clear(e->d);
2725 *e = p->arr[2];
2726 } else {
2727 evalue_set_si(e, 0, 1);
2728 r = 1;
2729 }
2730 free_evalue_refs(&(p->arr[1]));
2731 free_evalue_refs(&(p->arr[0]));
2732 free(p);
2733 value_clear(d);
2734 value_clear(min);
2735 value_clear(max);
2736 Matrix_Free(T);
2737 return r ? r : evalue_range_reduction_in_domain(e, D);
2738 } else if (bounded && equal) {
2739 /* Always zero */
2740 if (p->size == 3)
2741 free_evalue_refs(&(p->arr[2]));
2742 value_clear(e->d);
2743 *e = p->arr[1];
2744 free_evalue_refs(&(p->arr[0]));
2745 free(p);
2746 value_clear(d);
2747 value_clear(min);
2748 value_clear(max);
2749 Matrix_Free(T);
2750 return evalue_range_reduction_in_domain(e, D);
2751 } else if (bounded && value_eq(min, max)) {
2752 /* zero for a single value */
2753 Polyhedron *E;
2754 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2755 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2756 value_multiply(min, min, d);
2757 value_subtract(M->p[0][D->Dimension+1],
2758 M->p[0][D->Dimension+1], min);
2759 E = DomainAddConstraints(D, M, 0);
2760 value_clear(d);
2761 value_clear(min);
2762 value_clear(max);
2763 Matrix_Free(T);
2764 Matrix_Free(M);
2765 r = evalue_range_reduction_in_domain(&p->arr[1], E);
2766 if (p->size == 3)
2767 r |= evalue_range_reduction_in_domain(&p->arr[2], D);
2768 Domain_Free(E);
2769 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2770 return r;
2771 }
2772
2773 value_clear(d);
2774 value_clear(min);
2775 value_clear(max);
2776 Matrix_Free(T);
2777 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2778 }
2779
2780 i = p->type == relation ? 1 :
2781 p->type == fractional ? 1 : 0;
2782 for (; i<p->size; i++)
2783 r |= evalue_range_reduction_in_domain(&p->arr[i], D);
2784
2785 if (p->type != fractional) {
2786 if (r && p->type == polynomial) {
2787 evalue f;
2788 value_init(f.d);
2789 value_set_si(f.d, 0);
2790 f.x.p = new_enode(polynomial, 2, p->pos);
2791 evalue_set_si(&f.x.p->arr[0], 0, 1);
2792 evalue_set_si(&f.x.p->arr[1], 1, 1);
2793 reorder_terms_about(p, &f);
2794 value_clear(e->d);
2795 *e = p->arr[0];
2796 free(p);
2797 }
2798 return r;
2799 }
2800
2801 value_init(d);
2802 value_init(min);
2803 value_init(max);
2804 I = polynomial_projection(p, D, &d, &T, 1);
2805 Matrix_Free(T);
2806 bounded = line_minmax(I, &min, &max); /* frees I */
2807 mpz_fdiv_q(min, min, d);
2808 mpz_fdiv_q(max, max, d);
2809 value_subtract(d, max, min);
2810
2811 if (bounded && value_eq(min, max)) {
2812 evalue inc;
2813 value_init(inc.d);
2814 value_init(inc.x.n);
2815 value_set_si(inc.d, 1);
2816 value_oppose(inc.x.n, min);
2817 eadd(&inc, &p->arr[0]);
2818 reorder_terms_about(p, &p->arr[0]); /* frees arr[0] */
2819 value_clear(e->d);
2820 *e = p->arr[1];
2821 free(p);
2822 free_evalue_refs(&inc);
2823 r = 1;
2824 } else if (bounded && value_one_p(d) && p->size > 3) {
2825 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
2826 * See pages 199-200 of PhD thesis.
2827 */
2828 evalue rem;
2829 evalue inc;
2830 evalue t;
2831 evalue f;
2832 evalue factor;
2833 value_init(rem.d);
2834 value_set_si(rem.d, 0);
2835 rem.x.p = new_enode(fractional, 3, -1);
2836 evalue_copy(&rem.x.p->arr[0], &p->arr[0]);
2837 value_clear(rem.x.p->arr[1].d);
2838 value_clear(rem.x.p->arr[2].d);
2839 rem.x.p->arr[1] = p->arr[1];
2840 rem.x.p->arr[2] = p->arr[2];
2841 for (i = 3; i < p->size; ++i)
2842 p->arr[i-2] = p->arr[i];
2843 p->size -= 2;
2844
2845 value_init(inc.d);
2846 value_init(inc.x.n);
2847 value_set_si(inc.d, 1);
2848 value_oppose(inc.x.n, min);
2849
2850 value_init(t.d);
2851 evalue_copy(&t, &p->arr[0]);
2852 eadd(&inc, &t);
2853
2854 value_init(f.d);
2855 value_set_si(f.d, 0);
2856 f.x.p = new_enode(fractional, 3, -1);
2857 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2858 evalue_set_si(&f.x.p->arr[1], 1, 1);
2859 evalue_set_si(&f.x.p->arr[2], 2, 1);
2860
2861 value_init(factor.d);
2862 evalue_set_si(&factor, -1, 1);
2863 emul(&t, &factor);
2864
2865 eadd(&f, &factor);
2866 emul(&t, &factor);
2867
2868 value_clear(f.x.p->arr[1].x.n);
2869 value_clear(f.x.p->arr[2].x.n);
2870 evalue_set_si(&f.x.p->arr[1], 0, 1);
2871 evalue_set_si(&f.x.p->arr[2], -1, 1);
2872 eadd(&f, &factor);
2873
2874 if (r) {
2875 reorder_terms(&rem);
2876 reorder_terms(e);
2877 }
2878
2879 emul(&factor, e);
2880 eadd(&rem, e);
2881
2882 free_evalue_refs(&inc);
2883 free_evalue_refs(&t);
2884 free_evalue_refs(&f);
2885 free_evalue_refs(&factor);
2886 free_evalue_refs(&rem);
2887
2888 evalue_range_reduction_in_domain(e, D);
2889
2890 r = 1;
2891 } else {
2892 _reduce_evalue(&p->arr[0], 0, 1);
2893 if (r)
2894 reorder_terms(e);
2895 }
2896
2897 value_clear(d);
2898 value_clear(min);
2899 value_clear(max);
2900
2901 return r;
2902 }
2903
2904 void evalue_range_reduction(evalue *e)
2905 {
2906 int i;
2907 if (value_notzero_p(e->d) || e->x.p->type != partition)
2908 return;
2909
2910 for (i = 0; i < e->x.p->size/2; ++i)
2911 if (evalue_range_reduction_in_domain(&e->x.p->arr[2*i+1],
2912 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
2913 reduce_evalue(&e->x.p->arr[2*i+1]);
2914
2915 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2916 free_evalue_refs(&e->x.p->arr[2*i+1]);
2917 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
2918 value_clear(e->x.p->arr[2*i].d);
2919 e->x.p->size -= 2;
2920 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2921 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2922 --i;
2923 }
2924 }
2925 }
2926
2927 /* Frees EP
2928 */
2929 Enumeration* partition2enumeration(evalue *EP)
2930 {
2931 int i;
2932 Enumeration *en, *res = NULL;
2933
2934 if (EVALUE_IS_ZERO(*EP)) {
2935 free(EP);
2936 return res;
2937 }
2938
2939 for (i = 0; i < EP->x.p->size/2; ++i) {
2940 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[2*i])->Dimension);
2941 en = (Enumeration *)malloc(sizeof(Enumeration));
2942 en->next = res;
2943 res = en;
2944 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2945 value_clear(EP->x.p->arr[2*i].d);
2946 res->EP = EP->x.p->arr[2*i+1];
2947 }
2948 free(EP->x.p);
2949 value_clear(EP->d);
2950 free(EP);
2951 return res;
2952 }
2953
2954 int evalue_frac2floor_in_domain3(evalue *e, Polyhedron *D, int shift)
2955 {
2956 enode *p;
2957 int r = 0;
2958 int i;
2959 Polyhedron *I;
2960 Matrix *T;
2961 Value d, min;
2962 evalue fl;
2963
2964 if (value_notzero_p(e->d))
2965 return r;
2966
2967 p = e->x.p;
2968
2969 i = p->type == relation ? 1 :
2970 p->type == fractional ? 1 : 0;
2971 for (; i<p->size; i++)
2972 r |= evalue_frac2floor_in_domain3(&p->arr[i], D, shift);
2973
2974 if (p->type != fractional) {
2975 if (r && p->type == polynomial) {
2976 evalue f;
2977 value_init(f.d);
2978 value_set_si(f.d, 0);
2979 f.x.p = new_enode(polynomial, 2, p->pos);
2980 evalue_set_si(&f.x.p->arr[0], 0, 1);
2981 evalue_set_si(&f.x.p->arr[1], 1, 1);
2982 reorder_terms_about(p, &f);
2983 value_clear(e->d);
2984 *e = p->arr[0];
2985 free(p);
2986 }
2987 return r;
2988 }
2989
2990 if (shift) {
2991 value_init(d);
2992 I = polynomial_projection(p, D, &d, &T, 0);
2993
2994 /*
2995 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2996 */
2997
2998 assert(I->NbEq == 0); /* Should have been reduced */
2999
3000 /* Find minimum */
3001 for (i = 0; i < I->NbConstraints; ++i)
3002 if (value_pos_p(I->Constraint[i][1]))
3003 break;
3004
3005 if (i < I->NbConstraints) {
3006 value_init(min);
3007 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
3008 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
3009 if (value_neg_p(min)) {
3010 evalue offset;
3011 mpz_fdiv_q(min, min, d);
3012 value_init(offset.d);
3013 value_set_si(offset.d, 1);
3014 value_init(offset.x.n);
3015 value_oppose(offset.x.n, min);
3016 eadd(&offset, &p->arr[0]);
3017 free_evalue_refs(&offset);
3018 }
3019 value_clear(min);
3020 }
3021
3022 Polyhedron_Free(I);
3023 Matrix_Free(T);
3024 value_clear(d);
3025 }
3026
3027 value_init(fl.d);
3028 value_set_si(fl.d, 0);
3029 fl.x.p = new_enode(flooring, 3, -1);
3030 evalue_set_si(&fl.x.p->arr[1], 0, 1);
3031 evalue_set_si(&fl.x.p->arr[2], -1, 1);
3032 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
3033
3034 eadd(&fl, &p->arr[0]);
3035 reorder_terms_about(p, &p->arr[0]);
3036 value_clear(e->d);
3037 *e = p->arr[1];
3038 free(p);
3039 free_evalue_refs(&fl);
3040
3041 return 1;
3042 }
3043
3044 int evalue_frac2floor_in_domain(evalue *e, Polyhedron *D)
3045 {
3046 return evalue_frac2floor_in_domain3(e, D, 1);
3047 }
3048
3049 void evalue_frac2floor2(evalue *e, int shift)
3050 {
3051 int i;
3052 if (value_notzero_p(e->d) || e->x.p->type != partition)
3053 return;
3054
3055 for (i = 0; i < e->x.p->size/2; ++i)
3056 if (evalue_frac2floor_in_domain3(&e->x.p->arr[2*i+1],
3057 EVALUE_DOMAIN(e->x.p->arr[2*i]), shift))
3058 reduce_evalue(&e->x.p->arr[2*i+1]);
3059 }
3060
3061 void evalue_frac2floor(evalue *e)
3062 {
3063 evalue_frac2floor2(e, 1);
3064 }
3065
3066 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
3067 Vector *row)
3068 {
3069 int nr, nc;
3070 int i;
3071 int nparam = D->Dimension - nvar;
3072
3073 if (C == 0) {
3074 nr = D->NbConstraints + 2;
3075 nc = D->Dimension + 2 + 1;
3076 C = Matrix_Alloc(nr, nc);
3077 for (i = 0; i < D->NbConstraints; ++i) {
3078 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
3079 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3080 D->Dimension + 1 - nvar);
3081 }
3082 } else {
3083 Matrix *oldC = C;
3084 nr = C->NbRows + 2;
3085 nc = C->NbColumns + 1;
3086 C = Matrix_Alloc(nr, nc);
3087 for (i = 0; i < oldC->NbRows; ++i) {
3088 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
3089 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3090 oldC->NbColumns - 1 - nvar);
3091 }
3092 }
3093 value_set_si(C->p[nr-2][0], 1);
3094 value_set_si(C->p[nr-2][1 + nvar], 1);
3095 value_set_si(C->p[nr-2][nc - 1], -1);
3096
3097 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
3098 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
3099 1 + nparam);
3100
3101 return C;
3102 }
3103
3104 static void floor2frac_r(evalue *e, int nvar)
3105 {
3106 enode *p;
3107 int i;
3108 evalue f;
3109 evalue *pp;
3110
3111 if (value_notzero_p(e->d))
3112 return;
3113
3114 p = e->x.p;
3115
3116 assert(p->type == flooring);
3117 for (i = 1; i < p->size; i++)
3118 floor2frac_r(&p->arr[i], nvar);
3119
3120 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
3121 assert(pp->x.p->type == polynomial);
3122 pp->x.p->pos -= nvar;
3123 }
3124
3125 value_init(f.d);
3126 value_set_si(f.d, 0);
3127 f.x.p = new_enode(fractional, 3, -1);
3128 evalue_set_si(&f.x.p->arr[1], 0, 1);
3129 evalue_set_si(&f.x.p->arr[2], -1, 1);
3130 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3131
3132 eadd(&f, &p->arr[0]);
3133 reorder_terms_about(p, &p->arr[0]);
3134 value_clear(e->d);
3135 *e = p->arr[1];
3136 free(p);
3137 free_evalue_refs(&f);
3138 }
3139
3140 /* Convert flooring back to fractional and shift position
3141 * of the parameters by nvar
3142 */
3143 static void floor2frac(evalue *e, int nvar)
3144 {
3145 floor2frac_r(e, nvar);
3146 reduce_evalue(e);
3147 }
3148
3149 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
3150 {
3151 evalue *t;
3152 int nparam = D->Dimension - nvar;
3153
3154 if (C != 0) {
3155 C = Matrix_Copy(C);
3156 D = Constraints2Polyhedron(C, 0);
3157 Matrix_Free(C);
3158 }
3159
3160 t = barvinok_enumerate_e(D, 0, nparam, 0);
3161
3162 /* Double check that D was not unbounded. */
3163 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3164
3165 if (C != 0)
3166 Polyhedron_Free(D);
3167
3168 return t;
3169 }
3170
3171 evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3172 Matrix *C)
3173 {
3174 Vector *row = NULL;
3175 int i;
3176 evalue *res;
3177 Matrix *origC;
3178 evalue *factor = NULL;
3179 evalue cum;
3180
3181 if (EVALUE_IS_ZERO(*e))
3182 return 0;
3183
3184 if (D->next) {
3185 Polyhedron *DD = Disjoint_Domain(D, 0, 0);
3186 Polyhedron *Q;
3187
3188 Q = DD;
3189 DD = Q->next;
3190 Q->next = 0;
3191
3192 res = esum_over_domain(e, nvar, Q, C);
3193 Polyhedron_Free(Q);
3194
3195 for (Q = DD; Q; Q = DD) {
3196 evalue *t;
3197
3198 DD = Q->next;
3199 Q->next = 0;
3200
3201 t = esum_over_domain(e, nvar, Q, C);
3202 Polyhedron_Free(Q);
3203
3204 if (!res)
3205 res = t;
3206 else if (t) {
3207 eadd(t, res);
3208 free_evalue_refs(t);
3209 free(t);
3210 }
3211 }
3212 return res;
3213 }
3214
3215 if (value_notzero_p(e->d)) {
3216 evalue *t;
3217
3218 t = esum_over_domain_cst(nvar, D, C);
3219
3220 if (!EVALUE_IS_ONE(*e))
3221 emul(e, t);
3222
3223 return t;
3224 }
3225
3226 switch (e->x.p->type) {
3227 case flooring: {
3228 evalue *pp = &e->x.p->arr[0];
3229
3230 if (pp->x.p->pos > nvar) {
3231 /* remainder is independent of the summated vars */
3232 evalue f;
3233 evalue *t;
3234
3235 value_init(f.d);
3236 evalue_copy(&f, e);
3237 floor2frac(&f, nvar);
3238
3239 t = esum_over_domain_cst(nvar, D, C);
3240
3241 emul(&f, t);
3242
3243 free_evalue_refs(&f);
3244
3245 return t;
3246 }
3247
3248 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3249 poly_denom(pp, &row->p[1 + nvar]);
3250 value_set_si(row->p[0], 1);
3251 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3252 pp = &pp->x.p->arr[0]) {
3253 int pos;
3254 assert(pp->x.p->type == polynomial);
3255 pos = pp->x.p->pos;
3256 if (pos >= 1 + nvar)
3257 ++pos;
3258 value_assign(row->p[pos], row->p[1+nvar]);
3259 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3260 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3261 }
3262 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3263 value_division(row->p[1 + D->Dimension + 1],
3264 row->p[1 + D->Dimension + 1],
3265 pp->d);
3266 value_multiply(row->p[1 + D->Dimension + 1],
3267 row->p[1 + D->Dimension + 1],
3268 pp->x.n);
3269 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3270 break;
3271 }
3272 case polynomial: {
3273 int pos = e->x.p->pos;
3274
3275 if (pos > nvar) {
3276 factor = ALLOC(evalue);
3277 value_init(factor->d);
3278 value_set_si(factor->d, 0);
3279 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3280 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3281 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3282 break;
3283 }
3284
3285 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3286 for (i = 0; i < D->NbRays; ++i)
3287 if (value_notzero_p(D->Ray[i][pos]))
3288 break;
3289 assert(i < D->NbRays);
3290 if (value_neg_p(D->Ray[i][pos])) {
3291 factor = ALLOC(evalue);
3292 value_init(factor->d);
3293 evalue_set_si(factor, -1, 1);
3294 }
3295 value_set_si(row->p[0], 1);
3296 value_set_si(row->p[pos], 1);
3297 value_set_si(row->p[1 + nvar], -1);
3298 break;
3299 }
3300 default:
3301 assert(0);
3302 }
3303
3304 i = type_offset(e->x.p);
3305
3306 res = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3307 ++i;
3308
3309 if (factor) {
3310 value_init(cum.d);
3311 evalue_copy(&cum, factor);
3312 }
3313
3314 origC = C;
3315 for (; i < e->x.p->size; ++i) {
3316 evalue *t;
3317 if (row) {
3318 Matrix *prevC = C;
3319 C = esum_add_constraint(nvar, D, C, row);
3320 if (prevC != origC)
3321 Matrix_Free(prevC);
3322 }
3323 /*
3324 if (row)
3325 Vector_Print(stderr, P_VALUE_FMT, row);
3326 if (C)
3327 Matrix_Print(stderr, P_VALUE_FMT, C);
3328 */
3329 t = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3330
3331 if (t && factor)
3332 emul(&cum, t);
3333
3334 if (!res)
3335 res = t;
3336 else if (t) {
3337 eadd(t, res);
3338 free_evalue_refs(t);
3339 free(t);
3340 }
3341 if (factor && i+1 < e->x.p->size)
3342 emul(factor, &cum);
3343 }
3344 if (C != origC)
3345 Matrix_Free(C);
3346
3347 if (factor) {
3348 free_evalue_refs(factor);
3349 free_evalue_refs(&cum);
3350 free(factor);
3351 }
3352
3353 if (row)
3354 Vector_Free(row);
3355
3356 reduce_evalue(res);
3357
3358 return res;
3359 }
3360
3361 evalue *esum(evalue *e, int nvar)
3362 {
3363 int i;
3364 evalue *res = ALLOC(evalue);
3365 value_init(res->d);
3366
3367 assert(nvar >= 0);
3368 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3369 evalue_copy(res, e);
3370 return res;
3371 }
3372
3373 evalue_set_si(res, 0, 1);
3374
3375 assert(value_zero_p(e->d));
3376 assert(e->x.p->type == partition);
3377
3378 for (i = 0; i < e->x.p->size/2; ++i) {
3379 evalue *t;
3380 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3381 EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
3382 eadd(t, res);
3383 free_evalue_refs(t);
3384 free(t);
3385 }
3386
3387 reduce_evalue(res);
3388
3389 return res;
3390 }
3391
3392 /* Initial silly implementation */
3393 void eor(evalue *e1, evalue *res)
3394 {
3395 evalue E;
3396 evalue mone;
3397 value_init(E.d);
3398 value_init(mone.d);
3399 evalue_set_si(&mone, -1, 1);
3400
3401 evalue_copy(&E, res);
3402 eadd(e1, &E);
3403 emul(e1, res);
3404 emul(&mone, res);
3405 eadd(&E, res);
3406
3407 free_evalue_refs(&E);
3408 free_evalue_refs(&mone);
3409 }
3410
3411 /* computes denominator of polynomial evalue
3412 * d should point to a value initialized to 1
3413 */
3414 void evalue_denom(const evalue *e, Value *d)
3415 {
3416 int i, offset;
3417
3418 if (value_notzero_p(e->d)) {
3419 value_lcm(*d, e->d, d);
3420 return;
3421 }
3422 assert(e->x.p->type == polynomial ||
3423 e->x.p->type == fractional ||
3424 e->x.p->type == flooring);
3425 offset = type_offset(e->x.p);
3426 for (i = e->x.p->size-1; i >= offset; --i)
3427 evalue_denom(&e->x.p->arr[i], d);
3428 }
3429
3430 /* Divides the evalue e by the integer n */
3431 void evalue_div(evalue * e, Value n)
3432 {
3433 int i, offset;
3434
3435 if (value_notzero_p(e->d)) {
3436 Value gc;
3437 value_init(gc);
3438 value_multiply(e->d, e->d, n);
3439 Gcd(e->x.n, e->d, &gc);
3440 if (value_notone_p(gc)) {
3441 value_division(e->d, e->d, gc);
3442 value_division(e->x.n, e->x.n, gc);
3443 }
3444 value_clear(gc);
3445 return;
3446 }
3447 if (e->x.p->type == partition) {
3448 for (i = 0; i < e->x.p->size/2; ++i)
3449 evalue_div(&e->x.p->arr[2*i+1], n);
3450 return;
3451 }
3452 offset = type_offset(e->x.p);
3453 for (i = e->x.p->size-1; i >= offset; --i)
3454 evalue_div(&e->x.p->arr[i], n);
3455 }
3456
3457 static void evalue_frac2polynomial_r(evalue *e, int *signs, int sign, int in_frac)
3458 {
3459 int i, offset;
3460 Value d;
3461 enode *p;
3462
3463 if (value_notzero_p(e->d)) {
3464 if (in_frac && sign * value_sign(e->x.n) < 0) {
3465 value_set_si(e->x.n, 0);
3466 value_set_si(e->d, 1);
3467 }
3468 return;
3469 }
3470
3471 if (e->x.p->type == polynomial) {
3472 sign *= signs[e->x.p->pos-1];
3473 }
3474 offset = type_offset(e->x.p);
3475 evalue_frac2polynomial_r(&e->x.p->arr[offset], signs, sign, in_frac);
3476 in_frac |= e->x.p->type == fractional;
3477 for (i = e->x.p->size-1; i > offset; --i)
3478 evalue_frac2polynomial_r(&e->x.p->arr[i], signs, sign, in_frac);
3479
3480 if (e->x.p->type != fractional)
3481 return;
3482
3483 /* replace { a/m } by (m-1)/m */
3484 value_init(d);
3485 value_set_si(d, 1);
3486 evalue_denom(&e->x.p->arr[0], &d);
3487 free_evalue_refs(&e->x.p->arr[0]);
3488 value_init(e->x.p->arr[0].d);
3489 value_init(e->x.p->arr[0].x.n);
3490 value_assign(e->x.p->arr[0].d, d);
3491 value_decrement(e->x.p->arr[0].x.n, d);
3492 value_clear(d);
3493
3494 p = e->x.p;
3495 reorder_terms_about(p, &p->arr[0]);
3496 value_clear(e->d);
3497 *e = p->arr[1];
3498 free(p);
3499 }
3500
3501 /* Approximate the evalue in fractional representation by a polynomial.
3502 * If sign > 0, the result is an upper bound;
3503 * if sign < 0, the resutl is a lower bound.
3504 */
3505 void evalue_frac2polynomial(evalue *e, int sign, unsigned MaxRays)
3506 {
3507 int i, j, k, dim;
3508 int *signs;
3509
3510 if (value_notzero_p(e->d))
3511 return;
3512 assert(e->x.p->type == partition);
3513 /* make sure all variables in the domains have a fixed sign */
3514 evalue_split_domains_into_orthants(e, MaxRays);
3515
3516 assert(e->x.p->size >= 2);
3517 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3518
3519 signs = alloca(sizeof(int) * dim);
3520
3521 for (i = 0; i < e->x.p->size/2; ++i) {
3522 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
3523 POL_ENSURE_VERTICES(D);
3524 for (j = 0; j < dim; ++j) {
3525 signs[j] = 0;
3526 for (k = 0; k < D->NbRays; ++k) {
3527 signs[j] = value_sign(D->Ray[k][1+j]);
3528 if (signs[j])
3529 break;
3530 }
3531 }
3532 evalue_frac2polynomial_r(&e->x.p->arr[2*i+1], signs, sign, 0);
3533 }
3534 }
3535
3536 /* Split the domains of e (which is assumed to be a partition)
3537 * such that each resulting domain lies entirely in one orthant.
3538 */
3539 void evalue_split_domains_into_orthants(evalue *e, unsigned MaxRays)
3540 {
3541 int i, dim;
3542 assert(value_zero_p(e->d));
3543 assert(e->x.p->type == partition);
3544 assert(e->x.p->size >= 2);
3545 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3546
3547 for (i = 0; i < dim; ++i) {
3548 evalue split;
3549 Matrix *C, *C2;
3550 C = Matrix_Alloc(1, 1 + dim + 1);
3551 value_set_si(C->p[0][0], 1);
3552 value_init(split.d);
3553 value_set_si(split.d, 0);
3554 split.x.p = new_enode(partition, 4, dim);
3555 value_set_si(C->p[0][1+i], 1);
3556 C2 = Matrix_Copy(C);
3557 EVALUE_SET_DOMAIN(split.x.p->arr[0], Constraints2Polyhedron(C2, MaxRays));
3558 Matrix_Free(C2);
3559 evalue_set_si(&split.x.p->arr[1], 1, 1);
3560 value_set_si(C->p[0][1+i], -1);
3561 value_set_si(C->p[0][1+dim], -1);
3562 EVALUE_SET_DOMAIN(split.x.p->arr[2], Constraints2Polyhedron(C, MaxRays));
3563 evalue_set_si(&split.x.p->arr[3], 1, 1);
3564 emul(&split, e);
3565 free_evalue_refs(&split);
3566 Matrix_Free(C);
3567 }
3568 }
lattice point enumeration library