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