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