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