?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
force installation of (possibly) new version
[barvinok.git] / ev_operations.c
blob6f27d7be744a3d3790bc6b22bbbd5f4b4604eebe
1 #include <assert.h>
2 #include <math.h>
3 #include <stdlib.h>
4
5 #include "ev_operations.h"
6 #include "barvinok.h"
7 #include "util.h"
8
9 #ifndef value_pmodulus
10 #define value_pmodulus(ref,val1,val2) (mpz_fdiv_r((ref),(val1),(val2)))
11 #endif
12
13 #define ALLOC(p) p = (typeof(p))malloc(sizeof(*p))
14 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
15
16 void evalue_set_si(evalue *ev, int n, int d) {
17 value_set_si(ev->d, d);
18 value_init(ev->x.n);
19 value_set_si(ev->x.n, n);
20 }
21
22 void evalue_set(evalue *ev, Value n, Value d) {
23 value_assign(ev->d, d);
24 value_init(ev->x.n);
25 value_assign(ev->x.n, n);
26 }
27
28 void aep_evalue(evalue *e, int *ref) {
29
30 enode *p;
31 int i;
32
33 if (value_notzero_p(e->d))
34 return; /* a rational number, its already reduced */
35 if(!(p = e->x.p))
36 return; /* hum... an overflow probably occured */
37
38 /* First check the components of p */
39 for (i=0;i<p->size;i++)
40 aep_evalue(&p->arr[i],ref);
41
42 /* Then p itself */
43 switch (p->type) {
44 case polynomial:
45 case periodic:
46 case evector:
47 p->pos = ref[p->pos-1]+1;
48 }
49 return;
50 } /* aep_evalue */
51
52 /** Comments */
53 void addeliminatedparams_evalue(evalue *e,Matrix *CT) {
54
55 enode *p;
56 int i, j;
57 int *ref;
58
59 if (value_notzero_p(e->d))
60 return; /* a rational number, its already reduced */
61 if(!(p = e->x.p))
62 return; /* hum... an overflow probably occured */
63
64 /* Compute ref */
65 ref = (int *)malloc(sizeof(int)*(CT->NbRows-1));
66 for(i=0;i<CT->NbRows-1;i++)
67 for(j=0;j<CT->NbColumns;j++)
68 if(value_notzero_p(CT->p[i][j])) {
69 ref[i] = j;
70 break;
71 }
72
73 /* Transform the references in e, using ref */
74 aep_evalue(e,ref);
75 free( ref );
76 return;
77 } /* addeliminatedparams_evalue */
78
79 void addeliminatedparams_enum(evalue *e, Matrix *CT, Polyhedron *CEq,
80 unsigned MaxRays, unsigned nparam)
81 {
82 enode *p;
83 int i;
84
85 if (CT->NbRows == CT->NbColumns)
86 return;
87
88 if (EVALUE_IS_ZERO(*e))
89 return;
90
91 if (value_notzero_p(e->d)) {
92 evalue res;
93 value_init(res.d);
94 value_set_si(res.d, 0);
95 res.x.p = new_enode(partition, 2, nparam);
96 EVALUE_SET_DOMAIN(res.x.p->arr[0],
97 DomainConstraintSimplify(Polyhedron_Copy(CEq), MaxRays));
98 value_clear(res.x.p->arr[1].d);
99 res.x.p->arr[1] = *e;
100 *e = res;
101 return;
102 }
103
104 p = e->x.p;
105 assert(p);
106 assert(p->type == partition);
107 p->pos = nparam;
108
109 for (i=0; i<p->size/2; i++) {
110 Polyhedron *D = EVALUE_DOMAIN(p->arr[2*i]);
111 Polyhedron *T = DomainPreimage(D, CT, MaxRays);
112 Domain_Free(D);
113 D = T;
114 T = DomainIntersection(D, CEq, MaxRays);
115 Domain_Free(D);
116 EVALUE_SET_DOMAIN(p->arr[2*i], T);
117 addeliminatedparams_evalue(&p->arr[2*i+1], CT);
118 }
119 }
120
121 static int mod_rational_smaller(evalue *e1, evalue *e2)
122 {
123 int r;
124 Value m;
125 Value m2;
126 value_init(m);
127 value_init(m2);
128
129 assert(value_notzero_p(e1->d));
130 assert(value_notzero_p(e2->d));
131 value_multiply(m, e1->x.n, e2->d);
132 value_multiply(m2, e2->x.n, e1->d);
133 if (value_lt(m, m2))
134 r = 1;
135 else if (value_gt(m, m2))
136 r = 0;
137 else
138 r = -1;
139 value_clear(m);
140 value_clear(m2);
141
142 return r;
143 }
144
145 static int mod_term_smaller_r(evalue *e1, evalue *e2)
146 {
147 if (value_notzero_p(e1->d)) {
148 int r;
149 if (value_zero_p(e2->d))
150 return 1;
151 r = mod_rational_smaller(e1, e2);
152 return r == -1 ? 0 : r;
153 }
154 if (value_notzero_p(e2->d))
155 return 0;
156 if (e1->x.p->pos < e2->x.p->pos)
157 return 1;
158 else if (e1->x.p->pos > e2->x.p->pos)
159 return 0;
160 else {
161 int r = mod_rational_smaller(&e1->x.p->arr[1], &e2->x.p->arr[1]);
162 return r == -1
163 ? mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0])
164 : r;
165 }
166 }
167
168 static int mod_term_smaller(evalue *e1, evalue *e2)
169 {
170 assert(value_zero_p(e1->d));
171 assert(value_zero_p(e2->d));
172 assert(e1->x.p->type == fractional || e1->x.p->type == flooring);
173 assert(e2->x.p->type == fractional || e2->x.p->type == flooring);
174 return mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0]);
175 }
176
177 /* Negative pos means inequality */
178 /* s is negative of substitution if m is not zero */
179 struct fixed_param {
180 int pos;
181 evalue s;
182 Value d;
183 Value m;
184 };
185
186 struct subst {
187 struct fixed_param *fixed;
188 int n;
189 int max;
190 };
191
192 static int relations_depth(evalue *e)
193 {
194 int d;
195
196 for (d = 0;
197 value_zero_p(e->d) && e->x.p->type == relation;
198 e = &e->x.p->arr[1], ++d);
199 return d;
200 }
201
202 static void poly_denom_not_constant(evalue **pp, Value *d)
203 {
204 evalue *p = *pp;
205 value_set_si(*d, 1);
206
207 while (value_zero_p(p->d)) {
208 assert(p->x.p->type == polynomial);
209 assert(p->x.p->size == 2);
210 assert(value_notzero_p(p->x.p->arr[1].d));
211 value_lcm(*d, p->x.p->arr[1].d, d);
212 p = &p->x.p->arr[0];
213 }
214 *pp = p;
215 }
216
217 static void poly_denom(evalue *p, Value *d)
218 {
219 poly_denom_not_constant(&p, d);
220 value_lcm(*d, p->d, d);
221 }
222
223 #define EVALUE_IS_ONE(ev) (value_one_p((ev).d) && value_one_p((ev).x.n))
224
225 static void realloc_substitution(struct subst *s, int d)
226 {
227 struct fixed_param *n;
228 int i;
229 NALLOC(n, s->max+d);
230 for (i = 0; i < s->n; ++i)
231 n[i] = s->fixed[i];
232 free(s->fixed);
233 s->fixed = n;
234 s->max += d;
235 }
236
237 static int add_modulo_substitution(struct subst *s, evalue *r)
238 {
239 evalue *p;
240 evalue *f;
241 evalue *m;
242
243 assert(value_zero_p(r->d) && r->x.p->type == relation);
244 m = &r->x.p->arr[0];
245
246 /* May have been reduced already */
247 if (value_notzero_p(m->d))
248 return 0;
249
250 assert(value_zero_p(m->d) && m->x.p->type == fractional);
251 assert(m->x.p->size == 3);
252
253 /* fractional was inverted during reduction
254 * invert it back and move constant in
255 */
256 if (!EVALUE_IS_ONE(m->x.p->arr[2])) {
257 assert(value_pos_p(m->x.p->arr[2].d));
258 assert(value_mone_p(m->x.p->arr[2].x.n));
259 value_set_si(m->x.p->arr[2].x.n, 1);
260 value_increment(m->x.p->arr[1].x.n, m->x.p->arr[1].x.n);
261 assert(value_eq(m->x.p->arr[1].x.n, m->x.p->arr[1].d));
262 value_set_si(m->x.p->arr[1].x.n, 1);
263 eadd(&m->x.p->arr[1], &m->x.p->arr[0]);
264 value_set_si(m->x.p->arr[1].x.n, 0);
265 value_set_si(m->x.p->arr[1].d, 1);
266 }
267
268 /* Oops. Nested identical relations. */
269 if (!EVALUE_IS_ZERO(m->x.p->arr[1]))
270 return 0;
271
272 if (s->n >= s->max) {
273 int d = relations_depth(r);
274 realloc_substitution(s, d);
275 }
276
277 p = &m->x.p->arr[0];
278 assert(value_zero_p(p->d) && p->x.p->type == polynomial);
279 assert(p->x.p->size == 2);
280 f = &p->x.p->arr[1];
281
282 assert(value_pos_p(f->x.n));
283
284 value_init(s->fixed[s->n].m);
285 value_assign(s->fixed[s->n].m, f->d);
286 s->fixed[s->n].pos = p->x.p->pos;
287 value_init(s->fixed[s->n].d);
288 value_assign(s->fixed[s->n].d, f->x.n);
289 value_init(s->fixed[s->n].s.d);
290 evalue_copy(&s->fixed[s->n].s, &p->x.p->arr[0]);
291 ++s->n;
292
293 return 1;
294 }
295
296 static int type_offset(enode *p)
297 {
298 return p->type == fractional ? 1 :
299 p->type == flooring ? 1 : 0;
300 }
301
302 static void reorder_terms(enode *p, evalue *v)
303 {
304 int i;
305 int offset = type_offset(p);
306
307 for (i = p->size-1; i >= offset+1; i--) {
308 emul(v, &p->arr[i]);
309 eadd(&p->arr[i], &p->arr[i-1]);
310 free_evalue_refs(&(p->arr[i]));
311 }
312 p->size = offset+1;
313 free_evalue_refs(v);
314 }
315
316 void _reduce_evalue (evalue *e, struct subst *s, int fract) {
317
318 enode *p;
319 int i, j, k;
320 int add;
321
322 if (value_notzero_p(e->d)) {
323 if (fract)
324 mpz_fdiv_r(e->x.n, e->x.n, e->d);
325 return; /* a rational number, its already reduced */
326 }
327
328 if(!(p = e->x.p))
329 return; /* hum... an overflow probably occured */
330
331 /* First reduce the components of p */
332 add = p->type == relation;
333 for (i=0; i<p->size; i++) {
334 if (add && i == 1)
335 add = add_modulo_substitution(s, e);
336
337 if (i == 0 && p->type==fractional)
338 _reduce_evalue(&p->arr[i], s, 1);
339 else
340 _reduce_evalue(&p->arr[i], s, fract);
341
342 if (add && i == p->size-1) {
343 --s->n;
344 value_clear(s->fixed[s->n].m);
345 value_clear(s->fixed[s->n].d);
346 free_evalue_refs(&s->fixed[s->n].s);
347 } else if (add && i == 1)
348 s->fixed[s->n-1].pos *= -1;
349 }
350
351 if (p->type==periodic) {
352
353 /* Try to reduce the period */
354 for (i=1; i<=(p->size)/2; i++) {
355 if ((p->size % i)==0) {
356
357 /* Can we reduce the size to i ? */
358 for (j=0; j<i; j++)
359 for (k=j+i; k<e->x.p->size; k+=i)
360 if (!eequal(&p->arr[j], &p->arr[k])) goto you_lose;
361
362 /* OK, lets do it */
363 for (j=i; j<p->size; j++) free_evalue_refs(&p->arr[j]);
364 p->size = i;
365 break;
366
367 you_lose: /* OK, lets not do it */
368 continue;
369 }
370 }
371
372 /* Try to reduce its strength */
373 if (p->size == 1) {
374 value_clear(e->d);
375 memcpy(e,&p->arr[0],sizeof(evalue));
376 free(p);
377 }
378 }
379 else if (p->type==polynomial) {
380 for (k = 0; s && k < s->n; ++k) {
381 if (s->fixed[k].pos == p->pos) {
382 int divide = value_notone_p(s->fixed[k].d);
383 evalue d;
384
385 if (value_notzero_p(s->fixed[k].m)) {
386 if (!fract)
387 continue;
388 assert(p->size == 2);
389 if (divide && value_ne(s->fixed[k].d, p->arr[1].x.n))
390 continue;
391 if (!mpz_divisible_p(s->fixed[k].m, p->arr[1].d))
392 continue;
393 divide = 1;
394 }
395
396 if (divide) {
397 value_init(d.d);
398 value_assign(d.d, s->fixed[k].d);
399 value_init(d.x.n);
400 if (value_notzero_p(s->fixed[k].m))
401 value_oppose(d.x.n, s->fixed[k].m);
402 else
403 value_set_si(d.x.n, 1);
404 }
405
406 for (i=p->size-1;i>=1;i--) {
407 emul(&s->fixed[k].s, &p->arr[i]);
408 if (divide)
409 emul(&d, &p->arr[i]);
410 eadd(&p->arr[i], &p->arr[i-1]);
411 free_evalue_refs(&(p->arr[i]));
412 }
413 p->size = 1;
414 _reduce_evalue(&p->arr[0], s, fract);
415
416 if (divide)
417 free_evalue_refs(&d);
418
419 break;
420 }
421 }
422
423 /* Try to reduce the degree */
424 for (i=p->size-1;i>=1;i--) {
425 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
426 break;
427 /* Zero coefficient */
428 free_evalue_refs(&(p->arr[i]));
429 }
430 if (i+1<p->size)
431 p->size = i+1;
432
433 /* Try to reduce its strength */
434 if (p->size == 1) {
435 value_clear(e->d);
436 memcpy(e,&p->arr[0],sizeof(evalue));
437 free(p);
438 }
439 }
440 else if (p->type==fractional) {
441 int reorder = 0;
442 evalue v;
443
444 if (value_notzero_p(p->arr[0].d)) {
445 value_init(v.d);
446 value_assign(v.d, p->arr[0].d);
447 value_init(v.x.n);
448 mpz_fdiv_r(v.x.n, p->arr[0].x.n, p->arr[0].d);
449
450 reorder = 1;
451 } else {
452 evalue *f, *base;
453 evalue *pp = &p->arr[0];
454 assert(value_zero_p(pp->d) && pp->x.p->type == polynomial);
455 assert(pp->x.p->size == 2);
456
457 /* search for exact duplicate among the modulo inequalities */
458 do {
459 f = &pp->x.p->arr[1];
460 for (k = 0; s && k < s->n; ++k) {
461 if (-s->fixed[k].pos == pp->x.p->pos &&
462 value_eq(s->fixed[k].d, f->x.n) &&
463 value_eq(s->fixed[k].m, f->d) &&
464 eequal(&s->fixed[k].s, &pp->x.p->arr[0]))
465 break;
466 }
467 if (k < s->n) {
468 Value g;
469 value_init(g);
470
471 /* replace { E/m } by { (E-1)/m } + 1/m */
472 poly_denom(pp, &g);
473 if (reorder) {
474 evalue extra;
475 value_init(extra.d);
476 evalue_set_si(&extra, 1, 1);
477 value_assign(extra.d, g);
478 eadd(&extra, &v.x.p->arr[1]);
479 free_evalue_refs(&extra);
480
481 /* We've been going in circles; stop now */
482 if (value_ge(v.x.p->arr[1].x.n, v.x.p->arr[1].d)) {
483 free_evalue_refs(&v);
484 value_init(v.d);
485 evalue_set_si(&v, 0, 1);
486 break;
487 }
488 } else {
489 value_init(v.d);
490 value_set_si(v.d, 0);
491 v.x.p = new_enode(fractional, 3, -1);
492 evalue_set_si(&v.x.p->arr[1], 1, 1);
493 value_assign(v.x.p->arr[1].d, g);
494 evalue_set_si(&v.x.p->arr[2], 1, 1);
495 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
496 }
497
498 for (f = &v.x.p->arr[0]; value_zero_p(f->d);
499 f = &f->x.p->arr[0])
500 ;
501 value_division(f->d, g, f->d);
502 value_multiply(f->x.n, f->x.n, f->d);
503 value_assign(f->d, g);
504 value_decrement(f->x.n, f->x.n);
505 mpz_fdiv_r(f->x.n, f->x.n, f->d);
506
507 Gcd(f->d, f->x.n, &g);
508 value_division(f->d, f->d, g);
509 value_division(f->x.n, f->x.n, g);
510
511 value_clear(g);
512 pp = &v.x.p->arr[0];
513
514 reorder = 1;
515 }
516 } while (k < s->n);
517
518 /* reduction may have made this fractional arg smaller */
519 i = reorder ? p->size : 1;
520 for ( ; i < p->size; ++i)
521 if (value_zero_p(p->arr[i].d) &&
522 p->arr[i].x.p->type == fractional &&
523 !mod_term_smaller(e, &p->arr[i]))
524 break;
525 if (i < p->size) {
526 value_init(v.d);
527 value_set_si(v.d, 0);
528 v.x.p = new_enode(fractional, 3, -1);
529 evalue_set_si(&v.x.p->arr[1], 0, 1);
530 evalue_set_si(&v.x.p->arr[2], 1, 1);
531 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
532
533 reorder = 1;
534 }
535
536 if (!reorder) {
537 Value m;
538 Value r;
539 value_init(m);
540 value_init(r);
541 evalue *pp = &p->arr[0];
542 poly_denom_not_constant(&pp, &m);
543 mpz_fdiv_r(r, m, pp->d);
544 if (value_notzero_p(r)) {
545 value_init(v.d);
546 value_set_si(v.d, 0);
547 v.x.p = new_enode(fractional, 3, -1);
548
549 value_multiply(r, m, pp->x.n);
550 value_multiply(v.x.p->arr[1].d, m, pp->d);
551 value_init(v.x.p->arr[1].x.n);
552 mpz_fdiv_r(v.x.p->arr[1].x.n, r, pp->d);
553 mpz_fdiv_q(r, r, pp->d);
554
555 evalue_set_si(&v.x.p->arr[2], 1, 1);
556 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
557 pp = &v.x.p->arr[0];
558 while (value_zero_p(pp->d))
559 pp = &pp->x.p->arr[0];
560
561 value_assign(pp->d, m);
562 value_assign(pp->x.n, r);
563
564 Gcd(pp->d, pp->x.n, &r);
565 value_division(pp->d, pp->d, r);
566 value_division(pp->x.n, pp->x.n, r);
567
568 reorder = 1;
569 }
570 value_clear(m);
571 value_clear(r);
572 }
573
574 if (!reorder) {
575 int invert = 0;
576 Value twice;
577 value_init(twice);
578
579 for (pp = &p->arr[0]; value_zero_p(pp->d);
580 pp = &pp->x.p->arr[0]) {
581 f = &pp->x.p->arr[1];
582 assert(value_pos_p(f->d));
583 mpz_mul_ui(twice, f->x.n, 2);
584 if (value_lt(twice, f->d))
585 break;
586 if (value_eq(twice, f->d))
587 continue;
588 invert = 1;
589 break;
590 }
591
592 if (invert) {
593 value_init(v.d);
594 value_set_si(v.d, 0);
595 v.x.p = new_enode(fractional, 3, -1);
596 evalue_set_si(&v.x.p->arr[1], 0, 1);
597 poly_denom(&p->arr[0], &twice);
598 value_assign(v.x.p->arr[1].d, twice);
599 value_decrement(v.x.p->arr[1].x.n, twice);
600 evalue_set_si(&v.x.p->arr[2], -1, 1);
601 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
602
603 for (pp = &v.x.p->arr[0]; value_zero_p(pp->d);
604 pp = &pp->x.p->arr[0]) {
605 f = &pp->x.p->arr[1];
606 value_oppose(f->x.n, f->x.n);
607 mpz_fdiv_r(f->x.n, f->x.n, f->d);
608 }
609 value_division(pp->d, twice, pp->d);
610 value_multiply(pp->x.n, pp->x.n, pp->d);
611 value_assign(pp->d, twice);
612 value_oppose(pp->x.n, pp->x.n);
613 value_decrement(pp->x.n, pp->x.n);
614 mpz_fdiv_r(pp->x.n, pp->x.n, pp->d);
615
616 reorder = 1;
617 }
618
619 value_clear(twice);
620 }
621 }
622
623 if (reorder) {
624 reorder_terms(p, &v);
625 _reduce_evalue(&p->arr[1], s, fract);
626 }
627
628 /* Try to reduce the degree */
629 for (i=p->size-1;i>=2;i--) {
630 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
631 break;
632 /* Zero coefficient */
633 free_evalue_refs(&(p->arr[i]));
634 }
635 if (i+1<p->size)
636 p->size = i+1;
637
638 /* Try to reduce its strength */
639 if (p->size == 2) {
640 value_clear(e->d);
641 memcpy(e,&p->arr[1],sizeof(evalue));
642 free_evalue_refs(&(p->arr[0]));
643 free(p);
644 }
645 }
646 else if (p->type == flooring) {
647 /* Try to reduce the degree */
648 for (i=p->size-1;i>=2;i--) {
649 if (!EVALUE_IS_ZERO(p->arr[i]))
650 break;
651 /* Zero coefficient */
652 free_evalue_refs(&(p->arr[i]));
653 }
654 if (i+1<p->size)
655 p->size = i+1;
656
657 /* Try to reduce its strength */
658 if (p->size == 2) {
659 value_clear(e->d);
660 memcpy(e,&p->arr[1],sizeof(evalue));
661 free_evalue_refs(&(p->arr[0]));
662 free(p);
663 }
664 }
665 else if (p->type == relation) {
666 if (p->size == 3 && eequal(&p->arr[1], &p->arr[2])) {
667 free_evalue_refs(&(p->arr[2]));
668 free_evalue_refs(&(p->arr[0]));
669 p->size = 2;
670 value_clear(e->d);
671 *e = p->arr[1];
672 free(p);
673 return;
674 }
675 if (p->size == 3 && EVALUE_IS_ZERO(p->arr[2])) {
676 free_evalue_refs(&(p->arr[2]));
677 p->size = 2;
678 }
679 if (p->size == 2 && EVALUE_IS_ZERO(p->arr[1])) {
680 free_evalue_refs(&(p->arr[1]));
681 free_evalue_refs(&(p->arr[0]));
682 evalue_set_si(e, 0, 1);
683 free(p);
684 } else {
685 int reduced = 0;
686 evalue *m;
687 m = &p->arr[0];
688
689 /* Relation was reduced by means of an identical
690 * inequality => remove
691 */
692 if (value_zero_p(m->d) && !EVALUE_IS_ZERO(m->x.p->arr[1]))
693 reduced = 1;
694
695 if (reduced || value_notzero_p(p->arr[0].d)) {
696 if (!reduced && value_zero_p(p->arr[0].x.n)) {
697 value_clear(e->d);
698 memcpy(e,&p->arr[1],sizeof(evalue));
699 if (p->size == 3)
700 free_evalue_refs(&(p->arr[2]));
701 } else {
702 if (p->size == 3) {
703 value_clear(e->d);
704 memcpy(e,&p->arr[2],sizeof(evalue));
705 } else
706 evalue_set_si(e, 0, 1);
707 free_evalue_refs(&(p->arr[1]));
708 }
709 free_evalue_refs(&(p->arr[0]));
710 free(p);
711 }
712 }
713 }
714 return;
715 } /* reduce_evalue */
716
717 static void add_substitution(struct subst *s, Value *row, unsigned dim)
718 {
719 int k, l;
720 evalue *r;
721
722 for (k = 0; k < dim; ++k)
723 if (value_notzero_p(row[k+1]))
724 break;
725
726 Vector_Normalize_Positive(row+1, dim+1, k);
727 assert(s->n < s->max);
728 value_init(s->fixed[s->n].d);
729 value_init(s->fixed[s->n].m);
730 value_assign(s->fixed[s->n].d, row[k+1]);
731 s->fixed[s->n].pos = k+1;
732 value_set_si(s->fixed[s->n].m, 0);
733 r = &s->fixed[s->n].s;
734 value_init(r->d);
735 for (l = k+1; l < dim; ++l)
736 if (value_notzero_p(row[l+1])) {
737 value_set_si(r->d, 0);
738 r->x.p = new_enode(polynomial, 2, l + 1);
739 value_init(r->x.p->arr[1].x.n);
740 value_oppose(r->x.p->arr[1].x.n, row[l+1]);
741 value_set_si(r->x.p->arr[1].d, 1);
742 r = &r->x.p->arr[0];
743 }
744 value_init(r->x.n);
745 value_oppose(r->x.n, row[dim+1]);
746 value_set_si(r->d, 1);
747 ++s->n;
748 }
749
750 void reduce_evalue (evalue *e) {
751 struct subst s = { NULL, 0, 0 };
752
753 if (value_notzero_p(e->d))
754 return; /* a rational number, its already reduced */
755
756 if (e->x.p->type == partition) {
757 int i;
758 unsigned dim = -1;
759 for (i = 0; i < e->x.p->size/2; ++i) {
760 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
761 s.n = 0;
762 /* This shouldn't really happen;
763 * Empty domains should not be added.
764 */
765 if (emptyQ(D))
766 goto discard;
767
768 dim = D->Dimension;
769 if (D->next)
770 D = DomainConvex(D, 0);
771 if (!D->next && D->NbEq) {
772 int j, k;
773 if (s.max < dim) {
774 if (s.max != 0)
775 realloc_substitution(&s, dim);
776 else {
777 int d = relations_depth(&e->x.p->arr[2*i+1]);
778 s.max = dim+d;
779 NALLOC(s.fixed, s.max);
780 }
781 }
782 for (j = 0; j < D->NbEq; ++j)
783 add_substitution(&s, D->Constraint[j], dim);
784 }
785 if (D != EVALUE_DOMAIN(e->x.p->arr[2*i]))
786 Domain_Free(D);
787 _reduce_evalue(&e->x.p->arr[2*i+1], &s, 0);
788 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
789 discard:
790 free_evalue_refs(&e->x.p->arr[2*i+1]);
791 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
792 value_clear(e->x.p->arr[2*i].d);
793 e->x.p->size -= 2;
794 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
795 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
796 --i;
797 }
798 if (s.n != 0) {
799 int j;
800 for (j = 0; j < s.n; ++j) {
801 value_clear(s.fixed[j].d);
802 value_clear(s.fixed[j].m);
803 free_evalue_refs(&s.fixed[j].s);
804 }
805 }
806 }
807 if (e->x.p->size == 0) {
808 free(e->x.p);
809 evalue_set_si(e, 0, 1);
810 }
811 } else
812 _reduce_evalue(e, &s, 0);
813 if (s.max != 0)
814 free(s.fixed);
815 }
816
817 void print_evalue(FILE *DST,evalue *e,char **pname) {
818
819 if(value_notzero_p(e->d)) {
820 if(value_notone_p(e->d)) {
821 value_print(DST,VALUE_FMT,e->x.n);
822 fprintf(DST,"/");
823 value_print(DST,VALUE_FMT,e->d);
824 }
825 else {
826 value_print(DST,VALUE_FMT,e->x.n);
827 }
828 }
829 else
830 print_enode(DST,e->x.p,pname);
831 return;
832 } /* print_evalue */
833
834 void print_enode(FILE *DST,enode *p,char **pname) {
835
836 int i;
837
838 if (!p) {
839 fprintf(DST, "NULL");
840 return;
841 }
842 switch (p->type) {
843 case evector:
844 fprintf(DST, "{ ");
845 for (i=0; i<p->size; i++) {
846 print_evalue(DST, &p->arr[i], pname);
847 if (i!=(p->size-1))
848 fprintf(DST, ", ");
849 }
850 fprintf(DST, " }\n");
851 break;
852 case polynomial:
853 fprintf(DST, "( ");
854 for (i=p->size-1; i>=0; i--) {
855 print_evalue(DST, &p->arr[i], pname);
856 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
857 else if (i>1)
858 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
859 }
860 fprintf(DST, " )\n");
861 break;
862 case periodic:
863 fprintf(DST, "[ ");
864 for (i=0; i<p->size; i++) {
865 print_evalue(DST, &p->arr[i], pname);
866 if (i!=(p->size-1)) fprintf(DST, ", ");
867 }
868 fprintf(DST," ]_%s", pname[p->pos-1]);
869 break;
870 case flooring:
871 case fractional:
872 fprintf(DST, "( ");
873 for (i=p->size-1; i>=1; i--) {
874 print_evalue(DST, &p->arr[i], pname);
875 if (i >= 2) {
876 fprintf(DST, " * ");
877 fprintf(DST, p->type == flooring ? "[" : "{");
878 print_evalue(DST, &p->arr[0], pname);
879 fprintf(DST, p->type == flooring ? "]" : "}");
880 if (i>2)
881 fprintf(DST, "^%d + ", i-1);
882 else
883 fprintf(DST, " + ");
884 }
885 }
886 fprintf(DST, " )\n");
887 break;
888 case relation:
889 fprintf(DST, "[ ");
890 print_evalue(DST, &p->arr[0], pname);
891 fprintf(DST, "= 0 ] * \n");
892 print_evalue(DST, &p->arr[1], pname);
893 if (p->size > 2) {
894 fprintf(DST, " +\n [ ");
895 print_evalue(DST, &p->arr[0], pname);
896 fprintf(DST, "!= 0 ] * \n");
897 print_evalue(DST, &p->arr[2], pname);
898 }
899 break;
900 case partition: {
901 char **names = pname;
902 int maxdim = EVALUE_DOMAIN(p->arr[0])->Dimension;
903 if (p->pos < maxdim) {
904 NALLOC(names, maxdim);
905 for (i = 0; i < p->pos; ++i)
906 names[i] = pname[i];
907 for ( ; i < maxdim; ++i) {
908 NALLOC(names[i], 10);
909 snprintf(names[i], 10, "_p%d", i);
910 }
911 }
912
913 for (i=0; i<p->size/2; i++) {
914 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), names);
915 print_evalue(DST, &p->arr[2*i+1], names);
916 }
917
918 if (p->pos < maxdim) {
919 for (i = p->pos ; i < maxdim; ++i)
920 free(names[i]);
921 free(names);
922 }
923
924 break;
925 }
926 default:
927 assert(0);
928 }
929 return;
930 } /* print_enode */
931
932 static void eadd_rev(evalue *e1, evalue *res)
933 {
934 evalue ev;
935 value_init(ev.d);
936 evalue_copy(&ev, e1);
937 eadd(res, &ev);
938 free_evalue_refs(res);
939 *res = ev;
940 }
941
942 static void eadd_rev_cst (evalue *e1, evalue *res)
943 {
944 evalue ev;
945 value_init(ev.d);
946 evalue_copy(&ev, e1);
947 eadd(res, &ev.x.p->arr[type_offset(ev.x.p)]);
948 free_evalue_refs(res);
949 *res = ev;
950 }
951
952 struct section { Polyhedron * D; evalue E; };
953
954 void eadd_partitions (evalue *e1,evalue *res)
955 {
956 int n, i, j;
957 Polyhedron *d, *fd;
958 struct section *s;
959 s = (struct section *)
960 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
961 sizeof(struct section));
962 assert(s);
963 assert(e1->x.p->pos == res->x.p->pos);
964 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
965 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
966
967 n = 0;
968 for (j = 0; j < e1->x.p->size/2; ++j) {
969 assert(res->x.p->size >= 2);
970 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
971 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
972 if (!emptyQ(fd))
973 for (i = 1; i < res->x.p->size/2; ++i) {
974 Polyhedron *t = fd;
975 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
976 Domain_Free(t);
977 if (emptyQ(fd))
978 break;
979 }
980 if (emptyQ(fd)) {
981 Domain_Free(fd);
982 continue;
983 }
984 value_init(s[n].E.d);
985 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
986 s[n].D = fd;
987 ++n;
988 }
989 for (i = 0; i < res->x.p->size/2; ++i) {
990 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
991 for (j = 0; j < e1->x.p->size/2; ++j) {
992 Polyhedron *t;
993 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
994 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
995 if (emptyQ(d)) {
996 Domain_Free(d);
997 continue;
998 }
999 t = fd;
1000 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
1001 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
1002 Domain_Free(t);
1003 value_init(s[n].E.d);
1004 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1005 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
1006 s[n].D = d;
1007 ++n;
1008 }
1009 if (!emptyQ(fd)) {
1010 s[n].E = res->x.p->arr[2*i+1];
1011 s[n].D = fd;
1012 ++n;
1013 } else {
1014 free_evalue_refs(&res->x.p->arr[2*i+1]);
1015 Domain_Free(fd);
1016 }
1017 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
1018 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1019 value_clear(res->x.p->arr[2*i].d);
1020 }
1021
1022 free(res->x.p);
1023 assert(n > 0);
1024 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1025 for (j = 0; j < n; ++j) {
1026 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1027 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1028 value_clear(res->x.p->arr[2*j+1].d);
1029 res->x.p->arr[2*j+1] = s[j].E;
1030 }
1031
1032 free(s);
1033 }
1034
1035 static void explicit_complement(evalue *res)
1036 {
1037 enode *rel = new_enode(relation, 3, 0);
1038 assert(rel);
1039 value_clear(rel->arr[0].d);
1040 rel->arr[0] = res->x.p->arr[0];
1041 value_clear(rel->arr[1].d);
1042 rel->arr[1] = res->x.p->arr[1];
1043 value_set_si(rel->arr[2].d, 1);
1044 value_init(rel->arr[2].x.n);
1045 value_set_si(rel->arr[2].x.n, 0);
1046 free(res->x.p);
1047 res->x.p = rel;
1048 }
1049
1050 void eadd(evalue *e1,evalue *res) {
1051
1052 int i;
1053 if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
1054 /* Add two rational numbers */
1055 Value g,m1,m2;
1056 value_init(g);
1057 value_init(m1);
1058 value_init(m2);
1059
1060 value_multiply(m1,e1->x.n,res->d);
1061 value_multiply(m2,res->x.n,e1->d);
1062 value_addto(res->x.n,m1,m2);
1063 value_multiply(res->d,e1->d,res->d);
1064 Gcd(res->x.n,res->d,&g);
1065 if (value_notone_p(g)) {
1066 value_division(res->d,res->d,g);
1067 value_division(res->x.n,res->x.n,g);
1068 }
1069 value_clear(g); value_clear(m1); value_clear(m2);
1070 return ;
1071 }
1072 else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
1073 switch (res->x.p->type) {
1074 case polynomial:
1075 /* Add the constant to the constant term of a polynomial*/
1076 eadd(e1, &res->x.p->arr[0]);
1077 return ;
1078 case periodic:
1079 /* Add the constant to all elements of a periodic number */
1080 for (i=0; i<res->x.p->size; i++) {
1081 eadd(e1, &res->x.p->arr[i]);
1082 }
1083 return ;
1084 case evector:
1085 fprintf(stderr, "eadd: cannot add const with vector\n");
1086 return;
1087 case flooring:
1088 case fractional:
1089 eadd(e1, &res->x.p->arr[1]);
1090 return ;
1091 case partition:
1092 assert(EVALUE_IS_ZERO(*e1));
1093 break; /* Do nothing */
1094 case relation:
1095 /* Create (zero) complement if needed */
1096 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
1097 explicit_complement(res);
1098 for (i = 1; i < res->x.p->size; ++i)
1099 eadd(e1, &res->x.p->arr[i]);
1100 break;
1101 default:
1102 assert(0);
1103 }
1104 }
1105 /* add polynomial or periodic to constant
1106 * you have to exchange e1 and res, before doing addition */
1107
1108 else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
1109 eadd_rev(e1, res);
1110 return;
1111 }
1112 else { // ((e1->d==0) && (res->d==0))
1113 assert(!((e1->x.p->type == partition) ^
1114 (res->x.p->type == partition)));
1115 if (e1->x.p->type == partition) {
1116 eadd_partitions(e1, res);
1117 return;
1118 }
1119 if (e1->x.p->type == relation &&
1120 (res->x.p->type != relation ||
1121 mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1122 eadd_rev(e1, res);
1123 return;
1124 }
1125 if (res->x.p->type == relation) {
1126 if (e1->x.p->type == relation &&
1127 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1128 if (res->x.p->size < 3 && e1->x.p->size == 3)
1129 explicit_complement(res);
1130 if (e1->x.p->size < 3 && res->x.p->size == 3)
1131 explicit_complement(e1);
1132 for (i = 1; i < res->x.p->size; ++i)
1133 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1134 return;
1135 }
1136 if (res->x.p->size < 3)
1137 explicit_complement(res);
1138 for (i = 1; i < res->x.p->size; ++i)
1139 eadd(e1, &res->x.p->arr[i]);
1140 return;
1141 }
1142 if ((e1->x.p->type != res->x.p->type) ) {
1143 /* adding to evalues of different type. two cases are possible
1144 * res is periodic and e1 is polynomial, you have to exchange
1145 * e1 and res then to add e1 to the constant term of res */
1146 if (e1->x.p->type == polynomial) {
1147 eadd_rev_cst(e1, res);
1148 }
1149 else if (res->x.p->type == polynomial) {
1150 /* res is polynomial and e1 is periodic,
1151 add e1 to the constant term of res */
1152
1153 eadd(e1,&res->x.p->arr[0]);
1154 } else
1155 assert(0);
1156
1157 return;
1158 }
1159 else if (e1->x.p->pos != res->x.p->pos ||
1160 ((res->x.p->type == fractional ||
1161 res->x.p->type == flooring) &&
1162 !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1163 /* adding evalues of different position (i.e function of different unknowns
1164 * to case are possible */
1165
1166 switch (res->x.p->type) {
1167 case flooring:
1168 case fractional:
1169 if(mod_term_smaller(res, e1))
1170 eadd(e1,&res->x.p->arr[1]);
1171 else
1172 eadd_rev_cst(e1, res);
1173 return;
1174 case polynomial: // res and e1 are polynomials
1175 // add e1 to the constant term of res
1176
1177 if(res->x.p->pos < e1->x.p->pos)
1178 eadd(e1,&res->x.p->arr[0]);
1179 else
1180 eadd_rev_cst(e1, res);
1181 // value_clear(g); value_clear(m1); value_clear(m2);
1182 return;
1183 case periodic: // res and e1 are pointers to periodic numbers
1184 //add e1 to all elements of res
1185
1186 if(res->x.p->pos < e1->x.p->pos)
1187 for (i=0;i<res->x.p->size;i++) {
1188 eadd(e1,&res->x.p->arr[i]);
1189 }
1190 else
1191 eadd_rev(e1, res);
1192 return;
1193 default:
1194 assert(0);
1195 }
1196 }
1197
1198
1199 //same type , same pos and same size
1200 if (e1->x.p->size == res->x.p->size) {
1201 // add any element in e1 to the corresponding element in res
1202 i = type_offset(res->x.p);
1203 if (i == 1)
1204 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1205 for (; i<res->x.p->size; i++) {
1206 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1207 }
1208 return ;
1209 }
1210
1211 /* Sizes are different */
1212 switch(res->x.p->type) {
1213 case polynomial:
1214 case flooring:
1215 case fractional:
1216 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1217 /* new enode and add res to that new node. If you do not do */
1218 /* that, you lose the the upper weight part of e1 ! */
1219
1220 if(e1->x.p->size > res->x.p->size)
1221 eadd_rev(e1, res);
1222 else {
1223 i = type_offset(res->x.p);
1224 if (i == 1)
1225 assert(eequal(&e1->x.p->arr[0],
1226 &res->x.p->arr[0]));
1227 for (; i<e1->x.p->size ; i++) {
1228 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1229 }
1230
1231 return ;
1232 }
1233 break;
1234
1235 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1236 case periodic:
1237 {
1238 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1239 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1240 to add periodicaly elements of e1 to elements of ne, and finaly to
1241 return ne. */
1242 int x,y,p;
1243 Value ex, ey ,ep;
1244 evalue *ne;
1245 value_init(ex); value_init(ey);value_init(ep);
1246 x=e1->x.p->size;
1247 y= res->x.p->size;
1248 value_set_si(ex,e1->x.p->size);
1249 value_set_si(ey,res->x.p->size);
1250 value_assign (ep,*Lcm(ex,ey));
1251 p=(int)mpz_get_si(ep);
1252 ne= (evalue *) malloc (sizeof(evalue));
1253 value_init(ne->d);
1254 value_set_si( ne->d,0);
1255
1256 ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
1257 for(i=0;i<p;i++) {
1258 evalue_copy(&ne->x.p->arr[i], &res->x.p->arr[i%y]);
1259 }
1260 for(i=0;i<p;i++) {
1261 eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
1262 }
1263
1264 value_assign(res->d, ne->d);
1265 res->x.p=ne->x.p;
1266
1267 return ;
1268 }
1269 case evector:
1270 fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
1271 return ;
1272 default:
1273 assert(0);
1274 }
1275 }
1276 return ;
1277 }/* eadd */
1278
1279 static void emul_rev (evalue *e1, evalue *res)
1280 {
1281 evalue ev;
1282 value_init(ev.d);
1283 evalue_copy(&ev, e1);
1284 emul(res, &ev);
1285 free_evalue_refs(res);
1286 *res = ev;
1287 }
1288
1289 static void emul_poly (evalue *e1, evalue *res)
1290 {
1291 int i, j, o = type_offset(res->x.p);
1292 evalue tmp;
1293 int size=(e1->x.p->size + res->x.p->size - o - 1);
1294 value_init(tmp.d);
1295 value_set_si(tmp.d,0);
1296 tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
1297 if (o)
1298 evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
1299 for (i=o; i < e1->x.p->size; i++) {
1300 evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
1301 emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
1302 }
1303 for (; i<size; i++)
1304 evalue_set_si(&tmp.x.p->arr[i], 0, 1);
1305 for (i=o+1; i<res->x.p->size; i++)
1306 for (j=o; j<e1->x.p->size; j++) {
1307 evalue ev;
1308 value_init(ev.d);
1309 evalue_copy(&ev, &e1->x.p->arr[j]);
1310 emul(&res->x.p->arr[i], &ev);
1311 eadd(&ev, &tmp.x.p->arr[i+j-o]);
1312 free_evalue_refs(&ev);
1313 }
1314 free_evalue_refs(res);
1315 *res = tmp;
1316 }
1317
1318 void emul_partitions (evalue *e1,evalue *res)
1319 {
1320 int n, i, j, k;
1321 Polyhedron *d;
1322 struct section *s;
1323 s = (struct section *)
1324 malloc((e1->x.p->size/2) * (res->x.p->size/2) *
1325 sizeof(struct section));
1326 assert(s);
1327 assert(e1->x.p->pos == res->x.p->pos);
1328 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
1329 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1330
1331 n = 0;
1332 for (i = 0; i < res->x.p->size/2; ++i) {
1333 for (j = 0; j < e1->x.p->size/2; ++j) {
1334 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1335 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1336 if (emptyQ(d)) {
1337 Domain_Free(d);
1338 continue;
1339 }
1340
1341 /* This code is only needed because the partitions
1342 are not true partitions.
1343 */
1344 for (k = 0; k < n; ++k) {
1345 if (DomainIncludes(s[k].D, d))
1346 break;
1347 if (DomainIncludes(d, s[k].D)) {
1348 Domain_Free(s[k].D);
1349 free_evalue_refs(&s[k].E);
1350 if (n > k)
1351 s[k] = s[--n];
1352 --k;
1353 }
1354 }
1355 if (k < n) {
1356 Domain_Free(d);
1357 continue;
1358 }
1359
1360 value_init(s[n].E.d);
1361 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1362 emul(&e1->x.p->arr[2*j+1], &s[n].E);
1363 s[n].D = d;
1364 ++n;
1365 }
1366 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1367 value_clear(res->x.p->arr[2*i].d);
1368 free_evalue_refs(&res->x.p->arr[2*i+1]);
1369 }
1370
1371 free(res->x.p);
1372 if (n == 0)
1373 evalue_set_si(res, 0, 1);
1374 else {
1375 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1376 for (j = 0; j < n; ++j) {
1377 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1378 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1379 value_clear(res->x.p->arr[2*j+1].d);
1380 res->x.p->arr[2*j+1] = s[j].E;
1381 }
1382 }
1383
1384 free(s);
1385 }
1386
1387 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1388
1389 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1390 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1391 * evalues is not treated here */
1392
1393 void emul (evalue *e1, evalue *res ){
1394 int i,j;
1395
1396 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
1397 fprintf(stderr, "emul: do not proced on evector type !\n");
1398 return;
1399 }
1400
1401 if (EVALUE_IS_ZERO(*res))
1402 return;
1403
1404 if (value_zero_p(e1->d) && e1->x.p->type == partition) {
1405 if (value_zero_p(res->d) && res->x.p->type == partition)
1406 emul_partitions(e1, res);
1407 else
1408 emul_rev(e1, res);
1409 } else if (value_zero_p(res->d) && res->x.p->type == partition) {
1410 for (i = 0; i < res->x.p->size/2; ++i)
1411 emul(e1, &res->x.p->arr[2*i+1]);
1412 } else
1413 if (value_zero_p(res->d) && res->x.p->type == relation) {
1414 if (value_zero_p(e1->d) && e1->x.p->type == relation &&
1415 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1416 if (res->x.p->size < 3 && e1->x.p->size == 3)
1417 explicit_complement(res);
1418 if (e1->x.p->size < 3 && res->x.p->size == 3)
1419 explicit_complement(e1);
1420 for (i = 1; i < res->x.p->size; ++i)
1421 emul(&e1->x.p->arr[i], &res->x.p->arr[i]);
1422 return;
1423 }
1424 for (i = 1; i < res->x.p->size; ++i)
1425 emul(e1, &res->x.p->arr[i]);
1426 } else
1427 if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
1428 switch(e1->x.p->type) {
1429 case polynomial:
1430 switch(res->x.p->type) {
1431 case polynomial:
1432 if(e1->x.p->pos == res->x.p->pos) {
1433 /* Product of two polynomials of the same variable */
1434 emul_poly(e1, res);
1435 return;
1436 }
1437 else {
1438 /* Product of two polynomials of different variables */
1439
1440 if(res->x.p->pos < e1->x.p->pos)
1441 for( i=0; i<res->x.p->size ; i++)
1442 emul(e1, &res->x.p->arr[i]);
1443 else
1444 emul_rev(e1, res);
1445
1446 return ;
1447 }
1448 case periodic:
1449 case flooring:
1450 case fractional:
1451 /* Product of a polynomial and a periodic or fractional */
1452 emul_rev(e1, res);
1453 return;
1454 default:
1455 assert(0);
1456 }
1457 case periodic:
1458 switch(res->x.p->type) {
1459 case periodic:
1460 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
1461 /* Product of two periodics of the same parameter and period */
1462
1463 for(i=0; i<res->x.p->size;i++)
1464 emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
1465
1466 return;
1467 }
1468 else{
1469 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
1470 /* Product of two periodics of the same parameter and different periods */
1471 evalue *newp;
1472 Value x,y,z;
1473 int ix,iy,lcm;
1474 value_init(x); value_init(y);value_init(z);
1475 ix=e1->x.p->size;
1476 iy=res->x.p->size;
1477 value_set_si(x,e1->x.p->size);
1478 value_set_si(y,res->x.p->size);
1479 value_assign (z,*Lcm(x,y));
1480 lcm=(int)mpz_get_si(z);
1481 newp= (evalue *) malloc (sizeof(evalue));
1482 value_init(newp->d);
1483 value_set_si( newp->d,0);
1484 newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
1485 for(i=0;i<lcm;i++) {
1486 evalue_copy(&newp->x.p->arr[i],
1487 &res->x.p->arr[i%iy]);
1488 }
1489 for(i=0;i<lcm;i++)
1490 emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
1491
1492 value_assign(res->d,newp->d);
1493 res->x.p=newp->x.p;
1494
1495 value_clear(x); value_clear(y);value_clear(z);
1496 return ;
1497 }
1498 else {
1499 /* Product of two periodics of different parameters */
1500
1501 if(res->x.p->pos < e1->x.p->pos)
1502 for(i=0; i<res->x.p->size; i++)
1503 emul(e1, &(res->x.p->arr[i]));
1504 else
1505 emul_rev(e1, res);
1506
1507 return;
1508 }
1509 }
1510 case polynomial:
1511 /* Product of a periodic and a polynomial */
1512
1513 for(i=0; i<res->x.p->size ; i++)
1514 emul(e1, &(res->x.p->arr[i]));
1515
1516 return;
1517
1518 }
1519 case flooring:
1520 case fractional:
1521 switch(res->x.p->type) {
1522 case polynomial:
1523 for(i=0; i<res->x.p->size ; i++)
1524 emul(e1, &(res->x.p->arr[i]));
1525 return;
1526 default:
1527 case periodic:
1528 assert(0);
1529 case flooring:
1530 case fractional:
1531 assert(e1->x.p->type == res->x.p->type);
1532 if (e1->x.p->pos == res->x.p->pos &&
1533 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1534 evalue d;
1535 value_init(d.d);
1536 poly_denom(&e1->x.p->arr[0], &d.d);
1537 if (e1->x.p->type != fractional || !value_two_p(d.d))
1538 emul_poly(e1, res);
1539 else {
1540 evalue tmp;
1541 value_init(d.x.n);
1542 value_set_si(d.x.n, 1);
1543 /* { x }^2 == { x }/2 */
1544 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1545 assert(e1->x.p->size == 3);
1546 assert(res->x.p->size == 3);
1547 value_init(tmp.d);
1548 evalue_copy(&tmp, &res->x.p->arr[2]);
1549 emul(&d, &tmp);
1550 eadd(&res->x.p->arr[1], &tmp);
1551 emul(&e1->x.p->arr[2], &tmp);
1552 emul(&e1->x.p->arr[1], res);
1553 eadd(&tmp, &res->x.p->arr[2]);
1554 free_evalue_refs(&tmp);
1555 value_clear(d.x.n);
1556 }
1557 value_clear(d.d);
1558 } else {
1559 if(mod_term_smaller(res, e1))
1560 for(i=1; i<res->x.p->size ; i++)
1561 emul(e1, &(res->x.p->arr[i]));
1562 else
1563 emul_rev(e1, res);
1564 return;
1565 }
1566 }
1567 break;
1568 case relation:
1569 emul_rev(e1, res);
1570 break;
1571 default:
1572 assert(0);
1573 }
1574 }
1575 else {
1576 if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
1577 /* Product of two rational numbers */
1578
1579 Value g;
1580 value_init(g);
1581 value_multiply(res->d,e1->d,res->d);
1582 value_multiply(res->x.n,e1->x.n,res->x.n );
1583 Gcd(res->x.n, res->d,&g);
1584 if (value_notone_p(g)) {
1585 value_division(res->d,res->d,g);
1586 value_division(res->x.n,res->x.n,g);
1587 }
1588 value_clear(g);
1589 return ;
1590 }
1591 else {
1592 if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
1593 /* Product of an expression (polynomial or peririodic) and a rational number */
1594
1595 emul_rev(e1, res);
1596 return ;
1597 }
1598 else {
1599 /* Product of a rationel number and an expression (polynomial or peririodic) */
1600
1601 i = type_offset(res->x.p);
1602 for (; i<res->x.p->size; i++)
1603 emul(e1, &res->x.p->arr[i]);
1604
1605 return ;
1606 }
1607 }
1608 }
1609
1610 return ;
1611 }
1612
1613 /* Frees mask content ! */
1614 void emask(evalue *mask, evalue *res) {
1615 int n, i, j;
1616 Polyhedron *d, *fd;
1617 struct section *s;
1618 evalue mone;
1619 int pos;
1620
1621 if (EVALUE_IS_ZERO(*res)) {
1622 free_evalue_refs(mask);
1623 return;
1624 }
1625
1626 assert(value_zero_p(mask->d));
1627 assert(mask->x.p->type == partition);
1628 assert(value_zero_p(res->d));
1629 assert(res->x.p->type == partition);
1630 assert(mask->x.p->pos == res->x.p->pos);
1631 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1632 assert(mask->x.p->pos == EVALUE_DOMAIN(mask->x.p->arr[0])->Dimension);
1633 pos = res->x.p->pos;
1634
1635 s = (struct section *)
1636 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1637 sizeof(struct section));
1638 assert(s);
1639
1640 value_init(mone.d);
1641 evalue_set_si(&mone, -1, 1);
1642
1643 n = 0;
1644 for (j = 0; j < res->x.p->size/2; ++j) {
1645 assert(mask->x.p->size >= 2);
1646 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1647 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1648 if (!emptyQ(fd))
1649 for (i = 1; i < mask->x.p->size/2; ++i) {
1650 Polyhedron *t = fd;
1651 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1652 Domain_Free(t);
1653 if (emptyQ(fd))
1654 break;
1655 }
1656 if (emptyQ(fd)) {
1657 Domain_Free(fd);
1658 continue;
1659 }
1660 value_init(s[n].E.d);
1661 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1662 s[n].D = fd;
1663 ++n;
1664 }
1665 for (i = 0; i < mask->x.p->size/2; ++i) {
1666 if (EVALUE_IS_ONE(mask->x.p->arr[2*i+1]))
1667 continue;
1668
1669 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1670 eadd(&mone, &mask->x.p->arr[2*i+1]);
1671 emul(&mone, &mask->x.p->arr[2*i+1]);
1672 for (j = 0; j < res->x.p->size/2; ++j) {
1673 Polyhedron *t;
1674 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1675 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1676 if (emptyQ(d)) {
1677 Domain_Free(d);
1678 continue;
1679 }
1680 t = fd;
1681 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1682 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1683 Domain_Free(t);
1684 value_init(s[n].E.d);
1685 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1686 emul(&mask->x.p->arr[2*i+1], &s[n].E);
1687 s[n].D = d;
1688 ++n;
1689 }
1690
1691 if (!emptyQ(fd)) {
1692 /* Just ignore; this may have been previously masked off */
1693 }
1694 if (fd != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1695 Domain_Free(fd);
1696 }
1697
1698 free_evalue_refs(&mone);
1699 free_evalue_refs(mask);
1700 free_evalue_refs(res);
1701 value_init(res->d);
1702 if (n == 0)
1703 evalue_set_si(res, 0, 1);
1704 else {
1705 res->x.p = new_enode(partition, 2*n, pos);
1706 for (j = 0; j < n; ++j) {
1707 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1708 value_clear(res->x.p->arr[2*j+1].d);
1709 res->x.p->arr[2*j+1] = s[j].E;
1710 }
1711 }
1712
1713 free(s);
1714 }
1715
1716 void evalue_copy(evalue *dst, evalue *src)
1717 {
1718 value_assign(dst->d, src->d);
1719 if(value_notzero_p(src->d)) {
1720 value_init(dst->x.n);
1721 value_assign(dst->x.n, src->x.n);
1722 } else
1723 dst->x.p = ecopy(src->x.p);
1724 }
1725
1726 enode *new_enode(enode_type type,int size,int pos) {
1727
1728 enode *res;
1729 int i;
1730
1731 if(size == 0) {
1732 fprintf(stderr, "Allocating enode of size 0 !\n" );
1733 return NULL;
1734 }
1735 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1736 res->type = type;
1737 res->size = size;
1738 res->pos = pos;
1739 for(i=0; i<size; i++) {
1740 value_init(res->arr[i].d);
1741 value_set_si(res->arr[i].d,0);
1742 res->arr[i].x.p = 0;
1743 }
1744 return res;
1745 } /* new_enode */
1746
1747 enode *ecopy(enode *e) {
1748
1749 enode *res;
1750 int i;
1751
1752 res = new_enode(e->type,e->size,e->pos);
1753 for(i=0;i<e->size;++i) {
1754 value_assign(res->arr[i].d,e->arr[i].d);
1755 if(value_zero_p(res->arr[i].d))
1756 res->arr[i].x.p = ecopy(e->arr[i].x.p);
1757 else if (EVALUE_IS_DOMAIN(res->arr[i]))
1758 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1759 else {
1760 value_init(res->arr[i].x.n);
1761 value_assign(res->arr[i].x.n,e->arr[i].x.n);
1762 }
1763 }
1764 return(res);
1765 } /* ecopy */
1766
1767 int ecmp(const evalue *e1, const evalue *e2)
1768 {
1769 enode *p1, *p2;
1770 int i;
1771 int r;
1772
1773 if (value_notzero_p(e1->d) && value_notzero_p(e2->d)) {
1774 Value m, m2;
1775 value_init(m);
1776 value_init(m2);
1777 value_multiply(m, e1->x.n, e2->d);
1778 value_multiply(m2, e2->x.n, e1->d);
1779
1780 if (value_lt(m, m2))
1781 r = -1;
1782 else if (value_gt(m, m2))
1783 r = 1;
1784 else
1785 r = 0;
1786
1787 value_clear(m);
1788 value_clear(m2);
1789
1790 return r;
1791 }
1792 if (value_notzero_p(e1->d))
1793 return -1;
1794 if (value_notzero_p(e2->d))
1795 return 1;
1796
1797 p1 = e1->x.p;
1798 p2 = e2->x.p;
1799
1800 if (p1->type != p2->type)
1801 return p1->type - p2->type;
1802 if (p1->pos != p2->pos)
1803 return p1->pos - p2->pos;
1804 if (p1->size != p2->size)
1805 return p1->size - p2->size;
1806
1807 for (i = p1->size-1; i >= 0; --i)
1808 if ((r = ecmp(&p1->arr[i], &p2->arr[i])) != 0)
1809 return r;
1810
1811 return 0;
1812 }
1813
1814 int eequal(evalue *e1,evalue *e2) {
1815
1816 int i;
1817 enode *p1, *p2;
1818
1819 if (value_ne(e1->d,e2->d))
1820 return 0;
1821
1822 /* e1->d == e2->d */
1823 if (value_notzero_p(e1->d)) {
1824 if (value_ne(e1->x.n,e2->x.n))
1825 return 0;
1826
1827 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1828 return 1;
1829 }
1830
1831 /* e1->d == e2->d == 0 */
1832 p1 = e1->x.p;
1833 p2 = e2->x.p;
1834 if (p1->type != p2->type) return 0;
1835 if (p1->size != p2->size) return 0;
1836 if (p1->pos != p2->pos) return 0;
1837 for (i=0; i<p1->size; i++)
1838 if (!eequal(&p1->arr[i], &p2->arr[i]) )
1839 return 0;
1840 return 1;
1841 } /* eequal */
1842
1843 void free_evalue_refs(evalue *e) {
1844
1845 enode *p;
1846 int i;
1847
1848 if (EVALUE_IS_DOMAIN(*e)) {
1849 Domain_Free(EVALUE_DOMAIN(*e));
1850 value_clear(e->d);
1851 return;
1852 } else if (value_pos_p(e->d)) {
1853
1854 /* 'e' stores a constant */
1855 value_clear(e->d);
1856 value_clear(e->x.n);
1857 return;
1858 }
1859 assert(value_zero_p(e->d));
1860 value_clear(e->d);
1861 p = e->x.p;
1862 if (!p) return; /* null pointer */
1863 for (i=0; i<p->size; i++) {
1864 free_evalue_refs(&(p->arr[i]));
1865 }
1866 free(p);
1867 return;
1868 } /* free_evalue_refs */
1869
1870 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
1871 Vector * val, evalue *res)
1872 {
1873 unsigned nparam = periods->Size;
1874
1875 if (p == nparam) {
1876 double d = compute_evalue(e, val->p);
1877 d *= VALUE_TO_DOUBLE(m);
1878 if (d > 0)
1879 d += .25;
1880 else
1881 d -= .25;
1882 value_assign(res->d, m);
1883 value_init(res->x.n);
1884 value_set_double(res->x.n, d);
1885 mpz_fdiv_r(res->x.n, res->x.n, m);
1886 return;
1887 }
1888 if (value_one_p(periods->p[p]))
1889 mod2table_r(e, periods, m, p+1, val, res);
1890 else {
1891 Value tmp;
1892 value_init(tmp);
1893
1894 value_assign(tmp, periods->p[p]);
1895 value_set_si(res->d, 0);
1896 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
1897 do {
1898 value_decrement(tmp, tmp);
1899 value_assign(val->p[p], tmp);
1900 mod2table_r(e, periods, m, p+1, val,
1901 &res->x.p->arr[VALUE_TO_INT(tmp)]);
1902 } while (value_pos_p(tmp));
1903
1904 value_clear(tmp);
1905 }
1906 }
1907
1908 static void rel2table(evalue *e, int zero)
1909 {
1910 if (value_pos_p(e->d)) {
1911 if (value_zero_p(e->x.n) == zero)
1912 value_set_si(e->x.n, 1);
1913 else
1914 value_set_si(e->x.n, 0);
1915 value_set_si(e->d, 1);
1916 } else {
1917 int i;
1918 for (i = 0; i < e->x.p->size; ++i)
1919 rel2table(&e->x.p->arr[i], zero);
1920 }
1921 }
1922
1923 void evalue_mod2table(evalue *e, int nparam)
1924 {
1925 enode *p;
1926 int i;
1927
1928 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
1929 return;
1930 p = e->x.p;
1931 for (i=0; i<p->size; i++) {
1932 evalue_mod2table(&(p->arr[i]), nparam);
1933 }
1934 if (p->type == relation) {
1935 evalue copy;
1936
1937 if (p->size > 2) {
1938 value_init(copy.d);
1939 evalue_copy(&copy, &p->arr[0]);
1940 }
1941 rel2table(&p->arr[0], 1);
1942 emul(&p->arr[0], &p->arr[1]);
1943 if (p->size > 2) {
1944 rel2table(&copy, 0);
1945 emul(&copy, &p->arr[2]);
1946 eadd(&p->arr[2], &p->arr[1]);
1947 free_evalue_refs(&p->arr[2]);
1948 free_evalue_refs(&copy);
1949 }
1950 free_evalue_refs(&p->arr[0]);
1951 value_clear(e->d);
1952 *e = p->arr[1];
1953 free(p);
1954 } else if (p->type == fractional) {
1955 Vector *periods = Vector_Alloc(nparam);
1956 Vector *val = Vector_Alloc(nparam);
1957 Value tmp;
1958 evalue *ev;
1959 evalue EP, res;
1960
1961 value_init(tmp);
1962 value_set_si(tmp, 1);
1963 Vector_Set(periods->p, 1, nparam);
1964 Vector_Set(val->p, 0, nparam);
1965 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
1966 enode *p = ev->x.p;
1967
1968 assert(p->type == polynomial);
1969 assert(p->size == 2);
1970 value_assign(periods->p[p->pos-1], p->arr[1].d);
1971 value_lcm(tmp, p->arr[1].d, &tmp);
1972 }
1973 value_lcm(tmp, ev->d, &tmp);
1974 value_init(EP.d);
1975 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
1976
1977 value_init(res.d);
1978 evalue_set_si(&res, 0, 1);
1979 /* Compute the polynomial using Horner's rule */
1980 for (i=p->size-1;i>1;i--) {
1981 eadd(&p->arr[i], &res);
1982 emul(&EP, &res);
1983 }
1984 eadd(&p->arr[1], &res);
1985
1986 free_evalue_refs(e);
1987 free_evalue_refs(&EP);
1988 *e = res;
1989
1990 value_clear(tmp);
1991 Vector_Free(val);
1992 Vector_Free(periods);
1993 }
1994 } /* evalue_mod2table */
1995
1996 /********************************************************/
1997 /* function in domain */
1998 /* check if the parameters in list_args */
1999 /* verifies the constraints of Domain P */
2000 /********************************************************/
2001 int in_domain(Polyhedron *P, Value *list_args) {
2002
2003 int col,row;
2004 Value v; /* value of the constraint of a row when
2005 parameters are instanciated*/
2006 Value tmp;
2007
2008 value_init(v);
2009 value_init(tmp);
2010
2011 /*P->Constraint constraint matrice of polyhedron P*/
2012 for(row=0;row<P->NbConstraints;row++) {
2013 value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
2014 for(col=1;col<P->Dimension+1;col++) {
2015 value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
2016 value_addto(v,v,tmp);
2017 }
2018 if (value_notzero_p(P->Constraint[row][0])) {
2019
2020 /*if v is not >=0 then this constraint is not respected */
2021 if (value_neg_p(v)) {
2022 next:
2023 value_clear(v);
2024 value_clear(tmp);
2025 return P->next ? in_domain(P->next, list_args) : 0;
2026 }
2027 }
2028 else {
2029
2030 /*if v is not = 0 then this constraint is not respected */
2031 if (value_notzero_p(v))
2032 goto next;
2033 }
2034 }
2035
2036 /*if not return before this point => all
2037 the constraints are respected */
2038 value_clear(v);
2039 value_clear(tmp);
2040 return 1;
2041 } /* in_domain */
2042
2043 /****************************************************/
2044 /* function compute enode */
2045 /* compute the value of enode p with parameters */
2046 /* list "list_args */
2047 /* compute the polynomial or the periodic */
2048 /****************************************************/
2049
2050 static double compute_enode(enode *p, Value *list_args) {
2051
2052 int i;
2053 Value m, param;
2054 double res=0.0;
2055
2056 if (!p)
2057 return(0.);
2058
2059 value_init(m);
2060 value_init(param);
2061
2062 if (p->type == polynomial) {
2063 if (p->size > 1)
2064 value_assign(param,list_args[p->pos-1]);
2065
2066 /* Compute the polynomial using Horner's rule */
2067 for (i=p->size-1;i>0;i--) {
2068 res +=compute_evalue(&p->arr[i],list_args);
2069 res *=VALUE_TO_DOUBLE(param);
2070 }
2071 res +=compute_evalue(&p->arr[0],list_args);
2072 }
2073 else if (p->type == fractional) {
2074 double d = compute_evalue(&p->arr[0], list_args);
2075 d -= floor(d+1e-10);
2076
2077 /* Compute the polynomial using Horner's rule */
2078 for (i=p->size-1;i>1;i--) {
2079 res +=compute_evalue(&p->arr[i],list_args);
2080 res *=d;
2081 }
2082 res +=compute_evalue(&p->arr[1],list_args);
2083 }
2084 else if (p->type == flooring) {
2085 double d = compute_evalue(&p->arr[0], list_args);
2086 d = floor(d+1e-10);
2087
2088 /* Compute the polynomial using Horner's rule */
2089 for (i=p->size-1;i>1;i--) {
2090 res +=compute_evalue(&p->arr[i],list_args);
2091 res *=d;
2092 }
2093 res +=compute_evalue(&p->arr[1],list_args);
2094 }
2095 else if (p->type == periodic) {
2096 value_assign(m,list_args[p->pos-1]);
2097
2098 /* Choose the right element of the periodic */
2099 value_set_si(param,p->size);
2100 value_pmodulus(m,m,param);
2101 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
2102 }
2103 else if (p->type == relation) {
2104 if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
2105 res = compute_evalue(&p->arr[1], list_args);
2106 else if (p->size > 2)
2107 res = compute_evalue(&p->arr[2], list_args);
2108 }
2109 else if (p->type == partition) {
2110 int dim = EVALUE_DOMAIN(p->arr[0])->Dimension;
2111 Value *vals = list_args;
2112 if (p->pos < dim) {
2113 NALLOC(vals, dim);
2114 for (i = 0; i < dim; ++i) {
2115 value_init(vals[i]);
2116 if (i < p->pos)
2117 value_assign(vals[i], list_args[i]);
2118 }
2119 }
2120 for (i = 0; i < p->size/2; ++i)
2121 if (DomainContains(EVALUE_DOMAIN(p->arr[2*i]), vals, p->pos, 0, 1)) {
2122 res = compute_evalue(&p->arr[2*i+1], vals);
2123 break;
2124 }
2125 if (p->pos < dim) {
2126 for (i = 0; i < dim; ++i)
2127 value_clear(vals[i]);
2128 free(vals);
2129 }
2130 }
2131 else
2132 assert(0);
2133 value_clear(m);
2134 value_clear(param);
2135 return res;
2136 } /* compute_enode */
2137
2138 /*************************************************/
2139 /* return the value of Ehrhart Polynomial */
2140 /* It returns a double, because since it is */
2141 /* a recursive function, some intermediate value */
2142 /* might not be integral */
2143 /*************************************************/
2144
2145 double compute_evalue(evalue *e,Value *list_args) {
2146
2147 double res;
2148
2149 if (value_notzero_p(e->d)) {
2150 if (value_notone_p(e->d))
2151 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
2152 else
2153 res = VALUE_TO_DOUBLE(e->x.n);
2154 }
2155 else
2156 res = compute_enode(e->x.p,list_args);
2157 return res;
2158 } /* compute_evalue */
2159
2160
2161 /****************************************************/
2162 /* function compute_poly : */
2163 /* Check for the good validity domain */
2164 /* return the number of point in the Polyhedron */
2165 /* in allocated memory */
2166 /* Using the Ehrhart pseudo-polynomial */
2167 /****************************************************/
2168 Value *compute_poly(Enumeration *en,Value *list_args) {
2169
2170 Value *tmp;
2171 /* double d; int i; */
2172
2173 tmp = (Value *) malloc (sizeof(Value));
2174 assert(tmp != NULL);
2175 value_init(*tmp);
2176 value_set_si(*tmp,0);
2177
2178 if(!en)
2179 return(tmp); /* no ehrhart polynomial */
2180 if(en->ValidityDomain) {
2181 if(!en->ValidityDomain->Dimension) { /* no parameters */
2182 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2183 return(tmp);
2184 }
2185 }
2186 else
2187 return(tmp); /* no Validity Domain */
2188 while(en) {
2189 if(in_domain(en->ValidityDomain,list_args)) {
2190
2191 #ifdef EVAL_EHRHART_DEBUG
2192 Print_Domain(stdout,en->ValidityDomain);
2193 print_evalue(stdout,&en->EP);
2194 #endif
2195
2196 /* d = compute_evalue(&en->EP,list_args);
2197 i = d;
2198 printf("(double)%lf = %d\n", d, i ); */
2199 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2200 return(tmp);
2201 }
2202 else
2203 en=en->next;
2204 }
2205 value_set_si(*tmp,0);
2206 return(tmp); /* no compatible domain with the arguments */
2207 } /* compute_poly */
2208
2209 size_t value_size(Value v) {
2210 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2211 * sizeof(v[0]._mp_d[0]);
2212 }
2213
2214 size_t domain_size(Polyhedron *D)
2215 {
2216 int i, j;
2217 size_t s = sizeof(*D);
2218
2219 for (i = 0; i < D->NbConstraints; ++i)
2220 for (j = 0; j < D->Dimension+2; ++j)
2221 s += value_size(D->Constraint[i][j]);
2222
2223 /*
2224 for (i = 0; i < D->NbRays; ++i)
2225 for (j = 0; j < D->Dimension+2; ++j)
2226 s += value_size(D->Ray[i][j]);
2227 */
2228
2229 return D->next ? s+domain_size(D->next) : s;
2230 }
2231
2232 size_t enode_size(enode *p) {
2233 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2234 int i;
2235
2236 if (p->type == partition)
2237 for (i = 0; i < p->size/2; ++i) {
2238 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2239 s += evalue_size(&p->arr[2*i+1]);
2240 }
2241 else
2242 for (i = 0; i < p->size; ++i) {
2243 s += evalue_size(&p->arr[i]);
2244 }
2245 return s;
2246 }
2247
2248 size_t evalue_size(evalue *e)
2249 {
2250 size_t s = sizeof(*e);
2251 s += value_size(e->d);
2252 if (value_notzero_p(e->d))
2253 s += value_size(e->x.n);
2254 else
2255 s += enode_size(e->x.p);
2256 return s;
2257 }
2258
2259 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2260 {
2261 evalue *found = NULL;
2262 evalue offset;
2263 evalue copy;
2264 int i;
2265
2266 if (value_pos_p(e->d) || e->x.p->type != fractional)
2267 return NULL;
2268
2269 value_init(offset.d);
2270 value_init(offset.x.n);
2271 poly_denom(&e->x.p->arr[0], &offset.d);
2272 value_lcm(m, offset.d, &offset.d);
2273 value_set_si(offset.x.n, 1);
2274
2275 value_init(copy.d);
2276 evalue_copy(&copy, cst);
2277
2278 eadd(&offset, cst);
2279 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2280
2281 if (eequal(base, &e->x.p->arr[0]))
2282 found = &e->x.p->arr[0];
2283 else {
2284 value_set_si(offset.x.n, -2);
2285
2286 eadd(&offset, cst);
2287 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2288
2289 if (eequal(base, &e->x.p->arr[0]))
2290 found = base;
2291 }
2292 free_evalue_refs(cst);
2293 free_evalue_refs(&offset);
2294 *cst = copy;
2295
2296 for (i = 1; !found && i < e->x.p->size; ++i)
2297 found = find_second(base, cst, &e->x.p->arr[i], m);
2298
2299 return found;
2300 }
2301
2302 static evalue *find_relation_pair(evalue *e)
2303 {
2304 int i;
2305 evalue *found = NULL;
2306
2307 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2308 return NULL;
2309
2310 if (e->x.p->type == fractional) {
2311 Value m;
2312 evalue *cst;
2313
2314 value_init(m);
2315 poly_denom(&e->x.p->arr[0], &m);
2316
2317 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2318 cst = &cst->x.p->arr[0])
2319 ;
2320
2321 for (i = 1; !found && i < e->x.p->size; ++i)
2322 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2323
2324 value_clear(m);
2325 }
2326
2327 i = e->x.p->type == relation;
2328 for (; !found && i < e->x.p->size; ++i)
2329 found = find_relation_pair(&e->x.p->arr[i]);
2330
2331 return found;
2332 }
2333
2334 void evalue_mod2relation(evalue *e) {
2335 evalue *d;
2336
2337 if (value_zero_p(e->d) && e->x.p->type == partition) {
2338 int i;
2339
2340 for (i = 0; i < e->x.p->size/2; ++i) {
2341 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2342 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2343 value_clear(e->x.p->arr[2*i].d);
2344 free_evalue_refs(&e->x.p->arr[2*i+1]);
2345 e->x.p->size -= 2;
2346 if (2*i < e->x.p->size) {
2347 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2348 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2349 }
2350 --i;
2351 }
2352 }
2353 if (e->x.p->size == 0) {
2354 free(e->x.p);
2355 evalue_set_si(e, 0, 1);
2356 }
2357
2358 return;
2359 }
2360
2361 while ((d = find_relation_pair(e)) != NULL) {
2362 evalue split;
2363 evalue *ev;
2364
2365 value_init(split.d);
2366 value_set_si(split.d, 0);
2367 split.x.p = new_enode(relation, 3, 0);
2368 evalue_set_si(&split.x.p->arr[1], 1, 1);
2369 evalue_set_si(&split.x.p->arr[2], 1, 1);
2370
2371 ev = &split.x.p->arr[0];
2372 value_set_si(ev->d, 0);
2373 ev->x.p = new_enode(fractional, 3, -1);
2374 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2375 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2376 evalue_copy(&ev->x.p->arr[0], d);
2377
2378 emul(&split, e);
2379
2380 reduce_evalue(e);
2381
2382 free_evalue_refs(&split);
2383 }
2384 }
2385
2386 static int evalue_comp(const void * a, const void * b)
2387 {
2388 const evalue *e1 = *(const evalue **)a;
2389 const evalue *e2 = *(const evalue **)b;
2390 return ecmp(e1, e2);
2391 }
2392
2393 void evalue_combine(evalue *e)
2394 {
2395 evalue **evs;
2396 int i, k;
2397 enode *p;
2398 evalue tmp;
2399
2400 if (value_notzero_p(e->d) || e->x.p->type != partition)
2401 return;
2402
2403 NALLOC(evs, e->x.p->size/2);
2404 for (i = 0; i < e->x.p->size/2; ++i)
2405 evs[i] = &e->x.p->arr[2*i+1];
2406 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2407 p = new_enode(partition, e->x.p->size, e->x.p->pos);
2408 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2409 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2410 value_clear(p->arr[2*k].d);
2411 value_clear(p->arr[2*k+1].d);
2412 p->arr[2*k] = *(evs[i]-1);
2413 p->arr[2*k+1] = *(evs[i]);
2414 ++k;
2415 } else {
2416 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2417 Polyhedron *L = D;
2418
2419 value_clear((evs[i]-1)->d);
2420
2421 while (L->next)
2422 L = L->next;
2423 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2424 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2425 free_evalue_refs(evs[i]);
2426 }
2427 }
2428
2429 for (i = 2*k ; i < p->size; ++i)
2430 value_clear(p->arr[i].d);
2431
2432 free(evs);
2433 free(e->x.p);
2434 p->size = 2*k;
2435 e->x.p = p;
2436
2437 for (i = 0; i < e->x.p->size/2; ++i) {
2438 Polyhedron *H;
2439 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2440 continue;
2441 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2442 if (H == NULL)
2443 continue;
2444 for (k = 0; k < e->x.p->size/2; ++k) {
2445 Polyhedron *D, *N, **P;
2446 if (i == k)
2447 continue;
2448 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2449 D = *P;
2450 if (D == NULL)
2451 continue;
2452 for (; D; D = N) {
2453 *P = D;
2454 N = D->next;
2455 if (D->NbEq <= H->NbEq) {
2456 P = &D->next;
2457 continue;
2458 }
2459
2460 value_init(tmp.d);
2461 tmp.x.p = new_enode(partition, 2, e->x.p->pos);
2462 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2463 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2464 reduce_evalue(&tmp);
2465 if (value_notzero_p(tmp.d) ||
2466 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2467 P = &D->next;
2468 else {
2469 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2470 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2471 *P = NULL;
2472 }
2473 free_evalue_refs(&tmp);
2474 }
2475 }
2476 Polyhedron_Free(H);
2477 }
2478
2479 for (i = 0; i < e->x.p->size/2; ++i) {
2480 Polyhedron *H, *E;
2481 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2482 if (!D) {
2483 value_clear(e->x.p->arr[2*i].d);
2484 free_evalue_refs(&e->x.p->arr[2*i+1]);
2485 e->x.p->size -= 2;
2486 if (2*i < e->x.p->size) {
2487 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2488 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2489 }
2490 --i;
2491 continue;
2492 }
2493 if (!D->next)
2494 continue;
2495 H = DomainConvex(D, 0);
2496 E = DomainDifference(H, D, 0);
2497 Domain_Free(D);
2498 D = DomainDifference(H, E, 0);
2499 Domain_Free(H);
2500 Domain_Free(E);
2501 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2502 }
2503 }
2504
2505 /* May change coefficients to become non-standard if fiddle is set
2506 * => reduce p afterwards to correct
2507 */
2508 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2509 Matrix **R, int fiddle)
2510 {
2511 Polyhedron *I, *H;
2512 evalue *pp;
2513 unsigned dim = D->Dimension;
2514 Matrix *T = Matrix_Alloc(2, dim+1);
2515 Value twice;
2516 value_init(twice);
2517 assert(T);
2518
2519 assert(p->type == fractional);
2520 pp = &p->arr[0];
2521 value_set_si(T->p[1][dim], 1);
2522 poly_denom(pp, d);
2523 while (value_zero_p(pp->d)) {
2524 assert(pp->x.p->type == polynomial);
2525 assert(pp->x.p->size == 2);
2526 assert(value_notzero_p(pp->x.p->arr[1].d));
2527 value_division(T->p[0][pp->x.p->pos-1], *d, pp->x.p->arr[1].d);
2528 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2529 if (fiddle && value_gt(twice, pp->x.p->arr[1].d))
2530 value_substract(pp->x.p->arr[1].x.n,
2531 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2532 value_multiply(T->p[0][pp->x.p->pos-1],
2533 T->p[0][pp->x.p->pos-1], pp->x.p->arr[1].x.n);
2534 pp = &pp->x.p->arr[0];
2535 }
2536 value_division(T->p[0][dim], *d, pp->d);
2537 value_multiply(T->p[0][dim], T->p[0][dim], pp->x.n);
2538 I = DomainImage(D, T, 0);
2539 H = DomainConvex(I, 0);
2540 Domain_Free(I);
2541 *R = T;
2542
2543 value_clear(twice);
2544
2545 return H;
2546 }
2547
2548 static int reduce_in_domain(evalue *e, Polyhedron *D)
2549 {
2550 int i;
2551 enode *p;
2552 Value d, min, max;
2553 int r = 0;
2554 Polyhedron *I;
2555 Matrix *T;
2556 int bounded;
2557
2558 if (value_notzero_p(e->d))
2559 return r;
2560
2561 p = e->x.p;
2562
2563 if (p->type == relation) {
2564 int equal;
2565 value_init(d);
2566 value_init(min);
2567 value_init(max);
2568
2569 I = polynomial_projection(p->arr[0].x.p, D, &d, &T, 1);
2570 bounded = line_minmax(I, &min, &max); /* frees I */
2571 equal = value_eq(min, max);
2572 mpz_cdiv_q(min, min, d);
2573 mpz_fdiv_q(max, max, d);
2574
2575 if (bounded && value_gt(min, max)) {
2576 /* Never zero */
2577 if (p->size == 3) {
2578 value_clear(e->d);
2579 *e = p->arr[2];
2580 } else {
2581 evalue_set_si(e, 0, 1);
2582 r = 1;
2583 }
2584 free_evalue_refs(&(p->arr[1]));
2585 free_evalue_refs(&(p->arr[0]));
2586 free(p);
2587 value_clear(d);
2588 value_clear(min);
2589 value_clear(max);
2590 Matrix_Free(T);
2591 return r ? r : reduce_in_domain(e, D);
2592 } else if (bounded && equal) {
2593 /* Always zero */
2594 if (p->size == 3)
2595 free_evalue_refs(&(p->arr[2]));
2596 value_clear(e->d);
2597 *e = p->arr[1];
2598 free_evalue_refs(&(p->arr[0]));
2599 free(p);
2600 value_clear(d);
2601 value_clear(min);
2602 value_clear(max);
2603 Matrix_Free(T);
2604 return reduce_in_domain(e, D);
2605 } else if (bounded && value_eq(min, max)) {
2606 /* zero for a single value */
2607 Polyhedron *E;
2608 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2609 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2610 value_multiply(min, min, d);
2611 value_substract(M->p[0][D->Dimension+1],
2612 M->p[0][D->Dimension+1], min);
2613 E = DomainAddConstraints(D, M, 0);
2614 value_clear(d);
2615 value_clear(min);
2616 value_clear(max);
2617 Matrix_Free(T);
2618 Matrix_Free(M);
2619 r = reduce_in_domain(&p->arr[1], E);
2620 if (p->size == 3)
2621 r |= reduce_in_domain(&p->arr[2], D);
2622 Domain_Free(E);
2623 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2624 return r;
2625 }
2626
2627 value_clear(d);
2628 value_clear(min);
2629 value_clear(max);
2630 Matrix_Free(T);
2631 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2632 }
2633
2634 i = p->type == relation ? 1 :
2635 p->type == fractional ? 1 : 0;
2636 for (; i<p->size; i++)
2637 r |= reduce_in_domain(&p->arr[i], D);
2638
2639 if (p->type != fractional) {
2640 if (r && p->type == polynomial) {
2641 evalue f;
2642 value_init(f.d);
2643 value_set_si(f.d, 0);
2644 f.x.p = new_enode(polynomial, 2, p->pos);
2645 evalue_set_si(&f.x.p->arr[0], 0, 1);
2646 evalue_set_si(&f.x.p->arr[1], 1, 1);
2647 reorder_terms(p, &f);
2648 value_clear(e->d);
2649 *e = p->arr[0];
2650 free(p);
2651 }
2652 return r;
2653 }
2654
2655 value_init(d);
2656 value_init(min);
2657 value_init(max);
2658 I = polynomial_projection(p, D, &d, &T, 1);
2659 Matrix_Free(T);
2660 bounded = line_minmax(I, &min, &max); /* frees I */
2661 mpz_fdiv_q(min, min, d);
2662 mpz_fdiv_q(max, max, d);
2663 value_substract(d, max, min);
2664
2665 if (bounded && value_eq(min, max)) {
2666 evalue inc;
2667 value_init(inc.d);
2668 value_init(inc.x.n);
2669 value_set_si(inc.d, 1);
2670 value_oppose(inc.x.n, min);
2671 eadd(&inc, &p->arr[0]);
2672 reorder_terms(p, &p->arr[0]); /* frees arr[0] */
2673 value_clear(e->d);
2674 *e = p->arr[1];
2675 free(p);
2676 free_evalue_refs(&inc);
2677 r = 1;
2678 } else if (bounded && value_one_p(d) && p->size > 3) {
2679 evalue rem;
2680 evalue inc;
2681 evalue t;
2682 evalue f;
2683 evalue factor;
2684 value_init(rem.d);
2685 value_set_si(rem.d, 0);
2686 rem.x.p = new_enode(fractional, 3, -1);
2687 evalue_copy(&rem.x.p->arr[0], &p->arr[0]);
2688 rem.x.p->arr[1] = p->arr[1];
2689 rem.x.p->arr[2] = p->arr[2];
2690 for (i = 3; i < p->size; ++i)
2691 p->arr[i-2] = p->arr[i];
2692 p->size -= 2;
2693
2694 value_init(inc.d);
2695 value_init(inc.x.n);
2696 value_set_si(inc.d, 1);
2697 value_oppose(inc.x.n, min);
2698
2699 value_init(t.d);
2700 evalue_copy(&t, &p->arr[0]);
2701 eadd(&inc, &t);
2702
2703 value_init(f.d);
2704 value_set_si(f.d, 0);
2705 f.x.p = new_enode(fractional, 3, -1);
2706 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2707 evalue_set_si(&f.x.p->arr[1], 1, 1);
2708 evalue_set_si(&f.x.p->arr[2], 2, 1);
2709
2710 value_init(factor.d);
2711 evalue_set_si(&factor, -1, 1);
2712 emul(&t, &factor);
2713
2714 eadd(&f, &factor);
2715 emul(&t, &factor);
2716
2717 evalue_set_si(&f.x.p->arr[1], 0, 1);
2718 evalue_set_si(&f.x.p->arr[2], -1, 1);
2719 eadd(&f, &factor);
2720
2721 emul(&factor, e);
2722 eadd(&rem, e);
2723
2724 free_evalue_refs(&inc);
2725 free_evalue_refs(&t);
2726 free_evalue_refs(&f);
2727 free_evalue_refs(&factor);
2728 free_evalue_refs(&rem);
2729
2730 reduce_in_domain(e, D);
2731
2732 r = 1;
2733 } else {
2734 _reduce_evalue(&p->arr[0], 0, 1);
2735 if (r == 1) {
2736 evalue f;
2737 value_init(f.d);
2738 value_set_si(f.d, 0);
2739 f.x.p = new_enode(fractional, 3, -1);
2740 value_clear(f.x.p->arr[0].d);
2741 f.x.p->arr[0] = p->arr[0];
2742 evalue_set_si(&f.x.p->arr[1], 0, 1);
2743 evalue_set_si(&f.x.p->arr[2], 1, 1);
2744 reorder_terms(p, &f);
2745 value_clear(e->d);
2746 *e = p->arr[1];
2747 free(p);
2748 }
2749 }
2750
2751 value_clear(d);
2752 value_clear(min);
2753 value_clear(max);
2754
2755 return r;
2756 }
2757
2758 void evalue_range_reduction(evalue *e)
2759 {
2760 int i;
2761 if (value_notzero_p(e->d) || e->x.p->type != partition)
2762 return;
2763
2764 for (i = 0; i < e->x.p->size/2; ++i)
2765 if (reduce_in_domain(&e->x.p->arr[2*i+1],
2766 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
2767 reduce_evalue(&e->x.p->arr[2*i+1]);
2768
2769 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2770 free_evalue_refs(&e->x.p->arr[2*i+1]);
2771 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
2772 value_clear(e->x.p->arr[2*i].d);
2773 e->x.p->size -= 2;
2774 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2775 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2776 --i;
2777 }
2778 }
2779 }
2780
2781 /* Frees EP
2782 */
2783 Enumeration* partition2enumeration(evalue *EP)
2784 {
2785 int i;
2786 Enumeration *en, *res = NULL;
2787
2788 if (EVALUE_IS_ZERO(*EP)) {
2789 free(EP);
2790 return res;
2791 }
2792
2793 for (i = 0; i < EP->x.p->size/2; ++i) {
2794 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[2*i])->Dimension);
2795 en = (Enumeration *)malloc(sizeof(Enumeration));
2796 en->next = res;
2797 res = en;
2798 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2799 value_clear(EP->x.p->arr[2*i].d);
2800 res->EP = EP->x.p->arr[2*i+1];
2801 }
2802 free(EP->x.p);
2803 value_clear(EP->d);
2804 free(EP);
2805 return res;
2806 }
2807
2808 static int frac2floor_in_domain(evalue *e, Polyhedron *D)
2809 {
2810 enode *p;
2811 int r = 0;
2812 int i;
2813 Polyhedron *I;
2814 Matrix *T;
2815 Value d, min;
2816 evalue fl;
2817
2818 if (value_notzero_p(e->d))
2819 return r;
2820
2821 p = e->x.p;
2822
2823 i = p->type == relation ? 1 :
2824 p->type == fractional ? 1 : 0;
2825 for (; i<p->size; i++)
2826 r |= frac2floor_in_domain(&p->arr[i], D);
2827
2828 if (p->type != fractional) {
2829 if (r && p->type == polynomial) {
2830 evalue f;
2831 value_init(f.d);
2832 value_set_si(f.d, 0);
2833 f.x.p = new_enode(polynomial, 2, p->pos);
2834 evalue_set_si(&f.x.p->arr[0], 0, 1);
2835 evalue_set_si(&f.x.p->arr[1], 1, 1);
2836 reorder_terms(p, &f);
2837 value_clear(e->d);
2838 *e = p->arr[0];
2839 free(p);
2840 }
2841 return r;
2842 }
2843
2844 value_init(d);
2845 I = polynomial_projection(p, D, &d, &T, 0);
2846
2847 /*
2848 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2849 */
2850
2851 assert(I->NbEq == 0); /* Should have been reduced */
2852
2853 /* Find minimum */
2854 for (i = 0; i < I->NbConstraints; ++i)
2855 if (value_pos_p(I->Constraint[i][1]))
2856 break;
2857
2858 assert(i < I->NbConstraints);
2859 value_init(min);
2860 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
2861 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
2862 if (value_neg_p(min)) {
2863 evalue offset;
2864 mpz_fdiv_q(min, min, d);
2865 value_init(offset.d);
2866 value_set_si(offset.d, 1);
2867 value_init(offset.x.n);
2868 value_oppose(offset.x.n, min);
2869 eadd(&offset, &p->arr[0]);
2870 free_evalue_refs(&offset);
2871 }
2872
2873 Polyhedron_Free(I);
2874 Matrix_Free(T);
2875 value_clear(min);
2876 value_clear(d);
2877
2878 value_init(fl.d);
2879 value_set_si(fl.d, 0);
2880 fl.x.p = new_enode(flooring, 3, -1);
2881 evalue_set_si(&fl.x.p->arr[1], 0, 1);
2882 evalue_set_si(&fl.x.p->arr[2], -1, 1);
2883 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
2884
2885 eadd(&fl, &p->arr[0]);
2886 reorder_terms(p, &p->arr[0]);
2887 *e = p->arr[1];
2888 free(p);
2889 free_evalue_refs(&fl);
2890
2891 return 1;
2892 }
2893
2894 void evalue_frac2floor(evalue *e)
2895 {
2896 int i;
2897 if (value_notzero_p(e->d) || e->x.p->type != partition)
2898 return;
2899
2900 for (i = 0; i < e->x.p->size/2; ++i)
2901 if (frac2floor_in_domain(&e->x.p->arr[2*i+1],
2902 EVALUE_DOMAIN(e->x.p->arr[2*i])))
2903 reduce_evalue(&e->x.p->arr[2*i+1]);
2904 }
2905
2906 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
2907 Vector *row)
2908 {
2909 int nr, nc;
2910 int i;
2911 int nparam = D->Dimension - nvar;
2912
2913 if (C == 0) {
2914 nr = D->NbConstraints + 2;
2915 nc = D->Dimension + 2 + 1;
2916 C = Matrix_Alloc(nr, nc);
2917 for (i = 0; i < D->NbConstraints; ++i) {
2918 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
2919 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2920 D->Dimension + 1 - nvar);
2921 }
2922 } else {
2923 Matrix *oldC = C;
2924 nr = C->NbRows + 2;
2925 nc = C->NbColumns + 1;
2926 C = Matrix_Alloc(nr, nc);
2927 for (i = 0; i < oldC->NbRows; ++i) {
2928 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
2929 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2930 oldC->NbColumns - 1 - nvar);
2931 }
2932 }
2933 value_set_si(C->p[nr-2][0], 1);
2934 value_set_si(C->p[nr-2][1 + nvar], 1);
2935 value_set_si(C->p[nr-2][nc - 1], -1);
2936
2937 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
2938 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
2939 1 + nparam);
2940
2941 return C;
2942 }
2943
2944 static void floor2frac_r(evalue *e, int nvar)
2945 {
2946 enode *p;
2947 int i;
2948 evalue f;
2949 evalue *pp;
2950
2951 if (value_notzero_p(e->d))
2952 return;
2953
2954 p = e->x.p;
2955
2956 assert(p->type == flooring);
2957 for (i = 1; i < p->size; i++)
2958 floor2frac_r(&p->arr[i], nvar);
2959
2960 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
2961 assert(pp->x.p->type == polynomial);
2962 pp->x.p->pos -= nvar;
2963 }
2964
2965 value_init(f.d);
2966 value_set_si(f.d, 0);
2967 f.x.p = new_enode(fractional, 3, -1);
2968 evalue_set_si(&f.x.p->arr[1], 0, 1);
2969 evalue_set_si(&f.x.p->arr[2], -1, 1);
2970 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2971
2972 eadd(&f, &p->arr[0]);
2973 reorder_terms(p, &p->arr[0]);
2974 *e = p->arr[1];
2975 free(p);
2976 free_evalue_refs(&f);
2977 }
2978
2979 /* Convert flooring back to fractional and shift position
2980 * of the parameters by nvar
2981 */
2982 static void floor2frac(evalue *e, int nvar)
2983 {
2984 floor2frac_r(e, nvar);
2985 reduce_evalue(e);
2986 }
2987
2988 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
2989 {
2990 evalue *t;
2991 int nparam = D->Dimension - nvar;
2992
2993 if (C != 0) {
2994 C = Matrix_Copy(C);
2995 D = Constraints2Polyhedron(C, 0);
2996 Matrix_Free(C);
2997 }
2998
2999 t = barvinok_enumerate_e(D, 0, nparam, 0);
3000
3001 /* Double check that D was not unbounded. */
3002 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3003
3004 if (C != 0)
3005 Polyhedron_Free(D);
3006
3007 return t;
3008 }
3009
3010 evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3011 Matrix *C)
3012 {
3013 Vector *row = NULL;
3014 int i;
3015 evalue *res;
3016 Matrix *origC;
3017 evalue *factor = NULL;
3018 evalue cum;
3019
3020 if (EVALUE_IS_ZERO(*e))
3021 return 0;
3022
3023 if (D->next) {
3024 Polyhedron *DD = Disjoint_Domain(D, 0, 0);
3025 Polyhedron *Q;
3026
3027 Q = DD;
3028 DD = Q->next;
3029 Q->next = 0;
3030
3031 res = esum_over_domain(e, nvar, Q, C);
3032 Polyhedron_Free(Q);
3033
3034 for (Q = DD; Q; Q = DD) {
3035 evalue *t;
3036
3037 DD = Q->next;
3038 Q->next = 0;
3039
3040 t = esum_over_domain(e, nvar, Q, C);
3041 Polyhedron_Free(Q);
3042
3043 if (!res)
3044 res = t;
3045 else if (t) {
3046 eadd(t, res);
3047 free_evalue_refs(t);
3048 free(t);
3049 }
3050 }
3051 return res;
3052 }
3053
3054 if (value_notzero_p(e->d)) {
3055 evalue *t;
3056
3057 t = esum_over_domain_cst(nvar, D, C);
3058
3059 if (!EVALUE_IS_ONE(*e))
3060 emul(e, t);
3061
3062 return t;
3063 }
3064
3065 switch (e->x.p->type) {
3066 case flooring: {
3067 evalue *pp = &e->x.p->arr[0];
3068
3069 if (pp->x.p->pos > nvar) {
3070 /* remainder is independent of the summated vars */
3071 evalue f;
3072 evalue *t;
3073
3074 value_init(f.d);
3075 evalue_copy(&f, e);
3076 floor2frac(&f, nvar);
3077
3078 t = esum_over_domain_cst(nvar, D, C);
3079
3080 emul(&f, t);
3081
3082 free_evalue_refs(&f);
3083
3084 return t;
3085 }
3086
3087 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3088 poly_denom(pp, &row->p[1 + nvar]);
3089 value_set_si(row->p[0], 1);
3090 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3091 pp = &pp->x.p->arr[0]) {
3092 int pos;
3093 assert(pp->x.p->type == polynomial);
3094 pos = pp->x.p->pos;
3095 if (pos >= 1 + nvar)
3096 ++pos;
3097 value_assign(row->p[pos], row->p[1+nvar]);
3098 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3099 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3100 }
3101 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3102 value_division(row->p[1 + D->Dimension + 1],
3103 row->p[1 + D->Dimension + 1],
3104 pp->d);
3105 value_multiply(row->p[1 + D->Dimension + 1],
3106 row->p[1 + D->Dimension + 1],
3107 pp->x.n);
3108 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3109 break;
3110 }
3111 case polynomial: {
3112 int pos = e->x.p->pos;
3113
3114 if (pos > nvar) {
3115 ALLOC(factor);
3116 value_init(factor->d);
3117 value_set_si(factor->d, 0);
3118 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3119 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3120 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3121 break;
3122 }
3123
3124 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3125 for (i = 0; i < D->NbRays; ++i)
3126 if (value_notzero_p(D->Ray[i][pos]))
3127 break;
3128 assert(i < D->NbRays);
3129 if (value_neg_p(D->Ray[i][pos])) {
3130 ALLOC(factor);
3131 value_init(factor->d);
3132 evalue_set_si(factor, -1, 1);
3133 }
3134 value_set_si(row->p[0], 1);
3135 value_set_si(row->p[pos], 1);
3136 value_set_si(row->p[1 + nvar], -1);
3137 break;
3138 }
3139 default:
3140 assert(0);
3141 }
3142
3143 i = type_offset(e->x.p);
3144
3145 res = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3146 ++i;
3147
3148 if (factor) {
3149 value_init(cum.d);
3150 evalue_copy(&cum, factor);
3151 }
3152
3153 origC = C;
3154 for (; i < e->x.p->size; ++i) {
3155 evalue *t;
3156 if (row) {
3157 Matrix *prevC = C;
3158 C = esum_add_constraint(nvar, D, C, row);
3159 if (prevC != origC)
3160 Matrix_Free(prevC);
3161 }
3162 /*
3163 if (row)
3164 Vector_Print(stderr, P_VALUE_FMT, row);
3165 if (C)
3166 Matrix_Print(stderr, P_VALUE_FMT, C);
3167 */
3168 t = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3169
3170 if (t && factor)
3171 emul(&cum, t);
3172
3173 if (!res)
3174 res = t;
3175 else if (t) {
3176 eadd(t, res);
3177 free_evalue_refs(t);
3178 free(t);
3179 }
3180 if (factor && i+1 < e->x.p->size)
3181 emul(factor, &cum);
3182 }
3183 if (C != origC)
3184 Matrix_Free(C);
3185
3186 if (factor) {
3187 free_evalue_refs(factor);
3188 free_evalue_refs(&cum);
3189 free(factor);
3190 }
3191
3192 if (row)
3193 Vector_Free(row);
3194
3195 reduce_evalue(res);
3196
3197 return res;
3198 }
3199
3200 evalue *esum(evalue *e, int nvar)
3201 {
3202 int i;
3203 evalue *res;
3204 ALLOC(res);
3205 value_init(res->d);
3206
3207 assert(nvar >= 0);
3208 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3209 evalue_copy(res, e);
3210 return res;
3211 }
3212
3213 evalue_set_si(res, 0, 1);
3214
3215 assert(value_zero_p(e->d));
3216 assert(e->x.p->type == partition);
3217
3218 for (i = 0; i < e->x.p->size/2; ++i) {
3219 evalue *t;
3220 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3221 EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
3222 eadd(t, res);
3223 free_evalue_refs(t);
3224 free(t);
3225 }
3226
3227 reduce_evalue(res);
3228
3229 return res;
3230 }
3231
3232 /* Initial silly implementation */
3233 void eor(evalue *e1, evalue *res)
3234 {
3235 evalue E;
3236 evalue mone;
3237 value_init(E.d);
3238 value_init(mone.d);
3239 evalue_set_si(&mone, -1, 1);
3240
3241 evalue_copy(&E, res);
3242 eadd(e1, &E);
3243 emul(e1, res);
3244 emul(&mone, res);
3245 eadd(&E, res);
3246
3247 free_evalue_refs(&E);
3248 free_evalue_refs(&mone);
3249 }
lattice point enumeration library