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