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