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