?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
corrected a bug in reading matrices on MacOS (fgets) and cleaned a bit the code
[polylib.git] / source / kernel / Zpolyhedron.c
blob531d58e6cea3cdbc189dce587f9c6291f21a4b1c
1 /*
2 This file is part of PolyLib.
3
4 PolyLib is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 PolyLib is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with PolyLib. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 #include <polylib/polylib.h>
19 #include <stdlib.h>
20
21 static ZPolyhedron * ZPolyhedronIntersection(ZPolyhedron *, ZPolyhedron *);
22 static ZPolyhedron *ZPolyhedron_Copy(ZPolyhedron *A);
23 static void ZPolyhedron_Free(ZPolyhedron *Zpol);
24 static ZPolyhedron * ZPolyhedronDifference(ZPolyhedron *, ZPolyhedron *);
25 static ZPolyhedron * ZPolyhedronImage(ZPolyhedron *, Matrix *);
26 static ZPolyhedron * ZPolyhedronPreimage(ZPolyhedron *, Matrix *);
27 static ZPolyhedron *AddZPolytoZDomain(ZPolyhedron *A, ZPolyhedron *Head);
28 static void ZPolyhedronPrint(FILE *fp, const char *format, ZPolyhedron *A);
29
30 typedef struct forsimplify {
31 Polyhedron *Pol;
32 LatticeUnion *LatUni;
33 struct forsimplify *next;
34 } ForSimplify;
35
36
37 /*
38 * Returns True if 'Zpol' is empty, otherwise returns False
39 */
40 Bool isEmptyZPolyhedron (ZPolyhedron *Zpol) {
41
42 if(Zpol == NULL)
43 return True;
44 if((isEmptyLattice (Zpol->Lat)) || (emptyQ(Zpol->P)))
45 return True;
46 return False;
47 } /* isEmptyZPolyhedron */
48
49 /*
50 * Given Lattice 'Lat' and a Polyhderon 'Poly', allocate space, and return
51 * the Z-polyhderon corresponding to the image of the polyhderon 'Poly' by the
52 * lattice 'Lat'. If the input lattice 'Lat' is not integeral, it integralises
53 * it, i.e. the lattice of the Z-polyhderon returned is integeral.
54 */
55 ZPolyhedron *ZPolyhedron_Alloc(Lattice *Lat, Polyhedron *Poly) {
56
57 ZPolyhedron *A;
58
59 POL_ENSURE_FACETS(Poly);
60 POL_ENSURE_VERTICES(Poly);
61
62 if(Lat->NbRows != Poly->Dimension+1) {
63 fprintf(stderr,"\nInZPolyAlloc - The Lattice and the Polyhedron");
64 fprintf(stderr," are not compatible to form a ZPolyhedra\n");
65 return NULL;
66 }
67 if((!(isEmptyLattice(Lat))) && (!isfulldim (Lat))) {
68 fprintf(stderr,"\nZPolAlloc: Lattice not Full Dimensional\n");
69 return NULL;
70 }
71 A = (ZPolyhedron *)malloc(sizeof(ZPolyhedron));
72 if (!A) {
73 fprintf(stderr,"ZPolAlloc : Out of Memory\n");
74 return NULL;
75 }
76 A->next = NULL;
77 A->P = Domain_Copy(Poly);
78 A->Lat = Matrix_Copy(Lat);
79
80 if(IsLattice(Lat) == False) {
81 ZPolyhedron *Res;
82
83 Res = IntegraliseLattice (A);
84 ZPolyhedron_Free (A);
85 return Res;
86 }
87 return A;
88 } /* ZPolyhedron_Alloc */
89
90 /*
91 * Free the memory used by the Z-domain 'Head'
92 */
93 void ZDomain_Free (ZPolyhedron *Head) {
94
95 if (Head == NULL)
96 return;
97 if (Head->next != NULL)
98 ZDomain_Free(Head->next);
99 ZPolyhedron_Free(Head);
100 } /* ZDomain_Free */
101
102 /*
103 * Free the memory used by the Z-polyhderon 'Zpol'
104 */
105 static void ZPolyhedron_Free (ZPolyhedron *Zpol) {
106
107 if (Zpol == NULL)
108 return;
109 Matrix_Free((Matrix *) Zpol->Lat);
110 Domain_Free(Zpol->P);
111 free(Zpol);
112 return;
113 } /* ZPolyhderon_Free */
114
115 /*
116 * Return a copy of the Z-domain 'Head'
117 */
118 ZPolyhedron *ZDomain_Copy(ZPolyhedron *Head) {
119
120 ZPolyhedron *Zpol;
121 Zpol = ZPolyhedron_Copy(Head);
122
123 if (Head->next != NULL)
124 Zpol->next = ZDomain_Copy(Head->next);
125 return Zpol;
126 } /* ZDomain_Copy */
127
128 /*
129 * Return a copy of the Z-polyhderon 'A'
130 */
131 static ZPolyhedron *ZPolyhedron_Copy(ZPolyhedron *A) {
132
133 ZPolyhedron *Zpol;
134
135 Zpol = ZPolyhedron_Alloc(A->Lat, A->P);
136 return Zpol;
137 } /* ZPolyhderon_Copy */
138
139 /*
140 * Add the ZPolyhedron 'Zpol' to the Z-domain 'Result' and return a pointer
141 * to the new Z-domain.
142 */
143 static ZPolyhedron *AddZPoly2ZDomain(ZPolyhedron *Zpol, ZPolyhedron *Result) {
144
145 ZPolyhedron *A;
146
147 if (isEmptyZPolyhedron(Zpol))
148 return Result;
149 A = ZPolyhedron_Copy(Zpol);
150 A->next = NULL;
151
152 if (isEmptyZPolyhedron (Result)) {
153 ZDomain_Free (Result);
154 return A;
155 }
156 A->next = Result;
157 return A;
158 } /* AddZPoly2ZDomain */
159
160 /*
161 * Given a Z-polyhderon 'A' and a Z-domain 'Head', return a new Z-domain with
162 * 'A' added to it. If the new Z-polyhedron 'A', is already included in the
163 * Z-domain 'Head', it is not added in the list. Othewise, the function checks
164 * if the new Z-polyhedron 'A' to be added to the Z-domain 'Head' has a common
165 * lattice with some other Z-polyhderon already present in the Z-domain. If it
166 * is so, it takes the union of the underlying polyhdera; domains and returns.
167 * The function tries to make sure that the added Z-polyhedron 'A' is in the
168 * canonical form.
169 */
170 static ZPolyhedron *AddZPolytoZDomain(ZPolyhedron *A, ZPolyhedron *Head) {
171
172 ZPolyhedron *Zpol, *temp, *temp1;
173 Polyhedron *i;
174 Bool Added;
175
176 if ((A == NULL) || (isEmptyZPolyhedron(A)))
177 return Head;
178
179 /* For each "underlying" Pol, find the Cnf and add Zpol in Cnf*/
180 for(i=A->P; i!= NULL; i=i->next) {
181 ZPolyhedron *Z, *Z1;
182 Polyhedron *Image;
183 Matrix *H, *U;
184 Lattice *Lat ;
185
186 Added = False;
187 Image = Domain_Copy(i);
188 Domain_Free(Image->next);
189 Image->next = NULL;
190 Z1 = ZPolyhedron_Alloc(A->Lat,Image);
191 Domain_Free(Image);
192 CanonicalForm(Z1,&Z,&H);
193 ZDomain_Free(Z1);
194 Lat = (Lattice *)Matrix_Alloc(H->NbRows,Z->Lat->NbColumns);
195 Matrix_Product(H,Z->Lat,(Matrix *)Lat);
196 Matrix_Free(H);
197 AffineHermite(Lat,(Lattice **)&H,&U);
198 Image = DomainImage(Z->P,U,MAXNOOFRAYS);
199 ZDomain_Free(Z);
200
201 Zpol=ZPolyhedron_Alloc((Lattice *)H,Image);
202 Domain_Free(Image);
203 Matrix_Free((Matrix *)Lat);
204 Matrix_Free(H);
205 Matrix_Free(U);
206
207 if ((Head == NULL) || (isEmptyZPolyhedron (Head))) {
208 Head = Zpol;
209 continue;
210 }
211 temp1 = temp = Head;
212
213 /* Check if the curr pol is included in the zpol or vice versa. */
214 for(; temp != NULL; temp = temp->next) {
215 if (ZPolyhedronIncludes(Zpol, temp) == True) {
216 ZPolyhedron_Free (Zpol);
217 Added = True;
218 break;
219 }
220 else if (ZPolyhedronIncludes(temp, Zpol) == True) {
221 if (temp == Head) {
222 Zpol->next = temp->next;
223 Head = Zpol;
224 ZPolyhedron_Free (temp);
225 Added = True;
226 break;
227 }
228 temp1->next = Zpol;
229 Zpol->next = temp->next;
230 ZPolyhedron_Free (temp);
231 Added = True;
232 break ;
233 }
234 temp1 = temp ;
235 }
236 if(Added == True)
237 continue ;
238 for(temp = Head; temp != NULL; temp = temp->next) {
239 if(sameLattice(temp->Lat, Zpol->Lat) == True) {
240 Polyhedron *Union;
241
242 Union = DomainUnion (temp->P,Zpol->P,MAXNOOFRAYS);
243 if (!Union)
244 fprintf (stderr,"\n In AddZPolytoZDomain: Out of memory\n");
245 else {
246 Domain_Free(temp->P);
247 temp->P = Union;
248 Added = True;
249 ZPolyhedron_Free(Zpol);
250 }
251 break ;
252 }
253 temp1 = temp;
254 }
255 if (Added == False)
256 temp1->next = Zpol;
257 }
258 return Head ;
259 } /* AddZPolytoZDomain */
260
261 /*
262 * Return the empty Z-polyhedron of dimension 'dimension'
263 */
264 ZPolyhedron *EmptyZPolyhedron(int dimension) {
265
266 ZPolyhedron *Zpol;
267 Lattice *E ;
268 Polyhedron *P;
269
270 #ifdef DOMDEBUG
271 FILE *fp;
272 fp = fopen ("_debug", "a");
273 fprintf (fp, "\nEntered EMPTYZPOLYHEDRON\n");
274 fclose (fp);
275 #endif
276
277 E = EmptyLattice(dimension+1);
278 P = Empty_Polyhedron(dimension);
279 Zpol = ZPolyhedron_Alloc(E,P);
280 Matrix_Free((Matrix *) E);
281 Domain_Free(P);
282 return Zpol;
283 } /* EmptyZPolyhedron */
284
285 /*
286 * Given Z-domains 'A' and 'B', return True if A is included in 'B', otherwise
287 * return False.
288 */
289 Bool ZDomainIncludes(ZPolyhedron *A, ZPolyhedron *B) {
290
291 ZPolyhedron *Diff;
292 Bool ret = False;
293
294 Diff = ZDomainDifference(A,B);
295 if (isEmptyZPolyhedron(Diff))
296 ret = True;
297
298 ZDomain_Free(Diff);
299 return ret;
300 } /* ZDomainIncludes */
301
302 /*
303 * Given Z-polyhedra 'A' and 'B', return True if 'A' is included in 'B',
304 * otherwise return False
305 */
306 Bool ZPolyhedronIncludes(ZPolyhedron *A, ZPolyhedron *B) {
307
308 Polyhedron *Diff = NULL ;
309 Bool retval = False;
310
311 #ifdef DOMDEBUG
312 FILE *fp;
313 fp = fopen("_debug","a");
314 fprintf(fp,"\nEntered ZPOLYHEDRONINCLUDES\n");
315 fclose(fp);
316 #endif
317
318 if (LatticeIncludes(A->Lat, B->Lat) == True) {
319 Polyhedron *ImageA, *ImageB ;
320
321 ImageA = DomainImage(A->P,A->Lat,MAXNOOFRAYS);
322 ImageB = DomainImage(B->P,B->Lat,MAXNOOFRAYS);
323
324 Diff = DomainDifference(ImageA, ImageB, MAXNOOFRAYS);
325 if(emptyQ (Diff))
326 retval = True ;
327
328 Domain_Free (ImageA);
329 Domain_Free (ImageB);
330 Domain_Free (Diff);
331 }
332 return retval;
333 } /* ZPolyhedronIncludes */
334
335 /*
336 * Print the contents of a Z-domain 'A'
337 */
338 void ZDomainPrint(FILE *fp, const char *format, ZPolyhedron *A)
339 {
340 #ifdef DOMDEBUG
341 FILE *fp1;
342 fp1 = fopen("_debug", "a");
343 fprintf(fp1,"\nEntered ZDOMAINPRINT\n");
344 fclose(fp1);
345 #endif
346
347 ZPolyhedronPrint(fp,format,A);
348 if (A->next != NULL) {
349 fprintf(fp,"\nUNIONED with\n");
350 ZDomainPrint(fp,format,A->next);
351 }
352 return;
353 } /* ZDomainPrint */
354
355 /*
356 * Print the contents of a ZPolyhderon 'A'
357 */
358 static void ZPolyhedronPrint (FILE *fp, const char *format, ZPolyhedron *A)
359 {
360 if (A == NULL)
361 return ;
362 fprintf(fp,"\nZPOLYHEDRON: Dimension %d \n",A->Lat->NbRows-1);
363 fprintf(fp, "\nLATTICE: \n");
364 Matrix_Print(fp,format,(Matrix *)A->Lat);
365 Polyhedron_Print(fp,format,A->P);
366 return;
367 } /* ZPolyhedronPrint */
368
369 /*
370 * Return the Z-domain union of the Z-domain 'A' and 'B'. The dimensions of the
371 * Z-domains 'A' and 'B' must be equal. All the Z-polyhedra of the resulting
372 * union are expected to be in Canonical forms.
373 */
374 ZPolyhedron *ZDomainUnion (ZPolyhedron *A, ZPolyhedron *B) {
375
376 ZPolyhedron *Result = NULL, *temp;
377
378 #ifdef DOMDEBUG
379 FILE *fp;
380 fp = fopen("_debug", "a");
381 fprintf(fp,"\nEntered ZDOMAINUNION\n");
382 fclose(fp);
383 #endif
384
385 for(temp = A; temp != NULL; temp = temp->next)
386 Result = AddZPolytoZDomain(temp, Result);
387 for(temp = B; temp != NULL; temp = temp->next )
388 Result = AddZPolytoZDomain(temp, Result);
389 return Result;
390 } /* ZDomainUnion */
391
392 /*
393 * Return the Z-domain intersection of the Z-domains 'A' and 'B'.The dimensions
394 * of domains 'A' and 'B' must be equal.
395 */
396 ZPolyhedron *ZDomainIntersection (ZPolyhedron *A, ZPolyhedron *B) {
397
398 ZPolyhedron *Result = NULL, *tempA = NULL, *tempB = NULL;
399
400 #ifdef DOMDEBUG
401 FILE *fp;
402 fp = fopen("_debug", "a");
403 fprintf(fp,"\nEntered ZDOMAININTERSECTION\n");
404 fclose(fp);
405 #endif
406
407 for(tempA = A; tempA != NULL; tempA = tempA->next)
408 for(tempB = B; tempB != NULL; tempB = tempB->next) {
409 ZPolyhedron *Zpol;
410 Zpol = ZPolyhedronIntersection(tempA, tempB);
411 Result = AddZPolytoZDomain(Zpol, Result );
412 ZPolyhedron_Free (Zpol);
413 }
414 if (Result == NULL)
415 return EmptyZPolyhedron (A->Lat->NbRows-1);
416 return Result;
417 } /* ZDomainIntersection */
418
419 /*
420 * Return the Z-domain difference of the domains 'A' and 'B'. The dimensions of
421 * the Z-domains 'A' and 'B' must be equal. Note that the difference of two
422 * Z-polyhedra is a Union of Z-polyhedra. The algorithms is as given below :-
423 * Algorithm: (Given Z-domains A and B)
424 * Result <-- NULL
425 * for every Z-polyhderon Zpoly of A {
426 * temp <-- Zpoly;
427 * for every Z-polyhderon Z1 of B
428 * temp = temp - Z1;
429 * }
430 * Add temp to Result;
431 * return;
432 */
433 ZPolyhedron *ZDomainDifference(ZPolyhedron *A, ZPolyhedron *B) {
434
435 ZPolyhedron *Result = NULL, *tempA = NULL, *tempB = NULL;
436 ZPolyhedron *templist, *res, *i, *j;
437
438 #ifdef DOMDEBUG
439 FILE *fp;
440 fp = fopen("_debug", "a");
441 fprintf(fp,"\nEntered ZDOMAINDIFFERENCE\n");
442 fclose(fp);
443 #endif
444
445 if (A->Lat->NbRows != B->Lat->NbRows) {
446 fprintf(stderr, "In ZDomainDifference : the input ZDomains ");
447 fprintf(stderr, "do not have compatible dimensions\n");
448 fprintf(stderr, "ZDomainDifference not performed\n");
449 return NULL;
450 }
451
452 for(tempA = A; tempA != NULL; tempA = tempA->next) {
453 ZPolyhedron *temp = NULL;
454 temp = ZPolyhedron_Copy(tempA);
455
456 for(tempB = B; tempB != NULL; tempB = tempB->next) {
457 templist = NULL; res = NULL;
458 for(i = temp; i != NULL; i = i->next) {
459 res = ZPolyhedronDifference(i,tempB);
460 for (j = res; j != NULL; j = j->next )
461 templist = AddZPoly2ZDomain(j,templist);
462 ZDomain_Free(res);
463 }
464 ZDomain_Free (temp);
465 temp = NULL;
466 for(i = templist; i != NULL; i = i->next)
467 temp = AddZPoly2ZDomain(i, temp);
468 ZDomain_Free (templist);
469 }
470 for(i = temp; i != NULL; i = i->next)
471 Result = AddZPolytoZDomain(i, Result);
472 ZDomain_Free(temp);
473 }
474 if (Result==NULL)
475 return (EmptyZPolyhedron(A->Lat->NbRows-1));
476 return Result;
477 } /* ZDomainDifference */
478
479 /*
480 * Return the image of the Z-domain 'A' under the invertible, affine, rational
481 * transformation function 'Func'. The matrix representing the function 'Func'
482 * must be non-singular and the number of rows of the function must be equal
483 * to the number of rows in the matrix representing the lattice of 'A'.
484 * Note:: Image((Z1 U Z2),F) = Image(Z1,F) U Image(Z2 U F).
485 */
486 ZPolyhedron *ZDomainImage (ZPolyhedron *A, Matrix *Func) {
487
488 ZPolyhedron *Result = NULL, *temp;
489
490 #ifdef DOMDEBUG
491 FILE *fp;
492 fp = fopen ("_debug", "a");
493 fprintf (fp, "\nEntered ZDOMAINIMAGE\n");
494 fclose (fp);
495 #endif
496
497 for(temp = A; temp != NULL; temp = temp->next) {
498 ZPolyhedron *Zpol;
499 Zpol = ZPolyhedronImage (temp, Func);
500 Result = AddZPolytoZDomain (Zpol, Result);
501 ZPolyhedron_Free (Zpol);
502 }
503 if(Result == NULL)
504 return EmptyZPolyhedron(A->Lat->NbRows-1);
505 return Result;
506 } /* ZDomainImage */
507
508 /*
509 * Return the preimage of the Z-domain 'A' under the invertible, affine, ratio-
510 * nal transformation 'Func'. The number of rows of the matrix representing
511 * the function 'Func' must be equal to the number of rows of the matrix repr-
512 * senting the lattice of 'A'.
513 */
514 ZPolyhedron *ZDomainPreimage (ZPolyhedron *A, Matrix *Func) {
515
516 ZPolyhedron *Result = NULL, *temp ;
517
518 #ifdef DOMDEBUG
519 FILE *fp;
520 fp = fopen("_debug", "a");
521 fprintf(fp,"\nEntered ZDOMAINPREIMAGE\n");
522 fclose(fp);
523 #endif
524
525 if (A->Lat->NbRows != Func->NbRows) {
526 fprintf(stderr,"\nError : In ZDomainPreimage, ");
527 fprintf(stderr,"Incompatible dimensions of ZPolyhedron ");
528 fprintf(stderr,"and the Function \n");
529 return(EmptyZPolyhedron(Func->NbColumns-1));
530 }
531 for(temp = A; temp != NULL; temp = temp->next) {
532 ZPolyhedron *Zpol;
533 Zpol = ZPolyhedronPreimage(temp, Func);
534 Result = AddZPolytoZDomain(Zpol, Result);
535 ZPolyhedron_Free(Zpol);
536 }
537 if (Result == NULL)
538 return(EmptyZPolyhedron(Func->NbColumns-1));
539 return Result;
540 } /* ZDomainPreimage */
541
542 /*
543 * Return the Z-polyhedron intersection of the Z-polyhedra 'A' and 'B'.
544 * Note: If Z1 = L1 (intersect) P1 and Z2 = L2 (intersect) P2, then
545 * Z1 (intersect) Z2 = (L1 (intersect) L2) (intersect) (P1 (intersect) P2)
546 */
547 static ZPolyhedron *ZPolyhedronIntersection(ZPolyhedron *A, ZPolyhedron *B) {
548
549 ZPolyhedron *Result = NULL;
550 Lattice *LInter;
551 Polyhedron *PInter, *ImageA, *ImageB, *PreImage;
552
553 #ifdef DOMDEBUG
554 FILE *fp;
555 fp = fopen("_debug","a");
556 fprintf(fp,"\nEntered ZPOLYHEDRONINTERSECTION\n");
557 fclose(fp);
558 #endif
559
560 LInter = LatticeIntersection(A->Lat,B->Lat);
561 if(isEmptyLattice(LInter) == True) {
562 ZPolyhedron_Free (Result);
563 Matrix_Free(LInter);
564 return (EmptyZPolyhedron(A->Lat->NbRows-1));
565 }
566 ImageA = DomainImage(A->P,A->Lat,MAXNOOFRAYS);
567 ImageB = DomainImage(B->P,B->Lat,MAXNOOFRAYS);
568 PInter = DomainIntersection(ImageA,ImageB,MAXNOOFRAYS);
569 if (emptyQ(PInter))
570 Result = EmptyZPolyhedron(LInter->NbRows-1);
571 else {
572 PreImage = DomainPreimage(PInter,(Matrix *)LInter,MAXNOOFRAYS);
573 Result = ZPolyhedron_Alloc(LInter, PreImage);
574 Domain_Free (PreImage);
575 }
576 Matrix_Free(LInter);
577 Domain_Free(PInter);
578 Domain_Free(ImageA);
579 Domain_Free(ImageB);
580 return Result ;
581 } /* ZPolyhedronIntersection */
582
583 /*
584 * Return the difference of the two Z-polyhedra 'A' and 'B'. Below is the
585 * procedure to find the difference of 'A' and 'B' :-
586 * Procedure:
587 * Let A = L1 (intersect) P1' and B = L2 (intersect) P2' where
588 * (P1' = DomImage(P1,L1) and P2' = DomImage(P2,L2)). Then
589 * A-B = L1 (intersect) (P1'-P2') Union
590 * (L1-L2) (intersect) (P1' (intersect) P2')
591 */
592 static ZPolyhedron *ZPolyhedronDifference(ZPolyhedron *A, ZPolyhedron *B) {
593
594 ZPolyhedron *Result = NULL ;
595 LatticeUnion *LatDiff, *temp;
596 Polyhedron *DomDiff, *DomInter, *PreImage, *ImageA, *ImageB;
597 Bool flag = False;
598
599 #ifdef DOMDEBUG
600 FILE *fp;
601 fp = fopen("_debug", "a");
602 fprintf(fp,"\nEntered ZPOLYHEDRONDIFFERENCE\n");
603 fclose(fp);
604 #endif
605
606 if(isEmptyZPolyhedron (A))
607 return NULL;
608 if(isEmptyZPolyhedron (B)) {
609 Result = ZDomain_Copy (A);
610 return Result;
611 }
612 ImageA = DomainImage(A->P,(Matrix *)A->Lat,MAXNOOFRAYS);
613 ImageB = DomainImage(B->P,(Matrix *)B->Lat,MAXNOOFRAYS);
614 DomDiff = DomainDifference(ImageA,ImageB,MAXNOOFRAYS);
615 if (emptyQ (DomDiff))
616 flag = True;
617 else {
618 ZPolyhedron *Z;
619 PreImage = DomainPreimage(DomDiff,A->Lat,MAXNOOFRAYS);
620 Z = ZPolyhedron_Alloc(A->Lat,PreImage);
621 Result = AddZPolytoZDomain(Z,Result);
622 }
623 if (flag == True) /* DomDiff = NULL; DomInter = A */
624 DomInter = Domain_Copy(ImageA);
625 else {
626 DomInter = DomainIntersection(ImageA,ImageB,MAXNOOFRAYS);
627 if (emptyQ(DomInter)) {
628 if (flag == True)
629 return (EmptyZPolyhedron(A->Lat->NbRows-1));
630 else
631 return Result;
632 }
633 }
634 LatDiff = LatticeDifference(A->Lat, B->Lat);
635 if(LatDiff == NULL)
636 if(flag == True )
637 return(EmptyZPolyhedron (A->Lat->NbRows-1));
638
639 while (LatDiff != NULL) {
640 ZPolyhedron *tempZ = NULL;
641
642 PreImage = DomainPreimage(DomInter, LatDiff->M, MAXNOOFRAYS);
643 tempZ = ZPolyhedron_Alloc(LatDiff->M, PreImage);
644 Domain_Free(PreImage);
645 Result = AddZPoly2ZDomain(tempZ,Result);
646 ZPolyhedron_Free(tempZ);
647 temp = LatDiff;
648 LatDiff = LatDiff->next;
649 Matrix_Free ((Matrix *) temp->M);
650 free (temp);
651 }
652 Domain_Free (DomInter);
653 Domain_Free (DomDiff);
654 return Result;
655 } /* ZPolyhedronDifference */
656
657 /*
658 * Return the image of the Z-polyhedron 'ZPol' under the invertible, affine,
659 * rational transformation function 'Func'. The matrix representing the funct-
660 * ion must be non-singular and the number of rows of the function must be
661 * equal to the number of rows in the matrix representing the lattice of 'ZPol'
662 * Algorithm:
663 * 1) Let ZPol = L (intersect) Q
664 * 2) L1 = LatticeImage(L,F)
665 * 3) Q1 = DomainImage(Q,F)
666 * 4) Z1 = L1(Inverse(L1)*Q1)
667 * 5) Return Z1
668 */
669 static ZPolyhedron *ZPolyhedronImage(ZPolyhedron *ZPol,Matrix *Func) {
670
671 ZPolyhedron *Result = NULL ;
672 Matrix *LatIm ;
673 Polyhedron *Pol, *PolImage ;
674
675 #ifdef DOMDEBUG
676 FILE *fp;
677 fp = fopen("_debug", "a");
678 fprintf(fp,"\nEntered ZPOLYHEDRONIMAGE\n");
679 fclose(fp);
680 #endif
681
682 if ((Func->NbRows != ZPol->Lat->NbRows) || (Func->NbColumns != ZPol->Lat->NbColumns)) {
683 fprintf (stderr, "In ZPolImage - The Function, is not compatible with the ZPolyhedron\n");
684 return NULL;
685 }
686 LatIm = LatticeImage(ZPol->Lat,Func);
687 if (isEmptyLattice(LatIm)) {
688 Matrix_Free(LatIm);
689 return NULL;
690 }
691 Pol = DomainImage(ZPol->P,ZPol->Lat,MAXNOOFRAYS);
692 PolImage = DomainImage(Pol,Func,MAXNOOFRAYS);
693 Domain_Free(Pol);
694 if(emptyQ(PolImage)) {
695 Matrix_Free (LatIm);
696 Domain_Free (PolImage);
697 return NULL;
698 }
699 Pol = DomainPreimage(PolImage,LatIm,MAXNOOFRAYS);
700 Result = ZPolyhedron_Alloc(LatIm,Pol);
701 Domain_Free(Pol);
702 Domain_Free(PolImage);
703 Matrix_Free(LatIm);
704 return Result;
705 } /* ZPolyhedronImage */
706
707 /*
708 * Return the preimage of the Z-polyhedron 'Zpol' under an affine transformati-
709 * on function 'G'. The number of rows of matrix representing the function 'G',
710 * must be equal to the number of rows of the matrix representing the lattice
711 * of Z1.
712 * Algorithm:
713 * 1) Let Zpol = L (intersect) Q
714 * 2) L1 =LatticePreimage(L,F);
715 * 3) Q1 = DomainPreimage(Q,F);
716 * 4) Z1 = L1(Inverse(L1)*Q1);
717 * 5) Return Z1
718 */
719 static ZPolyhedron *ZPolyhedronPreimage(ZPolyhedron *Zpol, Matrix *G) {
720
721 Lattice *Latpreim;
722 Polyhedron *Qprime, *Q, *Polpreim;
723 ZPolyhedron *Result;
724
725 #ifdef DOMDEBUG
726 FILE *fp;
727 fp = fopen("_debug","a");
728 fprintf(fp,"\nEntered ZPOLYHEDRONPREIMAGE\n");
729 fclose(fp);
730 #endif
731
732 if(G->NbRows != Zpol->Lat->NbRows) {
733 fprintf(stderr,"\nIn ZPolyhedronPreimage: Error, The dimensions of the ");
734 fprintf(stderr,"function are not compatible with that of the Zpolyhedron");
735 return EmptyZPolyhedron(G->NbColumns-1);
736 }
737 Q = DomainImage(Zpol->P,Zpol->Lat,MAXNOOFRAYS);
738 Polpreim = DomainPreimage(Q,G,MAXNOOFRAYS);
739 if (emptyQ(Polpreim))
740 Result = NULL;
741 else {
742 Latpreim = LatticePreimage(Zpol->Lat,G);
743 if(isEmptyLattice(Latpreim))
744 Result = NULL;
745 else {
746 Qprime = DomainPreimage(Polpreim, Latpreim, MAXNOOFRAYS);
747 Result = ZPolyhedron_Alloc(Latpreim, Qprime);
748 Domain_Free(Qprime);
749 }
750 Matrix_Free(Latpreim);
751 }
752 Domain_Free(Q);
753 return Result;
754 } /* ZPolyhedronPreimage */
755
756 /*
757 * Return the Z-polyhderon 'Zpol' in canonical form: 'Result' (for the Z-poly-
758 * hedron in canonical form) and Basis 'Basis' (for the basis with respect to
759 * which 'Result' is in canonical form.
760 */
761 void CanonicalForm(ZPolyhedron *Zpol,ZPolyhedron **Result,Matrix **Basis) {
762
763 Matrix *B1 = NULL, *B2=NULL, *T1 , *B2inv;
764 int i, l1, l2;
765 Value tmp;
766 Polyhedron *Image, *ImageP;
767 Matrix *H, *U, *temp, *Hprime, *Uprime, *T2;
768
769 #ifdef DOMDEBUG
770 FILE *fp;
771 fp = fopen("_debug", "a");
772 fprintf(fp,"\nEntered CANONICALFORM\n");
773 fclose(fp);
774 #endif
775
776 if(isEmptyZPolyhedron (Zpol)) {
777 Basis[0] = Identity(Zpol->Lat->NbRows);
778 Result[0] = ZDomain_Copy (Zpol);
779 return ;
780 }
781 value_init(tmp);
782 l1 = FindHermiteBasisofDomain(Zpol->P,&B1);
783 Image = DomainImage (Zpol->P,(Matrix *)Zpol->Lat,MAXNOOFRAYS);
784 l2 = FindHermiteBasisofDomain(Image,&B2);
785
786 if (l1 != l2)
787 fprintf(stderr,"In CNF : Something wrong with the Input Zpolyhedra \n");
788
789 B2inv = Matrix_Alloc(B2->NbRows, B2->NbColumns);
790 temp = Matrix_Copy(B2);
791 Matrix_Inverse(temp,B2inv);
792 Matrix_Free(temp);
793
794 temp = Matrix_Alloc(B2inv->NbRows,Zpol->Lat->NbColumns);
795 T1 = Matrix_Alloc(temp->NbRows,B1->NbColumns);
796 Matrix_Product(B2inv,(Matrix *)Zpol->Lat,temp);
797 Matrix_Product(temp,B1,T1);
798 Matrix_Free(temp);
799
800 T2 = ChangeLatticeDimension(T1,l1);
801 temp = ChangeLatticeDimension(T2,T2->NbRows+1);
802
803 /* Adding the affine part */
804 for(i = 0; i < l1; i ++)
805 value_assign(temp->p[i][temp->NbColumns-1],T1->p[i][T1->NbColumns-1]);
806
807 AffineHermite(temp,&H,&U);
808 Hprime = ChangeLatticeDimension(H,Zpol->Lat->NbRows);
809
810 /* Exchanging the Affine part */
811 for(i = 0; i < l1; i ++) {
812 value_assign(tmp,Hprime->p[i][Hprime->NbColumns-1]);
813 value_assign(Hprime->p[i][Hprime->NbColumns-1],Hprime->p[i][H->NbColumns-1]);
814 value_assign(Hprime->p[i][H->NbColumns-1],tmp);
815 }
816 Uprime = ChangeLatticeDimension(U,Zpol->Lat->NbRows);
817
818 /* Exchanging the Affine part */
819 for (i = 0;i < l1; i++) {
820 value_assign(tmp,Uprime->p[i][Uprime->NbColumns-1]);
821 value_assign(Uprime->p[i][Uprime->NbColumns-1],Uprime->p[i][U->NbColumns-1]);
822 value_assign(Uprime->p[i][U->NbColumns-1],tmp);
823 }
824 Polyhedron_Free (Image);
825 Matrix_Free (B2inv);
826 B2inv = Matrix_Alloc(B1->NbRows, B1->NbColumns);
827 Matrix_Inverse(B1,B2inv);
828 ImageP = DomainImage(Zpol->P, B2inv, MAXNOOFRAYS);
829 Matrix_Free(B2inv);
830 Image = DomainImage(ImageP, Uprime, MAXNOOFRAYS);
831 Domain_Free(ImageP);
832 Result[0] = ZPolyhedron_Alloc(Hprime, Image);
833 Basis[0] = Matrix_Copy(B2);
834
835 /* Free the variables */
836 Polyhedron_Free (Image);
837 Matrix_Free (B1);
838 Matrix_Free (B2);
839 Matrix_Free (temp);
840 Matrix_Free (T1);
841 Matrix_Free (T2);
842 Matrix_Free (H);
843 Matrix_Free (U);
844 Matrix_Free (Hprime);
845 Matrix_Free (Uprime);
846 value_clear(tmp);
847 return;
848 } /* CanonicalForm */
849
850 /*
851 * Given a Z-polyhedron 'A' in which the Lattice is not integral, return the
852 * Z-polyhedron which contains all the integral points in the input lattice.
853 */
854 ZPolyhedron *IntegraliseLattice(ZPolyhedron *A) {
855
856 ZPolyhedron *Result;
857 Lattice *M = NULL, *Id;
858 Polyhedron *Im = NULL, *Preim = NULL;
859
860 #ifdef DOMDEBUG
861 FILE *fp;
862 fp = fopen("_debug","a");
863 fprintf(fp,"\nEntered INTEGRALISELATTICE\n");
864 fclose(fp);
865 #endif
866
867 Im = DomainImage(A->P,A->Lat,MAXNOOFRAYS);
868 Id = Identity(A->Lat->NbRows);
869 M = LatticeImage(Id, A->Lat);
870 if (isEmptyLattice(M))
871 Result = EmptyZPolyhedron(A->Lat->NbRows-1);
872 else {
873 Preim = DomainPreimage(Im,M,MAXNOOFRAYS);
874 Result = ZPolyhedron_Alloc(M,Preim);
875 }
876 Matrix_Free(M);
877 Domain_Free(Im);
878 Domain_Free(Preim);
879 return Result;
880 } /* IntegraliseLattice */
881
882 /*
883 * Return the simplified representation of the Z-domain 'ZDom'. It attempts to
884 * convexize unions of polyhedra when they correspond to the same lattices and
885 * to simplify union of lattices when they correspond to the same polyhdera.
886 */
887 ZPolyhedron *ZDomainSimplify(ZPolyhedron *ZDom) {
888
889 ZPolyhedron *Ztmp, *Result;
890 ForSimplify *Head, *Prev, *Curr;
891 ZPolyhedron *ZDomHead, *Emp;
892
893 if (ZDom == NULL) {
894 fprintf(stderr,"\nError in ZDomainSimplify - ZDomHead = NULL\n");
895 return NULL;
896 }
897 if (ZDom->next == NULL)
898 return (ZPolyhedron_Copy (ZDom));
899 Emp = EmptyZPolyhedron(ZDom->Lat->NbRows-1);
900 ZDomHead = ZDomainUnion(ZDom, Emp);
901 ZPolyhedron_Free(Emp);
902 Head = NULL;
903 Ztmp = ZDomHead;
904 do {
905 Polyhedron *Img;
906 Img = DomainImage(Ztmp->P,Ztmp->Lat,MAXNOOFRAYS);
907 for(Curr = Head; Curr != NULL; Curr = Curr->next) {
908 Polyhedron *Diff1;
909 Bool flag = False;
910
911 Diff1 = DomainDifference(Img,Curr->Pol,MAXNOOFRAYS);
912 if (emptyQ(Diff1)) {
913 Polyhedron *Diff2;
914
915 Diff2 = DomainDifference(Curr->Pol,Img,MAXNOOFRAYS);
916 if (emptyQ(Diff2))
917 flag = True;
918 Domain_Free(Diff2);
919 }
920 Domain_Free (Diff1);
921 if (flag == True) {
922 LatticeUnion *temp;
923
924 temp = (LatticeUnion *)malloc(sizeof(LatticeUnion));
925 temp->M = (Lattice *)Matrix_Copy((Matrix *)Ztmp->Lat);
926 temp->next = Curr->LatUni;
927 Curr->LatUni = temp;
928 break;
929 }
930 }
931 if(Curr == NULL) {
932 Curr = (ForSimplify *)malloc(sizeof(ForSimplify));
933 Curr->Pol = Domain_Copy(Img);
934 Curr->LatUni = (LatticeUnion *)malloc(sizeof(LatticeUnion));
935 Curr->LatUni->M = (Lattice *)Matrix_Copy((Matrix *)Ztmp->Lat);
936 Curr->LatUni->next = NULL;
937 Curr->next = Head;
938 Head = Curr;
939 }
940 Domain_Free (Img);
941 Ztmp = Ztmp->next;
942 } while(Ztmp != NULL);
943
944 for (Curr = Head; Curr != NULL; Curr = Curr->next)
945 Curr->LatUni = LatticeSimplify(Curr->LatUni);
946 Result = NULL;
947 for(Curr = Head; Curr != NULL; Curr = Curr->next) {
948 LatticeUnion *L;
949 for(L = Curr->LatUni; L != NULL; L = L->next) {
950 Polyhedron *Preim;
951 ZPolyhedron *Zpol;
952
953 Preim = DomainPreimage(Curr->Pol,L->M,MAXNOOFRAYS);
954 Zpol = ZPolyhedron_Alloc(L->M, Preim);
955 Zpol->next = Result;
956 Result = Zpol;
957 Domain_Free(Preim);
958 }
959 }
960 Curr = Head;
961 while (Curr != NULL) {
962 Prev = Curr;
963 Curr = Curr->next;
964 LatticeUnion_Free(Prev->LatUni);
965 Domain_Free(Prev->Pol);
966 free(Prev);
967 }
968 return Result;
969 } /* ZDomainSimplify */
970
971 ZPolyhedron *SplitZpolyhedron(ZPolyhedron *ZPol,Lattice *B) {
972
973 Lattice *Intersection = NULL;
974 Lattice *B1 = NULL, *B2 = NULL, *newB1 = NULL, *newB2 = NULL;
975 Matrix *U = NULL,*M1 = NULL, *M2 = NULL, *M1Inverse = NULL,*MtProduct = NULL;
976 Matrix *Vinv, *V , *temp, *DiagMatrix ;
977 Matrix *H , *U1 , *X, *Y ;
978 ZPolyhedron *zpnew, *Result;
979 LatticeUnion *Head = NULL, *tempHead = NULL;
980 int i;
981 Value k;
982
983 #ifdef DOMDEBUG
984 FILE *fp;
985 fp = fopen("_debug", "a");
986 fprintf(fp,"\nEntered SplitZpolyhedron \n");
987 fclose(fp);
988 #endif
989
990
991 if (B->NbRows != B->NbColumns) {
992 fprintf(stderr,"\n SplitZpolyhedron : The Input Matrix B is not a proper Lattice \n");
993 return NULL;
994 }
995
996 if (ZPol->Lat->NbRows != B->NbRows) {
997 fprintf(stderr,"\nSplitZpolyhedron : The Lattice in Zpolyhedron and B have ");
998 fprintf(stderr,"incompatible dimensions \n");
999 return NULL;
1000 }
1001
1002 if (isinHnf (ZPol->Lat) != True) {
1003 AffineHermite(ZPol->Lat,&H,&U1);
1004 X = Matrix_Copy(H);
1005 Matrix_Free(U1);
1006 Matrix_Free(H);
1007 }
1008 else
1009 X = Matrix_Copy(ZPol->Lat);
1010
1011 if (isinHnf(B) != True) {
1012 AffineHermite(B,&H,&U1);
1013 Y = Matrix_Copy(H);
1014 Matrix_Free(H);
1015 Matrix_Free(U1);
1016 }
1017 else
1018 Y = Matrix_Copy(B);
1019 if (isEmptyLattice(X)) {
1020 return NULL;
1021 }
1022
1023 Head=Lattice2LatticeUnion(X,Y);
1024
1025 /* If the spliting operation can't be done the result is the original Zplyhedron. */
1026
1027 if (Head == NULL) {
1028 Matrix_Free(X);
1029 Matrix_Free(Y);
1030 return ZPol;
1031 }
1032
1033
1034 Result=NULL;
1035
1036 if (Head)
1037 while(Head)
1038 {
1039 tempHead = Head;
1040 Head = Head->next;
1041 zpnew=ZPolyhedron_Alloc(tempHead->M,ZPol->P);
1042 Result=AddZPoly2ZDomain(zpnew,Result);
1043 ZPolyhedron_Free(zpnew);
1044 tempHead->next = NULL;
1045 free(tempHead);
1046 }
1047
1048 return Result;
1049 }
1050
1051
1052
Polyhedral library