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