?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Make laurent based summation the default
[barvinok.git] / evalue.c
blob4e823515eb2d731a28beafb5d11e89f633425bfa
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 evalue* evalue_zero()
47 {
48 evalue *EP = ALLOC(evalue);
49 value_init(EP->d);
50 evalue_set_si(EP, 0, 1);
51 return EP;
52 }
53
54 evalue *evalue_nan()
55 {
56 evalue *EP = ALLOC(evalue);
57 value_init(EP->d);
58 value_set_si(EP->d, -2);
59 EP->x.p = NULL;
60 return EP;
61 }
62
63 /* returns an evalue that corresponds to
64 *
65 * x_var
66 */
67 evalue *evalue_var(int var)
68 {
69 evalue *EP = ALLOC(evalue);
70 value_init(EP->d);
71 value_set_si(EP->d,0);
72 EP->x.p = new_enode(polynomial, 2, var + 1);
73 evalue_set_si(&EP->x.p->arr[0], 0, 1);
74 evalue_set_si(&EP->x.p->arr[1], 1, 1);
75 return EP;
76 }
77
78 void aep_evalue(evalue *e, int *ref) {
79
80 enode *p;
81 int i;
82
83 if (value_notzero_p(e->d))
84 return; /* a rational number, its already reduced */
85 if(!(p = e->x.p))
86 return; /* hum... an overflow probably occured */
87
88 /* First check the components of p */
89 for (i=0;i<p->size;i++)
90 aep_evalue(&p->arr[i],ref);
91
92 /* Then p itself */
93 switch (p->type) {
94 case polynomial:
95 case periodic:
96 case evector:
97 p->pos = ref[p->pos-1]+1;
98 }
99 return;
100 } /* aep_evalue */
101
102 /** Comments */
103 void addeliminatedparams_evalue(evalue *e,Matrix *CT) {
104
105 enode *p;
106 int i, j;
107 int *ref;
108
109 if (value_notzero_p(e->d))
110 return; /* a rational number, its already reduced */
111 if(!(p = e->x.p))
112 return; /* hum... an overflow probably occured */
113
114 /* Compute ref */
115 ref = (int *)malloc(sizeof(int)*(CT->NbRows-1));
116 for(i=0;i<CT->NbRows-1;i++)
117 for(j=0;j<CT->NbColumns;j++)
118 if(value_notzero_p(CT->p[i][j])) {
119 ref[i] = j;
120 break;
121 }
122
123 /* Transform the references in e, using ref */
124 aep_evalue(e,ref);
125 free( ref );
126 return;
127 } /* addeliminatedparams_evalue */
128
129 static void addeliminatedparams_partition(enode *p, Matrix *CT, Polyhedron *CEq,
130 unsigned nparam, unsigned MaxRays)
131 {
132 int i;
133 assert(p->type == partition);
134 p->pos = nparam;
135
136 for (i = 0; i < p->size/2; i++) {
137 Polyhedron *D = EVALUE_DOMAIN(p->arr[2*i]);
138 Polyhedron *T = DomainPreimage(D, CT, MaxRays);
139 Domain_Free(D);
140 if (CEq) {
141 D = T;
142 T = DomainIntersection(D, CEq, MaxRays);
143 Domain_Free(D);
144 }
145 EVALUE_SET_DOMAIN(p->arr[2*i], T);
146 }
147 }
148
149 void addeliminatedparams_enum(evalue *e, Matrix *CT, Polyhedron *CEq,
150 unsigned MaxRays, unsigned nparam)
151 {
152 enode *p;
153 int i;
154
155 if (CT->NbRows == CT->NbColumns)
156 return;
157
158 if (EVALUE_IS_ZERO(*e))
159 return;
160
161 if (value_notzero_p(e->d)) {
162 evalue res;
163 value_init(res.d);
164 value_set_si(res.d, 0);
165 res.x.p = new_enode(partition, 2, nparam);
166 EVALUE_SET_DOMAIN(res.x.p->arr[0],
167 DomainConstraintSimplify(Polyhedron_Copy(CEq), MaxRays));
168 value_clear(res.x.p->arr[1].d);
169 res.x.p->arr[1] = *e;
170 *e = res;
171 return;
172 }
173
174 p = e->x.p;
175 assert(p);
176
177 addeliminatedparams_partition(p, CT, CEq, nparam, MaxRays);
178 for (i = 0; i < p->size/2; i++)
179 addeliminatedparams_evalue(&p->arr[2*i+1], CT);
180 }
181
182 static int mod_rational_smaller(evalue *e1, evalue *e2)
183 {
184 int r;
185 Value m;
186 Value m2;
187 value_init(m);
188 value_init(m2);
189
190 assert(value_notzero_p(e1->d));
191 assert(value_notzero_p(e2->d));
192 value_multiply(m, e1->x.n, e2->d);
193 value_multiply(m2, e2->x.n, e1->d);
194 if (value_lt(m, m2))
195 r = 1;
196 else if (value_gt(m, m2))
197 r = 0;
198 else
199 r = -1;
200 value_clear(m);
201 value_clear(m2);
202
203 return r;
204 }
205
206 static int mod_term_smaller_r(evalue *e1, evalue *e2)
207 {
208 if (value_notzero_p(e1->d)) {
209 int r;
210 if (value_zero_p(e2->d))
211 return 1;
212 r = mod_rational_smaller(e1, e2);
213 return r == -1 ? 0 : r;
214 }
215 if (value_notzero_p(e2->d))
216 return 0;
217 if (e1->x.p->pos < e2->x.p->pos)
218 return 1;
219 else if (e1->x.p->pos > e2->x.p->pos)
220 return 0;
221 else {
222 int r = mod_rational_smaller(&e1->x.p->arr[1], &e2->x.p->arr[1]);
223 return r == -1
224 ? mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0])
225 : r;
226 }
227 }
228
229 static int mod_term_smaller(const evalue *e1, const evalue *e2)
230 {
231 assert(value_zero_p(e1->d));
232 assert(value_zero_p(e2->d));
233 assert(e1->x.p->type == fractional || e1->x.p->type == flooring);
234 assert(e2->x.p->type == fractional || e2->x.p->type == flooring);
235 return mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0]);
236 }
237
238 static void check_order(const evalue *e)
239 {
240 int i;
241 evalue *a;
242
243 if (value_notzero_p(e->d))
244 return;
245
246 switch (e->x.p->type) {
247 case partition:
248 for (i = 0; i < e->x.p->size/2; ++i)
249 check_order(&e->x.p->arr[2*i+1]);
250 break;
251 case relation:
252 for (i = 1; i < e->x.p->size; ++i) {
253 a = &e->x.p->arr[i];
254 if (value_notzero_p(a->d))
255 continue;
256 switch (a->x.p->type) {
257 case relation:
258 assert(mod_term_smaller(&e->x.p->arr[0], &a->x.p->arr[0]));
259 break;
260 case partition:
261 assert(0);
262 }
263 check_order(a);
264 }
265 break;
266 case polynomial:
267 for (i = 0; i < e->x.p->size; ++i) {
268 a = &e->x.p->arr[i];
269 if (value_notzero_p(a->d))
270 continue;
271 switch (a->x.p->type) {
272 case polynomial:
273 assert(e->x.p->pos < a->x.p->pos);
274 break;
275 case relation:
276 case partition:
277 assert(0);
278 }
279 check_order(a);
280 }
281 break;
282 case fractional:
283 case flooring:
284 for (i = 1; i < e->x.p->size; ++i) {
285 a = &e->x.p->arr[i];
286 if (value_notzero_p(a->d))
287 continue;
288 switch (a->x.p->type) {
289 case polynomial:
290 case relation:
291 case partition:
292 assert(0);
293 }
294 }
295 break;
296 }
297 }
298
299 /* Negative pos means inequality */
300 /* s is negative of substitution if m is not zero */
301 struct fixed_param {
302 int pos;
303 evalue s;
304 Value d;
305 Value m;
306 };
307
308 struct subst {
309 struct fixed_param *fixed;
310 int n;
311 int max;
312 };
313
314 static int relations_depth(evalue *e)
315 {
316 int d;
317
318 for (d = 0;
319 value_zero_p(e->d) && e->x.p->type == relation;
320 e = &e->x.p->arr[1], ++d);
321 return d;
322 }
323
324 static void poly_denom_not_constant(evalue **pp, Value *d)
325 {
326 evalue *p = *pp;
327 value_set_si(*d, 1);
328
329 while (value_zero_p(p->d)) {
330 assert(p->x.p->type == polynomial);
331 assert(p->x.p->size == 2);
332 assert(value_notzero_p(p->x.p->arr[1].d));
333 value_lcm(*d, *d, p->x.p->arr[1].d);
334 p = &p->x.p->arr[0];
335 }
336 *pp = p;
337 }
338
339 static void poly_denom(evalue *p, Value *d)
340 {
341 poly_denom_not_constant(&p, d);
342 value_lcm(*d, *d, p->d);
343 }
344
345 static void realloc_substitution(struct subst *s, int d)
346 {
347 struct fixed_param *n;
348 int i;
349 NALLOC(n, s->max+d);
350 for (i = 0; i < s->n; ++i)
351 n[i] = s->fixed[i];
352 free(s->fixed);
353 s->fixed = n;
354 s->max += d;
355 }
356
357 static int add_modulo_substitution(struct subst *s, evalue *r)
358 {
359 evalue *p;
360 evalue *f;
361 evalue *m;
362
363 assert(value_zero_p(r->d) && r->x.p->type == relation);
364 m = &r->x.p->arr[0];
365
366 /* May have been reduced already */
367 if (value_notzero_p(m->d))
368 return 0;
369
370 assert(value_zero_p(m->d) && m->x.p->type == fractional);
371 assert(m->x.p->size == 3);
372
373 /* fractional was inverted during reduction
374 * invert it back and move constant in
375 */
376 if (!EVALUE_IS_ONE(m->x.p->arr[2])) {
377 assert(value_pos_p(m->x.p->arr[2].d));
378 assert(value_mone_p(m->x.p->arr[2].x.n));
379 value_set_si(m->x.p->arr[2].x.n, 1);
380 value_increment(m->x.p->arr[1].x.n, m->x.p->arr[1].x.n);
381 assert(value_eq(m->x.p->arr[1].x.n, m->x.p->arr[1].d));
382 value_set_si(m->x.p->arr[1].x.n, 1);
383 eadd(&m->x.p->arr[1], &m->x.p->arr[0]);
384 value_set_si(m->x.p->arr[1].x.n, 0);
385 value_set_si(m->x.p->arr[1].d, 1);
386 }
387
388 /* Oops. Nested identical relations. */
389 if (!EVALUE_IS_ZERO(m->x.p->arr[1]))
390 return 0;
391
392 if (s->n >= s->max) {
393 int d = relations_depth(r);
394 realloc_substitution(s, d);
395 }
396
397 p = &m->x.p->arr[0];
398 assert(value_zero_p(p->d) && p->x.p->type == polynomial);
399 assert(p->x.p->size == 2);
400 f = &p->x.p->arr[1];
401
402 assert(value_pos_p(f->x.n));
403
404 value_init(s->fixed[s->n].m);
405 value_assign(s->fixed[s->n].m, f->d);
406 s->fixed[s->n].pos = p->x.p->pos;
407 value_init(s->fixed[s->n].d);
408 value_assign(s->fixed[s->n].d, f->x.n);
409 value_init(s->fixed[s->n].s.d);
410 evalue_copy(&s->fixed[s->n].s, &p->x.p->arr[0]);
411 ++s->n;
412
413 return 1;
414 }
415
416 static int type_offset(enode *p)
417 {
418 return p->type == fractional ? 1 :
419 p->type == flooring ? 1 :
420 p->type == relation ? 1 : 0;
421 }
422
423 static void reorder_terms_about(enode *p, evalue *v)
424 {
425 int i;
426 int offset = type_offset(p);
427
428 for (i = p->size-1; i >= offset+1; i--) {
429 emul(v, &p->arr[i]);
430 eadd(&p->arr[i], &p->arr[i-1]);
431 free_evalue_refs(&(p->arr[i]));
432 }
433 p->size = offset+1;
434 free_evalue_refs(v);
435 }
436
437 static void reorder_terms(evalue *e)
438 {
439 enode *p;
440 evalue f;
441 int offset;
442
443 assert(value_zero_p(e->d));
444 p = e->x.p;
445 assert(p->type == fractional ||
446 p->type == flooring ||
447 p->type == polynomial); /* for now */
448
449 offset = type_offset(p);
450 value_init(f.d);
451 value_set_si(f.d, 0);
452 f.x.p = new_enode(p->type, offset+2, p->pos);
453 if (offset == 1) {
454 value_clear(f.x.p->arr[0].d);
455 f.x.p->arr[0] = p->arr[0];
456 }
457 evalue_set_si(&f.x.p->arr[offset], 0, 1);
458 evalue_set_si(&f.x.p->arr[offset+1], 1, 1);
459 reorder_terms_about(p, &f);
460 value_clear(e->d);
461 *e = p->arr[offset];
462 free(p);
463 }
464
465 static void evalue_reduce_size(evalue *e)
466 {
467 int i, offset;
468 enode *p;
469 assert(value_zero_p(e->d));
470
471 p = e->x.p;
472 offset = type_offset(p);
473
474 /* Try to reduce the degree */
475 for (i = p->size-1; i >= offset+1; i--) {
476 if (!EVALUE_IS_ZERO(p->arr[i]))
477 break;
478 free_evalue_refs(&p->arr[i]);
479 }
480 if (i+1 < p->size)
481 p->size = i+1;
482
483 /* Try to reduce its strength */
484 if (p->type == relation) {
485 if (p->size == 1) {
486 free_evalue_refs(&p->arr[0]);
487 evalue_set_si(e, 0, 1);
488 free(p);
489 }
490 } else if (p->size == offset+1) {
491 value_clear(e->d);
492 memcpy(e, &p->arr[offset], sizeof(evalue));
493 if (offset == 1)
494 free_evalue_refs(&p->arr[0]);
495 free(p);
496 }
497 }
498
499 void _reduce_evalue (evalue *e, struct subst *s, int fract) {
500
501 enode *p;
502 int i, j, k;
503 int add;
504
505 if (value_notzero_p(e->d)) {
506 if (fract)
507 mpz_fdiv_r(e->x.n, e->x.n, e->d);
508 return; /* a rational number, its already reduced */
509 }
510
511 if(!(p = e->x.p))
512 return; /* hum... an overflow probably occured */
513
514 /* First reduce the components of p */
515 add = p->type == relation;
516 for (i=0; i<p->size; i++) {
517 if (add && i == 1)
518 add = add_modulo_substitution(s, e);
519
520 if (i == 0 && p->type==fractional)
521 _reduce_evalue(&p->arr[i], s, 1);
522 else
523 _reduce_evalue(&p->arr[i], s, fract);
524
525 if (add && i == p->size-1) {
526 --s->n;
527 value_clear(s->fixed[s->n].m);
528 value_clear(s->fixed[s->n].d);
529 free_evalue_refs(&s->fixed[s->n].s);
530 } else if (add && i == 1)
531 s->fixed[s->n-1].pos *= -1;
532 }
533
534 if (p->type==periodic) {
535
536 /* Try to reduce the period */
537 for (i=1; i<=(p->size)/2; i++) {
538 if ((p->size % i)==0) {
539
540 /* Can we reduce the size to i ? */
541 for (j=0; j<i; j++)
542 for (k=j+i; k<e->x.p->size; k+=i)
543 if (!eequal(&p->arr[j], &p->arr[k])) goto you_lose;
544
545 /* OK, lets do it */
546 for (j=i; j<p->size; j++) free_evalue_refs(&p->arr[j]);
547 p->size = i;
548 break;
549
550 you_lose: /* OK, lets not do it */
551 continue;
552 }
553 }
554
555 /* Try to reduce its strength */
556 if (p->size == 1) {
557 value_clear(e->d);
558 memcpy(e,&p->arr[0],sizeof(evalue));
559 free(p);
560 }
561 }
562 else if (p->type==polynomial) {
563 for (k = 0; s && k < s->n; ++k) {
564 if (s->fixed[k].pos == p->pos) {
565 int divide = value_notone_p(s->fixed[k].d);
566 evalue d;
567
568 if (value_notzero_p(s->fixed[k].m)) {
569 if (!fract)
570 continue;
571 assert(p->size == 2);
572 if (divide && value_ne(s->fixed[k].d, p->arr[1].x.n))
573 continue;
574 if (!mpz_divisible_p(s->fixed[k].m, p->arr[1].d))
575 continue;
576 divide = 1;
577 }
578
579 if (divide) {
580 value_init(d.d);
581 value_assign(d.d, s->fixed[k].d);
582 value_init(d.x.n);
583 if (value_notzero_p(s->fixed[k].m))
584 value_oppose(d.x.n, s->fixed[k].m);
585 else
586 value_set_si(d.x.n, 1);
587 }
588
589 for (i=p->size-1;i>=1;i--) {
590 emul(&s->fixed[k].s, &p->arr[i]);
591 if (divide)
592 emul(&d, &p->arr[i]);
593 eadd(&p->arr[i], &p->arr[i-1]);
594 free_evalue_refs(&(p->arr[i]));
595 }
596 p->size = 1;
597 _reduce_evalue(&p->arr[0], s, fract);
598
599 if (divide)
600 free_evalue_refs(&d);
601
602 break;
603 }
604 }
605
606 evalue_reduce_size(e);
607 }
608 else if (p->type==fractional) {
609 int reorder = 0;
610 evalue v;
611
612 if (value_notzero_p(p->arr[0].d)) {
613 value_init(v.d);
614 value_assign(v.d, p->arr[0].d);
615 value_init(v.x.n);
616 mpz_fdiv_r(v.x.n, p->arr[0].x.n, p->arr[0].d);
617
618 reorder = 1;
619 } else {
620 evalue *f, *base;
621 evalue *pp = &p->arr[0];
622 assert(value_zero_p(pp->d) && pp->x.p->type == polynomial);
623 assert(pp->x.p->size == 2);
624
625 /* search for exact duplicate among the modulo inequalities */
626 do {
627 f = &pp->x.p->arr[1];
628 for (k = 0; s && k < s->n; ++k) {
629 if (-s->fixed[k].pos == pp->x.p->pos &&
630 value_eq(s->fixed[k].d, f->x.n) &&
631 value_eq(s->fixed[k].m, f->d) &&
632 eequal(&s->fixed[k].s, &pp->x.p->arr[0]))
633 break;
634 }
635 if (k < s->n) {
636 Value g;
637 value_init(g);
638
639 /* replace { E/m } by { (E-1)/m } + 1/m */
640 poly_denom(pp, &g);
641 if (reorder) {
642 evalue extra;
643 value_init(extra.d);
644 evalue_set_si(&extra, 1, 1);
645 value_assign(extra.d, g);
646 eadd(&extra, &v.x.p->arr[1]);
647 free_evalue_refs(&extra);
648
649 /* We've been going in circles; stop now */
650 if (value_ge(v.x.p->arr[1].x.n, v.x.p->arr[1].d)) {
651 free_evalue_refs(&v);
652 value_init(v.d);
653 evalue_set_si(&v, 0, 1);
654 break;
655 }
656 } else {
657 value_init(v.d);
658 value_set_si(v.d, 0);
659 v.x.p = new_enode(fractional, 3, -1);
660 evalue_set_si(&v.x.p->arr[1], 1, 1);
661 value_assign(v.x.p->arr[1].d, g);
662 evalue_set_si(&v.x.p->arr[2], 1, 1);
663 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
664 }
665
666 for (f = &v.x.p->arr[0]; value_zero_p(f->d);
667 f = &f->x.p->arr[0])
668 ;
669 value_division(f->d, g, f->d);
670 value_multiply(f->x.n, f->x.n, f->d);
671 value_assign(f->d, g);
672 value_decrement(f->x.n, f->x.n);
673 mpz_fdiv_r(f->x.n, f->x.n, f->d);
674
675 value_gcd(g, f->d, f->x.n);
676 value_division(f->d, f->d, g);
677 value_division(f->x.n, f->x.n, g);
678
679 value_clear(g);
680 pp = &v.x.p->arr[0];
681
682 reorder = 1;
683 }
684 } while (k < s->n);
685
686 /* reduction may have made this fractional arg smaller */
687 i = reorder ? p->size : 1;
688 for ( ; i < p->size; ++i)
689 if (value_zero_p(p->arr[i].d) &&
690 p->arr[i].x.p->type == fractional &&
691 !mod_term_smaller(e, &p->arr[i]))
692 break;
693 if (i < p->size) {
694 value_init(v.d);
695 value_set_si(v.d, 0);
696 v.x.p = new_enode(fractional, 3, -1);
697 evalue_set_si(&v.x.p->arr[1], 0, 1);
698 evalue_set_si(&v.x.p->arr[2], 1, 1);
699 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
700
701 reorder = 1;
702 }
703
704 if (!reorder) {
705 Value m;
706 Value r;
707 evalue *pp = &p->arr[0];
708 value_init(m);
709 value_init(r);
710 poly_denom_not_constant(&pp, &m);
711 mpz_fdiv_r(r, m, pp->d);
712 if (value_notzero_p(r)) {
713 value_init(v.d);
714 value_set_si(v.d, 0);
715 v.x.p = new_enode(fractional, 3, -1);
716
717 value_multiply(r, m, pp->x.n);
718 value_multiply(v.x.p->arr[1].d, m, pp->d);
719 value_init(v.x.p->arr[1].x.n);
720 mpz_fdiv_r(v.x.p->arr[1].x.n, r, pp->d);
721 mpz_fdiv_q(r, r, pp->d);
722
723 evalue_set_si(&v.x.p->arr[2], 1, 1);
724 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
725 pp = &v.x.p->arr[0];
726 while (value_zero_p(pp->d))
727 pp = &pp->x.p->arr[0];
728
729 value_assign(pp->d, m);
730 value_assign(pp->x.n, r);
731
732 value_gcd(r, pp->d, pp->x.n);
733 value_division(pp->d, pp->d, r);
734 value_division(pp->x.n, pp->x.n, r);
735
736 reorder = 1;
737 }
738 value_clear(m);
739 value_clear(r);
740 }
741
742 if (!reorder) {
743 int invert = 0;
744 Value twice;
745 value_init(twice);
746
747 for (pp = &p->arr[0]; value_zero_p(pp->d);
748 pp = &pp->x.p->arr[0]) {
749 f = &pp->x.p->arr[1];
750 assert(value_pos_p(f->d));
751 mpz_mul_ui(twice, f->x.n, 2);
752 if (value_lt(twice, f->d))
753 break;
754 if (value_eq(twice, f->d))
755 continue;
756 invert = 1;
757 break;
758 }
759
760 if (invert) {
761 value_init(v.d);
762 value_set_si(v.d, 0);
763 v.x.p = new_enode(fractional, 3, -1);
764 evalue_set_si(&v.x.p->arr[1], 0, 1);
765 poly_denom(&p->arr[0], &twice);
766 value_assign(v.x.p->arr[1].d, twice);
767 value_decrement(v.x.p->arr[1].x.n, twice);
768 evalue_set_si(&v.x.p->arr[2], -1, 1);
769 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
770
771 for (pp = &v.x.p->arr[0]; value_zero_p(pp->d);
772 pp = &pp->x.p->arr[0]) {
773 f = &pp->x.p->arr[1];
774 value_oppose(f->x.n, f->x.n);
775 mpz_fdiv_r(f->x.n, f->x.n, f->d);
776 }
777 value_division(pp->d, twice, pp->d);
778 value_multiply(pp->x.n, pp->x.n, pp->d);
779 value_assign(pp->d, twice);
780 value_oppose(pp->x.n, pp->x.n);
781 value_decrement(pp->x.n, pp->x.n);
782 mpz_fdiv_r(pp->x.n, pp->x.n, pp->d);
783
784 /* Maybe we should do this during reduction of
785 * the constant.
786 */
787 value_gcd(twice, pp->d, pp->x.n);
788 value_division(pp->d, pp->d, twice);
789 value_division(pp->x.n, pp->x.n, twice);
790
791 reorder = 1;
792 }
793
794 value_clear(twice);
795 }
796 }
797
798 if (reorder) {
799 reorder_terms_about(p, &v);
800 _reduce_evalue(&p->arr[1], s, fract);
801 }
802
803 evalue_reduce_size(e);
804 }
805 else if (p->type == flooring) {
806 /* Replace floor(constant) by its value */
807 if (value_notzero_p(p->arr[0].d)) {
808 evalue v;
809 value_init(v.d);
810 value_set_si(v.d, 1);
811 value_init(v.x.n);
812 mpz_fdiv_q(v.x.n, p->arr[0].x.n, p->arr[0].d);
813 reorder_terms_about(p, &v);
814 _reduce_evalue(&p->arr[1], s, fract);
815 }
816 evalue_reduce_size(e);
817 }
818 else if (p->type == relation) {
819 if (p->size == 3 && eequal(&p->arr[1], &p->arr[2])) {
820 free_evalue_refs(&(p->arr[2]));
821 free_evalue_refs(&(p->arr[0]));
822 p->size = 2;
823 value_clear(e->d);
824 *e = p->arr[1];
825 free(p);
826 return;
827 }
828 evalue_reduce_size(e);
829 if (value_notzero_p(e->d) || p != e->x.p)
830 return;
831 else {
832 int reduced = 0;
833 evalue *m;
834 m = &p->arr[0];
835
836 /* Relation was reduced by means of an identical
837 * inequality => remove
838 */
839 if (value_zero_p(m->d) && !EVALUE_IS_ZERO(m->x.p->arr[1]))
840 reduced = 1;
841
842 if (reduced || value_notzero_p(p->arr[0].d)) {
843 if (!reduced && value_zero_p(p->arr[0].x.n)) {
844 value_clear(e->d);
845 memcpy(e,&p->arr[1],sizeof(evalue));
846 if (p->size == 3)
847 free_evalue_refs(&(p->arr[2]));
848 } else {
849 if (p->size == 3) {
850 value_clear(e->d);
851 memcpy(e,&p->arr[2],sizeof(evalue));
852 } else
853 evalue_set_si(e, 0, 1);
854 free_evalue_refs(&(p->arr[1]));
855 }
856 free_evalue_refs(&(p->arr[0]));
857 free(p);
858 }
859 }
860 }
861 return;
862 } /* reduce_evalue */
863
864 static void add_substitution(struct subst *s, Value *row, unsigned dim)
865 {
866 int k, l;
867 evalue *r;
868
869 for (k = 0; k < dim; ++k)
870 if (value_notzero_p(row[k+1]))
871 break;
872
873 Vector_Normalize_Positive(row+1, dim+1, k);
874 assert(s->n < s->max);
875 value_init(s->fixed[s->n].d);
876 value_init(s->fixed[s->n].m);
877 value_assign(s->fixed[s->n].d, row[k+1]);
878 s->fixed[s->n].pos = k+1;
879 value_set_si(s->fixed[s->n].m, 0);
880 r = &s->fixed[s->n].s;
881 value_init(r->d);
882 for (l = k+1; l < dim; ++l)
883 if (value_notzero_p(row[l+1])) {
884 value_set_si(r->d, 0);
885 r->x.p = new_enode(polynomial, 2, l + 1);
886 value_init(r->x.p->arr[1].x.n);
887 value_oppose(r->x.p->arr[1].x.n, row[l+1]);
888 value_set_si(r->x.p->arr[1].d, 1);
889 r = &r->x.p->arr[0];
890 }
891 value_init(r->x.n);
892 value_oppose(r->x.n, row[dim+1]);
893 value_set_si(r->d, 1);
894 ++s->n;
895 }
896
897 static void _reduce_evalue_in_domain(evalue *e, Polyhedron *D, struct subst *s)
898 {
899 unsigned dim;
900 Polyhedron *orig = D;
901
902 s->n = 0;
903 dim = D->Dimension;
904 if (D->next)
905 D = DomainConvex(D, 0);
906 /* We don't perform any substitutions if the domain is a union.
907 * We may therefore miss out on some possible simplifications,
908 * e.g., if a variable is always even in the whole union,
909 * while there is a relation in the evalue that evaluates
910 * to zero for even values of the variable.
911 */
912 if (!D->next && D->NbEq) {
913 int j, k;
914 if (s->max < dim) {
915 if (s->max != 0)
916 realloc_substitution(s, dim);
917 else {
918 int d = relations_depth(e);
919 s->max = dim+d;
920 NALLOC(s->fixed, s->max);
921 }
922 }
923 for (j = 0; j < D->NbEq; ++j)
924 add_substitution(s, D->Constraint[j], dim);
925 }
926 if (D != orig)
927 Domain_Free(D);
928 _reduce_evalue(e, s, 0);
929 if (s->n != 0) {
930 int j;
931 for (j = 0; j < s->n; ++j) {
932 value_clear(s->fixed[j].d);
933 value_clear(s->fixed[j].m);
934 free_evalue_refs(&s->fixed[j].s);
935 }
936 }
937 }
938
939 void reduce_evalue_in_domain(evalue *e, Polyhedron *D)
940 {
941 struct subst s = { NULL, 0, 0 };
942 POL_ENSURE_VERTICES(D);
943 if (emptyQ(D)) {
944 if (EVALUE_IS_ZERO(*e))
945 return;
946 free_evalue_refs(e);
947 value_init(e->d);
948 evalue_set_si(e, 0, 1);
949 return;
950 }
951 _reduce_evalue_in_domain(e, D, &s);
952 if (s.max != 0)
953 free(s.fixed);
954 }
955
956 void reduce_evalue (evalue *e) {
957 struct subst s = { NULL, 0, 0 };
958
959 if (value_notzero_p(e->d))
960 return; /* a rational number, its already reduced */
961
962 if (e->x.p->type == partition) {
963 int i;
964 unsigned dim = -1;
965 for (i = 0; i < e->x.p->size/2; ++i) {
966 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
967
968 /* This shouldn't really happen;
969 * Empty domains should not be added.
970 */
971 POL_ENSURE_VERTICES(D);
972 if (!emptyQ(D))
973 _reduce_evalue_in_domain(&e->x.p->arr[2*i+1], D, &s);
974
975 if (emptyQ(D) || EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
976 free_evalue_refs(&e->x.p->arr[2*i+1]);
977 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
978 value_clear(e->x.p->arr[2*i].d);
979 e->x.p->size -= 2;
980 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
981 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
982 --i;
983 }
984 }
985 if (e->x.p->size == 0) {
986 free(e->x.p);
987 evalue_set_si(e, 0, 1);
988 }
989 } else
990 _reduce_evalue(e, &s, 0);
991 if (s.max != 0)
992 free(s.fixed);
993 }
994
995 static void print_evalue_r(FILE *DST, const evalue *e, const char *const *pname)
996 {
997 if (EVALUE_IS_NAN(*e)) {
998 fprintf(DST, "NaN");
999 return;
1000 }
1001
1002 if(value_notzero_p(e->d)) {
1003 if(value_notone_p(e->d)) {
1004 value_print(DST,VALUE_FMT,e->x.n);
1005 fprintf(DST,"/");
1006 value_print(DST,VALUE_FMT,e->d);
1007 }
1008 else {
1009 value_print(DST,VALUE_FMT,e->x.n);
1010 }
1011 }
1012 else
1013 print_enode(DST,e->x.p,pname);
1014 return;
1015 } /* print_evalue */
1016
1017 void print_evalue(FILE *DST, const evalue *e, const char * const *pname)
1018 {
1019 print_evalue_r(DST, e, pname);
1020 if (value_notzero_p(e->d))
1021 fprintf(DST, "\n");
1022 }
1023
1024 void print_enode(FILE *DST, enode *p, const char *const *pname)
1025 {
1026 int i;
1027
1028 if (!p) {
1029 fprintf(DST, "NULL");
1030 return;
1031 }
1032 switch (p->type) {
1033 case evector:
1034 fprintf(DST, "{ ");
1035 for (i=0; i<p->size; i++) {
1036 print_evalue_r(DST, &p->arr[i], pname);
1037 if (i!=(p->size-1))
1038 fprintf(DST, ", ");
1039 }
1040 fprintf(DST, " }\n");
1041 break;
1042 case polynomial:
1043 fprintf(DST, "( ");
1044 for (i=p->size-1; i>=0; i--) {
1045 print_evalue_r(DST, &p->arr[i], pname);
1046 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
1047 else if (i>1)
1048 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
1049 }
1050 fprintf(DST, " )\n");
1051 break;
1052 case periodic:
1053 fprintf(DST, "[ ");
1054 for (i=0; i<p->size; i++) {
1055 print_evalue_r(DST, &p->arr[i], pname);
1056 if (i!=(p->size-1)) fprintf(DST, ", ");
1057 }
1058 fprintf(DST," ]_%s", pname[p->pos-1]);
1059 break;
1060 case flooring:
1061 case fractional:
1062 fprintf(DST, "( ");
1063 for (i=p->size-1; i>=1; i--) {
1064 print_evalue_r(DST, &p->arr[i], pname);
1065 if (i >= 2) {
1066 fprintf(DST, " * ");
1067 fprintf(DST, p->type == flooring ? "[" : "{");
1068 print_evalue_r(DST, &p->arr[0], pname);
1069 fprintf(DST, p->type == flooring ? "]" : "}");
1070 if (i>2)
1071 fprintf(DST, "^%d + ", i-1);
1072 else
1073 fprintf(DST, " + ");
1074 }
1075 }
1076 fprintf(DST, " )\n");
1077 break;
1078 case relation:
1079 fprintf(DST, "[ ");
1080 print_evalue_r(DST, &p->arr[0], pname);
1081 fprintf(DST, "= 0 ] * \n");
1082 print_evalue_r(DST, &p->arr[1], pname);
1083 if (p->size > 2) {
1084 fprintf(DST, " +\n [ ");
1085 print_evalue_r(DST, &p->arr[0], pname);
1086 fprintf(DST, "!= 0 ] * \n");
1087 print_evalue_r(DST, &p->arr[2], pname);
1088 }
1089 break;
1090 case partition: {
1091 char **new_names = NULL;
1092 const char *const *names = pname;
1093 int maxdim = EVALUE_DOMAIN(p->arr[0])->Dimension;
1094 if (!pname || p->pos < maxdim) {
1095 new_names = ALLOCN(char *, maxdim);
1096 for (i = 0; i < p->pos; ++i) {
1097 if (pname)
1098 new_names[i] = (char *)pname[i];
1099 else {
1100 new_names[i] = ALLOCN(char, 10);
1101 snprintf(new_names[i], 10, "%c", 'P'+i);
1102 }
1103 }
1104 for ( ; i < maxdim; ++i) {
1105 new_names[i] = ALLOCN(char, 10);
1106 snprintf(new_names[i], 10, "_p%d", i);
1107 }
1108 names = (const char**)new_names;
1109 }
1110
1111 for (i=0; i<p->size/2; i++) {
1112 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), names);
1113 print_evalue_r(DST, &p->arr[2*i+1], names);
1114 if (value_notzero_p(p->arr[2*i+1].d))
1115 fprintf(DST, "\n");
1116 }
1117
1118 if (!pname || p->pos < maxdim) {
1119 for (i = pname ? p->pos : 0; i < maxdim; ++i)
1120 free(new_names[i]);
1121 free(new_names);
1122 }
1123
1124 break;
1125 }
1126 default:
1127 assert(0);
1128 }
1129 return;
1130 } /* print_enode */
1131
1132 /* Returns
1133 * 0 if toplevels of e1 and e2 are at the same level
1134 * <0 if toplevel of e1 should be outside of toplevel of e2
1135 * >0 if toplevel of e2 should be outside of toplevel of e1
1136 */
1137 static int evalue_level_cmp(const evalue *e1, const evalue *e2)
1138 {
1139 if (value_notzero_p(e1->d) && value_notzero_p(e2->d))
1140 return 0;
1141 if (value_notzero_p(e1->d))
1142 return 1;
1143 if (value_notzero_p(e2->d))
1144 return -1;
1145 if (e1->x.p->type == partition && e2->x.p->type == partition)
1146 return 0;
1147 if (e1->x.p->type == partition)
1148 return -1;
1149 if (e2->x.p->type == partition)
1150 return 1;
1151 if (e1->x.p->type == relation && e2->x.p->type == relation) {
1152 if (eequal(&e1->x.p->arr[0], &e2->x.p->arr[0]))
1153 return 0;
1154 if (mod_term_smaller(&e1->x.p->arr[0], &e2->x.p->arr[0]))
1155 return -1;
1156 else
1157 return 1;
1158 }
1159 if (e1->x.p->type == relation)
1160 return -1;
1161 if (e2->x.p->type == relation)
1162 return 1;
1163 if (e1->x.p->type == polynomial && e2->x.p->type == polynomial)
1164 return e1->x.p->pos - e2->x.p->pos;
1165 if (e1->x.p->type == polynomial)
1166 return -1;
1167 if (e2->x.p->type == polynomial)
1168 return 1;
1169 if (e1->x.p->type == periodic && e2->x.p->type == periodic)
1170 return e1->x.p->pos - e2->x.p->pos;
1171 assert(e1->x.p->type != periodic);
1172 assert(e2->x.p->type != periodic);
1173 assert(e1->x.p->type == e2->x.p->type);
1174 if (eequal(&e1->x.p->arr[0], &e2->x.p->arr[0]))
1175 return 0;
1176 if (mod_term_smaller(e1, e2))
1177 return -1;
1178 else
1179 return 1;
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 value_init(v);
2209
2210 for (row = 0; row < P->NbConstraints; row++) {
2211 Inner_Product(P->Constraint[row]+1, list_args, P->Dimension, &v);
2212 value_addto(v, v, P->Constraint[row][P->Dimension+1]); /*constant part*/
2213 if (value_neg_p(v) ||
2214 value_zero_p(P->Constraint[row][0]) && value_notzero_p(v)) {
2215 in = 0;
2216 break;
2217 }
2218 }
2219
2220 value_clear(v);
2221 return in || (P->next && in_domain(P->next, list_args));
2222 } /* in_domain */
2223
2224 /****************************************************/
2225 /* function compute enode */
2226 /* compute the value of enode p with parameters */
2227 /* list "list_args */
2228 /* compute the polynomial or the periodic */
2229 /****************************************************/
2230
2231 static double compute_enode(enode *p, Value *list_args) {
2232
2233 int i;
2234 Value m, param;
2235 double res=0.0;
2236
2237 if (!p)
2238 return(0.);
2239
2240 value_init(m);
2241 value_init(param);
2242
2243 if (p->type == polynomial) {
2244 if (p->size > 1)
2245 value_assign(param,list_args[p->pos-1]);
2246
2247 /* Compute the polynomial using Horner's rule */
2248 for (i=p->size-1;i>0;i--) {
2249 res +=compute_evalue(&p->arr[i],list_args);
2250 res *=VALUE_TO_DOUBLE(param);
2251 }
2252 res +=compute_evalue(&p->arr[0],list_args);
2253 }
2254 else if (p->type == fractional) {
2255 double d = compute_evalue(&p->arr[0], list_args);
2256 d -= floor(d+1e-10);
2257
2258 /* Compute the polynomial using Horner's rule */
2259 for (i=p->size-1;i>1;i--) {
2260 res +=compute_evalue(&p->arr[i],list_args);
2261 res *=d;
2262 }
2263 res +=compute_evalue(&p->arr[1],list_args);
2264 }
2265 else if (p->type == flooring) {
2266 double d = compute_evalue(&p->arr[0], list_args);
2267 d = floor(d+1e-10);
2268
2269 /* Compute the polynomial using Horner's rule */
2270 for (i=p->size-1;i>1;i--) {
2271 res +=compute_evalue(&p->arr[i],list_args);
2272 res *=d;
2273 }
2274 res +=compute_evalue(&p->arr[1],list_args);
2275 }
2276 else if (p->type == periodic) {
2277 value_assign(m,list_args[p->pos-1]);
2278
2279 /* Choose the right element of the periodic */
2280 value_set_si(param,p->size);
2281 value_pmodulus(m,m,param);
2282 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
2283 }
2284 else if (p->type == relation) {
2285 if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
2286 res = compute_evalue(&p->arr[1], list_args);
2287 else if (p->size > 2)
2288 res = compute_evalue(&p->arr[2], list_args);
2289 }
2290 else if (p->type == partition) {
2291 int dim = EVALUE_DOMAIN(p->arr[0])->Dimension;
2292 Value *vals = list_args;
2293 if (p->pos < dim) {
2294 NALLOC(vals, dim);
2295 for (i = 0; i < dim; ++i) {
2296 value_init(vals[i]);
2297 if (i < p->pos)
2298 value_assign(vals[i], list_args[i]);
2299 }
2300 }
2301 for (i = 0; i < p->size/2; ++i)
2302 if (DomainContains(EVALUE_DOMAIN(p->arr[2*i]), vals, p->pos, 0, 1)) {
2303 res = compute_evalue(&p->arr[2*i+1], vals);
2304 break;
2305 }
2306 if (p->pos < dim) {
2307 for (i = 0; i < dim; ++i)
2308 value_clear(vals[i]);
2309 free(vals);
2310 }
2311 }
2312 else
2313 assert(0);
2314 value_clear(m);
2315 value_clear(param);
2316 return res;
2317 } /* compute_enode */
2318
2319 /*************************************************/
2320 /* return the value of Ehrhart Polynomial */
2321 /* It returns a double, because since it is */
2322 /* a recursive function, some intermediate value */
2323 /* might not be integral */
2324 /*************************************************/
2325
2326 double compute_evalue(const evalue *e, Value *list_args)
2327 {
2328 double res;
2329
2330 if (value_notzero_p(e->d)) {
2331 if (value_notone_p(e->d))
2332 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
2333 else
2334 res = VALUE_TO_DOUBLE(e->x.n);
2335 }
2336 else
2337 res = compute_enode(e->x.p,list_args);
2338 return res;
2339 } /* compute_evalue */
2340
2341
2342 /****************************************************/
2343 /* function compute_poly : */
2344 /* Check for the good validity domain */
2345 /* return the number of point in the Polyhedron */
2346 /* in allocated memory */
2347 /* Using the Ehrhart pseudo-polynomial */
2348 /****************************************************/
2349 Value *compute_poly(Enumeration *en,Value *list_args) {
2350
2351 Value *tmp;
2352 /* double d; int i; */
2353
2354 tmp = (Value *) malloc (sizeof(Value));
2355 assert(tmp != NULL);
2356 value_init(*tmp);
2357 value_set_si(*tmp,0);
2358
2359 if(!en)
2360 return(tmp); /* no ehrhart polynomial */
2361 if(en->ValidityDomain) {
2362 if(!en->ValidityDomain->Dimension) { /* no parameters */
2363 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2364 return(tmp);
2365 }
2366 }
2367 else
2368 return(tmp); /* no Validity Domain */
2369 while(en) {
2370 if(in_domain(en->ValidityDomain,list_args)) {
2371
2372 #ifdef EVAL_EHRHART_DEBUG
2373 Print_Domain(stdout,en->ValidityDomain);
2374 print_evalue(stdout,&en->EP);
2375 #endif
2376
2377 /* d = compute_evalue(&en->EP,list_args);
2378 i = d;
2379 printf("(double)%lf = %d\n", d, i ); */
2380 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2381 return(tmp);
2382 }
2383 else
2384 en=en->next;
2385 }
2386 value_set_si(*tmp,0);
2387 return(tmp); /* no compatible domain with the arguments */
2388 } /* compute_poly */
2389
2390 static evalue *eval_polynomial(const enode *p, int offset,
2391 evalue *base, Value *values)
2392 {
2393 int i;
2394 evalue *res, *c;
2395
2396 res = evalue_zero();
2397 for (i = p->size-1; i > offset; --i) {
2398 c = evalue_eval(&p->arr[i], values);
2399 eadd(c, res);
2400 evalue_free(c);
2401 emul(base, res);
2402 }
2403 c = evalue_eval(&p->arr[offset], values);
2404 eadd(c, res);
2405 evalue_free(c);
2406
2407 return res;
2408 }
2409
2410 evalue *evalue_eval(const evalue *e, Value *values)
2411 {
2412 evalue *res = NULL;
2413 evalue param;
2414 evalue *param2;
2415 int i;
2416
2417 if (value_notzero_p(e->d)) {
2418 res = ALLOC(evalue);
2419 value_init(res->d);
2420 evalue_copy(res, e);
2421 return res;
2422 }
2423 switch (e->x.p->type) {
2424 case polynomial:
2425 value_init(param.x.n);
2426 value_assign(param.x.n, values[e->x.p->pos-1]);
2427 value_init(param.d);
2428 value_set_si(param.d, 1);
2429
2430 res = eval_polynomial(e->x.p, 0, &param, values);
2431 free_evalue_refs(&param);
2432 break;
2433 case fractional:
2434 param2 = evalue_eval(&e->x.p->arr[0], values);
2435 mpz_fdiv_r(param2->x.n, param2->x.n, param2->d);
2436
2437 res = eval_polynomial(e->x.p, 1, param2, values);
2438 evalue_free(param2);
2439 break;
2440 case flooring:
2441 param2 = evalue_eval(&e->x.p->arr[0], values);
2442 mpz_fdiv_q(param2->x.n, param2->x.n, param2->d);
2443 value_set_si(param2->d, 1);
2444
2445 res = eval_polynomial(e->x.p, 1, param2, values);
2446 evalue_free(param2);
2447 break;
2448 case relation:
2449 param2 = evalue_eval(&e->x.p->arr[0], values);
2450 if (value_zero_p(param2->x.n))
2451 res = evalue_eval(&e->x.p->arr[1], values);
2452 else if (e->x.p->size > 2)
2453 res = evalue_eval(&e->x.p->arr[2], values);
2454 else
2455 res = evalue_zero();
2456 evalue_free(param2);
2457 break;
2458 case partition:
2459 assert(e->x.p->pos == EVALUE_DOMAIN(e->x.p->arr[0])->Dimension);
2460 for (i = 0; i < e->x.p->size/2; ++i)
2461 if (in_domain(EVALUE_DOMAIN(e->x.p->arr[2*i]), values)) {
2462 res = evalue_eval(&e->x.p->arr[2*i+1], values);
2463 break;
2464 }
2465 if (!res)
2466 res = evalue_zero();
2467 break;
2468 default:
2469 assert(0);
2470 }
2471 return res;
2472 }
2473
2474 size_t value_size(Value v) {
2475 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2476 * sizeof(v[0]._mp_d[0]);
2477 }
2478
2479 size_t domain_size(Polyhedron *D)
2480 {
2481 int i, j;
2482 size_t s = sizeof(*D);
2483
2484 for (i = 0; i < D->NbConstraints; ++i)
2485 for (j = 0; j < D->Dimension+2; ++j)
2486 s += value_size(D->Constraint[i][j]);
2487
2488 /*
2489 for (i = 0; i < D->NbRays; ++i)
2490 for (j = 0; j < D->Dimension+2; ++j)
2491 s += value_size(D->Ray[i][j]);
2492 */
2493
2494 return D->next ? s+domain_size(D->next) : s;
2495 }
2496
2497 size_t enode_size(enode *p) {
2498 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2499 int i;
2500
2501 if (p->type == partition)
2502 for (i = 0; i < p->size/2; ++i) {
2503 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2504 s += evalue_size(&p->arr[2*i+1]);
2505 }
2506 else
2507 for (i = 0; i < p->size; ++i) {
2508 s += evalue_size(&p->arr[i]);
2509 }
2510 return s;
2511 }
2512
2513 size_t evalue_size(evalue *e)
2514 {
2515 size_t s = sizeof(*e);
2516 s += value_size(e->d);
2517 if (value_notzero_p(e->d))
2518 s += value_size(e->x.n);
2519 else
2520 s += enode_size(e->x.p);
2521 return s;
2522 }
2523
2524 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2525 {
2526 evalue *found = NULL;
2527 evalue offset;
2528 evalue copy;
2529 int i;
2530
2531 if (value_pos_p(e->d) || e->x.p->type != fractional)
2532 return NULL;
2533
2534 value_init(offset.d);
2535 value_init(offset.x.n);
2536 poly_denom(&e->x.p->arr[0], &offset.d);
2537 value_lcm(offset.d, m, offset.d);
2538 value_set_si(offset.x.n, 1);
2539
2540 value_init(copy.d);
2541 evalue_copy(&copy, cst);
2542
2543 eadd(&offset, cst);
2544 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2545
2546 if (eequal(base, &e->x.p->arr[0]))
2547 found = &e->x.p->arr[0];
2548 else {
2549 value_set_si(offset.x.n, -2);
2550
2551 eadd(&offset, cst);
2552 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2553
2554 if (eequal(base, &e->x.p->arr[0]))
2555 found = base;
2556 }
2557 free_evalue_refs(cst);
2558 free_evalue_refs(&offset);
2559 *cst = copy;
2560
2561 for (i = 1; !found && i < e->x.p->size; ++i)
2562 found = find_second(base, cst, &e->x.p->arr[i], m);
2563
2564 return found;
2565 }
2566
2567 static evalue *find_relation_pair(evalue *e)
2568 {
2569 int i;
2570 evalue *found = NULL;
2571
2572 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2573 return NULL;
2574
2575 if (e->x.p->type == fractional) {
2576 Value m;
2577 evalue *cst;
2578
2579 value_init(m);
2580 poly_denom(&e->x.p->arr[0], &m);
2581
2582 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2583 cst = &cst->x.p->arr[0])
2584 ;
2585
2586 for (i = 1; !found && i < e->x.p->size; ++i)
2587 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2588
2589 value_clear(m);
2590 }
2591
2592 i = e->x.p->type == relation;
2593 for (; !found && i < e->x.p->size; ++i)
2594 found = find_relation_pair(&e->x.p->arr[i]);
2595
2596 return found;
2597 }
2598
2599 void evalue_mod2relation(evalue *e) {
2600 evalue *d;
2601
2602 if (value_zero_p(e->d) && e->x.p->type == partition) {
2603 int i;
2604
2605 for (i = 0; i < e->x.p->size/2; ++i) {
2606 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2607 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2608 value_clear(e->x.p->arr[2*i].d);
2609 free_evalue_refs(&e->x.p->arr[2*i+1]);
2610 e->x.p->size -= 2;
2611 if (2*i < e->x.p->size) {
2612 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2613 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2614 }
2615 --i;
2616 }
2617 }
2618 if (e->x.p->size == 0) {
2619 free(e->x.p);
2620 evalue_set_si(e, 0, 1);
2621 }
2622
2623 return;
2624 }
2625
2626 while ((d = find_relation_pair(e)) != NULL) {
2627 evalue split;
2628 evalue *ev;
2629
2630 value_init(split.d);
2631 value_set_si(split.d, 0);
2632 split.x.p = new_enode(relation, 3, 0);
2633 evalue_set_si(&split.x.p->arr[1], 1, 1);
2634 evalue_set_si(&split.x.p->arr[2], 1, 1);
2635
2636 ev = &split.x.p->arr[0];
2637 value_set_si(ev->d, 0);
2638 ev->x.p = new_enode(fractional, 3, -1);
2639 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2640 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2641 evalue_copy(&ev->x.p->arr[0], d);
2642
2643 emul(&split, e);
2644
2645 reduce_evalue(e);
2646
2647 free_evalue_refs(&split);
2648 }
2649 }
2650
2651 static int evalue_comp(const void * a, const void * b)
2652 {
2653 const evalue *e1 = *(const evalue **)a;
2654 const evalue *e2 = *(const evalue **)b;
2655 return ecmp(e1, e2);
2656 }
2657
2658 void evalue_combine(evalue *e)
2659 {
2660 evalue **evs;
2661 int i, k;
2662 enode *p;
2663 evalue tmp;
2664
2665 if (value_notzero_p(e->d) || e->x.p->type != partition)
2666 return;
2667
2668 NALLOC(evs, e->x.p->size/2);
2669 for (i = 0; i < e->x.p->size/2; ++i)
2670 evs[i] = &e->x.p->arr[2*i+1];
2671 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2672 p = new_enode(partition, e->x.p->size, e->x.p->pos);
2673 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2674 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2675 value_clear(p->arr[2*k].d);
2676 value_clear(p->arr[2*k+1].d);
2677 p->arr[2*k] = *(evs[i]-1);
2678 p->arr[2*k+1] = *(evs[i]);
2679 ++k;
2680 } else {
2681 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2682 Polyhedron *L = D;
2683
2684 value_clear((evs[i]-1)->d);
2685
2686 while (L->next)
2687 L = L->next;
2688 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2689 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2690 free_evalue_refs(evs[i]);
2691 }
2692 }
2693
2694 for (i = 2*k ; i < p->size; ++i)
2695 value_clear(p->arr[i].d);
2696
2697 free(evs);
2698 free(e->x.p);
2699 p->size = 2*k;
2700 e->x.p = p;
2701
2702 for (i = 0; i < e->x.p->size/2; ++i) {
2703 Polyhedron *H;
2704 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2705 continue;
2706 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2707 if (H == NULL)
2708 continue;
2709 for (k = 0; k < e->x.p->size/2; ++k) {
2710 Polyhedron *D, *N, **P;
2711 if (i == k)
2712 continue;
2713 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2714 D = *P;
2715 if (D == NULL)
2716 continue;
2717 for (; D; D = N) {
2718 *P = D;
2719 N = D->next;
2720 if (D->NbEq <= H->NbEq) {
2721 P = &D->next;
2722 continue;
2723 }
2724
2725 value_init(tmp.d);
2726 tmp.x.p = new_enode(partition, 2, e->x.p->pos);
2727 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2728 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2729 reduce_evalue(&tmp);
2730 if (value_notzero_p(tmp.d) ||
2731 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2732 P = &D->next;
2733 else {
2734 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2735 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2736 *P = NULL;
2737 }
2738 free_evalue_refs(&tmp);
2739 }
2740 }
2741 Polyhedron_Free(H);
2742 }
2743
2744 for (i = 0; i < e->x.p->size/2; ++i) {
2745 Polyhedron *H, *E;
2746 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2747 if (!D) {
2748 value_clear(e->x.p->arr[2*i].d);
2749 free_evalue_refs(&e->x.p->arr[2*i+1]);
2750 e->x.p->size -= 2;
2751 if (2*i < e->x.p->size) {
2752 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2753 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2754 }
2755 --i;
2756 continue;
2757 }
2758 if (!D->next)
2759 continue;
2760 H = DomainConvex(D, 0);
2761 E = DomainDifference(H, D, 0);
2762 Domain_Free(D);
2763 D = DomainDifference(H, E, 0);
2764 Domain_Free(H);
2765 Domain_Free(E);
2766 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2767 }
2768 }
2769
2770 /* Use smallest representative for coefficients in affine form in
2771 * argument of fractional.
2772 * Since any change will make the argument non-standard,
2773 * the containing evalue will have to be reduced again afterward.
2774 */
2775 static void fractional_minimal_coefficients(enode *p)
2776 {
2777 evalue *pp;
2778 Value twice;
2779 value_init(twice);
2780
2781 assert(p->type == fractional);
2782 pp = &p->arr[0];
2783 while (value_zero_p(pp->d)) {
2784 assert(pp->x.p->type == polynomial);
2785 assert(pp->x.p->size == 2);
2786 assert(value_notzero_p(pp->x.p->arr[1].d));
2787 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2788 if (value_gt(twice, pp->x.p->arr[1].d))
2789 value_subtract(pp->x.p->arr[1].x.n,
2790 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2791 pp = &pp->x.p->arr[0];
2792 }
2793
2794 value_clear(twice);
2795 }
2796
2797 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2798 Matrix **R)
2799 {
2800 Polyhedron *I, *H;
2801 evalue *pp;
2802 unsigned dim = D->Dimension;
2803 Matrix *T = Matrix_Alloc(2, dim+1);
2804 assert(T);
2805
2806 assert(p->type == fractional || p->type == flooring);
2807 value_set_si(T->p[1][dim], 1);
2808 evalue_extract_affine(&p->arr[0], T->p[0], &T->p[0][dim], d);
2809 I = DomainImage(D, T, 0);
2810 H = DomainConvex(I, 0);
2811 Domain_Free(I);
2812 if (R)
2813 *R = T;
2814 else
2815 Matrix_Free(T);
2816
2817 return H;
2818 }
2819
2820 static void replace_by_affine(evalue *e, Value offset)
2821 {
2822 enode *p;
2823 evalue inc;
2824
2825 p = e->x.p;
2826 value_init(inc.d);
2827 value_init(inc.x.n);
2828 value_set_si(inc.d, 1);
2829 value_oppose(inc.x.n, offset);
2830 eadd(&inc, &p->arr[0]);
2831 reorder_terms_about(p, &p->arr[0]); /* frees arr[0] */
2832 value_clear(e->d);
2833 *e = p->arr[1];
2834 free(p);
2835 free_evalue_refs(&inc);
2836 }
2837
2838 int evalue_range_reduction_in_domain(evalue *e, Polyhedron *D)
2839 {
2840 int i;
2841 enode *p;
2842 Value d, min, max;
2843 int r = 0;
2844 Polyhedron *I;
2845 int bounded;
2846
2847 if (value_notzero_p(e->d))
2848 return r;
2849
2850 p = e->x.p;
2851
2852 if (p->type == relation) {
2853 Matrix *T;
2854 int equal;
2855 value_init(d);
2856 value_init(min);
2857 value_init(max);
2858
2859 fractional_minimal_coefficients(p->arr[0].x.p);
2860 I = polynomial_projection(p->arr[0].x.p, D, &d, &T);
2861 bounded = line_minmax(I, &min, &max); /* frees I */
2862 equal = value_eq(min, max);
2863 mpz_cdiv_q(min, min, d);
2864 mpz_fdiv_q(max, max, d);
2865
2866 if (bounded && value_gt(min, max)) {
2867 /* Never zero */
2868 if (p->size == 3) {
2869 value_clear(e->d);
2870 *e = p->arr[2];
2871 } else {
2872 evalue_set_si(e, 0, 1);
2873 r = 1;
2874 }
2875 free_evalue_refs(&(p->arr[1]));
2876 free_evalue_refs(&(p->arr[0]));
2877 free(p);
2878 value_clear(d);
2879 value_clear(min);
2880 value_clear(max);
2881 Matrix_Free(T);
2882 return r ? r : evalue_range_reduction_in_domain(e, D);
2883 } else if (bounded && equal) {
2884 /* Always zero */
2885 if (p->size == 3)
2886 free_evalue_refs(&(p->arr[2]));
2887 value_clear(e->d);
2888 *e = p->arr[1];
2889 free_evalue_refs(&(p->arr[0]));
2890 free(p);
2891 value_clear(d);
2892 value_clear(min);
2893 value_clear(max);
2894 Matrix_Free(T);
2895 return evalue_range_reduction_in_domain(e, D);
2896 } else if (bounded && value_eq(min, max)) {
2897 /* zero for a single value */
2898 Polyhedron *E;
2899 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2900 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2901 value_multiply(min, min, d);
2902 value_subtract(M->p[0][D->Dimension+1],
2903 M->p[0][D->Dimension+1], min);
2904 E = DomainAddConstraints(D, M, 0);
2905 value_clear(d);
2906 value_clear(min);
2907 value_clear(max);
2908 Matrix_Free(T);
2909 Matrix_Free(M);
2910 r = evalue_range_reduction_in_domain(&p->arr[1], E);
2911 if (p->size == 3)
2912 r |= evalue_range_reduction_in_domain(&p->arr[2], D);
2913 Domain_Free(E);
2914 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2915 return r;
2916 }
2917
2918 value_clear(d);
2919 value_clear(min);
2920 value_clear(max);
2921 Matrix_Free(T);
2922 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2923 }
2924
2925 i = p->type == relation ? 1 :
2926 p->type == fractional ? 1 : 0;
2927 for (; i<p->size; i++)
2928 r |= evalue_range_reduction_in_domain(&p->arr[i], D);
2929
2930 if (p->type != fractional) {
2931 if (r && p->type == polynomial) {
2932 evalue f;
2933 value_init(f.d);
2934 value_set_si(f.d, 0);
2935 f.x.p = new_enode(polynomial, 2, p->pos);
2936 evalue_set_si(&f.x.p->arr[0], 0, 1);
2937 evalue_set_si(&f.x.p->arr[1], 1, 1);
2938 reorder_terms_about(p, &f);
2939 value_clear(e->d);
2940 *e = p->arr[0];
2941 free(p);
2942 }
2943 return r;
2944 }
2945
2946 value_init(d);
2947 value_init(min);
2948 value_init(max);
2949 fractional_minimal_coefficients(p);
2950 I = polynomial_projection(p, D, &d, NULL);
2951 bounded = line_minmax(I, &min, &max); /* frees I */
2952 mpz_fdiv_q(min, min, d);
2953 mpz_fdiv_q(max, max, d);
2954 value_subtract(d, max, min);
2955
2956 if (bounded && value_eq(min, max)) {
2957 replace_by_affine(e, min);
2958 r = 1;
2959 } else if (bounded && value_one_p(d) && p->size > 3) {
2960 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
2961 * See pages 199-200 of PhD thesis.
2962 */
2963 evalue rem;
2964 evalue inc;
2965 evalue t;
2966 evalue f;
2967 evalue factor;
2968 value_init(rem.d);
2969 value_set_si(rem.d, 0);
2970 rem.x.p = new_enode(fractional, 3, -1);
2971 evalue_copy(&rem.x.p->arr[0], &p->arr[0]);
2972 value_clear(rem.x.p->arr[1].d);
2973 value_clear(rem.x.p->arr[2].d);
2974 rem.x.p->arr[1] = p->arr[1];
2975 rem.x.p->arr[2] = p->arr[2];
2976 for (i = 3; i < p->size; ++i)
2977 p->arr[i-2] = p->arr[i];
2978 p->size -= 2;
2979
2980 value_init(inc.d);
2981 value_init(inc.x.n);
2982 value_set_si(inc.d, 1);
2983 value_oppose(inc.x.n, min);
2984
2985 value_init(t.d);
2986 evalue_copy(&t, &p->arr[0]);
2987 eadd(&inc, &t);
2988
2989 value_init(f.d);
2990 value_set_si(f.d, 0);
2991 f.x.p = new_enode(fractional, 3, -1);
2992 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2993 evalue_set_si(&f.x.p->arr[1], 1, 1);
2994 evalue_set_si(&f.x.p->arr[2], 2, 1);
2995
2996 value_init(factor.d);
2997 evalue_set_si(&factor, -1, 1);
2998 emul(&t, &factor);
2999
3000 eadd(&f, &factor);
3001 emul(&t, &factor);
3002
3003 value_clear(f.x.p->arr[1].x.n);
3004 value_clear(f.x.p->arr[2].x.n);
3005 evalue_set_si(&f.x.p->arr[1], 0, 1);
3006 evalue_set_si(&f.x.p->arr[2], -1, 1);
3007 eadd(&f, &factor);
3008
3009 if (r) {
3010 reorder_terms(&rem);
3011 reorder_terms(e);
3012 }
3013
3014 emul(&factor, e);
3015 eadd(&rem, e);
3016
3017 free_evalue_refs(&inc);
3018 free_evalue_refs(&t);
3019 free_evalue_refs(&f);
3020 free_evalue_refs(&factor);
3021 free_evalue_refs(&rem);
3022
3023 evalue_range_reduction_in_domain(e, D);
3024
3025 r = 1;
3026 } else {
3027 _reduce_evalue(&p->arr[0], 0, 1);
3028 if (r)
3029 reorder_terms(e);
3030 }
3031
3032 value_clear(d);
3033 value_clear(min);
3034 value_clear(max);
3035
3036 return r;
3037 }
3038
3039 void evalue_range_reduction(evalue *e)
3040 {
3041 int i;
3042 if (value_notzero_p(e->d) || e->x.p->type != partition)
3043 return;
3044
3045 for (i = 0; i < e->x.p->size/2; ++i)
3046 if (evalue_range_reduction_in_domain(&e->x.p->arr[2*i+1],
3047 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
3048 reduce_evalue(&e->x.p->arr[2*i+1]);
3049
3050 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
3051 free_evalue_refs(&e->x.p->arr[2*i+1]);
3052 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
3053 value_clear(e->x.p->arr[2*i].d);
3054 e->x.p->size -= 2;
3055 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
3056 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
3057 --i;
3058 }
3059 }
3060 }
3061
3062 /* Frees EP
3063 */
3064 Enumeration* partition2enumeration(evalue *EP)
3065 {
3066 int i;
3067 Enumeration *en, *res = NULL;
3068
3069 if (EVALUE_IS_ZERO(*EP)) {
3070 free(EP);
3071 return res;
3072 }
3073
3074 for (i = 0; i < EP->x.p->size/2; ++i) {
3075 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[2*i])->Dimension);
3076 en = (Enumeration *)malloc(sizeof(Enumeration));
3077 en->next = res;
3078 res = en;
3079 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
3080 value_clear(EP->x.p->arr[2*i].d);
3081 res->EP = EP->x.p->arr[2*i+1];
3082 }
3083 free(EP->x.p);
3084 value_clear(EP->d);
3085 free(EP);
3086 return res;
3087 }
3088
3089 int evalue_frac2floor_in_domain3(evalue *e, Polyhedron *D, int shift)
3090 {
3091 enode *p;
3092 int r = 0;
3093 int i;
3094 Polyhedron *I;
3095 Value d, min;
3096 evalue fl;
3097
3098 if (value_notzero_p(e->d))
3099 return r;
3100
3101 p = e->x.p;
3102
3103 i = p->type == relation ? 1 :
3104 p->type == fractional ? 1 : 0;
3105 for (; i<p->size; i++)
3106 r |= evalue_frac2floor_in_domain3(&p->arr[i], D, shift);
3107
3108 if (p->type != fractional) {
3109 if (r && p->type == polynomial) {
3110 evalue f;
3111 value_init(f.d);
3112 value_set_si(f.d, 0);
3113 f.x.p = new_enode(polynomial, 2, p->pos);
3114 evalue_set_si(&f.x.p->arr[0], 0, 1);
3115 evalue_set_si(&f.x.p->arr[1], 1, 1);
3116 reorder_terms_about(p, &f);
3117 value_clear(e->d);
3118 *e = p->arr[0];
3119 free(p);
3120 }
3121 return r;
3122 }
3123
3124 if (shift) {
3125 value_init(d);
3126 I = polynomial_projection(p, D, &d, NULL);
3127
3128 /*
3129 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3130 */
3131
3132 assert(I->NbEq == 0); /* Should have been reduced */
3133
3134 /* Find minimum */
3135 for (i = 0; i < I->NbConstraints; ++i)
3136 if (value_pos_p(I->Constraint[i][1]))
3137 break;
3138
3139 if (i < I->NbConstraints) {
3140 value_init(min);
3141 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
3142 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
3143 if (value_neg_p(min)) {
3144 evalue offset;
3145 mpz_fdiv_q(min, min, d);
3146 value_init(offset.d);
3147 value_set_si(offset.d, 1);
3148 value_init(offset.x.n);
3149 value_oppose(offset.x.n, min);
3150 eadd(&offset, &p->arr[0]);
3151 free_evalue_refs(&offset);
3152 }
3153 value_clear(min);
3154 }
3155
3156 Polyhedron_Free(I);
3157 value_clear(d);
3158 }
3159
3160 value_init(fl.d);
3161 value_set_si(fl.d, 0);
3162 fl.x.p = new_enode(flooring, 3, -1);
3163 evalue_set_si(&fl.x.p->arr[1], 0, 1);
3164 evalue_set_si(&fl.x.p->arr[2], -1, 1);
3165 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
3166
3167 eadd(&fl, &p->arr[0]);
3168 reorder_terms_about(p, &p->arr[0]);
3169 value_clear(e->d);
3170 *e = p->arr[1];
3171 free(p);
3172 free_evalue_refs(&fl);
3173
3174 return 1;
3175 }
3176
3177 int evalue_frac2floor_in_domain(evalue *e, Polyhedron *D)
3178 {
3179 return evalue_frac2floor_in_domain3(e, D, 1);
3180 }
3181
3182 void evalue_frac2floor2(evalue *e, int shift)
3183 {
3184 int i;
3185 if (value_notzero_p(e->d) || e->x.p->type != partition) {
3186 if (!shift) {
3187 if (evalue_frac2floor_in_domain3(e, NULL, 0))
3188 reduce_evalue(e);
3189 }
3190 return;
3191 }
3192
3193 for (i = 0; i < e->x.p->size/2; ++i)
3194 if (evalue_frac2floor_in_domain3(&e->x.p->arr[2*i+1],
3195 EVALUE_DOMAIN(e->x.p->arr[2*i]), shift))
3196 reduce_evalue(&e->x.p->arr[2*i+1]);
3197 }
3198
3199 void evalue_frac2floor(evalue *e)
3200 {
3201 evalue_frac2floor2(e, 1);
3202 }
3203
3204 /* Add a new variable with lower bound 1 and upper bound specified
3205 * by row. If negative is true, then the new variable has upper
3206 * bound -1 and lower bound specified by row.
3207 */
3208 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
3209 Vector *row, int negative)
3210 {
3211 int nr, nc;
3212 int i;
3213 int nparam = D->Dimension - nvar;
3214
3215 if (C == 0) {
3216 nr = D->NbConstraints + 2;
3217 nc = D->Dimension + 2 + 1;
3218 C = Matrix_Alloc(nr, nc);
3219 for (i = 0; i < D->NbConstraints; ++i) {
3220 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
3221 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3222 D->Dimension + 1 - nvar);
3223 }
3224 } else {
3225 Matrix *oldC = C;
3226 nr = C->NbRows + 2;
3227 nc = C->NbColumns + 1;
3228 C = Matrix_Alloc(nr, nc);
3229 for (i = 0; i < oldC->NbRows; ++i) {
3230 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
3231 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3232 oldC->NbColumns - 1 - nvar);
3233 }
3234 }
3235 value_set_si(C->p[nr-2][0], 1);
3236 if (negative)
3237 value_set_si(C->p[nr-2][1 + nvar], -1);
3238 else
3239 value_set_si(C->p[nr-2][1 + nvar], 1);
3240 value_set_si(C->p[nr-2][nc - 1], -1);
3241
3242 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
3243 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
3244 1 + nparam);
3245
3246 return C;
3247 }
3248
3249 static void floor2frac_r(evalue *e, int nvar)
3250 {
3251 enode *p;
3252 int i;
3253 evalue f;
3254 evalue *pp;
3255
3256 if (value_notzero_p(e->d))
3257 return;
3258
3259 p = e->x.p;
3260
3261 assert(p->type == flooring);
3262 for (i = 1; i < p->size; i++)
3263 floor2frac_r(&p->arr[i], nvar);
3264
3265 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
3266 assert(pp->x.p->type == polynomial);
3267 pp->x.p->pos -= nvar;
3268 }
3269
3270 value_init(f.d);
3271 value_set_si(f.d, 0);
3272 f.x.p = new_enode(fractional, 3, -1);
3273 evalue_set_si(&f.x.p->arr[1], 0, 1);
3274 evalue_set_si(&f.x.p->arr[2], -1, 1);
3275 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3276
3277 eadd(&f, &p->arr[0]);
3278 reorder_terms_about(p, &p->arr[0]);
3279 value_clear(e->d);
3280 *e = p->arr[1];
3281 free(p);
3282 free_evalue_refs(&f);
3283 }
3284
3285 /* Convert flooring back to fractional and shift position
3286 * of the parameters by nvar
3287 */
3288 static void floor2frac(evalue *e, int nvar)
3289 {
3290 floor2frac_r(e, nvar);
3291 reduce_evalue(e);
3292 }
3293
3294 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
3295 {
3296 evalue *t;
3297 int nparam = D->Dimension - nvar;
3298
3299 if (C != 0) {
3300 C = Matrix_Copy(C);
3301 D = Constraints2Polyhedron(C, 0);
3302 Matrix_Free(C);
3303 }
3304
3305 t = barvinok_enumerate_e(D, 0, nparam, 0);
3306
3307 /* Double check that D was not unbounded. */
3308 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3309
3310 if (C != 0)
3311 Polyhedron_Free(D);
3312
3313 return t;
3314 }
3315
3316 static void domain_signs(Polyhedron *D, int *signs)
3317 {
3318 int j, k;
3319
3320 POL_ENSURE_VERTICES(D);
3321 for (j = 0; j < D->Dimension; ++j) {
3322 signs[j] = 0;
3323 for (k = 0; k < D->NbRays; ++k) {
3324 signs[j] = value_sign(D->Ray[k][1+j]);
3325 if (signs[j])
3326 break;
3327 }
3328 }
3329 }
3330
3331 static evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3332 int *signs, Matrix *C, unsigned MaxRays)
3333 {
3334 Vector *row = NULL;
3335 int i, offset;
3336 evalue *res;
3337 Matrix *origC;
3338 evalue *factor = NULL;
3339 evalue cum;
3340 int negative = 0;
3341
3342 if (EVALUE_IS_ZERO(*e))
3343 return 0;
3344
3345 if (D->next) {
3346 Polyhedron *DD = Disjoint_Domain(D, 0, MaxRays);
3347 Polyhedron *Q;
3348
3349 Q = DD;
3350 DD = Q->next;
3351 Q->next = 0;
3352
3353 res = esum_over_domain(e, nvar, Q, signs, C, MaxRays);
3354 Polyhedron_Free(Q);
3355
3356 for (Q = DD; Q; Q = DD) {
3357 evalue *t;
3358
3359 DD = Q->next;
3360 Q->next = 0;
3361
3362 t = esum_over_domain(e, nvar, Q, signs, C, MaxRays);
3363 Polyhedron_Free(Q);
3364
3365 if (!res)
3366 res = t;
3367 else if (t) {
3368 eadd(t, res);
3369 evalue_free(t);
3370 }
3371 }
3372 return res;
3373 }
3374
3375 if (value_notzero_p(e->d)) {
3376 evalue *t;
3377
3378 t = esum_over_domain_cst(nvar, D, C);
3379
3380 if (!EVALUE_IS_ONE(*e))
3381 emul(e, t);
3382
3383 return t;
3384 }
3385
3386 if (!signs) {
3387 signs = alloca(sizeof(int) * D->Dimension);
3388 domain_signs(D, signs);
3389 }
3390
3391 switch (e->x.p->type) {
3392 case flooring: {
3393 evalue *pp = &e->x.p->arr[0];
3394
3395 if (pp->x.p->pos > nvar) {
3396 /* remainder is independent of the summated vars */
3397 evalue f;
3398 evalue *t;
3399
3400 value_init(f.d);
3401 evalue_copy(&f, e);
3402 floor2frac(&f, nvar);
3403
3404 t = esum_over_domain_cst(nvar, D, C);
3405
3406 emul(&f, t);
3407
3408 free_evalue_refs(&f);
3409
3410 return t;
3411 }
3412
3413 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3414 poly_denom(pp, &row->p[1 + nvar]);
3415 value_set_si(row->p[0], 1);
3416 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3417 pp = &pp->x.p->arr[0]) {
3418 int pos;
3419 assert(pp->x.p->type == polynomial);
3420 pos = pp->x.p->pos;
3421 if (pos >= 1 + nvar)
3422 ++pos;
3423 value_assign(row->p[pos], row->p[1+nvar]);
3424 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3425 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3426 }
3427 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3428 value_division(row->p[1 + D->Dimension + 1],
3429 row->p[1 + D->Dimension + 1],
3430 pp->d);
3431 value_multiply(row->p[1 + D->Dimension + 1],
3432 row->p[1 + D->Dimension + 1],
3433 pp->x.n);
3434 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3435 break;
3436 }
3437 case polynomial: {
3438 int pos = e->x.p->pos;
3439
3440 if (pos > nvar) {
3441 factor = ALLOC(evalue);
3442 value_init(factor->d);
3443 value_set_si(factor->d, 0);
3444 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3445 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3446 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3447 break;
3448 }
3449
3450 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3451 negative = signs[pos-1] < 0;
3452 value_set_si(row->p[0], 1);
3453 if (negative) {
3454 value_set_si(row->p[pos], -1);
3455 value_set_si(row->p[1 + nvar], 1);
3456 } else {
3457 value_set_si(row->p[pos], 1);
3458 value_set_si(row->p[1 + nvar], -1);
3459 }
3460 break;
3461 }
3462 default:
3463 assert(0);
3464 }
3465
3466 offset = type_offset(e->x.p);
3467
3468 res = esum_over_domain(&e->x.p->arr[offset], nvar, D, signs, C, MaxRays);
3469
3470 if (factor) {
3471 value_init(cum.d);
3472 evalue_copy(&cum, factor);
3473 }
3474
3475 origC = C;
3476 for (i = 1; offset+i < e->x.p->size; ++i) {
3477 evalue *t;
3478 if (row) {
3479 Matrix *prevC = C;
3480 C = esum_add_constraint(nvar, D, C, row, negative);
3481 if (prevC != origC)
3482 Matrix_Free(prevC);
3483 }
3484 /*
3485 if (row)
3486 Vector_Print(stderr, P_VALUE_FMT, row);
3487 if (C)
3488 Matrix_Print(stderr, P_VALUE_FMT, C);
3489 */
3490 t = esum_over_domain(&e->x.p->arr[offset+i], nvar, D, signs, C, MaxRays);
3491
3492 if (t) {
3493 if (factor)
3494 emul(&cum, t);
3495 if (negative && (i % 2))
3496 evalue_negate(t);
3497 }
3498
3499 if (!res)
3500 res = t;
3501 else if (t) {
3502 eadd(t, res);
3503 evalue_free(t);
3504 }
3505 if (factor && offset+i+1 < e->x.p->size)
3506 emul(factor, &cum);
3507 }
3508 if (C != origC)
3509 Matrix_Free(C);
3510
3511 if (factor) {
3512 free_evalue_refs(&cum);
3513 evalue_free(factor);
3514 }
3515
3516 if (row)
3517 Vector_Free(row);
3518
3519 reduce_evalue(res);
3520
3521 return res;
3522 }
3523
3524 static Polyhedron_Insert(Polyhedron ***next, Polyhedron *Q)
3525 {
3526 if (emptyQ(Q))
3527 Polyhedron_Free(Q);
3528 else {
3529 **next = Q;
3530 *next = &(Q->next);
3531 }
3532 }
3533
3534 static Polyhedron *Polyhedron_Split_Into_Orthants(Polyhedron *P,
3535 unsigned MaxRays)
3536 {
3537 int i = 0;
3538 Polyhedron *D = P;
3539 Vector *c = Vector_Alloc(1 + P->Dimension + 1);
3540 value_set_si(c->p[0], 1);
3541
3542 if (P->Dimension == 0)
3543 return Polyhedron_Copy(P);
3544
3545 for (i = 0; i < P->Dimension; ++i) {
3546 Polyhedron *L = NULL;
3547 Polyhedron **next = &L;
3548 Polyhedron *I;
3549
3550 for (I = D; I; I = I->next) {
3551 Polyhedron *Q;
3552 value_set_si(c->p[1+i], 1);
3553 value_set_si(c->p[1+P->Dimension], 0);
3554 Q = AddConstraints(c->p, 1, I, MaxRays);
3555 Polyhedron_Insert(&next, Q);
3556 value_set_si(c->p[1+i], -1);
3557 value_set_si(c->p[1+P->Dimension], -1);
3558 Q = AddConstraints(c->p, 1, I, MaxRays);
3559 Polyhedron_Insert(&next, Q);
3560 value_set_si(c->p[1+i], 0);
3561 }
3562 if (D != P)
3563 Domain_Free(D);
3564 D = L;
3565 }
3566 Vector_Free(c);
3567 return D;
3568 }
3569
3570 /* Make arguments of all floors non-negative */
3571 static void shift_floor_in_domain(evalue *e, Polyhedron *D)
3572 {
3573 Value d, m;
3574 Polyhedron *I;
3575 int i;
3576 enode *p;
3577
3578 if (value_notzero_p(e->d))
3579 return;
3580
3581 p = e->x.p;
3582
3583 for (i = type_offset(p); i < p->size; ++i)
3584 shift_floor_in_domain(&p->arr[i], D);
3585
3586 if (p->type != flooring)
3587 return;
3588
3589 value_init(d);
3590 value_init(m);
3591
3592 I = polynomial_projection(p, D, &d, NULL);
3593 assert(I->NbEq == 0); /* Should have been reduced */
3594
3595 for (i = 0; i < I->NbConstraints; ++i)
3596 if (value_pos_p(I->Constraint[i][1]))
3597 break;
3598 assert(i < I->NbConstraints);
3599 if (i < I->NbConstraints) {
3600 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
3601 mpz_fdiv_q(m, I->Constraint[i][2], I->Constraint[i][1]);
3602 if (value_neg_p(m)) {
3603 /* replace [e] by [e-m]+m such that e-m >= 0 */
3604 evalue f;
3605
3606 value_init(f.d);
3607 value_init(f.x.n);
3608 value_set_si(f.d, 1);
3609 value_oppose(f.x.n, m);
3610 eadd(&f, &p->arr[0]);
3611 value_clear(f.x.n);
3612
3613 value_set_si(f.d, 0);
3614 f.x.p = new_enode(flooring, 3, -1);
3615 value_clear(f.x.p->arr[0].d);
3616 f.x.p->arr[0] = p->arr[0];
3617 evalue_set_si(&f.x.p->arr[2], 1, 1);
3618 value_set_si(f.x.p->arr[1].d, 1);
3619 value_init(f.x.p->arr[1].x.n);
3620 value_assign(f.x.p->arr[1].x.n, m);
3621 reorder_terms_about(p, &f);
3622 value_clear(e->d);
3623 *e = p->arr[1];
3624 free(p);
3625 }
3626 }
3627 Polyhedron_Free(I);
3628 value_clear(d);
3629 value_clear(m);
3630 }
3631
3632 evalue *box_summate(Polyhedron *P, evalue *E, unsigned nvar, unsigned MaxRays)
3633 {
3634 evalue *sum = evalue_zero();
3635 Polyhedron *D = Polyhedron_Split_Into_Orthants(P, MaxRays);
3636 for (P = D; P; P = P->next) {
3637 evalue *t;
3638 evalue *fe = evalue_dup(E);
3639 Polyhedron *next = P->next;
3640 P->next = NULL;
3641 reduce_evalue_in_domain(fe, P);
3642 evalue_frac2floor2(fe, 0);
3643 shift_floor_in_domain(fe, P);
3644 t = esum_over_domain(fe, nvar, P, NULL, NULL, MaxRays);
3645 if (t) {
3646 eadd(t, sum);
3647 evalue_free(t);
3648 }
3649 evalue_free(fe);
3650 P->next = next;
3651 }
3652 Domain_Free(D);
3653 return sum;
3654 }
3655
3656 /* Initial silly implementation */
3657 void eor(evalue *e1, evalue *res)
3658 {
3659 evalue E;
3660 evalue mone;
3661 value_init(E.d);
3662 value_init(mone.d);
3663 evalue_set_si(&mone, -1, 1);
3664
3665 evalue_copy(&E, res);
3666 eadd(e1, &E);
3667 emul(e1, res);
3668 emul(&mone, res);
3669 eadd(&E, res);
3670
3671 free_evalue_refs(&E);
3672 free_evalue_refs(&mone);
3673 }
3674
3675 /* computes denominator of polynomial evalue
3676 * d should point to a value initialized to 1
3677 */
3678 void evalue_denom(const evalue *e, Value *d)
3679 {
3680 int i, offset;
3681
3682 if (value_notzero_p(e->d)) {
3683 value_lcm(*d, *d, e->d);
3684 return;
3685 }
3686 assert(e->x.p->type == polynomial ||
3687 e->x.p->type == fractional ||
3688 e->x.p->type == flooring);
3689 offset = type_offset(e->x.p);
3690 for (i = e->x.p->size-1; i >= offset; --i)
3691 evalue_denom(&e->x.p->arr[i], d);
3692 }
3693
3694 /* Divides the evalue e by the integer n */
3695 void evalue_div(evalue *e, Value n)
3696 {
3697 int i, offset;
3698
3699 if (value_one_p(n) || EVALUE_IS_ZERO(*e))
3700 return;
3701
3702 if (value_notzero_p(e->d)) {
3703 Value gc;
3704 value_init(gc);
3705 value_multiply(e->d, e->d, n);
3706 value_gcd(gc, e->x.n, e->d);
3707 if (value_notone_p(gc)) {
3708 value_division(e->d, e->d, gc);
3709 value_division(e->x.n, e->x.n, gc);
3710 }
3711 value_clear(gc);
3712 return;
3713 }
3714 if (e->x.p->type == partition) {
3715 for (i = 0; i < e->x.p->size/2; ++i)
3716 evalue_div(&e->x.p->arr[2*i+1], n);
3717 return;
3718 }
3719 offset = type_offset(e->x.p);
3720 for (i = e->x.p->size-1; i >= offset; --i)
3721 evalue_div(&e->x.p->arr[i], n);
3722 }
3723
3724 /* Multiplies the evalue e by the integer n */
3725 void evalue_mul(evalue *e, Value n)
3726 {
3727 int i, offset;
3728
3729 if (value_one_p(n) || EVALUE_IS_ZERO(*e))
3730 return;
3731
3732 if (value_notzero_p(e->d)) {
3733 Value gc;
3734 value_init(gc);
3735 value_multiply(e->x.n, e->x.n, n);
3736 value_gcd(gc, e->x.n, e->d);
3737 if (value_notone_p(gc)) {
3738 value_division(e->d, e->d, gc);
3739 value_division(e->x.n, e->x.n, gc);
3740 }
3741 value_clear(gc);
3742 return;
3743 }
3744 if (e->x.p->type == partition) {
3745 for (i = 0; i < e->x.p->size/2; ++i)
3746 evalue_mul(&e->x.p->arr[2*i+1], n);
3747 return;
3748 }
3749 offset = type_offset(e->x.p);
3750 for (i = e->x.p->size-1; i >= offset; --i)
3751 evalue_mul(&e->x.p->arr[i], n);
3752 }
3753
3754 /* Multiplies the evalue e by the n/d */
3755 void evalue_mul_div(evalue *e, Value n, Value d)
3756 {
3757 int i, offset;
3758
3759 if ((value_one_p(n) && value_one_p(d)) || EVALUE_IS_ZERO(*e))
3760 return;
3761
3762 if (value_notzero_p(e->d)) {
3763 Value gc;
3764 value_init(gc);
3765 value_multiply(e->x.n, e->x.n, n);
3766 value_multiply(e->d, e->d, d);
3767 value_gcd(gc, e->x.n, e->d);
3768 if (value_notone_p(gc)) {
3769 value_division(e->d, e->d, gc);
3770 value_division(e->x.n, e->x.n, gc);
3771 }
3772 value_clear(gc);
3773 return;
3774 }
3775 if (e->x.p->type == partition) {
3776 for (i = 0; i < e->x.p->size/2; ++i)
3777 evalue_mul_div(&e->x.p->arr[2*i+1], n, d);
3778 return;
3779 }
3780 offset = type_offset(e->x.p);
3781 for (i = e->x.p->size-1; i >= offset; --i)
3782 evalue_mul_div(&e->x.p->arr[i], n, d);
3783 }
3784
3785 void evalue_negate(evalue *e)
3786 {
3787 int i, offset;
3788
3789 if (value_notzero_p(e->d)) {
3790 value_oppose(e->x.n, e->x.n);
3791 return;
3792 }
3793 if (e->x.p->type == partition) {
3794 for (i = 0; i < e->x.p->size/2; ++i)
3795 evalue_negate(&e->x.p->arr[2*i+1]);
3796 return;
3797 }
3798 offset = type_offset(e->x.p);
3799 for (i = e->x.p->size-1; i >= offset; --i)
3800 evalue_negate(&e->x.p->arr[i]);
3801 }
3802
3803 void evalue_add_constant(evalue *e, const Value cst)
3804 {
3805 int i;
3806
3807 if (value_zero_p(e->d)) {
3808 if (e->x.p->type == partition) {
3809 for (i = 0; i < e->x.p->size/2; ++i)
3810 evalue_add_constant(&e->x.p->arr[2*i+1], cst);
3811 return;
3812 }
3813 if (e->x.p->type == relation) {
3814 for (i = 1; i < e->x.p->size; ++i)
3815 evalue_add_constant(&e->x.p->arr[i], cst);
3816 return;
3817 }
3818 do {
3819 e = &e->x.p->arr[type_offset(e->x.p)];
3820 } while (value_zero_p(e->d));
3821 }
3822 value_addmul(e->x.n, cst, e->d);
3823 }
3824
3825 static void evalue_frac2polynomial_r(evalue *e, int *signs, int sign, int in_frac)
3826 {
3827 int i, offset;
3828 Value d;
3829 enode *p;
3830 int sign_odd = sign;
3831
3832 if (value_notzero_p(e->d)) {
3833 if (in_frac && sign * value_sign(e->x.n) < 0) {
3834 value_set_si(e->x.n, 0);
3835 value_set_si(e->d, 1);
3836 }
3837 return;
3838 }
3839
3840 if (e->x.p->type == relation) {
3841 for (i = e->x.p->size-1; i >= 1; --i)
3842 evalue_frac2polynomial_r(&e->x.p->arr[i], signs, sign, in_frac);
3843 return;
3844 }
3845
3846 if (e->x.p->type == polynomial)
3847 sign_odd *= signs[e->x.p->pos-1];
3848 offset = type_offset(e->x.p);
3849 evalue_frac2polynomial_r(&e->x.p->arr[offset], signs, sign, in_frac);
3850 in_frac |= e->x.p->type == fractional;
3851 for (i = e->x.p->size-1; i > offset; --i)
3852 evalue_frac2polynomial_r(&e->x.p->arr[i], signs,
3853 (i - offset) % 2 ? sign_odd : sign, in_frac);
3854
3855 if (e->x.p->type != fractional)
3856 return;
3857
3858 /* replace { a/m } by (m-1)/m if sign != 0
3859 * and by (m-1)/(2m) if sign == 0
3860 */
3861 value_init(d);
3862 value_set_si(d, 1);
3863 evalue_denom(&e->x.p->arr[0], &d);
3864 free_evalue_refs(&e->x.p->arr[0]);
3865 value_init(e->x.p->arr[0].d);
3866 value_init(e->x.p->arr[0].x.n);
3867 if (sign == 0)
3868 value_addto(e->x.p->arr[0].d, d, d);
3869 else
3870 value_assign(e->x.p->arr[0].d, d);
3871 value_decrement(e->x.p->arr[0].x.n, d);
3872 value_clear(d);
3873
3874 p = e->x.p;
3875 reorder_terms_about(p, &p->arr[0]);
3876 value_clear(e->d);
3877 *e = p->arr[1];
3878 free(p);
3879 }
3880
3881 /* Approximate the evalue in fractional representation by a polynomial.
3882 * If sign > 0, the result is an upper bound;
3883 * if sign < 0, the result is a lower bound;
3884 * if sign = 0, the result is an intermediate approximation.
3885 */
3886 void evalue_frac2polynomial(evalue *e, int sign, unsigned MaxRays)
3887 {
3888 int i, dim;
3889 int *signs;
3890
3891 if (value_notzero_p(e->d))
3892 return;
3893 assert(e->x.p->type == partition);
3894 /* make sure all variables in the domains have a fixed sign */
3895 if (sign) {
3896 evalue_split_domains_into_orthants(e, MaxRays);
3897 if (EVALUE_IS_ZERO(*e))
3898 return;
3899 }
3900
3901 assert(e->x.p->size >= 2);
3902 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3903
3904 signs = alloca(sizeof(int) * dim);
3905
3906 if (!sign)
3907 for (i = 0; i < dim; ++i)
3908 signs[i] = 0;
3909 for (i = 0; i < e->x.p->size/2; ++i) {
3910 if (sign)
3911 domain_signs(EVALUE_DOMAIN(e->x.p->arr[2*i]), signs);
3912 evalue_frac2polynomial_r(&e->x.p->arr[2*i+1], signs, sign, 0);
3913 }
3914 }
3915
3916 /* Split the domains of e (which is assumed to be a partition)
3917 * such that each resulting domain lies entirely in one orthant.
3918 */
3919 void evalue_split_domains_into_orthants(evalue *e, unsigned MaxRays)
3920 {
3921 int i, dim;
3922 assert(value_zero_p(e->d));
3923 assert(e->x.p->type == partition);
3924 assert(e->x.p->size >= 2);
3925 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3926
3927 for (i = 0; i < dim; ++i) {
3928 evalue split;
3929 Matrix *C, *C2;
3930 C = Matrix_Alloc(1, 1 + dim + 1);
3931 value_set_si(C->p[0][0], 1);
3932 value_init(split.d);
3933 value_set_si(split.d, 0);
3934 split.x.p = new_enode(partition, 4, dim);
3935 value_set_si(C->p[0][1+i], 1);
3936 C2 = Matrix_Copy(C);
3937 EVALUE_SET_DOMAIN(split.x.p->arr[0], Constraints2Polyhedron(C2, MaxRays));
3938 Matrix_Free(C2);
3939 evalue_set_si(&split.x.p->arr[1], 1, 1);
3940 value_set_si(C->p[0][1+i], -1);
3941 value_set_si(C->p[0][1+dim], -1);
3942 EVALUE_SET_DOMAIN(split.x.p->arr[2], Constraints2Polyhedron(C, MaxRays));
3943 evalue_set_si(&split.x.p->arr[3], 1, 1);
3944 emul(&split, e);
3945 free_evalue_refs(&split);
3946 Matrix_Free(C);
3947 }
3948 }
3949
3950 static evalue *find_fractional_with_max_periods(evalue *e, Polyhedron *D,
3951 int max_periods,
3952 Matrix **TT,
3953 Value *min, Value *max)
3954 {
3955 Matrix *T;
3956 evalue *f = NULL;
3957 Value d;
3958 int i;
3959
3960 if (value_notzero_p(e->d))
3961 return NULL;
3962
3963 if (e->x.p->type == fractional) {
3964 Polyhedron *I;
3965 int bounded;
3966
3967 value_init(d);
3968 I = polynomial_projection(e->x.p, D, &d, &T);
3969 bounded = line_minmax(I, min, max); /* frees I */
3970 if (bounded) {
3971 Value mp;
3972 value_init(mp);
3973 value_set_si(mp, max_periods);
3974 mpz_fdiv_q(*min, *min, d);
3975 mpz_fdiv_q(*max, *max, d);
3976 value_assign(T->p[1][D->Dimension], d);
3977 value_subtract(d, *max, *min);
3978 if (value_ge(d, mp))
3979 Matrix_Free(T);
3980 else
3981 f = evalue_dup(&e->x.p->arr[0]);
3982 value_clear(mp);
3983 } else
3984 Matrix_Free(T);
3985 value_clear(d);
3986 if (f) {
3987 *TT = T;
3988 return f;
3989 }
3990 }
3991
3992 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
3993 if ((f = find_fractional_with_max_periods(&e->x.p->arr[i], D, max_periods,
3994 TT, min, max)))
3995 return f;
3996
3997 return NULL;
3998 }
3999
4000 static void replace_fract_by_affine(evalue *e, evalue *f, Value val)
4001 {
4002 int i, offset;
4003
4004 if (value_notzero_p(e->d))
4005 return;
4006
4007 offset = type_offset(e->x.p);
4008 for (i = e->x.p->size-1; i >= offset; --i)
4009 replace_fract_by_affine(&e->x.p->arr[i], f, val);
4010
4011 if (e->x.p->type != fractional)
4012 return;
4013
4014 if (!eequal(&e->x.p->arr[0], f))
4015 return;
4016
4017 replace_by_affine(e, val);
4018 }
4019
4020 /* Look for fractional parts that can be removed by splitting the corresponding
4021 * domain into at most max_periods parts.
4022 * We use a very simply strategy that looks for the first fractional part
4023 * that satisfies the condition, performs the split and then continues
4024 * looking for other fractional parts in the split domains until no
4025 * such fractional part can be found anymore.
4026 */
4027 void evalue_split_periods(evalue *e, int max_periods, unsigned int MaxRays)
4028 {
4029 int i, j, n;
4030 Value min;
4031 Value max;
4032 Value d;
4033
4034 if (EVALUE_IS_ZERO(*e))
4035 return;
4036 if (value_notzero_p(e->d) || e->x.p->type != partition) {
4037 fprintf(stderr,
4038 "WARNING: evalue_split_periods called on incorrect evalue type\n");
4039 return;
4040 }
4041
4042 value_init(min);
4043 value_init(max);
4044 value_init(d);
4045
4046 for (i = 0; i < e->x.p->size/2; ++i) {
4047 enode *p;
4048 evalue *f;
4049 Matrix *T;
4050 Matrix *M;
4051 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
4052 Polyhedron *E;
4053 f = find_fractional_with_max_periods(&e->x.p->arr[2*i+1], D, max_periods,
4054 &T, &min, &max);
4055 if (!f)
4056 continue;
4057
4058 M = Matrix_Alloc(2, 2+D->Dimension);
4059
4060 value_subtract(d, max, min);
4061 n = VALUE_TO_INT(d)+1;
4062
4063 value_set_si(M->p[0][0], 1);
4064 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
4065 value_multiply(d, max, T->p[1][D->Dimension]);
4066 value_subtract(M->p[0][1+D->Dimension], M->p[0][1+D->Dimension], d);
4067 value_set_si(d, -1);
4068 value_set_si(M->p[1][0], 1);
4069 Vector_Scale(T->p[0], M->p[1]+1, d, D->Dimension+1);
4070 value_addmul(M->p[1][1+D->Dimension], max, T->p[1][D->Dimension]);
4071 value_addto(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension],
4072 T->p[1][D->Dimension]);
4073 value_decrement(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension]);
4074
4075 p = new_enode(partition, e->x.p->size + (n-1)*2, e->x.p->pos);
4076 for (j = 0; j < 2*i; ++j) {
4077 value_clear(p->arr[j].d);
4078 p->arr[j] = e->x.p->arr[j];
4079 }
4080 for (j = 2*i+2; j < e->x.p->size; ++j) {
4081 value_clear(p->arr[j+2*(n-1)].d);
4082 p->arr[j+2*(n-1)] = e->x.p->arr[j];
4083 }
4084 for (j = n-1; j >= 0; --j) {
4085 if (j == 0) {
4086 value_clear(p->arr[2*i+1].d);
4087 p->arr[2*i+1] = e->x.p->arr[2*i+1];
4088 } else
4089 evalue_copy(&p->arr[2*(i+j)+1], &e->x.p->arr[2*i+1]);
4090 if (j != n-1) {
4091 value_subtract(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension],
4092 T->p[1][D->Dimension]);
4093 value_addto(M->p[0][1+D->Dimension], M->p[0][1+D->Dimension],
4094 T->p[1][D->Dimension]);
4095 }
4096 replace_fract_by_affine(&p->arr[2*(i+j)+1], f, max);
4097 E = DomainAddConstraints(D, M, MaxRays);
4098 EVALUE_SET_DOMAIN(p->arr[2*(i+j)], E);
4099 if (evalue_range_reduction_in_domain(&p->arr[2*(i+j)+1], E))
4100 reduce_evalue(&p->arr[2*(i+j)+1]);
4101 value_decrement(max, max);
4102 }
4103 value_clear(e->x.p->arr[2*i].d);
4104 Domain_Free(D);
4105 Matrix_Free(M);
4106 Matrix_Free(T);
4107 evalue_free(f);
4108 free(e->x.p);
4109 e->x.p = p;
4110 --i;
4111 }
4112
4113 value_clear(d);
4114 value_clear(min);
4115 value_clear(max);
4116 }
4117
4118 void evalue_extract_affine(const evalue *e, Value *coeff, Value *cst, Value *d)
4119 {
4120 value_set_si(*d, 1);
4121 evalue_denom(e, d);
4122 for ( ; value_zero_p(e->d); e = &e->x.p->arr[0]) {
4123 evalue *c;
4124 assert(e->x.p->type == polynomial);
4125 assert(e->x.p->size == 2);
4126 c = &e->x.p->arr[1];
4127 value_multiply(coeff[e->x.p->pos-1], *d, c->x.n);
4128 value_division(coeff[e->x.p->pos-1], coeff[e->x.p->pos-1], c->d);
4129 }
4130 value_multiply(*cst, *d, e->x.n);
4131 value_division(*cst, *cst, e->d);
4132 }
4133
4134 /* returns an evalue that corresponds to
4135 *
4136 * c/den x_param
4137 */
4138 static evalue *term(int param, Value c, Value den)
4139 {
4140 evalue *EP = ALLOC(evalue);
4141 value_init(EP->d);
4142 value_set_si(EP->d,0);
4143 EP->x.p = new_enode(polynomial, 2, param + 1);
4144 evalue_set_si(&EP->x.p->arr[0], 0, 1);
4145 value_init(EP->x.p->arr[1].x.n);
4146 value_assign(EP->x.p->arr[1].d, den);
4147 value_assign(EP->x.p->arr[1].x.n, c);
4148 return EP;
4149 }
4150
4151 evalue *affine2evalue(Value *coeff, Value denom, int nvar)
4152 {
4153 int i;
4154 evalue *E = ALLOC(evalue);
4155 value_init(E->d);
4156 evalue_set(E, coeff[nvar], denom);
4157 for (i = 0; i < nvar; ++i) {
4158 evalue *t;
4159 if (value_zero_p(coeff[i]))
4160 continue;
4161 t = term(i, coeff[i], denom);
4162 eadd(t, E);
4163 evalue_free(t);
4164 }
4165 return E;
4166 }
4167
4168 void evalue_substitute(evalue *e, evalue **subs)
4169 {
4170 int i, offset;
4171 evalue *v;
4172 enode *p;
4173
4174 if (value_notzero_p(e->d))
4175 return;
4176
4177 p = e->x.p;
4178 assert(p->type != partition);
4179
4180 for (i = 0; i < p->size; ++i)
4181 evalue_substitute(&p->arr[i], subs);
4182
4183 if (p->type == relation) {
4184 /* For relation a ? b : c, compute (a' ? 1) * b' + (a' ? 0 : 1) * c' */
4185 if (p->size == 3) {
4186 v = ALLOC(evalue);
4187 value_init(v->d);
4188 value_set_si(v->d, 0);
4189 v->x.p = new_enode(relation, 3, 0);
4190 evalue_copy(&v->x.p->arr[0], &p->arr[0]);
4191 evalue_set_si(&v->x.p->arr[1], 0, 1);
4192 evalue_set_si(&v->x.p->arr[2], 1, 1);
4193 emul(v, &p->arr[2]);
4194 evalue_free(v);
4195 }
4196 v = ALLOC(evalue);
4197 value_init(v->d);
4198 value_set_si(v->d, 0);
4199 v->x.p = new_enode(relation, 2, 0);
4200 value_clear(v->x.p->arr[0].d);
4201 v->x.p->arr[0] = p->arr[0];
4202 evalue_set_si(&v->x.p->arr[1], 1, 1);
4203 emul(v, &p->arr[1]);
4204 evalue_free(v);
4205 if (p->size == 3) {
4206 eadd(&p->arr[2], &p->arr[1]);
4207 free_evalue_refs(&p->arr[2]);
4208 }
4209 value_clear(e->d);
4210 *e = p->arr[1];
4211 free(p);
4212 return;
4213 }
4214
4215 if (p->type == polynomial)
4216 v = subs[p->pos-1];
4217 else {
4218 v = ALLOC(evalue);
4219 value_init(v->d);
4220 value_set_si(v->d, 0);
4221 v->x.p = new_enode(p->type, 3, -1);
4222 value_clear(v->x.p->arr[0].d);
4223 v->x.p->arr[0] = p->arr[0];
4224 evalue_set_si(&v->x.p->arr[1], 0, 1);
4225 evalue_set_si(&v->x.p->arr[2], 1, 1);
4226 }
4227
4228 offset = type_offset(p);
4229
4230 for (i = p->size-1; i >= offset+1; i--) {
4231 emul(v, &p->arr[i]);
4232 eadd(&p->arr[i], &p->arr[i-1]);
4233 free_evalue_refs(&(p->arr[i]));
4234 }
4235
4236 if (p->type != polynomial)
4237 evalue_free(v);
4238
4239 value_clear(e->d);
4240 *e = p->arr[offset];
4241 free(p);
4242 }
4243
4244 /* evalue e is given in terms of "new" parameter; CP maps the new
4245 * parameters back to the old parameters.
4246 * Transforms e such that it refers back to the old parameters and
4247 * adds appropriate constraints to the domain.
4248 * In particular, if CP maps the new parameters onto an affine
4249 * subspace of the old parameters, then the corresponding equalities
4250 * are added to the domain.
4251 * Also, if any of the new parameters was a rational combination
4252 * of the old parameters $p' = (<a, p> + c)/m$, then modulo
4253 * constraints ${<a, p> + c)/m} = 0$ are added to ensure
4254 * the new evalue remains non-zero only for integer parameters
4255 * of the new parameters (which have been removed by the substitution).
4256 */
4257 void evalue_backsubstitute(evalue *e, Matrix *CP, unsigned MaxRays)
4258 {
4259 Matrix *eq;
4260 Matrix *inv;
4261 evalue **subs;
4262 enode *p;
4263 int i, j;
4264 unsigned nparam = CP->NbColumns-1;
4265 Polyhedron *CEq;
4266 Value gcd;
4267
4268 if (EVALUE_IS_ZERO(*e))
4269 return;
4270
4271 assert(value_zero_p(e->d));
4272 p = e->x.p;
4273 assert(p->type == partition);
4274
4275 inv = left_inverse(CP, &eq);
4276 subs = ALLOCN(evalue *, nparam);
4277 for (i = 0; i < nparam; ++i)
4278 subs[i] = affine2evalue(inv->p[i], inv->p[nparam][inv->NbColumns-1],
4279 inv->NbColumns-1);
4280
4281 CEq = Constraints2Polyhedron(eq, MaxRays);
4282 addeliminatedparams_partition(p, inv, CEq, inv->NbColumns-1, MaxRays);
4283 Polyhedron_Free(CEq);
4284
4285 for (i = 0; i < p->size/2; ++i)
4286 evalue_substitute(&p->arr[2*i+1], subs);
4287
4288 for (i = 0; i < nparam; ++i)
4289 evalue_free(subs[i]);
4290 free(subs);
4291
4292 value_init(gcd);
4293 for (i = 0; i < inv->NbRows-1; ++i) {
4294 Vector_Gcd(inv->p[i], inv->NbColumns, &gcd);
4295 value_gcd(gcd, gcd, inv->p[inv->NbRows-1][inv->NbColumns-1]);
4296 if (value_eq(gcd, inv->p[inv->NbRows-1][inv->NbColumns-1]))
4297 continue;
4298 Vector_AntiScale(inv->p[i], inv->p[i], gcd, inv->NbColumns);
4299 value_divexact(gcd, inv->p[inv->NbRows-1][inv->NbColumns-1], gcd);
4300
4301 for (j = 0; j < p->size/2; ++j) {
4302 evalue *arg = affine2evalue(inv->p[i], gcd, inv->NbColumns-1);
4303 evalue *ev;
4304 evalue rel;
4305
4306 value_init(rel.d);
4307 value_set_si(rel.d, 0);
4308 rel.x.p = new_enode(relation, 2, 0);
4309 value_clear(rel.x.p->arr[1].d);
4310 rel.x.p->arr[1] = p->arr[2*j+1];
4311 ev = &rel.x.p->arr[0];
4312 value_set_si(ev->d, 0);
4313 ev->x.p = new_enode(fractional, 3, -1);
4314 evalue_set_si(&ev->x.p->arr[1], 0, 1);
4315 evalue_set_si(&ev->x.p->arr[2], 1, 1);
4316 value_clear(ev->x.p->arr[0].d);
4317 ev->x.p->arr[0] = *arg;
4318 free(arg);
4319
4320 p->arr[2*j+1] = rel;
4321 }
4322 }
4323 value_clear(gcd);
4324
4325 Matrix_Free(eq);
4326 Matrix_Free(inv);
4327 }
4328
4329 /* Computes
4330 *
4331 * \sum_{i=0}^n c_i/d X^i
4332 *
4333 * where d is the last element in the vector c.
4334 */
4335 evalue *evalue_polynomial(Vector *c, const evalue* X)
4336 {
4337 unsigned dim = c->Size-2;
4338 evalue EC;
4339 evalue *EP = ALLOC(evalue);
4340 int i;
4341
4342 value_init(EP->d);
4343
4344 if (EVALUE_IS_ZERO(*X) || dim == 0) {
4345 evalue_set(EP, c->p[0], c->p[dim+1]);
4346 reduce_constant(EP);
4347 return EP;
4348 }
4349
4350 evalue_set(EP, c->p[dim], c->p[dim+1]);
4351
4352 value_init(EC.d);
4353 evalue_set(&EC, c->p[dim], c->p[dim+1]);
4354
4355 for (i = dim-1; i >= 0; --i) {
4356 emul(X, EP);
4357 value_assign(EC.x.n, c->p[i]);
4358 eadd(&EC, EP);
4359 }
4360 free_evalue_refs(&EC);
4361 return EP;
4362 }
4363
4364 /* Create an evalue from an array of pairs of domains and evalues. */
4365 evalue *evalue_from_section_array(struct evalue_section *s, int n)
4366 {
4367 int i;
4368 evalue *res;
4369
4370 res = ALLOC(evalue);
4371 value_init(res->d);
4372
4373 if (n == 0)
4374 evalue_set_si(res, 0, 1);
4375 else {
4376 value_set_si(res->d, 0);
4377 res->x.p = new_enode(partition, 2*n, s[0].D->Dimension);
4378 for (i = 0; i < n; ++i) {
4379 EVALUE_SET_DOMAIN(res->x.p->arr[2*i], s[i].D);
4380 value_clear(res->x.p->arr[2*i+1].d);
4381 res->x.p->arr[2*i+1] = *s[i].E;
4382 free(s[i].E);
4383 }
4384 }
4385 return res;
4386 }
4387
4388 /* shift variables (>= first, 0-based) in polynomial n up (may be negative) */
4389 void evalue_shift_variables(evalue *e, int first, int n)
4390 {
4391 int i;
4392 if (value_notzero_p(e->d))
4393 return;
4394 assert(e->x.p->type == polynomial ||
4395 e->x.p->type == flooring ||
4396 e->x.p->type == fractional);
4397 if (e->x.p->type == polynomial && e->x.p->pos >= first+1) {
4398 assert(e->x.p->pos + n >= 1);
4399 e->x.p->pos += n;
4400 }
4401 for (i = 0; i < e->x.p->size; ++i)
4402 evalue_shift_variables(&e->x.p->arr[i], first, n);
4403 }
4404
4405 static const evalue *outer_floor(evalue *e, const evalue *outer)
4406 {
4407 int i;
4408
4409 if (value_notzero_p(e->d))
4410 return outer;
4411 switch (e->x.p->type) {
4412 case flooring:
4413 if (!outer || evalue_level_cmp(outer, &e->x.p->arr[0]) > 0)
4414 return &e->x.p->arr[0];
4415 else
4416 return outer;
4417 case polynomial:
4418 case fractional:
4419 case relation:
4420 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
4421 outer = outer_floor(&e->x.p->arr[i], outer);
4422 return outer;
4423 case partition:
4424 for (i = 0; i < e->x.p->size/2; ++i)
4425 outer = outer_floor(&e->x.p->arr[2*i+1], outer);
4426 return outer;
4427 default:
4428 assert(0);
4429 }
4430 }
4431
4432 /* Find and return outermost floor argument or NULL if e has no floors */
4433 evalue *evalue_outer_floor(evalue *e)
4434 {
4435 const evalue *floor = outer_floor(e, NULL);
4436 return floor ? evalue_dup(floor): NULL;
4437 }
4438
4439 static void evalue_set_to_zero(evalue *e)
4440 {
4441 if (EVALUE_IS_ZERO(*e))
4442 return;
4443 if (value_zero_p(e->d)) {
4444 free_evalue_refs(e);
4445 value_init(e->d);
4446 value_init(e->x.n);
4447 }
4448 value_set_si(e->d, 1);
4449 value_set_si(e->x.n, 0);
4450 }
4451
4452 /* Replace (outer) floor with argument "floor" by variable "var" (0-based)
4453 * and drop terms not containing the floor.
4454 * Returns true if e contains the floor.
4455 */
4456 int evalue_replace_floor(evalue *e, const evalue *floor, int var)
4457 {
4458 int i;
4459 int contains = 0;
4460 int reorder = 0;
4461
4462 if (value_notzero_p(e->d))
4463 return 0;
4464 switch (e->x.p->type) {
4465 case flooring:
4466 if (!eequal(floor, &e->x.p->arr[0]))
4467 return 0;
4468 e->x.p->type = polynomial;
4469 e->x.p->pos = 1 + var;
4470 e->x.p->size--;
4471 free_evalue_refs(&e->x.p->arr[0]);
4472 for (i = 0; i < e->x.p->size; ++i)
4473 e->x.p->arr[i] = e->x.p->arr[i+1];
4474 evalue_set_to_zero(&e->x.p->arr[0]);
4475 return 1;
4476 case polynomial:
4477 case fractional:
4478 case relation:
4479 for (i = type_offset(e->x.p); i < e->x.p->size; ++i) {
4480 int c = evalue_replace_floor(&e->x.p->arr[i], floor, var);
4481 contains |= c;
4482 if (!c)
4483 evalue_set_to_zero(&e->x.p->arr[i]);
4484 if (c && !reorder && evalue_level_cmp(&e->x.p->arr[i], e) < 0)
4485 reorder = 1;
4486 }
4487 evalue_reduce_size(e);
4488 if (reorder)
4489 reorder_terms(e);
4490 return contains;
4491 case partition:
4492 default:
4493 assert(0);
4494 }
4495 }
4496
4497 /* Replace (outer) floor with argument "floor" by variable zero */
4498 void evalue_drop_floor(evalue *e, const evalue *floor)
4499 {
4500 int i;
4501 enode *p;
4502
4503 if (value_notzero_p(e->d))
4504 return;
4505 switch (e->x.p->type) {
4506 case flooring:
4507 if (!eequal(floor, &e->x.p->arr[0]))
4508 return;
4509 p = e->x.p;
4510 free_evalue_refs(&p->arr[0]);
4511 for (i = 2; i < p->size; ++i)
4512 free_evalue_refs(&p->arr[i]);
4513 value_clear(e->d);
4514 *e = p->arr[1];
4515 free(p);
4516 break;
4517 case polynomial:
4518 case fractional:
4519 case relation:
4520 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
4521 evalue_drop_floor(&e->x.p->arr[i], floor);
4522 evalue_reduce_size(e);
4523 break;
4524 case partition:
4525 default:
4526 assert(0);
4527 }
4528 }
lattice point enumeration library