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