?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
install headers
[barvinok.git] / ev_operations.c
blobba9c1c3f992e4c819cbf3de7ffd8a0b2f3d111ed
1 #include <assert.h>
2 #include <math.h>
3
4 #include "ev_operations.h"
5
6 void evalue_set_si(evalue *ev, int n, int d) {
7 value_set_si(ev->d, d);
8 value_init(ev->x.n);
9 value_set_si(ev->x.n, n);
10 }
11
12 void evalue_set(evalue *ev, Value n, Value d) {
13 value_assign(ev->d, d);
14 value_init(ev->x.n);
15 value_assign(ev->x.n, n);
16 }
17
18 void aep_evalue(evalue *e, int *ref) {
19
20 enode *p;
21 int i;
22
23 if (value_notzero_p(e->d))
24 return; /* a rational number, its already reduced */
25 if(!(p = e->x.p))
26 return; /* hum... an overflow probably occured */
27
28 /* First check the components of p */
29 for (i=0;i<p->size;i++)
30 aep_evalue(&p->arr[i],ref);
31
32 /* Then p itself */
33 if (p->type != modulo)
34 p->pos = ref[p->pos-1]+1;
35 return;
36 } /* aep_evalue */
37
38 /** Comments */
39 void addeliminatedparams_evalue(evalue *e,Matrix *CT) {
40
41 enode *p;
42 int i, j;
43 int *ref;
44
45 if (value_notzero_p(e->d))
46 return; /* a rational number, its already reduced */
47 if(!(p = e->x.p))
48 return; /* hum... an overflow probably occured */
49
50 /* Compute ref */
51 ref = (int *)malloc(sizeof(int)*(CT->NbRows-1));
52 for(i=0;i<CT->NbRows-1;i++)
53 for(j=0;j<CT->NbColumns;j++)
54 if(value_notzero_p(CT->p[i][j])) {
55 ref[i] = j;
56 break;
57 }
58
59 /* Transform the references in e, using ref */
60 aep_evalue(e,ref);
61 free( ref );
62 return;
63 } /* addeliminatedparams_evalue */
64
65 void reduce_evalue (evalue *e) {
66
67 enode *p;
68 int i, j, k;
69
70 if (value_notzero_p(e->d))
71 return; /* a rational number, its already reduced */
72 if(!(p = e->x.p))
73 return; /* hum... an overflow probably occured */
74
75 /* First reduce the components of p */
76 for (i=0; i<p->size; i++)
77 reduce_evalue(&p->arr[i]);
78
79 if (p->type==periodic) {
80
81 /* Try to reduce the period */
82 for (i=1; i<=(p->size)/2; i++) {
83 if ((p->size % i)==0) {
84
85 /* Can we reduce the size to i ? */
86 for (j=0; j<i; j++)
87 for (k=j+i; k<e->x.p->size; k+=i)
88 if (!eequal(&p->arr[j], &p->arr[k])) goto you_lose;
89
90 /* OK, lets do it */
91 for (j=i; j<p->size; j++) free_evalue_refs(&p->arr[j]);
92 p->size = i;
93 break;
94
95 you_lose: /* OK, lets not do it */
96 continue;
97 }
98 }
99
100 /* Try to reduce its strength */
101 if (p->size == 1) {
102 value_clear(e->d);
103 memcpy(e,&p->arr[0],sizeof(evalue));
104 free(p);
105 }
106 }
107 else if (p->type==polynomial) {
108
109 /* Try to reduce the degree */
110 for (i=p->size-1;i>=1;i--) {
111 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
112 break;
113 /* Zero coefficient */
114 free_evalue_refs(&(p->arr[i]));
115 }
116 if (i+1<p->size)
117 p->size = i+1;
118
119 /* Try to reduce its strength */
120 if (p->size == 1) {
121 value_clear(e->d);
122 memcpy(e,&p->arr[0],sizeof(evalue));
123 free(p);
124 }
125 }
126 else if (p->type==modulo) {
127
128 /* Try to reduce the degree */
129 for (i=p->size-1;i>=2;i--) {
130 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
131 break;
132 /* Zero coefficient */
133 free_evalue_refs(&(p->arr[i]));
134 }
135 if (i+1<p->size)
136 p->size = i+1;
137
138 /* Try to reduce its strength */
139 if (p->size == 2) {
140 value_clear(e->d);
141 memcpy(e,&p->arr[1],sizeof(evalue));
142 free_evalue_refs(&(p->arr[0]));
143 free(p);
144 }
145 }
146 } /* reduce_evalue */
147
148 void print_evalue(FILE *DST,evalue *e,char **pname) {
149
150 if(value_notzero_p(e->d)) {
151 if(value_notone_p(e->d)) {
152 value_print(DST,VALUE_FMT,e->x.n);
153 fprintf(DST,"/");
154 value_print(DST,VALUE_FMT,e->d);
155 }
156 else {
157 value_print(DST,VALUE_FMT,e->x.n);
158 }
159 }
160 else
161 print_enode(DST,e->x.p,pname);
162 return;
163 } /* print_evalue */
164
165 void print_enode(FILE *DST,enode *p,char **pname) {
166
167 int i;
168
169 if (!p) {
170 fprintf(DST, "NULL");
171 return;
172 }
173 switch (p->type) {
174 case evector:
175 fprintf(DST, "{ ");
176 for (i=0; i<p->size; i++) {
177 print_evalue(DST, &p->arr[i], pname);
178 if (i!=(p->size-1))
179 fprintf(DST, ", ");
180 }
181 fprintf(DST, " }\n");
182 break;
183 case polynomial:
184 fprintf(DST, "( ");
185 for (i=p->size-1; i>=0; i--) {
186 print_evalue(DST, &p->arr[i], pname);
187 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
188 else if (i>1)
189 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
190 }
191 fprintf(DST, " )\n");
192 break;
193 case periodic:
194 fprintf(DST, "[ ");
195 for (i=0; i<p->size; i++) {
196 print_evalue(DST, &p->arr[i], pname);
197 if (i!=(p->size-1)) fprintf(DST, ", ");
198 }
199 fprintf(DST," ]_%s", pname[p->pos-1]);
200 break;
201 case modulo:
202 fprintf(DST, "( ");
203 for (i=p->size-1; i>=1; i--) {
204 print_evalue(DST, &p->arr[i], pname);
205 if (i >= 2) {
206 fprintf(DST, " * ");
207 if (i > 2)
208 fprintf(DST, "(");
209 fprintf(DST, "(");
210 print_evalue(DST, &p->arr[0], pname);
211 fprintf(DST, ") mod %d", p->pos);
212 if (i>2)
213 fprintf(DST, ")^%d + ", i-1);
214 else
215 fprintf(DST, " + ", i-1);
216 }
217 }
218 fprintf(DST, " )\n");
219 break;
220 case relation:
221 fprintf(DST, "[ ");
222 print_evalue(DST, &p->arr[0], pname);
223 fprintf(DST, "= 0 ] * \n");
224 print_evalue(DST, &p->arr[1], pname);
225 break;
226 case partition:
227 for (i=0; i<p->size/2; i++) {
228 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), pname);
229 print_evalue(DST, &p->arr[2*i+1], pname);
230 }
231 break;
232 default:
233 assert(0);
234 }
235 return;
236 } /* print_enode */
237
238 static int mod_term_smaller(evalue *e1, evalue *e2)
239 {
240 if (value_notzero_p(e1->d)) {
241 if (value_zero_p(e2->d))
242 return 1;
243 return value_lt(e1->x.n, e2->x.n);
244 }
245 if (value_notzero_p(e2->d))
246 return 0;
247 if (e1->x.p->pos < e2->x.p->pos)
248 return 1;
249 else if (e1->x.p->pos > e2->x.p->pos)
250 return 0;
251 else
252 return mod_term_smaller(&e1->x.p->arr[0], &e2->x.p->arr[0]);
253 }
254
255 static void eadd_rev(evalue *e1, evalue *res)
256 {
257 evalue ev;
258 value_init(ev.d);
259 evalue_copy(&ev, e1);
260 eadd(res, &ev);
261 free_evalue_refs(res);
262 *res = ev;
263 }
264
265 static void eadd_rev_cst (evalue *e1, evalue *res)
266 {
267 evalue ev;
268 value_init(ev.d);
269 evalue_copy(&ev, e1);
270 eadd(res, &ev.x.p->arr[ev.x.p->type==modulo]);
271 free_evalue_refs(res);
272 *res = ev;
273 }
274
275 struct section { Polyhedron * D; evalue E; };
276
277 void eadd_partitions (evalue *e1,evalue *res)
278 {
279 int n, i, j;
280 Polyhedron *d, *fd;
281 struct section *s;
282 s = (struct section *)
283 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
284 sizeof(struct section));
285 assert(s);
286
287 n = 0;
288 for (j = 0; j < e1->x.p->size/2; ++j) {
289 assert(res->x.p->size >= 2);
290 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
291 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
292 if (!emptyQ(fd))
293 for (i = 1; i < res->x.p->size/2; ++i) {
294 Polyhedron *t = fd;
295 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
296 Domain_Free(t);
297 if (emptyQ(fd))
298 break;
299 }
300 if (emptyQ(fd)) {
301 Domain_Free(fd);
302 continue;
303 }
304 value_init(s[n].E.d);
305 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
306 s[n].D = fd;
307 ++n;
308 }
309 for (i = 0; i < res->x.p->size/2; ++i) {
310 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
311 for (j = 0; j < e1->x.p->size/2; ++j) {
312 Polyhedron *t;
313 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
314 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
315 if (emptyQ(d)) {
316 Domain_Free(d);
317 continue;
318 }
319 t = fd;
320 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
321 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
322 Domain_Free(t);
323 value_init(s[n].E.d);
324 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
325 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
326 s[n].D = d;
327 ++n;
328 }
329 if (!emptyQ(fd)) {
330 s[n].E = res->x.p->arr[2*i+1];
331 s[n].D = fd;
332 ++n;
333 } else
334 free_evalue_refs(&res->x.p->arr[2*i+1]);
335 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
336 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
337 value_clear(res->x.p->arr[2*i].d);
338 }
339
340 free(res->x.p);
341 res->x.p = new_enode(partition, 2*n, -1);
342 for (j = 0; j < n; ++j) {
343 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
344 value_clear(res->x.p->arr[2*j+1].d);
345 res->x.p->arr[2*j+1] = s[j].E;
346 }
347
348 free(s);
349 }
350
351 void eadd(evalue *e1,evalue *res) {
352
353 int i;
354 if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
355 /* Add two rational numbers */
356 Value g,m1,m2;
357 value_init(g);
358 value_init(m1);
359 value_init(m2);
360
361 value_multiply(m1,e1->x.n,res->d);
362 value_multiply(m2,res->x.n,e1->d);
363 value_addto(res->x.n,m1,m2);
364 value_multiply(res->d,e1->d,res->d);
365 Gcd(res->x.n,res->d,&g);
366 if (value_notone_p(g)) {
367 value_division(res->d,res->d,g);
368 value_division(res->x.n,res->x.n,g);
369 }
370 value_clear(g); value_clear(m1); value_clear(m2);
371 return ;
372 }
373 else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
374 switch (res->x.p->type) {
375 case polynomial:
376 /* Add the constant to the constant term of a polynomial*/
377 eadd(e1, &res->x.p->arr[0]);
378 return ;
379 case periodic:
380 /* Add the constant to all elements of a periodic number */
381 for (i=0; i<res->x.p->size; i++) {
382 eadd(e1, &res->x.p->arr[i]);
383 }
384 return ;
385 case evector:
386 fprintf(stderr, "eadd: cannot add const with vector\n");
387 return;
388 case modulo:
389 eadd(e1, &res->x.p->arr[1]);
390 return ;
391 case partition:
392 assert(EVALUE_IS_ZERO(*e1));
393 break; /* Do nothing */
394 case relation:
395 for (i = 1; i < res->x.p->size; ++i)
396 eadd(e1, &res->x.p->arr[i]);
397 break;
398 default:
399 assert(0);
400 }
401 }
402 /* add polynomial or periodic to constant
403 * you have to exchange e1 and res, before doing addition */
404
405 else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
406 eadd_rev(e1, res);
407 return;
408 }
409 else { // ((e1->d==0) && (res->d==0))
410 assert(!((e1->x.p->type == partition) ^
411 (res->x.p->type == partition)));
412 if (e1->x.p->type == partition) {
413 eadd_partitions(e1, res);
414 return;
415 }
416 if (e1->x.p->type == relation) {
417 eadd_rev(e1, res);
418 return;
419 }
420 if (res->x.p->type == relation) {
421 for (i = 1; i < res->x.p->size; ++i)
422 eadd(e1, &res->x.p->arr[i]);
423 return;
424 }
425 if ((e1->x.p->type != res->x.p->type) ) {
426 /* adding to evalues of different type. two cases are possible
427 * res is periodic and e1 is polynomial, you have to exchange
428 * e1 and res then to add e1 to the constant term of res */
429 if (e1->x.p->type == polynomial) {
430 eadd_rev_cst(e1, res);
431 }
432 else if (res->x.p->type == polynomial) {
433 /* res is polynomial and e1 is periodic,
434 add e1 to the constant term of res */
435
436 eadd(e1,&res->x.p->arr[0]);
437 } else
438 assert(0);
439
440 return;
441 }
442 else if (e1->x.p->pos != res->x.p->pos ||
443 (res->x.p->type == modulo &&
444 !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
445 /* adding evalues of different position (i.e function of different unknowns
446 * to case are possible */
447
448 switch (res->x.p->type) {
449 case modulo:
450 if(mod_term_smaller(res, e1))
451 eadd(e1,&res->x.p->arr[1]);
452 else
453 eadd_rev_cst(e1, res);
454 return;
455 case polynomial: // res and e1 are polynomials
456 // add e1 to the constant term of res
457
458 if(res->x.p->pos < e1->x.p->pos)
459 eadd(e1,&res->x.p->arr[0]);
460 else
461 eadd_rev_cst(e1, res);
462 // value_clear(g); value_clear(m1); value_clear(m2);
463 return;
464 case periodic: // res and e1 are pointers to periodic numbers
465 //add e1 to all elements of res
466
467 if(res->x.p->pos < e1->x.p->pos)
468 for (i=0;i<res->x.p->size;i++) {
469 eadd(e1,&res->x.p->arr[i]);
470 }
471 else
472 eadd_rev(e1, res);
473 return;
474 }
475 }
476
477
478 //same type , same pos and same size
479 if (e1->x.p->size == res->x.p->size) {
480 // add any element in e1 to the corresponding element in res
481 if (res->x.p->type == modulo)
482 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
483 i = res->x.p->type == modulo ? 1 : 0;
484 for (; i<res->x.p->size; i++) {
485 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
486 }
487 return ;
488 }
489
490 /* Sizes are different */
491 switch(res->x.p->type) {
492 case polynomial:
493 case modulo:
494 /* VIN100: if e1-size > res-size you have to copy e1 in a */
495 /* new enode and add res to that new node. If you do not do */
496 /* that, you lose the the upper weight part of e1 ! */
497
498 if(e1->x.p->size > res->x.p->size)
499 eadd_rev(e1, res);
500 else {
501
502 if (res->x.p->type == modulo)
503 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
504 i = res->x.p->type == modulo ? 1 : 0;
505 for (; i<e1->x.p->size ; i++) {
506 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
507 }
508
509 return ;
510 }
511 break;
512
513 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
514 case periodic:
515 {
516 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
517 of the sizes of e1 and res, then to copy res periodicaly in ne, after
518 to add periodicaly elements of e1 to elements of ne, and finaly to
519 return ne. */
520 int x,y,p;
521 Value ex, ey ,ep;
522 evalue *ne;
523 value_init(ex); value_init(ey);value_init(ep);
524 x=e1->x.p->size;
525 y= res->x.p->size;
526 value_set_si(ex,e1->x.p->size);
527 value_set_si(ey,res->x.p->size);
528 value_assign (ep,*Lcm(ex,ey));
529 p=(int)mpz_get_si(ep);
530 ne= (evalue *) malloc (sizeof(evalue));
531 value_init(ne->d);
532 value_set_si( ne->d,0);
533
534 ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
535 for(i=0;i<p;i++) {
536 value_assign(ne->x.p->arr[i].d, res->x.p->arr[i%y].d);
537 if (value_notzero_p(ne->x.p->arr[i].d)) {
538 value_init(ne->x.p->arr[i].x.n);
539 value_assign(ne->x.p->arr[i].x.n, res->x.p->arr[i%y].x.n);
540 }
541 else {
542 ne->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%y].x.p);
543 }
544 }
545 for(i=0;i<p;i++) {
546 eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
547 }
548
549 value_assign(res->d, ne->d);
550 res->x.p=ne->x.p;
551
552 return ;
553 }
554 case evector:
555 fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
556 return ;
557 }
558 }
559 return ;
560 }/* eadd */
561
562 static void emul_rev (evalue *e1, evalue *res)
563 {
564 evalue ev;
565 value_init(ev.d);
566 evalue_copy(&ev, e1);
567 emul(res, &ev);
568 free_evalue_refs(res);
569 *res = ev;
570 }
571
572 static void emul_poly (evalue *e1, evalue *res)
573 {
574 int i, j, o = res->x.p->type == modulo;
575 evalue tmp;
576 int size=(e1->x.p->size + res->x.p->size - o - 1);
577 value_init(tmp.d);
578 value_set_si(tmp.d,0);
579 tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
580 if (o)
581 evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
582 for (i=o; i < e1->x.p->size; i++) {
583 evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
584 emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
585 }
586 for (; i<size; i++)
587 evalue_set_si(&tmp.x.p->arr[i], 0, 1);
588 for (i=o+1; i<res->x.p->size; i++)
589 for (j=o; j<e1->x.p->size; j++) {
590 evalue ev;
591 value_init(ev.d);
592 evalue_copy(&ev, &e1->x.p->arr[j]);
593 emul(&res->x.p->arr[i], &ev);
594 eadd(&ev, &tmp.x.p->arr[i+j-o]);
595 free_evalue_refs(&ev);
596 }
597 free_evalue_refs(res);
598 *res = tmp;
599 }
600
601 void emul_partitions (evalue *e1,evalue *res)
602 {
603 int n, i, j, k;
604 Polyhedron *d;
605 struct section *s;
606 s = (struct section *)
607 malloc((e1->x.p->size/2) * (res->x.p->size/2) *
608 sizeof(struct section));
609 assert(s);
610
611 n = 0;
612 for (i = 0; i < res->x.p->size/2; ++i) {
613 for (j = 0; j < e1->x.p->size/2; ++j) {
614 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
615 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
616 if (emptyQ(d)) {
617 Domain_Free(d);
618 continue;
619 }
620
621 /* This code is only needed because the partitions
622 are not true partitions.
623 */
624 for (k = 0; k < n; ++k) {
625 if (DomainIncludes(s[k].D, d))
626 break;
627 if (DomainIncludes(d, s[k].D)) {
628 Domain_Free(s[k].D);
629 free_evalue_refs(&s[k].E);
630 if (n > k)
631 s[k] = s[--n];
632 --k;
633 }
634 }
635 if (k < n) {
636 Domain_Free(d);
637 continue;
638 }
639
640 value_init(s[n].E.d);
641 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
642 emul(&e1->x.p->arr[2*j+1], &s[n].E);
643 s[n].D = d;
644 ++n;
645 }
646 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
647 value_clear(res->x.p->arr[2*i].d);
648 free_evalue_refs(&res->x.p->arr[2*i+1]);
649 }
650
651 free(res->x.p);
652 res->x.p = new_enode(partition, 2*n, -1);
653 for (j = 0; j < n; ++j) {
654 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
655 value_clear(res->x.p->arr[2*j+1].d);
656 res->x.p->arr[2*j+1] = s[j].E;
657 }
658
659 free(s);
660 }
661
662 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
663 * do a copy of "res" befor calling this function if you nead it after. The vector type of
664 * evalues is not treated here */
665
666 void emul (evalue *e1, evalue *res ){
667 int i,j;
668
669 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
670 fprintf(stderr, "emul: do not proced on evector type !\n");
671 return;
672 }
673
674 if (value_zero_p(e1->d) && e1->x.p->type == partition) {
675 if (value_zero_p(res->d) && res->x.p->type == partition)
676 emul_partitions(e1, res);
677 else
678 emul_rev(e1, res);
679 } else if (value_zero_p(res->d) && res->x.p->type == partition) {
680 for (i = 0; i < res->x.p->size/2; ++i)
681 emul(e1, &res->x.p->arr[2*i+1]);
682 } else
683 if (value_zero_p(res->d) && res->x.p->type == relation)
684 emul(e1, &res->x.p->arr[1]);
685 else
686 if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
687 switch(e1->x.p->type) {
688 case polynomial:
689 switch(res->x.p->type) {
690 case polynomial:
691 if(e1->x.p->pos == res->x.p->pos) {
692 /* Product of two polynomials of the same variable */
693 emul_poly(e1, res);
694 return;
695 }
696 else {
697 /* Product of two polynomials of different variables */
698
699 if(res->x.p->pos < e1->x.p->pos)
700 for( i=0; i<res->x.p->size ; i++)
701 emul(e1, &res->x.p->arr[i]);
702 else
703 emul_rev(e1, res);
704
705 return ;
706 }
707 case periodic:
708 case modulo:
709 /* Product of a polynomial and a periodic or modulo */
710 emul_rev(e1, res);
711 return;
712 }
713 case periodic:
714 switch(res->x.p->type) {
715 case periodic:
716 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
717 /* Product of two periodics of the same parameter and period */
718
719 for(i=0; i<res->x.p->size;i++)
720 emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
721
722 return;
723 }
724 else{
725 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
726 /* Product of two periodics of the same parameter and different periods */
727 evalue *newp;
728 Value x,y,z;
729 int ix,iy,lcm;
730 value_init(x); value_init(y);value_init(z);
731 ix=e1->x.p->size;
732 iy=res->x.p->size;
733 value_set_si(x,e1->x.p->size);
734 value_set_si(y,res->x.p->size);
735 value_assign (z,*Lcm(x,y));
736 lcm=(int)mpz_get_si(z);
737 newp= (evalue *) malloc (sizeof(evalue));
738 value_init(newp->d);
739 value_set_si( newp->d,0);
740 newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
741 for(i=0;i<lcm;i++) {
742 value_assign(newp->x.p->arr[i].d, res->x.p->arr[i%iy].d);
743 if (value_notzero_p(newp->x.p->arr[i].d)) {
744 value_assign(newp->x.p->arr[i].x.n, res->x.p->arr[i%iy].x.n);
745 }
746 else {
747 newp->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%iy].x.p);
748 }
749
750 }
751 for(i=0;i<lcm;i++)
752 emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
753
754 value_assign(res->d,newp->d);
755 res->x.p=newp->x.p;
756
757 value_clear(x); value_clear(y);value_clear(z);
758 return ;
759 }
760 else {
761 /* Product of two periodics of different parameters */
762
763 for(i=0; i<res->x.p->size; i++)
764 emul(e1, &(res->x.p->arr[i]));
765
766 return;
767 }
768 }
769 case polynomial:
770 /* Product of a periodic and a polynomial */
771
772 for(i=0; i<res->x.p->size ; i++)
773 emul(e1, &(res->x.p->arr[i]));
774
775 return;
776
777 }
778 case modulo:
779 switch(res->x.p->type) {
780 case polynomial:
781 for(i=0; i<res->x.p->size ; i++)
782 emul(e1, &(res->x.p->arr[i]));
783 return;
784 case periodic:
785 assert(0);
786 case modulo:
787 if (e1->x.p->pos == res->x.p->pos &&
788 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
789 if (e1->x.p->pos != 2)
790 emul_poly(e1, res);
791 else {
792 evalue tmp;
793 /* x mod 2 == (x mod 2)^2 */
794 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1) (x mod 2) */
795 assert(e1->x.p->size == 3);
796 assert(res->x.p->size == 3);
797 value_init(tmp.d);
798 evalue_copy(&tmp, &res->x.p->arr[1]);
799 eadd(&res->x.p->arr[2], &tmp);
800 emul(&e1->x.p->arr[2], &tmp);
801 emul(&e1->x.p->arr[1], res);
802 eadd(&tmp, &res->x.p->arr[2]);
803 free_evalue_refs(&tmp);
804 }
805 } else {
806 if(mod_term_smaller(res, e1))
807 for(i=1; i<res->x.p->size ; i++)
808 emul(e1, &(res->x.p->arr[i]));
809 else
810 emul_rev(e1, res);
811 return;
812 }
813 }
814 break;
815 case relation:
816 emul_rev(e1, res);
817 break;
818 default:
819 assert(0);
820 }
821 }
822 else {
823 if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
824 /* Product of two rational numbers */
825
826 Value g;
827 value_init(g);
828 value_multiply(res->d,e1->d,res->d);
829 value_multiply(res->x.n,e1->x.n,res->x.n );
830 Gcd(res->x.n, res->d,&g);
831 if (value_notone_p(g)) {
832 value_division(res->d,res->d,g);
833 value_division(res->x.n,res->x.n,g);
834 }
835 value_clear(g);
836 return ;
837 }
838 else {
839 if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
840 /* Product of an expression (polynomial or peririodic) and a rational number */
841
842 emul_rev(e1, res);
843 return ;
844 }
845 else {
846 /* Product of a rationel number and an expression (polynomial or peririodic) */
847
848 i = res->x.p->type == modulo ? 1 : 0;
849 for (; i<res->x.p->size; i++)
850 emul(e1, &res->x.p->arr[i]);
851
852 return ;
853 }
854 }
855 }
856
857 return ;
858 }
859
860 /* Frees mask ! */
861 void emask(evalue *mask, evalue *res) {
862 int n, i, j;
863 Polyhedron *d, *fd;
864 struct section *s;
865 evalue mone;
866
867 assert(mask->x.p->type == partition);
868 assert(res->x.p->type == partition);
869
870 s = (struct section *)
871 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
872 sizeof(struct section));
873 assert(s);
874
875 value_init(mone.d);
876 evalue_set_si(&mone, -1, 1);
877
878 n = 0;
879 for (j = 0; j < res->x.p->size/2; ++j) {
880 assert(mask->x.p->size >= 2);
881 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
882 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
883 if (!emptyQ(fd))
884 for (i = 1; i < mask->x.p->size/2; ++i) {
885 Polyhedron *t = fd;
886 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
887 Domain_Free(t);
888 if (emptyQ(fd))
889 break;
890 }
891 if (emptyQ(fd)) {
892 Domain_Free(fd);
893 continue;
894 }
895 value_init(s[n].E.d);
896 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
897 s[n].D = fd;
898 ++n;
899 }
900 for (i = 0; i < mask->x.p->size/2; ++i) {
901 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
902 eadd(&mone, &mask->x.p->arr[2*i+1]);
903 emul(&mone, &mask->x.p->arr[2*i+1]);
904 for (j = 0; j < res->x.p->size/2; ++j) {
905 Polyhedron *t;
906 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
907 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
908 if (emptyQ(d)) {
909 Domain_Free(d);
910 continue;
911 }
912 t = fd;
913 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
914 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
915 Domain_Free(t);
916 value_init(s[n].E.d);
917 evalue_copy(&s[n].E, &mask->x.p->arr[2*i+1]);
918 emul(&res->x.p->arr[2*j+1], &s[n].E);
919 s[n].D = d;
920 ++n;
921 }
922 if (!emptyQ(fd)) {
923 assert(0); // We don't allow this.
924 }
925 }
926
927 free_evalue_refs(&mone);
928 free_evalue_refs(mask);
929 free_evalue_refs(res);
930 value_init(res->d);
931 res->x.p = new_enode(partition, 2*n, -1);
932 for (j = 0; j < n; ++j) {
933 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
934 value_clear(res->x.p->arr[2*j+1].d);
935 res->x.p->arr[2*j+1] = s[j].E;
936 }
937
938 free(s);
939 }
940
941 void evalue_copy(evalue *dst, evalue *src)
942 {
943 value_assign(dst->d, src->d);
944 if(value_notzero_p(src->d)) {
945 value_init(dst->x.n);
946 value_assign(dst->x.n, src->x.n);
947 } else
948 dst->x.p = ecopy(src->x.p);
949 }
950
951 enode *new_enode(enode_type type,int size,int pos) {
952
953 enode *res;
954 int i;
955
956 if(size == 0) {
957 fprintf(stderr, "Allocating enode of size 0 !\n" );
958 return NULL;
959 }
960 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
961 res->type = type;
962 res->size = size;
963 res->pos = pos;
964 for(i=0; i<size; i++) {
965 value_init(res->arr[i].d);
966 value_set_si(res->arr[i].d,0);
967 res->arr[i].x.p = 0;
968 }
969 return res;
970 } /* new_enode */
971
972 enode *ecopy(enode *e) {
973
974 enode *res;
975 int i;
976
977 res = new_enode(e->type,e->size,e->pos);
978 for(i=0;i<e->size;++i) {
979 value_assign(res->arr[i].d,e->arr[i].d);
980 if(value_zero_p(res->arr[i].d))
981 res->arr[i].x.p = ecopy(e->arr[i].x.p);
982 else if (EVALUE_IS_DOMAIN(res->arr[i]))
983 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
984 else {
985 value_init(res->arr[i].x.n);
986 value_assign(res->arr[i].x.n,e->arr[i].x.n);
987 }
988 }
989 return(res);
990 } /* ecopy */
991
992 int eequal(evalue *e1,evalue *e2) {
993
994 int i;
995 enode *p1, *p2;
996
997 if (value_ne(e1->d,e2->d))
998 return 0;
999
1000 /* e1->d == e2->d */
1001 if (value_notzero_p(e1->d)) {
1002 if (value_ne(e1->x.n,e2->x.n))
1003 return 0;
1004
1005 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1006 return 1;
1007 }
1008
1009 /* e1->d == e2->d == 0 */
1010 p1 = e1->x.p;
1011 p2 = e2->x.p;
1012 if (p1->type != p2->type) return 0;
1013 if (p1->size != p2->size) return 0;
1014 if (p1->pos != p2->pos) return 0;
1015 for (i=0; i<p1->size; i++)
1016 if (!eequal(&p1->arr[i], &p2->arr[i]) )
1017 return 0;
1018 return 1;
1019 } /* eequal */
1020
1021 void free_evalue_refs(evalue *e) {
1022
1023 enode *p;
1024 int i;
1025
1026 if (EVALUE_IS_DOMAIN(*e)) {
1027 Domain_Free(EVALUE_DOMAIN(*e));
1028 value_clear(e->d);
1029 return;
1030 } else if (value_pos_p(e->d)) {
1031
1032 /* 'e' stores a constant */
1033 value_clear(e->d);
1034 value_clear(e->x.n);
1035 return;
1036 }
1037 assert(value_zero_p(e->d));
1038 value_clear(e->d);
1039 p = e->x.p;
1040 if (!p) return; /* null pointer */
1041 for (i=0; i<p->size; i++) {
1042 free_evalue_refs(&(p->arr[i]));
1043 }
1044 free(p);
1045 return;
1046 } /* free_evalue_refs */
1047
1048 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
1049 Vector * val, evalue *res)
1050 {
1051 unsigned nparam = periods->Size;
1052
1053 if (p == nparam) {
1054 double d = compute_evalue(e, val->p);
1055 if (d > 0)
1056 d += .25;
1057 else
1058 d -= .25;
1059 value_set_si(res->d, 1);
1060 value_init(res->x.n);
1061 value_set_double(res->x.n, d);
1062 mpz_fdiv_r(res->x.n, res->x.n, m);
1063 return;
1064 }
1065 if (value_one_p(periods->p[p]))
1066 mod2table_r(e, periods, m, p+1, val, res);
1067 else {
1068 Value tmp;
1069 value_init(tmp);
1070
1071 value_assign(tmp, periods->p[p]);
1072 value_set_si(res->d, 0);
1073 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
1074 do {
1075 value_decrement(tmp, tmp);
1076 value_assign(val->p[p], tmp);
1077 mod2table_r(e, periods, m, p+1, val,
1078 &res->x.p->arr[VALUE_TO_INT(tmp)]);
1079 } while (value_pos_p(tmp));
1080
1081 value_clear(tmp);
1082 }
1083 }
1084
1085 void evalue_mod2table(evalue *e, int nparam)
1086 {
1087 enode *p;
1088 int i;
1089
1090 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
1091 return;
1092 p = e->x.p;
1093 for (i=0; i<p->size; i++) {
1094 evalue_mod2table(&(p->arr[i]), nparam);
1095 }
1096 if (p->type == modulo) {
1097 Vector *periods = Vector_Alloc(nparam);
1098 Vector *val = Vector_Alloc(nparam);
1099 Value tmp;
1100 evalue *ev;
1101 evalue EP, res;
1102
1103 value_init(tmp);
1104 value_set_si(tmp, p->pos);
1105 Vector_Set(periods->p, 1, nparam);
1106 Vector_Set(val->p, 0, nparam);
1107 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
1108 enode *p = ev->x.p;
1109
1110 assert(p->type == polynomial);
1111 assert(p->size == 2);
1112 assert(value_one_p(p->arr[1].d));
1113 Gcd(tmp, p->arr[1].x.n, &periods->p[p->pos-1]);
1114 value_division(periods->p[p->pos-1], tmp, periods->p[p->pos-1]);
1115 }
1116 value_set_si(tmp, p->pos);
1117 value_init(EP.d);
1118 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
1119
1120 value_init(res.d);
1121 evalue_set_si(&res, 0, 1);
1122 /* Compute the polynomial using Horner's rule */
1123 for (i=p->size-1;i>1;i--) {
1124 eadd(&p->arr[i], &res);
1125 emul(&EP, &res);
1126 }
1127 eadd(&p->arr[1], &res);
1128
1129 free_evalue_refs(e);
1130 free_evalue_refs(&EP);
1131 *e = res;
1132
1133 value_clear(tmp);
1134 Vector_Free(val);
1135 Vector_Free(periods);
1136 }
1137 } /* evalue_mod2table */
1138
1139 /********************************************************/
1140 /* function in domain */
1141 /* check if the parameters in list_args */
1142 /* verifies the constraints of Domain P */
1143 /********************************************************/
1144 int in_domain(Polyhedron *P, Value *list_args) {
1145
1146 int col,row;
1147 Value v; /* value of the constraint of a row when
1148 parameters are instanciated*/
1149 Value tmp;
1150
1151 value_init(v);
1152 value_init(tmp);
1153
1154 /*P->Constraint constraint matrice of polyhedron P*/
1155 for(row=0;row<P->NbConstraints;row++) {
1156 value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
1157 for(col=1;col<P->Dimension+1;col++) {
1158 value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
1159 value_addto(v,v,tmp);
1160 }
1161 if (value_notzero_p(P->Constraint[row][0])) {
1162
1163 /*if v is not >=0 then this constraint is not respected */
1164 if (value_neg_p(v)) {
1165 next:
1166 value_clear(v);
1167 value_clear(tmp);
1168 return P->next ? in_domain(P->next, list_args) : 0;
1169 }
1170 }
1171 else {
1172
1173 /*if v is not = 0 then this constraint is not respected */
1174 if (value_notzero_p(v))
1175 goto next;
1176 }
1177 }
1178
1179 /*if not return before this point => all
1180 the constraints are respected */
1181 value_clear(v);
1182 value_clear(tmp);
1183 return 1;
1184 } /* in_domain */
1185
1186 /****************************************************/
1187 /* function compute enode */
1188 /* compute the value of enode p with parameters */
1189 /* list "list_args */
1190 /* compute the polynomial or the periodic */
1191 /****************************************************/
1192
1193 static double compute_enode(enode *p, Value *list_args) {
1194
1195 int i;
1196 Value m, param;
1197 double res=0.0;
1198
1199 if (!p)
1200 return(0.);
1201
1202 value_init(m);
1203 value_init(param);
1204
1205 if (p->type == polynomial) {
1206 if (p->size > 1)
1207 value_assign(param,list_args[p->pos-1]);
1208
1209 /* Compute the polynomial using Horner's rule */
1210 for (i=p->size-1;i>0;i--) {
1211 res +=compute_evalue(&p->arr[i],list_args);
1212 res *=VALUE_TO_DOUBLE(param);
1213 }
1214 res +=compute_evalue(&p->arr[0],list_args);
1215 }
1216 else if (p->type == modulo) {
1217 double d = compute_evalue(&p->arr[0], list_args);
1218 if (d > 0)
1219 d += .25;
1220 else
1221 d -= .25;
1222 value_set_double(param, d);
1223 value_set_si(m, p->pos);
1224 mpz_fdiv_r(param, param, m);
1225
1226 /* Compute the polynomial using Horner's rule */
1227 for (i=p->size-1;i>1;i--) {
1228 res +=compute_evalue(&p->arr[i],list_args);
1229 res *=VALUE_TO_DOUBLE(param);
1230 }
1231 res +=compute_evalue(&p->arr[1],list_args);
1232 }
1233 else if (p->type == periodic) {
1234 value_assign(param,list_args[p->pos-1]);
1235
1236 /* Choose the right element of the periodic */
1237 value_absolute(m,param);
1238 value_set_si(param,p->size);
1239 value_modulus(m,m,param);
1240 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
1241 }
1242 else if (p->type == relation) {
1243 if (fabs(compute_evalue(&p->arr[0], list_args)) < 0.5)
1244 res = compute_evalue(&p->arr[1], list_args);
1245 }
1246 value_clear(m);
1247 value_clear(param);
1248 return res;
1249 } /* compute_enode */
1250
1251 /*************************************************/
1252 /* return the value of Ehrhart Polynomial */
1253 /* It returns a double, because since it is */
1254 /* a recursive function, some intermediate value */
1255 /* might not be integral */
1256 /*************************************************/
1257
1258 double compute_evalue(evalue *e,Value *list_args) {
1259
1260 double res;
1261
1262 if (value_notzero_p(e->d)) {
1263 if (value_notone_p(e->d))
1264 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
1265 else
1266 res = VALUE_TO_DOUBLE(e->x.n);
1267 }
1268 else
1269 res = compute_enode(e->x.p,list_args);
1270 return res;
1271 } /* compute_evalue */
1272
1273
1274 /****************************************************/
1275 /* function compute_poly : */
1276 /* Check for the good validity domain */
1277 /* return the number of point in the Polyhedron */
1278 /* in allocated memory */
1279 /* Using the Ehrhart pseudo-polynomial */
1280 /****************************************************/
1281 Value *compute_poly(Enumeration *en,Value *list_args) {
1282
1283 Value *tmp;
1284 /* double d; int i; */
1285
1286 tmp = (Value *) malloc (sizeof(Value));
1287 assert(tmp != NULL);
1288 value_init(*tmp);
1289 value_set_si(*tmp,0);
1290
1291 if(!en)
1292 return(tmp); /* no ehrhart polynomial */
1293 if(en->ValidityDomain) {
1294 if(!en->ValidityDomain->Dimension) { /* no parameters */
1295 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1296 return(tmp);
1297 }
1298 }
1299 else
1300 return(tmp); /* no Validity Domain */
1301 while(en) {
1302 if(in_domain(en->ValidityDomain,list_args)) {
1303
1304 #ifdef EVAL_EHRHART_DEBUG
1305 Print_Domain(stdout,en->ValidityDomain);
1306 print_evalue(stdout,&en->EP);
1307 #endif
1308
1309 /* d = compute_evalue(&en->EP,list_args);
1310 i = d;
1311 printf("(double)%lf = %d\n", d, i ); */
1312 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1313 return(tmp);
1314 }
1315 else
1316 en=en->next;
1317 }
1318 value_set_si(*tmp,0);
1319 return(tmp); /* no compatible domain with the arguments */
1320 } /* compute_poly */
1321
lattice point enumeration library