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