?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
evalue.c: reorder_terms: fix typo
[barvinok.git] / evalue.c
blobc955b055b9a6adb4b19c3a4c741b87e152324d62
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, char **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);
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 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
3233 Vector *row)
3234 {
3235 int nr, nc;
3236 int i;
3237 int nparam = D->Dimension - nvar;
3238
3239 if (C == 0) {
3240 nr = D->NbConstraints + 2;
3241 nc = D->Dimension + 2 + 1;
3242 C = Matrix_Alloc(nr, nc);
3243 for (i = 0; i < D->NbConstraints; ++i) {
3244 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
3245 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3246 D->Dimension + 1 - nvar);
3247 }
3248 } else {
3249 Matrix *oldC = C;
3250 nr = C->NbRows + 2;
3251 nc = C->NbColumns + 1;
3252 C = Matrix_Alloc(nr, nc);
3253 for (i = 0; i < oldC->NbRows; ++i) {
3254 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
3255 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3256 oldC->NbColumns - 1 - nvar);
3257 }
3258 }
3259 value_set_si(C->p[nr-2][0], 1);
3260 value_set_si(C->p[nr-2][1 + nvar], 1);
3261 value_set_si(C->p[nr-2][nc - 1], -1);
3262
3263 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
3264 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
3265 1 + nparam);
3266
3267 return C;
3268 }
3269
3270 static void floor2frac_r(evalue *e, int nvar)
3271 {
3272 enode *p;
3273 int i;
3274 evalue f;
3275 evalue *pp;
3276
3277 if (value_notzero_p(e->d))
3278 return;
3279
3280 p = e->x.p;
3281
3282 assert(p->type == flooring);
3283 for (i = 1; i < p->size; i++)
3284 floor2frac_r(&p->arr[i], nvar);
3285
3286 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
3287 assert(pp->x.p->type == polynomial);
3288 pp->x.p->pos -= nvar;
3289 }
3290
3291 value_init(f.d);
3292 value_set_si(f.d, 0);
3293 f.x.p = new_enode(fractional, 3, -1);
3294 evalue_set_si(&f.x.p->arr[1], 0, 1);
3295 evalue_set_si(&f.x.p->arr[2], -1, 1);
3296 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3297
3298 eadd(&f, &p->arr[0]);
3299 reorder_terms_about(p, &p->arr[0]);
3300 value_clear(e->d);
3301 *e = p->arr[1];
3302 free(p);
3303 free_evalue_refs(&f);
3304 }
3305
3306 /* Convert flooring back to fractional and shift position
3307 * of the parameters by nvar
3308 */
3309 static void floor2frac(evalue *e, int nvar)
3310 {
3311 floor2frac_r(e, nvar);
3312 reduce_evalue(e);
3313 }
3314
3315 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
3316 {
3317 evalue *t;
3318 int nparam = D->Dimension - nvar;
3319
3320 if (C != 0) {
3321 C = Matrix_Copy(C);
3322 D = Constraints2Polyhedron(C, 0);
3323 Matrix_Free(C);
3324 }
3325
3326 t = barvinok_enumerate_e(D, 0, nparam, 0);
3327
3328 /* Double check that D was not unbounded. */
3329 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3330
3331 if (C != 0)
3332 Polyhedron_Free(D);
3333
3334 return t;
3335 }
3336
3337 static evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3338 int *signs, Matrix *C, unsigned MaxRays)
3339 {
3340 Vector *row = NULL;
3341 int i, offset;
3342 evalue *res;
3343 Matrix *origC;
3344 evalue *factor = NULL;
3345 evalue cum;
3346 int negate_odd = 0;
3347
3348 if (EVALUE_IS_ZERO(*e))
3349 return 0;
3350
3351 if (D->next) {
3352 Polyhedron *DD = Disjoint_Domain(D, 0, MaxRays);
3353 Polyhedron *Q;
3354
3355 Q = DD;
3356 DD = Q->next;
3357 Q->next = 0;
3358
3359 res = esum_over_domain(e, nvar, Q, signs, C, MaxRays);
3360 Polyhedron_Free(Q);
3361
3362 for (Q = DD; Q; Q = DD) {
3363 evalue *t;
3364
3365 DD = Q->next;
3366 Q->next = 0;
3367
3368 t = esum_over_domain(e, nvar, Q, signs, C, MaxRays);
3369 Polyhedron_Free(Q);
3370
3371 if (!res)
3372 res = t;
3373 else if (t) {
3374 eadd(t, res);
3375 free_evalue_refs(t);
3376 free(t);
3377 }
3378 }
3379 return res;
3380 }
3381
3382 if (value_notzero_p(e->d)) {
3383 evalue *t;
3384
3385 t = esum_over_domain_cst(nvar, D, C);
3386
3387 if (!EVALUE_IS_ONE(*e))
3388 emul(e, t);
3389
3390 return t;
3391 }
3392
3393 switch (e->x.p->type) {
3394 case flooring: {
3395 evalue *pp = &e->x.p->arr[0];
3396
3397 if (pp->x.p->pos > nvar) {
3398 /* remainder is independent of the summated vars */
3399 evalue f;
3400 evalue *t;
3401
3402 value_init(f.d);
3403 evalue_copy(&f, e);
3404 floor2frac(&f, nvar);
3405
3406 t = esum_over_domain_cst(nvar, D, C);
3407
3408 emul(&f, t);
3409
3410 free_evalue_refs(&f);
3411
3412 return t;
3413 }
3414
3415 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3416 poly_denom(pp, &row->p[1 + nvar]);
3417 value_set_si(row->p[0], 1);
3418 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3419 pp = &pp->x.p->arr[0]) {
3420 int pos;
3421 assert(pp->x.p->type == polynomial);
3422 pos = pp->x.p->pos;
3423 if (pos >= 1 + nvar)
3424 ++pos;
3425 value_assign(row->p[pos], row->p[1+nvar]);
3426 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3427 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3428 }
3429 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3430 value_division(row->p[1 + D->Dimension + 1],
3431 row->p[1 + D->Dimension + 1],
3432 pp->d);
3433 value_multiply(row->p[1 + D->Dimension + 1],
3434 row->p[1 + D->Dimension + 1],
3435 pp->x.n);
3436 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3437 break;
3438 }
3439 case polynomial: {
3440 int pos = e->x.p->pos;
3441
3442 if (pos > nvar) {
3443 factor = ALLOC(evalue);
3444 value_init(factor->d);
3445 value_set_si(factor->d, 0);
3446 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3447 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3448 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3449 break;
3450 }
3451
3452 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3453 negate_odd = signs[pos-1] < 0;
3454 value_set_si(row->p[0], 1);
3455 value_set_si(row->p[pos], 1);
3456 value_set_si(row->p[1 + nvar], -1);
3457 break;
3458 }
3459 default:
3460 assert(0);
3461 }
3462
3463 offset = type_offset(e->x.p);
3464
3465 res = esum_over_domain(&e->x.p->arr[offset], nvar, D, signs, C, MaxRays);
3466
3467 if (factor) {
3468 value_init(cum.d);
3469 evalue_copy(&cum, factor);
3470 }
3471
3472 origC = C;
3473 for (i = 1; offset+i < e->x.p->size; ++i) {
3474 evalue *t;
3475 if (row) {
3476 Matrix *prevC = C;
3477 C = esum_add_constraint(nvar, D, C, row);
3478 if (prevC != origC)
3479 Matrix_Free(prevC);
3480 }
3481 /*
3482 if (row)
3483 Vector_Print(stderr, P_VALUE_FMT, row);
3484 if (C)
3485 Matrix_Print(stderr, P_VALUE_FMT, C);
3486 */
3487 t = esum_over_domain(&e->x.p->arr[offset+i], nvar, D, signs, C, MaxRays);
3488
3489 if (t) {
3490 if (factor)
3491 emul(&cum, t);
3492 if (negate_odd && (i % 2))
3493 evalue_negate(t);
3494 }
3495
3496 if (!res)
3497 res = t;
3498 else if (t) {
3499 eadd(t, res);
3500 free_evalue_refs(t);
3501 free(t);
3502 }
3503 if (factor && offset+i+1 < e->x.p->size)
3504 emul(factor, &cum);
3505 }
3506 if (C != origC)
3507 Matrix_Free(C);
3508
3509 if (factor) {
3510 free_evalue_refs(factor);
3511 free_evalue_refs(&cum);
3512 free(factor);
3513 }
3514
3515 if (row)
3516 Vector_Free(row);
3517
3518 reduce_evalue(res);
3519
3520 return res;
3521 }
3522
3523 static void domain_signs(Polyhedron *D, int *signs)
3524 {
3525 int j, k;
3526
3527 POL_ENSURE_VERTICES(D);
3528 for (j = 0; j < D->Dimension; ++j) {
3529 signs[j] = 0;
3530 for (k = 0; k < D->NbRays; ++k) {
3531 signs[j] = value_sign(D->Ray[k][1+j]);
3532 if (signs[j])
3533 break;
3534 }
3535 }
3536 }
3537
3538 evalue *evalue_sum(evalue *e, int nvar, unsigned MaxRays)
3539 {
3540 int i, dim;
3541 int *signs;
3542 evalue *res = ALLOC(evalue);
3543 value_init(res->d);
3544
3545 assert(nvar >= 0);
3546 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3547 evalue_copy(res, e);
3548 return res;
3549 }
3550
3551 evalue_set_si(res, 0, 1);
3552
3553 assert(value_zero_p(e->d));
3554 assert(e->x.p->type == partition);
3555 evalue_split_domains_into_orthants(e, MaxRays);
3556
3557 assert(e->x.p->size >= 2);
3558 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3559
3560 signs = alloca(sizeof(int) * dim);
3561
3562 for (i = 0; i < e->x.p->size/2; ++i) {
3563 evalue *t;
3564 domain_signs(EVALUE_DOMAIN(e->x.p->arr[2*i]), signs);
3565 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3566 EVALUE_DOMAIN(e->x.p->arr[2*i]), signs, 0,
3567 MaxRays);
3568 eadd(t, res);
3569 free_evalue_refs(t);
3570 free(t);
3571 }
3572
3573 reduce_evalue(res);
3574
3575 return res;
3576 }
3577
3578 evalue *esum(evalue *e, int nvar)
3579 {
3580 return evalue_sum(e, nvar, 0);
3581 }
3582
3583 /* Initial silly implementation */
3584 void eor(evalue *e1, evalue *res)
3585 {
3586 evalue E;
3587 evalue mone;
3588 value_init(E.d);
3589 value_init(mone.d);
3590 evalue_set_si(&mone, -1, 1);
3591
3592 evalue_copy(&E, res);
3593 eadd(e1, &E);
3594 emul(e1, res);
3595 emul(&mone, res);
3596 eadd(&E, res);
3597
3598 free_evalue_refs(&E);
3599 free_evalue_refs(&mone);
3600 }
3601
3602 /* computes denominator of polynomial evalue
3603 * d should point to a value initialized to 1
3604 */
3605 void evalue_denom(const evalue *e, Value *d)
3606 {
3607 int i, offset;
3608
3609 if (value_notzero_p(e->d)) {
3610 value_lcm(*d, e->d, d);
3611 return;
3612 }
3613 assert(e->x.p->type == polynomial ||
3614 e->x.p->type == fractional ||
3615 e->x.p->type == flooring);
3616 offset = type_offset(e->x.p);
3617 for (i = e->x.p->size-1; i >= offset; --i)
3618 evalue_denom(&e->x.p->arr[i], d);
3619 }
3620
3621 /* Divides the evalue e by the integer n */
3622 void evalue_div(evalue * e, Value n)
3623 {
3624 int i, offset;
3625
3626 if (value_notzero_p(e->d)) {
3627 Value gc;
3628 value_init(gc);
3629 value_multiply(e->d, e->d, n);
3630 Gcd(e->x.n, e->d, &gc);
3631 if (value_notone_p(gc)) {
3632 value_division(e->d, e->d, gc);
3633 value_division(e->x.n, e->x.n, gc);
3634 }
3635 value_clear(gc);
3636 return;
3637 }
3638 if (e->x.p->type == partition) {
3639 for (i = 0; i < e->x.p->size/2; ++i)
3640 evalue_div(&e->x.p->arr[2*i+1], n);
3641 return;
3642 }
3643 offset = type_offset(e->x.p);
3644 for (i = e->x.p->size-1; i >= offset; --i)
3645 evalue_div(&e->x.p->arr[i], n);
3646 }
3647
3648 void evalue_negate(evalue *e)
3649 {
3650 int i, offset;
3651
3652 if (value_notzero_p(e->d)) {
3653 value_oppose(e->x.n, e->x.n);
3654 return;
3655 }
3656 if (e->x.p->type == partition) {
3657 for (i = 0; i < e->x.p->size/2; ++i)
3658 evalue_negate(&e->x.p->arr[2*i+1]);
3659 return;
3660 }
3661 offset = type_offset(e->x.p);
3662 for (i = e->x.p->size-1; i >= offset; --i)
3663 evalue_negate(&e->x.p->arr[i]);
3664 }
3665
3666 void evalue_add_constant(evalue *e, const Value cst)
3667 {
3668 int i;
3669
3670 if (value_zero_p(e->d)) {
3671 if (e->x.p->type == partition) {
3672 for (i = 0; i < e->x.p->size/2; ++i)
3673 evalue_add_constant(&e->x.p->arr[2*i+1], cst);
3674 return;
3675 }
3676 if (e->x.p->type == relation) {
3677 for (i = 1; i < e->x.p->size; ++i)
3678 evalue_add_constant(&e->x.p->arr[i], cst);
3679 return;
3680 }
3681 do {
3682 e = &e->x.p->arr[type_offset(e->x.p)];
3683 } while (value_zero_p(e->d));
3684 }
3685 value_addmul(e->x.n, cst, e->d);
3686 }
3687
3688 static void evalue_frac2polynomial_r(evalue *e, int *signs, int sign, int in_frac)
3689 {
3690 int i, offset;
3691 Value d;
3692 enode *p;
3693 int sign_odd = sign;
3694
3695 if (value_notzero_p(e->d)) {
3696 if (in_frac && sign * value_sign(e->x.n) < 0) {
3697 value_set_si(e->x.n, 0);
3698 value_set_si(e->d, 1);
3699 }
3700 return;
3701 }
3702
3703 if (e->x.p->type == relation) {
3704 for (i = e->x.p->size-1; i >= 1; --i)
3705 evalue_frac2polynomial_r(&e->x.p->arr[i], signs, sign, in_frac);
3706 return;
3707 }
3708
3709 if (e->x.p->type == polynomial)
3710 sign_odd *= signs[e->x.p->pos-1];
3711 offset = type_offset(e->x.p);
3712 evalue_frac2polynomial_r(&e->x.p->arr[offset], signs, sign, in_frac);
3713 in_frac |= e->x.p->type == fractional;
3714 for (i = e->x.p->size-1; i > offset; --i)
3715 evalue_frac2polynomial_r(&e->x.p->arr[i], signs,
3716 (i - offset) % 2 ? sign_odd : sign, in_frac);
3717
3718 if (e->x.p->type != fractional)
3719 return;
3720
3721 /* replace { a/m } by (m-1)/m if sign != 0
3722 * and by (m-1)/(2m) if sign == 0
3723 */
3724 value_init(d);
3725 value_set_si(d, 1);
3726 evalue_denom(&e->x.p->arr[0], &d);
3727 free_evalue_refs(&e->x.p->arr[0]);
3728 value_init(e->x.p->arr[0].d);
3729 value_init(e->x.p->arr[0].x.n);
3730 if (sign == 0)
3731 value_addto(e->x.p->arr[0].d, d, d);
3732 else
3733 value_assign(e->x.p->arr[0].d, d);
3734 value_decrement(e->x.p->arr[0].x.n, d);
3735 value_clear(d);
3736
3737 p = e->x.p;
3738 reorder_terms_about(p, &p->arr[0]);
3739 value_clear(e->d);
3740 *e = p->arr[1];
3741 free(p);
3742 }
3743
3744 /* Approximate the evalue in fractional representation by a polynomial.
3745 * If sign > 0, the result is an upper bound;
3746 * if sign < 0, the result is a lower bound;
3747 * if sign = 0, the result is an intermediate approximation.
3748 */
3749 void evalue_frac2polynomial(evalue *e, int sign, unsigned MaxRays)
3750 {
3751 int i, dim;
3752 int *signs;
3753
3754 if (value_notzero_p(e->d))
3755 return;
3756 assert(e->x.p->type == partition);
3757 /* make sure all variables in the domains have a fixed sign */
3758 if (sign)
3759 evalue_split_domains_into_orthants(e, MaxRays);
3760
3761 assert(e->x.p->size >= 2);
3762 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3763
3764 signs = alloca(sizeof(int) * dim);
3765
3766 if (!sign)
3767 for (i = 0; i < dim; ++i)
3768 signs[i] = 0;
3769 for (i = 0; i < e->x.p->size/2; ++i) {
3770 if (sign)
3771 domain_signs(EVALUE_DOMAIN(e->x.p->arr[2*i]), signs);
3772 evalue_frac2polynomial_r(&e->x.p->arr[2*i+1], signs, sign, 0);
3773 }
3774 }
3775
3776 /* Split the domains of e (which is assumed to be a partition)
3777 * such that each resulting domain lies entirely in one orthant.
3778 */
3779 void evalue_split_domains_into_orthants(evalue *e, unsigned MaxRays)
3780 {
3781 int i, dim;
3782 assert(value_zero_p(e->d));
3783 assert(e->x.p->type == partition);
3784 assert(e->x.p->size >= 2);
3785 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3786
3787 for (i = 0; i < dim; ++i) {
3788 evalue split;
3789 Matrix *C, *C2;
3790 C = Matrix_Alloc(1, 1 + dim + 1);
3791 value_set_si(C->p[0][0], 1);
3792 value_init(split.d);
3793 value_set_si(split.d, 0);
3794 split.x.p = new_enode(partition, 4, dim);
3795 value_set_si(C->p[0][1+i], 1);
3796 C2 = Matrix_Copy(C);
3797 EVALUE_SET_DOMAIN(split.x.p->arr[0], Constraints2Polyhedron(C2, MaxRays));
3798 Matrix_Free(C2);
3799 evalue_set_si(&split.x.p->arr[1], 1, 1);
3800 value_set_si(C->p[0][1+i], -1);
3801 value_set_si(C->p[0][1+dim], -1);
3802 EVALUE_SET_DOMAIN(split.x.p->arr[2], Constraints2Polyhedron(C, MaxRays));
3803 evalue_set_si(&split.x.p->arr[3], 1, 1);
3804 emul(&split, e);
3805 free_evalue_refs(&split);
3806 Matrix_Free(C);
3807 }
3808 }
3809
3810 static evalue *find_fractional_with_max_periods(evalue *e, Polyhedron *D,
3811 int max_periods,
3812 Matrix **TT,
3813 Value *min, Value *max)
3814 {
3815 Matrix *T;
3816 evalue *f = NULL;
3817 Value d;
3818 int i;
3819
3820 if (value_notzero_p(e->d))
3821 return NULL;
3822
3823 if (e->x.p->type == fractional) {
3824 Polyhedron *I;
3825 int bounded;
3826
3827 value_init(d);
3828 I = polynomial_projection(e->x.p, D, &d, &T);
3829 bounded = line_minmax(I, min, max); /* frees I */
3830 if (bounded) {
3831 Value mp;
3832 value_init(mp);
3833 value_set_si(mp, max_periods);
3834 mpz_fdiv_q(*min, *min, d);
3835 mpz_fdiv_q(*max, *max, d);
3836 value_assign(T->p[1][D->Dimension], d);
3837 value_subtract(d, *max, *min);
3838 if (value_ge(d, mp))
3839 Matrix_Free(T);
3840 else
3841 f = evalue_dup(&e->x.p->arr[0]);
3842 value_clear(mp);
3843 } else
3844 Matrix_Free(T);
3845 value_clear(d);
3846 if (f) {
3847 *TT = T;
3848 return f;
3849 }
3850 }
3851
3852 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
3853 if ((f = find_fractional_with_max_periods(&e->x.p->arr[i], D, max_periods,
3854 TT, min, max)))
3855 return f;
3856
3857 return NULL;
3858 }
3859
3860 static void replace_fract_by_affine(evalue *e, evalue *f, Value val)
3861 {
3862 int i, offset;
3863
3864 if (value_notzero_p(e->d))
3865 return;
3866
3867 offset = type_offset(e->x.p);
3868 for (i = e->x.p->size-1; i >= offset; --i)
3869 replace_fract_by_affine(&e->x.p->arr[i], f, val);
3870
3871 if (e->x.p->type != fractional)
3872 return;
3873
3874 if (!eequal(&e->x.p->arr[0], f))
3875 return;
3876
3877 replace_by_affine(e, val);
3878 }
3879
3880 /* Look for fractional parts that can be removed by splitting the corresponding
3881 * domain into at most max_periods parts.
3882 * We use a very simply strategy that looks for the first fractional part
3883 * that satisfies the condition, performs the split and then continues
3884 * looking for other fractional parts in the split domains until no
3885 * such fractional part can be found anymore.
3886 */
3887 void evalue_split_periods(evalue *e, int max_periods, unsigned int MaxRays)
3888 {
3889 int i, j, n;
3890 Value min;
3891 Value max;
3892 Value d;
3893
3894 if (EVALUE_IS_ZERO(*e))
3895 return;
3896 if (value_notzero_p(e->d) || e->x.p->type != partition) {
3897 fprintf(stderr,
3898 "WARNING: evalue_split_periods called on incorrect evalue type\n");
3899 return;
3900 }
3901
3902 value_init(min);
3903 value_init(max);
3904 value_init(d);
3905
3906 for (i = 0; i < e->x.p->size/2; ++i) {
3907 enode *p;
3908 evalue *f;
3909 Matrix *T;
3910 Matrix *M;
3911 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
3912 Polyhedron *E;
3913 f = find_fractional_with_max_periods(&e->x.p->arr[2*i+1], D, max_periods,
3914 &T, &min, &max);
3915 if (!f)
3916 continue;
3917
3918 M = Matrix_Alloc(2, 2+D->Dimension);
3919
3920 value_subtract(d, max, min);
3921 n = VALUE_TO_INT(d)+1;
3922
3923 value_set_si(M->p[0][0], 1);
3924 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
3925 value_multiply(d, max, T->p[1][D->Dimension]);
3926 value_subtract(M->p[0][1+D->Dimension], M->p[0][1+D->Dimension], d);
3927 value_set_si(d, -1);
3928 value_set_si(M->p[1][0], 1);
3929 Vector_Scale(T->p[0], M->p[1]+1, d, D->Dimension+1);
3930 value_addmul(M->p[1][1+D->Dimension], max, T->p[1][D->Dimension]);
3931 value_addto(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension],
3932 T->p[1][D->Dimension]);
3933 value_decrement(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension]);
3934
3935 p = new_enode(partition, e->x.p->size + (n-1)*2, e->x.p->pos);
3936 for (j = 0; j < 2*i; ++j) {
3937 value_clear(p->arr[j].d);
3938 p->arr[j] = e->x.p->arr[j];
3939 }
3940 for (j = 2*i+2; j < e->x.p->size; ++j) {
3941 value_clear(p->arr[j+2*(n-1)].d);
3942 p->arr[j+2*(n-1)] = e->x.p->arr[j];
3943 }
3944 for (j = n-1; j >= 0; --j) {
3945 if (j == 0) {
3946 value_clear(p->arr[2*i+1].d);
3947 p->arr[2*i+1] = e->x.p->arr[2*i+1];
3948 } else
3949 evalue_copy(&p->arr[2*(i+j)+1], &e->x.p->arr[2*i+1]);
3950 if (j != n-1) {
3951 value_subtract(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension],
3952 T->p[1][D->Dimension]);
3953 value_addto(M->p[0][1+D->Dimension], M->p[0][1+D->Dimension],
3954 T->p[1][D->Dimension]);
3955 }
3956 replace_fract_by_affine(&p->arr[2*(i+j)+1], f, max);
3957 E = DomainAddConstraints(D, M, MaxRays);
3958 EVALUE_SET_DOMAIN(p->arr[2*(i+j)], E);
3959 if (evalue_range_reduction_in_domain(&p->arr[2*(i+j)+1], E))
3960 reduce_evalue(&p->arr[2*(i+j)+1]);
3961 value_decrement(max, max);
3962 }
3963 value_clear(e->x.p->arr[2*i].d);
3964 Domain_Free(D);
3965 Matrix_Free(M);
3966 Matrix_Free(T);
3967 free_evalue_refs(f);
3968 free(f);
3969 free(e->x.p);
3970 e->x.p = p;
3971 --i;
3972 }
3973
3974 value_clear(d);
3975 value_clear(min);
3976 value_clear(max);
3977 }
3978
3979 void evalue_extract_affine(const evalue *e, Value *coeff, Value *cst, Value *d)
3980 {
3981 value_set_si(*d, 1);
3982 evalue_denom(e, d);
3983 for ( ; value_zero_p(e->d); e = &e->x.p->arr[0]) {
3984 evalue *c;
3985 assert(e->x.p->type == polynomial);
3986 assert(e->x.p->size == 2);
3987 c = &e->x.p->arr[1];
3988 value_multiply(coeff[e->x.p->pos-1], *d, c->x.n);
3989 value_division(coeff[e->x.p->pos-1], coeff[e->x.p->pos-1], c->d);
3990 }
3991 value_multiply(*cst, *d, e->x.n);
3992 value_division(*cst, *cst, e->d);
3993 }
3994
3995 /* returns an evalue that corresponds to
3996 *
3997 * c/den x_param
3998 */
3999 static evalue *term(int param, Value c, Value den)
4000 {
4001 evalue *EP = ALLOC(evalue);
4002 value_init(EP->d);
4003 value_set_si(EP->d,0);
4004 EP->x.p = new_enode(polynomial, 2, param + 1);
4005 evalue_set_si(&EP->x.p->arr[0], 0, 1);
4006 value_init(EP->x.p->arr[1].x.n);
4007 value_assign(EP->x.p->arr[1].d, den);
4008 value_assign(EP->x.p->arr[1].x.n, c);
4009 return EP;
4010 }
4011
4012 evalue *affine2evalue(Value *coeff, Value denom, int nvar)
4013 {
4014 int i;
4015 evalue *E = ALLOC(evalue);
4016 value_init(E->d);
4017 evalue_set(E, coeff[nvar], denom);
4018 for (i = 0; i < nvar; ++i) {
4019 if (value_zero_p(coeff[i]))
4020 continue;
4021 evalue *t = term(i, coeff[i], denom);
4022 eadd(t, E);
4023 free_evalue_refs(t);
4024 free(t);
4025 }
4026 return E;
4027 }
4028
4029 void evalue_substitute(evalue *e, evalue **subs)
4030 {
4031 int i, offset;
4032 evalue *v;
4033 enode *p;
4034
4035 if (value_notzero_p(e->d))
4036 return;
4037
4038 p = e->x.p;
4039 assert(p->type != partition);
4040
4041 for (i = 0; i < p->size; ++i)
4042 evalue_substitute(&p->arr[i], subs);
4043
4044 if (p->type == polynomial)
4045 v = subs[p->pos-1];
4046 else {
4047 v = ALLOC(evalue);
4048 value_init(v->d);
4049 value_set_si(v->d, 0);
4050 v->x.p = new_enode(p->type, 3, -1);
4051 value_clear(v->x.p->arr[0].d);
4052 v->x.p->arr[0] = p->arr[0];
4053 evalue_set_si(&v->x.p->arr[1], 0, 1);
4054 evalue_set_si(&v->x.p->arr[2], 1, 1);
4055 }
4056
4057 offset = type_offset(p);
4058
4059 for (i = p->size-1; i >= offset+1; i--) {
4060 emul(v, &p->arr[i]);
4061 eadd(&p->arr[i], &p->arr[i-1]);
4062 free_evalue_refs(&(p->arr[i]));
4063 }
4064
4065 if (p->type != polynomial) {
4066 free_evalue_refs(v);
4067 free(v);
4068 }
4069
4070 value_clear(e->d);
4071 *e = p->arr[offset];
4072 free(p);
4073 }
4074
4075 /* evalue e is given in terms of "new" parameter; CP maps the new
4076 * parameters back to the old parameters.
4077 * Transforms e such that it refers back to the old parameters.
4078 */
4079 void evalue_backsubstitute(evalue *e, Matrix *CP, unsigned MaxRays)
4080 {
4081 Matrix *eq;
4082 Matrix *inv;
4083 evalue **subs;
4084 enode *p;
4085 int i;
4086 unsigned nparam = CP->NbColumns-1;
4087 Polyhedron *CEq;
4088
4089 if (EVALUE_IS_ZERO(*e))
4090 return;
4091
4092 assert(value_zero_p(e->d));
4093 p = e->x.p;
4094 assert(p->type == partition);
4095
4096 inv = left_inverse(CP, &eq);
4097 subs = ALLOCN(evalue *, nparam);
4098 for (i = 0; i < nparam; ++i)
4099 subs[i] = affine2evalue(inv->p[i], inv->p[nparam][inv->NbColumns-1],
4100 inv->NbColumns-1);
4101
4102 CEq = Constraints2Polyhedron(eq, MaxRays);
4103 addeliminatedparams_partition(p, inv, CEq, inv->NbColumns-1, MaxRays);
4104 Polyhedron_Free(CEq);
4105
4106 for (i = 0; i < p->size/2; ++i)
4107 evalue_substitute(&p->arr[2*i+1], subs);
4108
4109 for (i = 0; i < nparam; ++i) {
4110 free_evalue_refs(subs[i]);
4111 free(subs[i]);
4112 }
4113 free(subs);
4114 Matrix_Free(eq);
4115 Matrix_Free(inv);
4116 }
4117
4118 evalue *evalue_polynomial(Vector *c, evalue* X)
4119 {
4120 unsigned dim = c->Size-2;
4121 evalue EC;
4122 evalue *EP = ALLOC(evalue);
4123 int i;
4124
4125 value_init(EC.d);
4126 evalue_set(&EC, c->p[dim], c->p[dim+1]);
4127
4128 value_init(EP->d);
4129 evalue_set(EP, c->p[dim], c->p[dim+1]);
4130
4131 for (i = dim-1; i >= 0; --i) {
4132 emul(X, EP);
4133 value_assign(EC.x.n, c->p[i]);
4134 eadd(&EC, EP);
4135 }
4136 free_evalue_refs(&EC);
4137 return EP;
4138 }
lattice point enumeration library