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