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