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