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