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