?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
options.c: add argp parse
[barvinok.git] / evalue.c
bloba281b6e315e02cae2440a5fa1abb90a4b03fbb3b
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_domain3(evalue *e, Polyhedron *D, int shift)
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_domain3(&p->arr[i], D, shift);
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 if (shift) {
2912 value_init(d);
2913 I = polynomial_projection(p, D, &d, &T, 0);
2914
2915 /*
2916 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2917 */
2918
2919 assert(I->NbEq == 0); /* Should have been reduced */
2920
2921 /* Find minimum */
2922 for (i = 0; i < I->NbConstraints; ++i)
2923 if (value_pos_p(I->Constraint[i][1]))
2924 break;
2925
2926 if (i < I->NbConstraints) {
2927 value_init(min);
2928 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
2929 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
2930 if (value_neg_p(min)) {
2931 evalue offset;
2932 mpz_fdiv_q(min, min, d);
2933 value_init(offset.d);
2934 value_set_si(offset.d, 1);
2935 value_init(offset.x.n);
2936 value_oppose(offset.x.n, min);
2937 eadd(&offset, &p->arr[0]);
2938 free_evalue_refs(&offset);
2939 }
2940 value_clear(min);
2941 }
2942
2943 Polyhedron_Free(I);
2944 Matrix_Free(T);
2945 value_clear(d);
2946 }
2947
2948 value_init(fl.d);
2949 value_set_si(fl.d, 0);
2950 fl.x.p = new_enode(flooring, 3, -1);
2951 evalue_set_si(&fl.x.p->arr[1], 0, 1);
2952 evalue_set_si(&fl.x.p->arr[2], -1, 1);
2953 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
2954
2955 eadd(&fl, &p->arr[0]);
2956 reorder_terms(p, &p->arr[0]);
2957 value_clear(e->d);
2958 *e = p->arr[1];
2959 free(p);
2960 free_evalue_refs(&fl);
2961
2962 return 1;
2963 }
2964
2965 int evalue_frac2floor_in_domain(evalue *e, Polyhedron *D)
2966 {
2967 return evalue_frac2floor_in_domain3(e, D, 1);
2968 }
2969
2970 void evalue_frac2floor2(evalue *e, int shift)
2971 {
2972 int i;
2973 if (value_notzero_p(e->d) || e->x.p->type != partition)
2974 return;
2975
2976 for (i = 0; i < e->x.p->size/2; ++i)
2977 if (evalue_frac2floor_in_domain3(&e->x.p->arr[2*i+1],
2978 EVALUE_DOMAIN(e->x.p->arr[2*i]), shift))
2979 reduce_evalue(&e->x.p->arr[2*i+1]);
2980 }
2981
2982 void evalue_frac2floor(evalue *e)
2983 {
2984 evalue_frac2floor2(e, 1);
2985 }
2986
2987 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
2988 Vector *row)
2989 {
2990 int nr, nc;
2991 int i;
2992 int nparam = D->Dimension - nvar;
2993
2994 if (C == 0) {
2995 nr = D->NbConstraints + 2;
2996 nc = D->Dimension + 2 + 1;
2997 C = Matrix_Alloc(nr, nc);
2998 for (i = 0; i < D->NbConstraints; ++i) {
2999 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
3000 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3001 D->Dimension + 1 - nvar);
3002 }
3003 } else {
3004 Matrix *oldC = C;
3005 nr = C->NbRows + 2;
3006 nc = C->NbColumns + 1;
3007 C = Matrix_Alloc(nr, nc);
3008 for (i = 0; i < oldC->NbRows; ++i) {
3009 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
3010 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3011 oldC->NbColumns - 1 - nvar);
3012 }
3013 }
3014 value_set_si(C->p[nr-2][0], 1);
3015 value_set_si(C->p[nr-2][1 + nvar], 1);
3016 value_set_si(C->p[nr-2][nc - 1], -1);
3017
3018 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
3019 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
3020 1 + nparam);
3021
3022 return C;
3023 }
3024
3025 static void floor2frac_r(evalue *e, int nvar)
3026 {
3027 enode *p;
3028 int i;
3029 evalue f;
3030 evalue *pp;
3031
3032 if (value_notzero_p(e->d))
3033 return;
3034
3035 p = e->x.p;
3036
3037 assert(p->type == flooring);
3038 for (i = 1; i < p->size; i++)
3039 floor2frac_r(&p->arr[i], nvar);
3040
3041 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
3042 assert(pp->x.p->type == polynomial);
3043 pp->x.p->pos -= nvar;
3044 }
3045
3046 value_init(f.d);
3047 value_set_si(f.d, 0);
3048 f.x.p = new_enode(fractional, 3, -1);
3049 evalue_set_si(&f.x.p->arr[1], 0, 1);
3050 evalue_set_si(&f.x.p->arr[2], -1, 1);
3051 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3052
3053 eadd(&f, &p->arr[0]);
3054 reorder_terms(p, &p->arr[0]);
3055 value_clear(e->d);
3056 *e = p->arr[1];
3057 free(p);
3058 free_evalue_refs(&f);
3059 }
3060
3061 /* Convert flooring back to fractional and shift position
3062 * of the parameters by nvar
3063 */
3064 static void floor2frac(evalue *e, int nvar)
3065 {
3066 floor2frac_r(e, nvar);
3067 reduce_evalue(e);
3068 }
3069
3070 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
3071 {
3072 evalue *t;
3073 int nparam = D->Dimension - nvar;
3074
3075 if (C != 0) {
3076 C = Matrix_Copy(C);
3077 D = Constraints2Polyhedron(C, 0);
3078 Matrix_Free(C);
3079 }
3080
3081 t = barvinok_enumerate_e(D, 0, nparam, 0);
3082
3083 /* Double check that D was not unbounded. */
3084 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3085
3086 if (C != 0)
3087 Polyhedron_Free(D);
3088
3089 return t;
3090 }
3091
3092 evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3093 Matrix *C)
3094 {
3095 Vector *row = NULL;
3096 int i;
3097 evalue *res;
3098 Matrix *origC;
3099 evalue *factor = NULL;
3100 evalue cum;
3101
3102 if (EVALUE_IS_ZERO(*e))
3103 return 0;
3104
3105 if (D->next) {
3106 Polyhedron *DD = Disjoint_Domain(D, 0, 0);
3107 Polyhedron *Q;
3108
3109 Q = DD;
3110 DD = Q->next;
3111 Q->next = 0;
3112
3113 res = esum_over_domain(e, nvar, Q, C);
3114 Polyhedron_Free(Q);
3115
3116 for (Q = DD; Q; Q = DD) {
3117 evalue *t;
3118
3119 DD = Q->next;
3120 Q->next = 0;
3121
3122 t = esum_over_domain(e, nvar, Q, C);
3123 Polyhedron_Free(Q);
3124
3125 if (!res)
3126 res = t;
3127 else if (t) {
3128 eadd(t, res);
3129 free_evalue_refs(t);
3130 free(t);
3131 }
3132 }
3133 return res;
3134 }
3135
3136 if (value_notzero_p(e->d)) {
3137 evalue *t;
3138
3139 t = esum_over_domain_cst(nvar, D, C);
3140
3141 if (!EVALUE_IS_ONE(*e))
3142 emul(e, t);
3143
3144 return t;
3145 }
3146
3147 switch (e->x.p->type) {
3148 case flooring: {
3149 evalue *pp = &e->x.p->arr[0];
3150
3151 if (pp->x.p->pos > nvar) {
3152 /* remainder is independent of the summated vars */
3153 evalue f;
3154 evalue *t;
3155
3156 value_init(f.d);
3157 evalue_copy(&f, e);
3158 floor2frac(&f, nvar);
3159
3160 t = esum_over_domain_cst(nvar, D, C);
3161
3162 emul(&f, t);
3163
3164 free_evalue_refs(&f);
3165
3166 return t;
3167 }
3168
3169 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3170 poly_denom(pp, &row->p[1 + nvar]);
3171 value_set_si(row->p[0], 1);
3172 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3173 pp = &pp->x.p->arr[0]) {
3174 int pos;
3175 assert(pp->x.p->type == polynomial);
3176 pos = pp->x.p->pos;
3177 if (pos >= 1 + nvar)
3178 ++pos;
3179 value_assign(row->p[pos], row->p[1+nvar]);
3180 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3181 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3182 }
3183 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3184 value_division(row->p[1 + D->Dimension + 1],
3185 row->p[1 + D->Dimension + 1],
3186 pp->d);
3187 value_multiply(row->p[1 + D->Dimension + 1],
3188 row->p[1 + D->Dimension + 1],
3189 pp->x.n);
3190 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3191 break;
3192 }
3193 case polynomial: {
3194 int pos = e->x.p->pos;
3195
3196 if (pos > nvar) {
3197 factor = ALLOC(evalue);
3198 value_init(factor->d);
3199 value_set_si(factor->d, 0);
3200 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3201 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3202 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3203 break;
3204 }
3205
3206 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3207 for (i = 0; i < D->NbRays; ++i)
3208 if (value_notzero_p(D->Ray[i][pos]))
3209 break;
3210 assert(i < D->NbRays);
3211 if (value_neg_p(D->Ray[i][pos])) {
3212 factor = ALLOC(evalue);
3213 value_init(factor->d);
3214 evalue_set_si(factor, -1, 1);
3215 }
3216 value_set_si(row->p[0], 1);
3217 value_set_si(row->p[pos], 1);
3218 value_set_si(row->p[1 + nvar], -1);
3219 break;
3220 }
3221 default:
3222 assert(0);
3223 }
3224
3225 i = type_offset(e->x.p);
3226
3227 res = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3228 ++i;
3229
3230 if (factor) {
3231 value_init(cum.d);
3232 evalue_copy(&cum, factor);
3233 }
3234
3235 origC = C;
3236 for (; i < e->x.p->size; ++i) {
3237 evalue *t;
3238 if (row) {
3239 Matrix *prevC = C;
3240 C = esum_add_constraint(nvar, D, C, row);
3241 if (prevC != origC)
3242 Matrix_Free(prevC);
3243 }
3244 /*
3245 if (row)
3246 Vector_Print(stderr, P_VALUE_FMT, row);
3247 if (C)
3248 Matrix_Print(stderr, P_VALUE_FMT, C);
3249 */
3250 t = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3251
3252 if (t && factor)
3253 emul(&cum, t);
3254
3255 if (!res)
3256 res = t;
3257 else if (t) {
3258 eadd(t, res);
3259 free_evalue_refs(t);
3260 free(t);
3261 }
3262 if (factor && i+1 < e->x.p->size)
3263 emul(factor, &cum);
3264 }
3265 if (C != origC)
3266 Matrix_Free(C);
3267
3268 if (factor) {
3269 free_evalue_refs(factor);
3270 free_evalue_refs(&cum);
3271 free(factor);
3272 }
3273
3274 if (row)
3275 Vector_Free(row);
3276
3277 reduce_evalue(res);
3278
3279 return res;
3280 }
3281
3282 evalue *esum(evalue *e, int nvar)
3283 {
3284 int i;
3285 evalue *res = ALLOC(evalue);
3286 value_init(res->d);
3287
3288 assert(nvar >= 0);
3289 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3290 evalue_copy(res, e);
3291 return res;
3292 }
3293
3294 evalue_set_si(res, 0, 1);
3295
3296 assert(value_zero_p(e->d));
3297 assert(e->x.p->type == partition);
3298
3299 for (i = 0; i < e->x.p->size/2; ++i) {
3300 evalue *t;
3301 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3302 EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
3303 eadd(t, res);
3304 free_evalue_refs(t);
3305 free(t);
3306 }
3307
3308 reduce_evalue(res);
3309
3310 return res;
3311 }
3312
3313 /* Initial silly implementation */
3314 void eor(evalue *e1, evalue *res)
3315 {
3316 evalue E;
3317 evalue mone;
3318 value_init(E.d);
3319 value_init(mone.d);
3320 evalue_set_si(&mone, -1, 1);
3321
3322 evalue_copy(&E, res);
3323 eadd(e1, &E);
3324 emul(e1, res);
3325 emul(&mone, res);
3326 eadd(&E, res);
3327
3328 free_evalue_refs(&E);
3329 free_evalue_refs(&mone);
3330 }
3331
3332 /* computes denominator of polynomial evalue
3333 * d should point to a value initialized to 1
3334 */
3335 void evalue_denom(evalue *e, Value *d)
3336 {
3337 int i, offset;
3338
3339 if (value_notzero_p(e->d)) {
3340 value_lcm(*d, e->d, d);
3341 return;
3342 }
3343 assert(e->x.p->type == polynomial ||
3344 e->x.p->type == fractional ||
3345 e->x.p->type == flooring);
3346 offset = type_offset(e->x.p);
3347 for (i = e->x.p->size-1; i >= offset; --i)
3348 evalue_denom(&e->x.p->arr[i], d);
3349 }
lattice point enumeration library