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