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