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