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 /***********************************************************************/
19 /* Parametrized polyhedra V4.20 */
20 /* copyright 1995-2000 Vincent Loechner */
21 /* copyright 1996-1997, Doran Wilde */
22 /* Permission is granted to copy, use, and distribute */
23 /* for any commercial or noncommercial purpose under the terms */
24 /* of the GNU General Public license, version 2, June 1991 */
25 /* (see file : LICENSING). */
26 /***********************************************************************/
27
28 /********************* -----------USER #DEFS-------- ***********************/
29 /* These are mainly for debug purposes. You shouldn't need to change */
30 /* anything for daily usage... */
31 /***************************************************************************/
32
33 /* you may define each macro independently */
34 /* #define DEBUGPP */
35 /* #define DEBUGPP3 */ /* initialization of domain, context, ... */
36 /* #define DEBUGPP31 */ /* even more init-domains */
37 /* #define DEBUGPP32 */ /* even even more... (Elim_Columns) */
38 /* #define DEBUGPP4 */ /* m-faces scan */
39 /* #define DEBUGPP41 */ /* inverse Di in scan */
40 /* #define DEBUGPP5 */ /* Compute_PDomains */
41 /********************* ---------END USER #DEFS------ ***********************/
42
43 #include <stdio.h>
44 #include <stdlib.h>
45 #include <string.h>
46 #include <assert.h>
47 #ifdef DEBUGPP
48 #include <time.h>
49 #endif
50
51 #include <polylib/polylib.h>
52
53 static void traite_m_face(Polyhedron *, unsigned int *, unsigned int *);
54 static void scan_m_face(int,int,Polyhedron *,unsigned int *);
55
56 /*
57 * Return the intersection of two polyhedral domains 'Pol1' and 'Pol2' such
58 * that if the intersection is a polyhedron of lower dimension (a degenerate
59 * polyhedron) than the operands, it is discarded from the resulting polyhedra
60 * list.
61 */
PDomainIntersection(Polyhedron * Pol1,Polyhedron * Pol2,unsigned NbMaxRays)62 Polyhedron *PDomainIntersection(Polyhedron *Pol1,Polyhedron *Pol2,unsigned NbMaxRays) {
63
64 Polyhedron *p1, *p2, *p3, *d;
65
66 if (!Pol1 || !Pol2) return (Polyhedron*) 0;
67 if((Pol1->Dimension != Pol2->Dimension) || (Pol1->NbEq != Pol2->NbEq)) {
68 fprintf(stderr,
69 "? PDomainIntersection: operation on different dimensions\n");
70 return (Polyhedron*) 0;
71 }
72
73 POL_ENSURE_FACETS(Pol1);
74 POL_ENSURE_VERTICES(Pol1);
75 POL_ENSURE_FACETS(Pol2);
76 POL_ENSURE_VERTICES(Pol2);
77
78 d = (Polyhedron *)0;
79 for (p1=Pol1; p1; p1=p1->next) {
80 for (p2=Pol2; p2; p2=p2->next) {
81 p3 = AddConstraints(p2->Constraint[0],
82 p2->NbConstraints,p1,NbMaxRays);
83 if (!p3) continue;
84
85 /* If the new polyhedron 'p3' has lower dimension, discard it */
86 if (p3->NbEq!=Pol1->NbEq)
87 Polyhedron_Free(p3) ;
88
89 /* Otherwise add it to the new polyhderal domain 'd'. */
90 else
91 d = AddPolyToDomain(p3,d);
92 }
93 }
94 return d;
95 } /* PDomainIntersection */
96
97 /*
98 * Given polyhderal domains 'Pol1' and 'Pol2', return the difference of the
99 * two domains with a modification that the resulting polyhedra in the new
100 * domain don't have a 1 unit space around cut and the degenerate results
101 * (of smaller dimension) are discarded.
102 */
PDomainDifference(Polyhedron * Pol1,Polyhedron * Pol2,unsigned NbMaxRays)103 Polyhedron *PDomainDifference(Polyhedron *Pol1,Polyhedron *Pol2,unsigned NbMaxRays) {
104
105 Polyhedron *p1, *p2, *p3, *d;
106 int i;
107
108 if (!Pol1 || !Pol2)
109 return (Polyhedron*) 0;
110 if((Pol1->Dimension != Pol2->Dimension) || (Pol1->NbEq != Pol2->NbEq)) {
111 fprintf(stderr,
112 "? PDomainDifference: operation on different dimensions\n");
113 return (Polyhedron*) 0;
114 }
115
116 POL_ENSURE_FACETS(Pol1);
117 POL_ENSURE_VERTICES(Pol1);
118 POL_ENSURE_FACETS(Pol2);
119 POL_ENSURE_VERTICES(Pol2);
120
121 d = (Polyhedron *)0;
122 for (p2=Pol2; p2; p2=p2->next) {
123 for (p1=Pol1; p1; p1=p1->next) {
124 for (i=0; i<p2->NbConstraints; i++) {
125
126 /* Add the constraint (-p2->Constraint[i]) >= 0 in 'p1' */
127 p3 = SubConstraint(p2->Constraint[i],p1,NbMaxRays,2);
128 if (!p3) continue;
129
130 /* If the new polyhedron 'p3' is empty or is a polyhedron of lower */
131 /* dimension, discard it. */
132 if (emptyQ(p3) || p3->NbEq!=Pol1->NbEq)
133 Polyhedron_Free(p3);
134
135 /* Otherwise add 'p3' to the new polyhderal domain 'd' */
136 else
137 d = AddPolyToDomain(p3,d);
138 }
139 }
140 if (p2 != Pol2)
141 Domain_Free(Pol1);
142 Pol1 = d;
143 d = (Polyhedron *)0;
144 }
145 return Pol1;
146 } /* PDomainDifference */
147
148 /*
149 * Return 1 if matrix 'Mat' is full column ranked, otherwise return 0.
150 */
TestRank(Matrix * Mat)151 static int TestRank(Matrix *Mat) {
152
153 int i,j,k;
154 Value m1,m2,m3,gcd,tmp;
155
156 /* Initialize all the 'Value' variables */
157 value_init(m1); value_init(m2);
158 value_init(m3); value_init(gcd); value_init(tmp);
159
160 for(k=0;k<Mat->NbColumns;++k) {
161
162 /* If the digonal entry (k,k) is zero, search down the column(k) */
163 /* starting from row(k) to find a non-zero entry */
164 if(value_zero_p(Mat->p[k][k])) {
165 for(j=k+1;j<Mat->NbRows;++j) {
166
167 /* If a non-zero entry (j,k) is found */
168 if(value_notzero_p(Mat->p[j][k])) {
169
170 /* Exchange row(k) and row(j) */
171 for(i=k;i<Mat->NbColumns;++i) {
172 value_assign(tmp,Mat->p[j][i]);
173 value_assign(Mat->p[j][i],Mat->p[k][i]);
174 value_assign(Mat->p[k][i],tmp);
175 }
176 break;
177 }
178 }
179
180 /* If no non-zero entry is found then the matrix 'Mat' is not full */
181 /* ranked. Return zero. */
182 if(j>=Mat->NbRows) {
183
184 /* Clear all the 'Value' variables */
185 value_clear(m1); value_clear(m2);
186 value_clear(m3); value_clear(gcd); value_clear(tmp);
187 return 0;
188 }
189 }
190
191 /* Now Mat[k][k] is the pivot element */
192 for(j=k+1;j<Mat->NbRows;++j) {
193
194 /* Make every other entry (below row(k)) in column(k) zero */
195 value_gcd(gcd, Mat->p[j][k], Mat->p[k][k]);
196 for(i=k+1;i<Mat->NbColumns;++i) {
197
198 /* pour tous les indices i > k */
199 value_multiply(m1,Mat->p[j][i],Mat->p[k][k]);
200 value_multiply(m2,Mat->p[j][k],Mat->p[k][i]);
201 value_subtract(m3,m1,m2);
202 value_division(Mat->p[j][i],m3,gcd);
203 }
204 }
205 }
206
207 /* Clear all the 'Value' variables */
208 value_clear(m1); value_clear(m2);
209 value_clear(m3); value_clear(gcd); value_clear(tmp);
210
211 /* The matrix 'Mat' is full ranked, return 1 */
212 return 1;
213 } /* TestRank */
214
215 /*
216 * The Saturation matrix is defined to be an integer (int type) matrix. It is
217 * a boolean matrix which has a row for every constraint and a column for
218 * every line or ray. The bits in the binary format of each integer in the
219 * saturation matrix stores the information whether the corresponding constr-
220 * aint is saturated by ray(line) or not.
221 */
222 typedef struct {
223 unsigned int NbRows;
224 unsigned int NbColumns;
225 unsigned int **p;
226 unsigned int *p_init;
227 } SatMatrix;
228
SMAlloc(int rows,int cols)229 static SatMatrix *SMAlloc(int rows,int cols) {
230
231 unsigned int **q, *p;
232 int i;
233
234 SatMatrix *result = (SatMatrix *)malloc(sizeof(SatMatrix));
235 assert (result != NULL);
236
237 result->NbRows = rows;
238 result->NbColumns = cols;
239 result->p = q = (unsigned int **)malloc(rows * sizeof(unsigned int *));
240 assert (result->p != NULL);
241 result->p_init = p = (unsigned int *)malloc(rows * cols * sizeof(unsigned int));
242 assert (result->p_init != NULL);
243
244 for (i=0;i<rows;i++) {
245 *q++ = p;
246 p += cols;
247 }
248
249 return result;
250 } /* SMAlloc */
251
SMPrint(SatMatrix * matrix)252 static void SMPrint (SatMatrix *matrix) {
253
254 unsigned int *p;
255 int i, j;
256 unsigned NbRows, NbColumns;
257
258 fprintf(stderr,"%d %d\n",NbRows=matrix->NbRows, NbColumns=matrix->NbColumns);
259 for (i=0;i<NbRows;i++) {
260 p = *(matrix->p+i);
261 for (j=0;j<NbColumns;j++)
262 fprintf(stderr, " %10X ", *p++);
263 fprintf(stderr, "\n");
264 }
265 } /* SMPrint */
266
267
SMFree(SatMatrix * matrix)268 static void SMFree (SatMatrix *matrix) {
269
270 free ((char *) matrix->p_init);
271 free ((char *) matrix->p);
272 free ((char *) matrix);
273 return;
274 } /* SMFree */
275
276 /* -------------------------------------------------------------------------
277 * Shared Global Variables:
278 * Used by procedures: Find_m_face, scan_m_face, Poly2Sat, traite_m_face,
279 * count_sat
280 * -------------------------------------------------------------------------
281 */
282 static int m; /* number of parameters */
283 static int m_dim; /* dimension of m-face */
284 static int n; /* dimension (not including parameters) */
285 static int ws; /* Working Space size */
286 static int nr; /* (NbRays-1)/32 + 1 */
287
288 static Polyhedron *CEqualities;/* Equalities in the context */
289 static SatMatrix *Sat; /* Saturation Matrix (row=constraint, col=ray)*/
290 static unsigned int *egalite; /* Bool vector marking constraints in m-face */
291 static Matrix *Xi, *Pi; /* Xi and Pi */
292 static Matrix *PiTest; /* Matrix used to test if Pi is full ranked? */
293 static Matrix *CTest;
294 static Matrix *PiInv; /* Matrix inverse Pi, with the last col of */
295 /* each line = denominator of the line */
296 static Matrix *RaysDi; /* Constraint matrix for computing Di */
297
298 static int KD; /* Flag : keep the full domains in memory ? */
299 /* 1 = yes; 0 = no, keep constraints only */
300
301 static int nbPV; /* The number of parameterized vertices */
302 static Param_Vertices *PV_Result; /* List of parameterized vertices */
303 static Param_Domain *PDomains; /* List of domains. */
304
305 #ifdef DEBUGPP
306 static int nbfaces;
307 #endif
308
309 /*
310 * Add the constraints from the context polyhedron 'CEqualities' to the
311 * constraints of polyhedra in the polyhedral domain 'D' and return the new
312 * polyhedral domain. Polyhedral domain 'D' is subsequently deleted from memory
313 */
Add_CEqualities(Polyhedron * D)314 static Polyhedron *Add_CEqualities(Polyhedron *D) {
315
316 Polyhedron *d,*r,*tmp;
317
318 if(!CEqualities)
319 return D;
320 else {
321 if(!D || emptyQ(D)) {
322 if(D)
323 Domain_Free(D);
324 return(Polyhedron_Copy(CEqualities));
325 }
326 r = AddConstraints(D->Constraint[0],D->NbConstraints,
327 CEqualities,ws);
328 tmp = r;
329 for(d=D->next;d;d=d->next) {
330 tmp->next = AddConstraints(d->Constraint[0],d->NbConstraints,
331 CEqualities,ws);
332 tmp = tmp->next;
333 }
334 Domain_Free(D);
335 return(r);
336 }
337 } /* Add_CEqualities */
338
339 #define INT_BITS (sizeof(unsigned) * 8)
340
int_array2bit_vector(unsigned int * array,int n)341 unsigned int *int_array2bit_vector(unsigned int *array, int n)
342 {
343 int i, ix;
344 unsigned bx;
345 int words = (n+INT_BITS-1)/INT_BITS;
346 unsigned int *bv = (unsigned int *)calloc(words, sizeof(unsigned));
347 assert(bv);
348 for (i = 0, ix = 0, bx = MSB; i < n; ++i) {
349 if (array[i])
350 bv[ix] |= bx;
351 NEXT(ix, bx);
352 }
353 return bv;
354 }
355
356 /*----------------------------------------------------------------------*/
357 /* traite_m_face */
358 /* Given an m-face, compute the parameterized vertex */
359 /* D - The entire domain */
360 /* mf - Bit vector marking the lines/rays in the m-face */
361 /* egalite - boolean vector marking saturated constraints in m-face */
362 /*----------------------------------------------------------------------*/
traite_m_face(Polyhedron * D,unsigned int * mf,unsigned int * egalite)363 static void traite_m_face(Polyhedron *D, unsigned int *mf,
364 unsigned int *egalite)
365 {
366 Matrix *Si; /* Solution */
367 Polyhedron *PDi; /* polyhedron Di */
368 Param_Vertices *PV;
369 int j,k,c,r;
370 unsigned kx, bx;
371
372 #ifdef DEBUGPP
373 ++nbfaces;
374 #endif
375
376 /* Extract Xi, Pi, and RaysDi from D */
377 RaysDi->NbRows = 0;
378 for(k=0,c=0,kx=0,bx=MSB;k<D->NbRays;++k) {
379 if(mf[kx]&bx) { /* this ray is in the current m-face */
380 if(c<m+1) {
381 int i;
382
383 /* tester si cette nouvelle colonne est lin. indep. des autres */
384 /* i.e. si gauss ne donne pas de '0' sur la colonne Pi */
385 /* jusqu'a l'indice 'c' */
386
387 /* construit PiTest */
388 for(j=0;j<m+1;++j) {
389 for(i=0;i<c;++i)
390
391 /* les c premieres colonnes */
392 value_assign(PiTest->p[j][i],Pi->p[j][i]);
393
394 /* la nouvelle */
395 value_assign(PiTest->p[j][c],D->Ray[k][j+1+n]);
396 }
397 PiTest->NbColumns = c+1;
398 r = TestRank(PiTest);
399 if(r /* TestRank(PiTest) */) {
400
401 /* Ok, c'est lin. indep. */
402 for (j=0;j<n;j++)
403 value_assign(Xi->p[j][c],D->Ray[k][j+1]); /* Xi */
404 for (j=0;j<m;j++)
405 value_assign(Pi->p[j][c],D->Ray[k][j+1+n]); /* Pi */
406 value_assign(Xi->p[n][c],D->Ray[k][n+m+1]); /* const */
407 value_assign(Pi->p[m][c],D->Ray[k][n+m+1]); /* const */
408 c++;
409 }
410 }
411
412 /* Status bit */
413 value_assign(RaysDi->p[RaysDi->NbRows][0],D->Ray[k][0]);
414 Vector_Copy(&D->Ray[k][n+1],&RaysDi->p[RaysDi->NbRows][1],(m+1));
415 ++RaysDi->NbRows;
416 }
417 NEXT(kx,bx);
418 }
419
420 #ifdef DEBUGPP41
421 fprintf(stderr, "\nRaysDi=\n");
422 Matrix_Print(stderr,P_VALUE_FMT,RaysDi);
423 if(c < m+1)
424 fprintf(stderr, "Invalid ");
425 fprintf(stderr, "Pi=\n");
426 Matrix_Print(stderr,P_VALUE_FMT,Pi);
427 #endif
428
429 #ifdef DEBUGPP4
430 if(c < m+1)
431 fprintf(stderr,"Eliminated because of no vertex\n");
432 #endif
433
434 if(c < m+1)
435 return;
436
437 /* RaysDi->numColumns = m+2; */ /* stays the same */
438
439 /* Xi->NbColumns = m+1;*/ /* VIN100: stays the same. was 'c'! */
440 /* Xi->NbRows = n+1; */ /* stays the same */
441 /* Pi->NbColumns = m+1;*/ /* VIN100: stays the same. was 'c'! */
442 /* Pi->NbRows = m+1; */ /* stays the same */
443
444 #ifdef DEBUGPP4
445 fprintf(stderr,"Xi = ");
446 Matrix_Print(stderr,P_VALUE_FMT,Xi);
447 fprintf(stderr,"Pi = ");
448 Matrix_Print(stderr,P_VALUE_FMT,Pi);
449 #endif
450
451 /* (Right) invert Pi if POSSIBLE, if not then next m-face */
452 /* Pi is destroyed */
453 if(!MatInverse(Pi,PiInv)) {
454
455 #ifdef DEBUGPP4
456 fprintf(stderr, "Eliminated because of no inverse Pi\n");
457 #endif
458
459 return;
460 }
461
462 #ifdef DEBUGPP4
463 fprintf(stderr,"FACE GENERATED!\n");
464 fprintf(stderr,"PiInv = ");
465 Matrix_Print(stderr,P_VALUE_FMT,PiInv);
466 #endif
467
468 /* Compute Si (now called Ti in the paper) */
469 Si = Matrix_Alloc(Xi->NbRows,PiInv->NbColumns);
470 rat_prodmat(Si,Xi,PiInv);
471
472 #ifdef DEBUGPP4
473 fprintf(stderr,"Si = ");
474 Matrix_Print(stderr,P_VALUE_FMT,Si);
475 #endif
476
477 Si->NbRows--; /* throw out the last row = 0 ... 0 1 */
478
479 /* Copy all of that into the PV structure */
480 PV = (Param_Vertices *) malloc(sizeof(Param_Vertices));
481 PV->next = PV_Result;
482 PV->Vertex = Si;
483 PV->Domain = NULL;
484 PV->Facets = int_array2bit_vector(egalite, D->NbConstraints);
485 PV_Result = PV;
486 nbPV++; /* increment vertex count */
487
488 /* Ok... now compute the parameter domain */
489 PDi = Rays2Polyhedron(RaysDi,ws);
490
491 #ifdef DEBUGPP3
492 fprintf(stderr,"RaysDi = ");
493 Matrix_Print(stderr,P_VALUE_FMT,RaysDi);
494 fprintf(stderr,"PDi = ");
495 Polyhedron_Print(stderr,P_VALUE_FMT,PDi);
496 #endif
497
498 if(KD==0) {
499
500 /* Add the equalities again to the domain */
501 PDi = Add_CEqualities(PDi);
502 PV->Domain = Polyhedron2Constraints(PDi);
503 Polyhedron_Free(PDi);
504 }
505 else {
506 Param_Domain *PD;
507 PD = (Param_Domain *) malloc(sizeof(Param_Domain));
508 PD->Domain = PDi;
509 PD->F = NULL;
510 PD->next = PDomains;
511 PDomains = PD;
512 }
513 return;
514 } /* traite_m_face */
515
516 /*----------------------------------------------------------------------*/
517 /* count_sat */
518 /* count the number of saturated rays in the bit vector mf */
519 /* Uses nr from global area */
520 /*----------------------------------------------------------------------*/
521 int cntbit[256] = { /* counts for 8 bits */
522 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
523 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
524 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
525 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
526
527 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
528 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
529 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
530 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
531
532 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
533 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
534 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
535 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
536
537 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
538 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
539 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
540 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8 };
541
count_sat(unsigned int * mf)542 static int count_sat (unsigned int *mf) {
543
544 register unsigned int i, tmp, cnt=0;
545
546 for (i=0; i<nr; i++) {
547 tmp = mf[i];
548 cnt = cnt
549 + cntbit[ tmp & 0xff ]
550 + cntbit[ (tmp>>8) & 0xff ]
551 + cntbit[ (tmp>>16) & 0xff ]
552 + cntbit[ (tmp>>24) & 0xff ]
553 ;
554 }
555 return cnt;
556 } /* count_sat */
557
558 /* Returns true if all bits in part are also set in bv */
bit_vector_includes(unsigned int * bv,int len,unsigned int * part)559 static int bit_vector_includes(unsigned int *bv, int len, unsigned int *part)
560 {
561 int j;
562
563 for (j = 0; j < len; j++) {
564 #ifdef DEBUGPP4
565 fprintf(stderr, "mf=%08X Sat=%08X &=%08X\n",
566 part[j],bv[j], (part[j] & bv[j]));
567 #endif
568 if ((part[j] & bv[j]) != part[j])
569 return 0;
570 }
571 return 1;
572 }
573
574 /*----------------------------------------------------------------------*/
575 /* let D + E + L be the dimension of the polyhedron */
576 /* D = Dimension of polytope (ray space) */
577 /* L = Dimension of Lineality space (number of lines, bid) */
578 /* E = Dimension of Affine hull (number of equations) */
579 /* n = number of data indices */
580 /* m = number of parameters */
581 /* full domain: */
582 /* n + m = D + E + L */
583 /* projected domains: */
584 /* m = D_m + E_m + L_m */
585 /* n = D_n + E_n + L_n */
586 /* What dimension M-face, when projected onto parameter space, */
587 /* will give an L_m-face? */
588 /* What are the conditions? */
589 /* - at least one vertex */
590 /* - number of rays >= D_m+1 after removal of redundants */
591 /* */
592 /* dim of face nb saturated constraints nb saturated lines,rays */
593 /* ----------- ------------------------ ----------------------- */
594 /* (0+L)-face all E eqns + >=D ineq all L lines + 1 ray */
595 /* (M+L)-face all E eqns + >=(D-M) ineq all L lines + >=(M+1) rays */
596 /* (D+L)-face all E eqns + 0 ineq all L lines + >=(D+1) rays */
597 /*----------------------------------------------------------------------*/
598 /*----------------------------------------------------------------------*/
599 /* scan_m_face */
600 /* pos : the next candidate constraint position */
601 /* nb_un : number of saturated constraints needed to finish a face */
602 /* D : the source polyhedron (context included ) */
603 /* mf : bit-array marking rays which are saturated so far */
604 /* From Global area: */
605 /* ---------------- */
606 /* n : number of data indices */
607 /* m : number of parameters */
608 /* egalite : boolean vector marking saturated constraints in m-face */
609 /* Sat : Saturation Matrix (row=constraints, col=rays) */
610 /* ws : working space size */
611 /* nr : (NbRays-1)/32 + 1 */
612 /* */
613 /* Recursive function to find the rays and vertices of each m-face */
614 /*----------------------------------------------------------------------*/
scan_m_face(int pos,int nb_un,Polyhedron * D,unsigned int * mf)615 static void scan_m_face(int pos,int nb_un,Polyhedron *D,unsigned int *mf) {
616 /* pos - the next candidate constraint position */
617 /* nb_un - the number of constraints needed to finish a face */
618 /* D - the source polyhedron */
619 /* mf - (bit vector) marks rays that are saturated so far */
620
621 unsigned int *new_mf;
622
623 #ifdef DEBUGPP4
624 fprintf(stderr,"Start scan_m_face(pos=%d, nb_un=%d, n=%d, m=%d\n",
625 pos,nb_un,n,m);
626 fprintf(stderr,"mf = ");
627 {
628 int i;
629 for(i=0;i<nr;i++)
630 fprintf(stderr,"%08X", mf[i]);
631 fprintf(stderr,"\nequality = [");
632 for(i=0;i<D->NbConstraints;i++)
633 fprintf(stderr," %1d",egalite[i]);
634 fprintf(stderr,"]\n");
635 }
636 #endif
637
638 if(nb_un == 0) { /* Base case */
639 int i;
640
641 /*********** ELIMINATION OF REDUNDANT FACES ***********/
642 /* if all these vertices also verify a previous constraint */
643 /* which is NOT already selected, we eliminate this face */
644 /* This keeps the lexicographically greatest selection */
645 for(i=0;i<pos-1;i++)
646 {
647 if(egalite[i])
648 continue; /* already selected */
649
650 /* if Sat[i] & mf == mf then it's redundant */
651 #ifdef DEBUGPP4
652 fprintf(stderr, "Sat[%d]\n", i);
653 #endif
654 if (bit_vector_includes(Sat->p[i], nr, mf)) {
655 #ifdef DEBUGPP4
656 fprintf(stderr, "Redundant with constraint %d\n", i);
657 #endif
658 return; /* it is redundant */
659 }
660 }
661 /********* END OF ELIMINATION OF DEGENERATE FACES *********/
662 /* Now check for other constraints that are verified */
663 for (i = pos; i < D->NbConstraints; ++i) {
664 if (bit_vector_includes(Sat->p[i], nr, mf))
665 egalite[i] = 1;
666 }
667 /* if we haven't found a constraint verified by all */
668 /* the rays, its OK, it's a new face. */
669 traite_m_face(D, mf, egalite);
670 for (i = pos; i < D->NbConstraints; ++i)
671 egalite[i] = 0;
672 return;
673 }
674
675 /* See if there are enough constraints left to finish */
676 if((pos+nb_un)>D->NbConstraints) return;
677
678 /* Recurring part of the procedure */
679 /* Add the pos'th constraint, compute new saturation vector */
680 {
681 int k;
682 new_mf = (unsigned int *)malloc(nr*sizeof(unsigned int));
683 for (k=0; k<nr; k++)
684 new_mf[k] = mf[k] & Sat->p[pos][k];
685 }
686 #ifdef DEBUGPP4
687 fprintf(stderr,"new_mf = ");
688 {
689 int i;
690 for(i=0;i<nr;i++) {
691 fprintf(stderr,"%08X", new_mf[i]);
692 }
693 fprintf(stderr,"\ncount(new_mf) = %d\n",count_sat(new_mf));
694 }
695 #endif
696
697 {
698 int c;
699 c = count_sat(new_mf);
700 /* optimization : at least m_dim+1 rays must be saturated to add this constraint */
701 if (c>m_dim )
702 {
703 int redundant = 0;
704
705 egalite[pos]=1; /* Try it with the pos-th constraint */
706
707 /* If this constraint does not change anything,
708 * it is redundant with respect to the selected
709 * equalities and the remaining inequalities.
710 * Check whether it is redundant with respect
711 * to just the selected equalities.
712 */
713 if( c==count_sat(mf) ) {
714 int i, c, j;
715
716 for (i = 0, c = 0; i < D->NbConstraints; ++i) {
717 if (egalite[i] == 0 || egalite[i] == -1)
718 continue;
719 for (j = 0; j < D->Dimension+1; ++j)
720 value_assign(CTest->p[j][c],
721 D->Constraint[i][j+1]);
722 ++c;
723 }
724 CTest->NbColumns = c;
725 #ifdef DEBUGPP41
726 Matrix_Print(stderr,P_VALUE_FMT,CTest);
727 #endif
728 redundant = !TestRank(CTest);
729 }
730
731 /* Do not decrement nb_un if equality is redundant. */
732 if( redundant )
733 {
734 egalite[pos]=-1; /* Don't use in further redundance test
735 */
736 scan_m_face(pos+1,nb_un,D,new_mf);
737 }
738 else
739 {
740 scan_m_face(pos+1,nb_un-1,D,new_mf);
741 }
742 }
743 }
744 free(new_mf);
745 egalite[pos]=0; /* Try it without the pos-th constraint */
746 if ((pos+nb_un)>=D->NbConstraints) return;
747 scan_m_face(pos+1,nb_un,D,mf);
748 return;
749 } /* scan_m_face */
750
751 /*
752 * Create a saturation matrix with rows correspond to the constraints and
753 * columns correspond to the rays of the polyhedron 'Pol'. Global variable
754 * 'nr' is set in the function.
755 */
Poly2Sat(Polyhedron * Pol,unsigned int ** L)756 static SatMatrix *Poly2Sat(Polyhedron *Pol,unsigned int **L) {
757
758 SatMatrix *Sat;
759 int i, j, k, kx;
760 unsigned int *Temp;
761 Value *p1, *p2, p3,tmp;
762 unsigned Dimension, NbRay, NbCon, bx;
763
764 /* Initialize all the 'Value' variables */
765 value_init(p3); value_init(tmp);
766
767 NbRay = Pol->NbRays;
768 NbCon = Pol->NbConstraints;
769 Dimension = Pol->Dimension+1; /* Homogeneous Dimension */
770
771 /* Build the Sat matrix */
772 nr = (NbRay - 1)/(sizeof(int)*8) + 1; /* Set globally */
773 Sat = SMAlloc(NbCon,nr);
774 Temp = (unsigned int *)malloc(nr*sizeof(unsigned int));
775 memset(Sat->p_init,0,nr*NbCon*sizeof(int));
776 memset(Temp,0,nr*sizeof(unsigned int));
777 kx=0; bx=MSB;
778 for (k=0; k<NbRay; k++) {
779 for (i=0; i<NbCon; i++) {
780 p1 = &Pol->Constraint[i][1];
781 p2 = &Pol->Ray[k][1];
782 value_set_si(p3,0);
783 for (j=0;j<Dimension;j++) {
784 value_multiply(tmp,*p1,*p2);
785 value_addto(p3,p3,tmp);
786 p1++; p2++;
787 }
788 if (value_zero_p(p3))
789 Sat->p[i][kx]|=bx;
790 }
791 Temp[kx] |= bx;
792 NEXT(kx, bx);
793 }
794
795 /* Set 'L' to an array containing ones in every bit position of its */
796 /* elements. */
797 *L = Temp;
798
799 /* Clear all the 'Value' variables */
800 value_clear(p3); value_clear(tmp);
801
802 return Sat;
803 } /* Poly2Sat */
804
805 /*
806 * Create a parametrized polyhedron with zero parameters. This function was
807 * first written by Xavier Redon, and was later modified by others.
808 */
GenParamPolyhedron(Polyhedron * Pol,Matrix * Rays)809 Param_Polyhedron *GenParamPolyhedron(Polyhedron *Pol, Matrix *Rays)
810 {
811 Param_Polyhedron *result;
812 int nbRows, nbColumns;
813 int i, size, rays;
814
815 nbRows=Pol->NbRays;
816 nbColumns=Pol->Dimension+2;
817
818 /* Count the number of rays */
819 for(i=0, rays=0; i<nbRows; i++)
820 if(value_notone_p(Pol->Ray[i][0]) ||
821 value_zero_p(Pol->Ray[i][nbColumns-1]))
822 ++rays;
823
824 /* Initialize the result */
825 result=(Param_Polyhedron *)malloc(sizeof(Param_Polyhedron));
826 result->nbV=nbRows-rays;
827 result->V=NULL;
828 result->Constraints = Polyhedron2Constraints(Pol);
829 result->Rays = Rays;
830
831 /* Build the parametric vertices */
832 for(i=0;i<nbRows;i++) {
833 Matrix *vertex;
834 Param_Vertices *paramVertex;
835 int j;
836
837 if (value_notone_p(Pol->Ray[i][0]) ||
838 value_zero_p(Pol->Ray[i][nbColumns-1]))
839 continue;
840
841 vertex=Matrix_Alloc(nbColumns-2,2);
842 for(j=1;j<nbColumns-1;j++) {
843 value_assign(vertex->p[j-1][0],Pol->Ray[i][j]);
844 value_assign(vertex->p[j-1][1],Pol->Ray[i][nbColumns-1]);
845 }
846 paramVertex=(Param_Vertices *)malloc(sizeof(Param_Vertices));
847 paramVertex->Vertex=vertex;
848
849 /* There is one validity domain : universe of dimension 0 */
850 paramVertex->Domain=Matrix_Alloc(1,2);
851 value_set_si(paramVertex->Domain->p[0][0],1);
852 value_set_si(paramVertex->Domain->p[0][1],1);
853 paramVertex->Facets = NULL;
854 paramVertex->next=result->V;
855 result->V=paramVertex;
856 }
857
858 /* Build the parametric domains (only one here) */
859 if (nbRows > 1)
860 size=(nbRows-1)/(8*sizeof(int))+1;
861 else
862 size = 1;
863 result->D=(Param_Domain *)malloc(sizeof(Param_Domain));
864 result->D->next=NULL;
865 result->D->Domain=Universe_Polyhedron(0);
866 result->D->F=(unsigned int *)malloc(size*sizeof(int));
867 memset(&result->D->F[0],0xFF,size*sizeof(int));
868
869 return result;
870 } /* GenParamPolyhedron */
871
872
873 /*----------------------------------------------------------------------*/
874 /* PreElim_Columns */
875 /* function being called before Elim_Columns */
876 /* Equalities in E are analysed to initialize ref and p. */
877 /* These two vectors are used to construct the new constraint matrix */
878 /* PreElim_Columns returns the transformation matrix to re-convert the */
879 /* resulting domains in the same format (by adding empty columns) */
880 /* in the parameter space */
881 /*----------------------------------------------------------------------*/
PreElim_Columns(Polyhedron * E,int * p,int * ref,int m)882 Matrix *PreElim_Columns(Polyhedron *E,int *p,int *ref,int m) {
883
884 int i,j,l;
885 Matrix *T;
886
887 /* find which columns to eliminate */
888 /* p contains, for each line in E, the column to eliminate */
889 /* (i.e. the corresponding parameter index to eliminate) */
890 /* 0 <= i <= E->NbEq, and 1 <= p[i] <= E->Dimension */
891 memset(p,0,sizeof(int)*(E->NbEq));
892
893 #ifdef DEBUGPP32
894 fprintf(stderr,"\n\nPreElim_Columns starting\n");
895 fprintf(stderr,"E =\n");
896 Polyhedron_Print(stderr,P_VALUE_FMT,E);
897 #endif
898
899 for(l=0;l<E->NbEq;++l) {
900 if(value_notzero_p(E->Constraint[l][0])) {
901 fprintf(stderr,"Internal error: Elim_Columns (polyparam.c), expecting equalities first in E.\n");
902 free(p);
903 return(NULL);
904 }
905 for(i=1;value_zero_p(E->Constraint[l][i]);++i) {
906 if(i==E->Dimension+1) {
907 fprintf(stderr,"Internal error: Elim_Columns (polyparam.c), expecting non-empty constraint in E.\n");
908 free(p);
909 return( NULL );
910 }
911 }
912 p[l] = i;
913
914 #ifdef DEBUGPP32
915 fprintf(stderr,"p[%d] = %d,",l,p[l]);
916 #endif
917 }
918
919 /* Reference vector: column ref[i] in A corresponds to column i in M */
920 for(i=0;i<E->Dimension+2-E->NbEq;++i) {
921 ref[i]=i;
922 for(j=0;j<E->NbEq;++j)
923 if(p[j]<=ref[i])
924 ref[i]++;
925
926 #ifdef DEBUGPP32
927 fprintf(stderr,"ref[%d] = %d,",i,ref[i]);
928 #endif
929 }
930
931 /* Size of T : phdim-nbEq * phdim */
932 T = Matrix_Alloc(m+1-E->NbEq,m+1);
933 for(i=0;i<T->NbColumns;i++)
934 for(j=0;j<T->NbRows;j++)
935 if(ref[E->Dimension-m+j+1] == E->Dimension-m+i+1)
936 value_set_si(T->p[j][i],1);
937 else
938 value_set_si(T->p[j][i],0);
939
940 #ifdef DEBUGPP32
941 fprintf(stderr,"\nTransMatrix =\n");
942 Matrix_Print(stderr,P_VALUE_FMT,T);
943 #endif
944
945 return(T);
946
947 } /* PreElim_Columns */
948
949 /*----------------------------------------------------------------------*/
950 /* Elim_Columns */
951 /* Eliminate columns from A, using the equalities in E. */
952 /* ref and p are precalculated by PreElim_Columns, using E; */
953 /* these two vectors are used to construct the new constraint matrix */
954 /*----------------------------------------------------------------------*/
Elim_Columns(Polyhedron * A,Polyhedron * E,int * p,int * ref)955 Polyhedron *Elim_Columns(Polyhedron *A,Polyhedron *E,int *p,int *ref) {
956
957 int i,l,c;
958 Matrix *M, *C;
959 Polyhedron *R;
960 Value tmp1,tmp2;
961
962 /* Initialize all the 'Value' variables */
963 value_init(tmp1); value_init(tmp2);
964
965 #ifdef DEBUGPP32
966 fprintf(stderr,"\nElim_Columns starting\n");
967 fprintf(stderr,"A =\n" );
968 Polyhedron_Print(stderr,P_VALUE_FMT,A);
969 #endif
970
971 /* Builds M = constraint matrix of A, useless columns zeroed */
972 M = Polyhedron2Constraints(A);
973 for(l=0;l<E->NbEq;++l) {
974 for(c=0;c<M->NbRows;++c) {
975 if(value_notzero_p(M->p[c][p[l]])) {
976
977 /* A parameter to eliminate here ! */
978 for(i=1;i<M->NbColumns;++i) {
979 value_multiply(tmp1,M->p[c][i],E->Constraint[l][p[l]]);
980 value_multiply(tmp2,E->Constraint[l][i],M->p[c][p[l]]);
981 value_subtract(M->p[c][i],tmp1,tmp2);
982 }
983 }
984 }
985 }
986
987 #ifdef DEBUGPP32
988 fprintf(stderr,"\nElim_Columns after zeroing columns of A.\n");
989 fprintf(stderr,"M =\n");
990 Matrix_Print(stderr,P_VALUE_FMT,M);
991 #endif
992
993 /* Builds C = constraint matrix, useless columns eliminated */
994 C = Matrix_Alloc(M->NbRows,M->NbColumns - E->NbEq);
995 for(l=0;l<C->NbRows;++l)
996 for(c=0;c<C->NbColumns;++c) {
997 value_assign(C->p[l][c],M->p[l][ref[c]]);
998 }
999
1000 #ifdef DEBUGPP32
1001 fprintf(stderr,"\nElim_Columns after eliminating columns of A.\n");
1002 fprintf(stderr,"C =\n");
1003 Matrix_Print(stderr,P_VALUE_FMT,C);
1004 #endif
1005
1006 R = Constraints2Polyhedron(C,ws);
1007 Matrix_Free(C);
1008 Matrix_Free(M);
1009 value_clear(tmp1); value_clear(tmp2);
1010 return(R);
1011 } /* Elim_Columns */
1012
1013
Recession_Cone(Polyhedron * P,unsigned nvar,unsigned MaxRays)1014 static Polyhedron *Recession_Cone(Polyhedron *P, unsigned nvar, unsigned MaxRays)
1015 {
1016 int i;
1017 Matrix *M = Matrix_Alloc(P->NbConstraints, 1 + nvar + 1);
1018 Polyhedron *R;
1019 for (i = 0; i < P->NbConstraints; ++i)
1020 Vector_Copy(P->Constraint[i], M->p[i], 1+nvar);
1021 R = Constraints2Polyhedron(M, MaxRays);
1022 Matrix_Free(M);
1023 return R;
1024 }
1025
1026 /* Compute lines/unidirectional rays of the (non parametric) polyhedron */
1027 /* Input :
1028 * D1 : combined polyhedron,
1029 * Output :
1030 * Rays : non parametric ray matrix
1031 * return value : number of lines
1032 */
ComputeNPLinesRays(int n,Polyhedron * D1,Matrix ** Rays)1033 static int ComputeNPLinesRays(int n, Polyhedron *D1, Matrix **Rays)
1034 {
1035 int i, j, nbr, dimfaces;
1036 Polyhedron *RC; /* Recession Cone */
1037
1038 RC = Recession_Cone(D1, n, ws);
1039
1040 /* get the rays/lines from RC */
1041 nbr = 0;
1042 for (i = 0; i < RC->NbRays ;i++)
1043 if (value_zero_p(RC->Ray[i][n+1]))
1044 nbr++;
1045 *Rays=Matrix_Alloc(nbr, n+2);
1046 for (i = 0, j = 0; j < nbr ;i++)
1047 if (value_zero_p(RC->Ray[i][n+1]))
1048 Vector_Copy(RC->Ray[i], (*Rays)->p[j++], n+2);
1049
1050 dimfaces = RC->NbBid;
1051 Polyhedron_Free(RC);
1052
1053 #ifdef DEBUGPP31
1054 fprintf(stderr, "Rays = ");
1055 Matrix_Print(stderr, P_VALUE_FMT, *Rays);
1056 fprintf(stderr, "dimfaces = %d\n", dimfaces);
1057 #endif
1058
1059 return dimfaces;
1060 }
1061
1062 /*
1063 * Given a polyhedron 'Di' in combined data and parameter space and a context
1064 * polyhedron 'C' representing the constraints on the parameter space, create
1065 * a list of parameterized vertices and assign values to global variables:
1066 * m,n,ws,Sat,egalite,mf,Xi,Pi,PiInv,RaysDi,CEqualities.
1067 */
Find_m_faces(Polyhedron ** Di,Polyhedron * C,int keep_dom,int working_space,Polyhedron ** CEq,Matrix ** CT)1068 Param_Polyhedron *Find_m_faces(Polyhedron **Di,Polyhedron *C,int keep_dom,int working_space,Polyhedron **CEq,Matrix **CT) {
1069
1070 unsigned int *mf;
1071 int i, j, dimfaces;
1072 Polyhedron *D=*Di,
1073 *C1, /* true context */
1074 *D1; /* the combined polyhedron, including context C */
1075 Matrix *Rays, /* lines/rays (non parametric) */
1076 *M;
1077 Param_Polyhedron *res;
1078 int *p, *ref;
1079
1080 CEqualities = NULL;
1081
1082 if(CT) {
1083 *CEq = NULL;
1084 *CT = NULL;
1085 }
1086
1087 if(!D || !C)
1088 return (Param_Polyhedron *) 0;
1089
1090 ws = working_space;
1091 m = C->Dimension;
1092 n = D->Dimension - m;
1093 if(n<0) {
1094 fprintf(stderr,
1095 "Find_m_faces: ?%d parameters of a %d-polyhedron !\n",m,n);
1096 return (Param_Polyhedron *) 0;
1097 }
1098 if (m==0)
1099 return GenParamPolyhedron(D,Matrix_Alloc(0,2));
1100
1101 /* Add constraints from Context to D -> result in D1 */
1102 C1 = align_context(C,D->Dimension,ws);
1103
1104 #ifdef DEBUGPP31
1105 fprintf(stderr,"m = %d\n",m);
1106 fprintf(stderr, "D = ");
1107 Polyhedron_Print(stderr,P_VALUE_FMT,D);
1108 fprintf(stderr,"C1 = ");
1109 Polyhedron_Print(stderr,P_VALUE_FMT,C1);
1110 #endif
1111
1112 D1 = DomainIntersection(D,C1,ws);
1113
1114 #ifdef DEBUGPP31
1115 fprintf(stderr,"D1 = ");
1116 Polyhedron_Print(stderr,P_VALUE_FMT,D1);
1117 #endif
1118
1119 Domain_Free(C1);
1120
1121 if (!D1)
1122 return(NULL);
1123 if (emptyQ(D1)) {
1124 Polyhedron_Free(D1);
1125 return(NULL);
1126 }
1127
1128 /* Compute the true context C1 */
1129 /* M : lines in the direction of the first n indices (index space) */
1130 M = Matrix_Alloc(n, D1->Dimension+2);
1131 for (i=0; i<n; i++)
1132 value_set_si(M->p[i][i+1],1);
1133 C1 = DomainAddRays(D1,M,ws);
1134 Matrix_Free(M);
1135
1136 #ifdef DEBUGPP31
1137 fprintf(stderr,"True context C1 = ");
1138 Polyhedron_Print(stderr,P_VALUE_FMT,C1);
1139 #endif
1140
1141 dimfaces = ComputeNPLinesRays(n, D1, &Rays);
1142
1143 /* CEqualities contains the constraints (to be added again later) */
1144 /* *CT is the transformation matrix to add the removed parameters */
1145 if(!CT) {
1146 if (C1->NbEq == 0) {
1147 Polyhedron_Free(C1);
1148 C1 = NULL;
1149 } else {
1150 Polyhedron *CEq1, /* CEqualities, in homogeneous dim */
1151 *D2; /* D1, (temporary) simplified */
1152
1153 /* Remove equalities from true context C1 and from D1 */
1154 /* Compute CEqualities = matrix of equalities in C1, projected in */
1155 /* the parameter space */
1156 M = Matrix_Alloc(C1->NbEq,m+2);
1157 for(j=0,i=0;i<C1->NbEq;++i,++j) {
1158 while(value_notzero_p(C1->Constraint[j][0]))
1159 ++j;
1160 value_assign(M->p[i][0],C1->Constraint[j][0]);
1161 Vector_Copy(&C1->Constraint[j][D->Dimension-m+1],&M->p[i][1],(m+1));
1162 }
1163 CEqualities = Constraints2Polyhedron(M,ws);
1164 Matrix_Free(M);
1165 CEq1 = align_context(CEqualities,D->Dimension,ws);
1166
1167 /* Simplify D1 and C1 (remove the equalities) */
1168 D2 = DomainSimplify(D1,CEq1,ws);
1169 Polyhedron_Free(D1);
1170 Polyhedron_Free(C1);
1171 Polyhedron_Free(CEq1);
1172 D1 = D2;
1173 C1 = NULL;
1174 }
1175 }
1176 else { /* if( CT ) */
1177 Polyhedron *CEq1, /* CEqualities */
1178 *D2; /* D1, (temporary) simplified */
1179
1180 /* Suppress all useless constraints in parameter domain */
1181 /* when CT is not NULL (ehrhart) */
1182 /* Vin100, march 01 */
1183 CEq1 = C1;
1184 M = Matrix_Alloc(C1->NbConstraints,m+2);
1185 for(i=0;i<C1->NbConstraints;++i) {
1186 value_assign(M->p[i][0],C1->Constraint[i][0]);
1187 Vector_Copy(&C1->Constraint[i][D->Dimension-m+1],&M->p[i][1],(m+1));
1188 }
1189 CEqualities = Constraints2Polyhedron( M, ws );
1190 Matrix_Free(M);
1191
1192 D2 = DomainSimplify(D1,CEq1,ws);
1193 Polyhedron_Free(D1);
1194 D1 = D2;
1195 C1 = Universe_Polyhedron(D2->Dimension);
1196
1197 /* if CT is not NULL, the constraints are eliminated */
1198 /* *CT will contain the transformation matrix to come back to the */
1199 /* original dimension (for a polyhedron, in the parameter space) */
1200 if( CEq1->NbEq )
1201 {
1202 m -= CEq1->NbEq;
1203 p = (int *)malloc(sizeof(int)*(CEq1->NbEq));
1204 }
1205 else
1206 p = NULL;
1207 ref = (int*) malloc(sizeof(int)*
1208 (CEq1->Dimension+2-CEq1->NbEq));
1209 *CT = PreElim_Columns(CEq1,p,ref,CEqualities->Dimension);
1210 D2 = Elim_Columns(D1,CEq1,p,ref);
1211 if (p)
1212 free(p);
1213 free(ref);
1214
1215 #ifdef DEBUGPP3
1216 fprintf(stderr,"D2\t Dim = %3d\tNbEq = %3d\tLines = %3d\n",
1217 D2->Dimension,D2->NbEq,D2->NbBid);
1218 C2 = Elim_Columns(C1,CEq1,p,ref);
1219 fprintf(stderr,"C2\t Dim = %3d\tNbEq = %3d\tLines = %3d\n",
1220 C2->Dimension,C2->NbEq,C2->NbBid);
1221 Polyhedron_Free(C2);
1222 #endif
1223
1224 Polyhedron_Free(D1);
1225 Polyhedron_Free(C1);
1226 D1 = D2;
1227 C1 = NULL;
1228 *CEq = CEqualities;
1229
1230 #ifdef DEBUGPP3
1231 fprintf(stderr,"Polyhedron CEq = ");
1232 Polyhedron_Print(stderr,P_VALUE_FMT,*CEq);
1233 fprintf(stderr,"Matrix CT = ");
1234 Matrix_Print(stderr,P_VALUE_FMT,*CT);
1235 #endif
1236
1237 Polyhedron_Free(CEq1);
1238 CEqualities = NULL; /* don't simplify ! */
1239
1240 /* m changed !!! */
1241 if(m==0) {
1242 /* return the new D1 too */
1243 *Di = D1;
1244
1245 return GenParamPolyhedron(D1, Rays);
1246 }
1247 }
1248
1249 #ifdef DEBUGPP3
1250 fprintf(stderr,"Polyhedron D1 (D AND C) = ");
1251 Polyhedron_Print(stderr,P_VALUE_FMT, D1);
1252 fprintf(stderr,"Polyhedron CEqualities = ");
1253 if(CEqualities) Polyhedron_Print(stderr,P_VALUE_FMT, CEqualities);
1254 else fprintf(stderr,"NULL\n");
1255 #endif
1256
1257 KD = keep_dom;
1258 PDomains = NULL;
1259 PV_Result = NULL;
1260 nbPV = 0;
1261
1262 if (emptyQ(D1)) {
1263 Polyhedron_Free(D1);
1264 Matrix_Free(Rays);
1265 return NULL;
1266 }
1267
1268 /* mf : a bit array indicating which rays are part of the m-face */
1269 /* Poly2Sat initializes mf to all ones */
1270 /* set global variable nr to size (number of words) of mf */
1271 Sat = Poly2Sat(D1,&mf);
1272
1273 #ifdef DEBUGPP4
1274 fprintf(stderr,"Sat = ");
1275 SMPrint(Sat);
1276
1277 fprintf(stderr,"mf = ");
1278 for (i=0; i<nr; i++)
1279 fprintf(stderr,"%08X", mf[i]);
1280 fprintf(stderr, "\n");
1281 #endif
1282
1283 /* A boolean array saying which constraints are part of the m-face */
1284 egalite = (unsigned int *)malloc(sizeof(int)*D1->NbConstraints);
1285 memset(egalite,0, sizeof(int)*D1->NbConstraints);
1286
1287 for (i=0; i<D1->NbEq; i++)
1288 egalite[i] = 1;
1289
1290 Xi = Matrix_Alloc(n+1,m+1);
1291 Pi = Matrix_Alloc(m+1,m+1);
1292 PiTest = Matrix_Alloc(m+1,m+1);
1293 CTest = Matrix_Alloc(D->Dimension+1,D->NbConstraints);
1294 PiInv = Matrix_Alloc(m+1,m+2);
1295 RaysDi = Matrix_Alloc(D1->NbRays,m+2);
1296 m_dim = m;
1297
1298 /* m_dim has to be increased by the dimension of the smallest faces
1299 * of the (non parametric) polyhedron
1300 */
1301 m_dim += dimfaces;
1302
1303 /* if the smallest face is of smaller dimension than m_dim,
1304 * then increase m_dim (I think this should never happen --Vincent)
1305 */
1306 #ifdef DEBUGPP3
1307 if (m_dim < D1->NbBid)
1308 fprintf(stderr, "m_dim (%d) < D1->NbBid (%d)\n", m_dim, D1->NbBid );
1309 #endif
1310 if (m_dim < D1->NbBid)
1311 m_dim = D1->NbBid;
1312
1313 #ifdef DEBUGPP
1314 nbfaces=0;
1315 #endif
1316 #ifdef DEBUGPP3
1317 fprintf(stderr, "m_dim = %d\n", m_dim);
1318 fprintf(stderr,
1319 "Target: find faces that saturate %d constraints and %d rays/lines\n",
1320 D1->Dimension - m_dim,m_dim+1);
1321 #endif
1322
1323 /* D1->NbEq constraints already saturated ! */
1324 scan_m_face(D1->NbEq,(D1->Dimension - m_dim - D1->NbEq),D1,mf);
1325
1326 /* pos, number of constraints needed */
1327
1328 #ifdef DEBUGPP
1329 fprintf( stderr, "Number of m-faces: %d\n", nbfaces );
1330 #endif
1331
1332 Matrix_Free(RaysDi);
1333 Matrix_Free(PiInv);
1334 Matrix_Free(PiTest);
1335 Matrix_Free(CTest);
1336 Matrix_Free(Pi);
1337 Matrix_Free(Xi);
1338 free(egalite);
1339 free(mf);
1340 SMFree(Sat);
1341
1342 /* if(CEqualities && keep_dom==0) {
1343 Domain_Free(CEqualities);
1344 }
1345 */
1346
1347 res = (Param_Polyhedron *) malloc (sizeof(Param_Polyhedron));
1348 res->nbV = nbPV;
1349 res->V = PV_Result;
1350 res->D = PDomains;
1351 res->Constraints = Polyhedron2Constraints(D1);
1352 res->Rays = Rays;
1353
1354 if(CT) /* return the new D1 too ! */
1355 *Di = D1;
1356 else
1357 Domain_Free(D1);
1358
1359 return(res);
1360 } /* Find_m_faces */
1361
1362 /*
1363 * Given parametric domain 'PD' and number of parametric vertices 'nb_domains',
1364 * find the vertices that belong to distinct sub-domains.
1365 */
Compute_PDomains(Param_Domain * PD,int nb_domains,int working_space)1366 void Compute_PDomains(Param_Domain *PD,int nb_domains,int working_space) {
1367
1368 unsigned bx;
1369 int i, ix, nv;
1370 Polyhedron *dx, *d1, *d2;
1371 Param_Domain *p1, *p2, *p2prev, *PDNew;
1372
1373 if (nb_domains==0) {
1374
1375 #ifdef DEBUGPP5
1376 fprintf(stderr,"No domains\n");
1377 #endif
1378
1379 return;
1380 }
1381
1382 /* Already filled out by GenParamPolyhedron */
1383 if (!PD->next && PD->F)
1384 return;
1385
1386 /* Initialization */
1387 nv = (nb_domains - 1)/(8*sizeof(int)) + 1;
1388
1389 #ifdef DEBUGPP5
1390 fprintf(stderr,"nv = %d\n",nv);
1391 #endif
1392
1393 for(p1=PD,i=0,ix=0,bx=MSB;p1;p1=p1->next,i++) {
1394
1395 /* Assign a bit array 'p1->F' of suitable size to include the vertices */
1396 p1->F = (unsigned *) malloc (nv * sizeof(unsigned));
1397
1398 /* Set the bit array to zeros */
1399 memset(p1->F,0,nv * sizeof(unsigned));
1400 p1->F[ix] |= bx; /* Set i'th bit to one */
1401 NEXT(ix, bx);
1402 }
1403
1404 #ifdef DEBUGPP5
1405 fprintf(stderr,"nb of vertices=%d\n",i);
1406 #endif
1407
1408 /* Walk the PD list with two pointers */
1409 ix = 0; bx=MSB;
1410 for (p1=PD;p1;p1=p1->next) {
1411 for (p2prev=p1,p2=p1->next;p2;p2prev=p2,p2=p2->next) {
1412
1413 /* Find intersection */
1414 dx = PDomainIntersection(p1->Domain,p2->Domain,working_space);
1415
1416 if (!dx || emptyQ(dx)) {
1417
1418 #ifdef DEBUGPP5
1419 fprintf( stderr, "Empty dx (p1 inter p2). Continuing\n");
1420 #endif
1421 if(dx)
1422 Domain_Free(dx);
1423 continue;
1424 }
1425
1426 #ifdef DEBUGPP5
1427 fprintf(stderr,"Begin PDomainDifference\n");
1428 fprintf(stderr, "p1=");
1429 Polyhedron_Print(stderr,P_VALUE_FMT,p1->Domain);
1430 fprintf(stderr,"p2=");
1431 Polyhedron_Print(stderr,P_VALUE_FMT,p2->Domain);
1432 #endif
1433 d1 = PDomainDifference(p1->Domain,p2->Domain,working_space);
1434 d2 = PDomainDifference(p2->Domain,p1->Domain,working_space);
1435
1436 #ifdef DEBUGPP5
1437 fprintf(stderr,"p1\\p2=");
1438 Polyhedron_Print(stderr,P_VALUE_FMT,d1);
1439 fprintf(stderr,"p2\\p1=");
1440 Polyhedron_Print(stderr,P_VALUE_FMT,d2);
1441 fprintf(stderr,"END PDomainDifference\n\n");
1442 #endif
1443 if (!d1 || emptyQ(d1) || d1->NbEq!=0) {
1444
1445 #ifdef DEBUGPP5
1446 fprintf(stderr,"Empty d1\n");
1447 #endif
1448 if (d1)
1449 Domain_Free(d1);
1450 Domain_Free(dx);
1451
1452 if (!d2 || emptyQ(d2) || d2->NbEq!=0) {
1453
1454 #ifdef DEBUGPP5
1455 fprintf( stderr, "Empty d2 (deleting)\n");
1456 #endif
1457 /* dx = p1->Domain = p2->Domain */
1458 if (d2) Domain_Free(d2);
1459
1460 /* Update p1 */
1461 for (i=0;i<nv;i++)
1462 p1->F[i] |= p2->F[i];
1463
1464 /* Delete p2 */
1465 p2prev->next = p2->next;
1466 Domain_Free(p2->Domain);
1467 free(p2->F);
1468 free(p2);
1469 p2 = p2prev;
1470 }
1471 else { /* d2 is not empty --> dx==p1->domain */
1472
1473 #ifdef DEBUGPP5
1474 fprintf( stderr, "p2 replaced by d2\n");
1475 #endif
1476 /* Update p1 */
1477 for(i=0;i<nv;i++)
1478 p1->F[i] |= p2->F[i];
1479
1480 /* Replace p2 with d2 */
1481 Domain_Free( p2->Domain );
1482 p2->Domain = d2;
1483 }
1484 }
1485 else { /* d1 is not empty */
1486 if (!d2 || emptyQ(d2) || d2->NbEq!=0) {
1487
1488 #ifdef DEBUGPP5
1489 fprintf( stderr, "p1 replaced by d1\n");
1490 #endif
1491 if (d2) Domain_Free(d2);
1492
1493 /* dx = p2->domain */
1494 Domain_Free(dx);
1495
1496 /* Update p2 */
1497 for(i=0;i<nv;i++)
1498 p2->F[i] |= p1->F[i];
1499
1500 /* Replace p1 with d1 */
1501 Domain_Free(p1->Domain);
1502 p1->Domain = d1;
1503 }
1504 else { /*d2 is not empty-->d1,d2,dx are distinct */
1505
1506 #ifdef DEBUGPP5
1507 fprintf(stderr,"Non-empty d1 and d2\nNew node created\n");
1508 #endif
1509 /* Create a new node for dx */
1510 PDNew = (Param_Domain *) malloc( sizeof(Param_Domain) );
1511 PDNew->F = (unsigned int *)malloc( nv*sizeof(int) );
1512 memset(PDNew->F,0,nv*sizeof(int));
1513 PDNew->Domain = dx;
1514
1515 for (i=0;i<nv;i++)
1516 PDNew->F[i] = p1->F[i] | p2->F[i];
1517
1518 /* Replace p1 with d1 */
1519 Domain_Free( p1->Domain );
1520 p1->Domain = d1;
1521
1522 /* Replace p2 with d2 */
1523 Domain_Free( p2->Domain );
1524 p2->Domain = d2;
1525
1526 /* Insert new node after p1 */
1527 PDNew->next = p1->next;
1528 p1->next = PDNew;
1529 }
1530 }
1531 } /* end of p2 scan */
1532 if (p1->Domain->next) {
1533 Polyhedron *C = DomainConvex(p1->Domain, working_space);
1534 Domain_Free(p1->Domain);
1535 p1->Domain = C;
1536 }
1537 } /* end of p1 scan */
1538 } /* Compute_PDomains */
1539
1540 /*
1541 * Given a polyhedron 'Din' in combined data and parametre space, a context
1542 * polyhedron 'Cin' representing the constraints on the parameter space and
1543 * a working space size 'working_space', return a parametric polyhedron with
1544 * a list of parametric vertices and their defining domains.
1545 */
Polyhedron2Param_Vertices(Polyhedron * Din,Polyhedron * Cin,int working_space)1546 Param_Polyhedron *Polyhedron2Param_Vertices(Polyhedron *Din,Polyhedron *Cin,int working_space) {
1547
1548 Param_Polyhedron *result;
1549
1550 POL_ENSURE_FACETS(Din);
1551 POL_ENSURE_VERTICES(Din);
1552 POL_ENSURE_FACETS(Cin);
1553 POL_ENSURE_VERTICES(Cin);
1554
1555 #ifdef DEBUGPP
1556 fprintf(stderr,"Polyhedron2Param_Vertices algorithm starting at : %.2fs\n",
1557 (float)clock()/CLOCKS_PER_SEC);
1558 #endif
1559
1560 /***************** Scan the m-faces ****************/
1561 result = Find_m_faces(&Din,Cin,0,working_space,NULL,NULL);
1562
1563 #ifdef DEBUGPP
1564 fprintf(stderr, "nb of points : %d\n",result->nbV);
1565 #endif
1566
1567 #ifdef DEBUGPP
1568 fprintf(stderr, "end main loop : %.2fs\n", (float)clock()/CLOCKS_PER_SEC);
1569 #endif
1570
1571 return(result);
1572 } /* Polyhedron2Param_Vertices */
1573
1574 /*
1575 * Free the memory allocated to a list of parametrized vertices
1576 */
Param_Vertices_Free(Param_Vertices * PV)1577 void Param_Vertices_Free(Param_Vertices *PV) {
1578
1579 Param_Vertices *next_pv;
1580
1581 while(PV) {
1582 next_pv = PV->next;
1583 if (PV->Vertex) Matrix_Free(PV->Vertex);
1584 if (PV->Domain) Matrix_Free(PV->Domain);
1585 if (PV->Facets) free(PV->Facets);
1586 free(PV);
1587 PV = next_pv;
1588 }
1589 } /* Param_Vertices_Free */
1590
1591 /*
1592 * Print a list of parametrized vertices *
1593 */
Print_Vertex(FILE * DST,Matrix * V,const char ** param_names)1594 void Print_Vertex(FILE *DST, Matrix *V, const char **param_names)
1595 {
1596 int l, v;
1597 int first;
1598 Value gcd,tmp;
1599
1600 value_init(gcd); value_init(tmp);
1601
1602 fprintf(DST, "[" );
1603 for(l=0;l<V->NbRows;++l){
1604
1605 /* Variables */
1606 first=1;
1607 fprintf(DST, " " );
1608 for(v=0;v < V->NbColumns-2;++v) {
1609 if(value_notzero_p(V->p[l][v])) {
1610 value_gcd(gcd, V->p[l][v], V->p[l][V->NbColumns-1]);
1611 value_divexact(tmp, V->p[l][v], gcd);
1612 if(value_posz_p(tmp)) {
1613 if(!first)
1614 fprintf(DST, "+");
1615 if(value_notone_p(tmp)) {
1616 value_print(DST,VALUE_FMT,tmp);
1617 }
1618 }
1619 else { /* V->p[l][v]/gcd<0 */
1620 if(value_mone_p(tmp))
1621 fprintf(DST, "-" );
1622 else {
1623 value_print(DST,VALUE_FMT,tmp);
1624 }
1625 }
1626 value_divexact(tmp, V->p[l][V->NbColumns-1], gcd);
1627 if(value_notone_p(tmp)) {
1628 fprintf(DST, "%s/", param_names[v]);
1629 value_print(DST,VALUE_FMT,tmp);
1630 }
1631 else
1632 fprintf(DST, "%s", param_names[v]);
1633 first=0;
1634 }
1635 }
1636
1637 /* Constant */
1638 if(value_notzero_p(V->p[l][v]) || first) {
1639 if(value_posz_p(V->p[l][v]) && !first)
1640 fprintf(DST,"+");
1641 value_gcd(gcd, V->p[l][v], V->p[l][V->NbColumns-1]);
1642 value_divexact(tmp, V->p[l][v], gcd);
1643 value_print(DST,VALUE_FMT,tmp);
1644 value_divexact(tmp, V->p[l][V->NbColumns-1], gcd);
1645 if(value_notone_p(tmp)) {
1646 fprintf(DST,"/");
1647 value_print(DST,VALUE_FMT,tmp);
1648 fprintf(DST," ");
1649 }
1650 }
1651 if (l<V->NbRows-1)
1652 fprintf(DST, ", ");
1653 }
1654 fprintf(DST, " ]");
1655 value_clear(gcd); value_clear(tmp);
1656 return;
1657 } /* Print_Vertex */
1658
1659 /*----------------------------------------------------------------------*/
1660 /* VertexCT */
1661 /* convert a paramvertex from reduced space to normal m-space */
1662 /*----------------------------------------------------------------------*/
VertexCT(Matrix * V,Matrix * CT)1663 Matrix *VertexCT(Matrix *V,Matrix *CT) {
1664
1665 Matrix *Vt;
1666 int i,j,k;
1667
1668 if(CT) {
1669
1670 /* Have to transform the vertices to original dimension */
1671 Vt = Matrix_Alloc(V->NbRows,CT->NbColumns+1);
1672 for(i=0;i<V->NbRows;++i) {
1673 value_assign(Vt->p[i][CT->NbColumns],V->p[i][V->NbColumns-1]);
1674 for(j=0;j<CT->NbColumns;j++) {
1675 for(k=0;k<CT->NbRows;k++)
1676 if(value_notzero_p(CT->p[k][j]))
1677 break;
1678 if(k<CT->NbRows)
1679 value_assign(Vt->p[i][j],V->p[i][k]);
1680 else
1681 value_set_si(Vt->p[i][j],0);
1682 }
1683 }
1684 return(Vt);
1685 }
1686 else
1687 return(NULL);
1688 } /* VertexCT */
1689
1690 /*
1691 * Print the validity Domain 'D' of a parametric polyhedron
1692 */
Print_Domain(FILE * DST,Polyhedron * D,const char ** pname)1693 void Print_Domain(FILE *DST, Polyhedron *D, const char **pname)
1694 {
1695 int l, v;
1696 int first;
1697
1698 POL_ENSURE_FACETS(D);
1699 POL_ENSURE_VERTICES(D);
1700
1701 for(l=0;l<D->NbConstraints;++l) {
1702 fprintf(DST, " ");
1703 first = 1;
1704 for(v=1;v<=D->Dimension;++v) {
1705 if(value_notzero_p(D->Constraint[l][v])) {
1706 if(value_one_p(D->Constraint[l][v])) {
1707 if(first)
1708 fprintf(DST, "%s ", pname[v-1]);
1709 else
1710 fprintf(DST, "+ %s ", pname[v-1] );
1711 }
1712 else if(value_mone_p(D->Constraint[l][v]))
1713 fprintf(DST, "- %s ", pname[v-1] );
1714 else {
1715 if(value_pos_p(D->Constraint[l][v]) && !first )
1716 fprintf(DST, "+ " );
1717 value_print(DST,VALUE_FMT,D->Constraint[l][v]);
1718 fprintf(DST,"%s ",pname[v-1]);
1719 }
1720 first = 0;
1721 }
1722 }
1723 if(value_notzero_p(D->Constraint[l][v])) {
1724 if(value_pos_p(D->Constraint[l][v]) && !first)
1725 fprintf(DST,"+");
1726 fprintf(DST," ");
1727 value_print(DST,VALUE_FMT,D->Constraint[l][v]);
1728 }
1729 fprintf(DST,(value_notzero_p(D->Constraint[l][0])) ?" >= 0":" = 0");
1730 fprintf(DST, "\n" );
1731 }
1732 fprintf(DST, "\n");
1733
1734 if( D->next )
1735 {
1736 fprintf( DST, "UNION\n" );
1737 Print_Domain( DST, D->next, pname );
1738 }
1739 return;
1740 } /* Print_Domain */
1741
1742 /*
1743 * Given a list of parametrized vertices and an array of parameter names, Print
1744 * a list of parametrized vertices in a comprehensible format.
1745 */
Param_Vertices_Print(FILE * DST,Param_Vertices * PV,const char ** param_names)1746 void Param_Vertices_Print(FILE *DST, Param_Vertices *PV, const char **param_names)
1747 {
1748 Polyhedron *poly;
1749
1750 while(PV) {
1751 fprintf(DST, "Vertex :\n" );
1752 Print_Vertex(DST,PV->Vertex,param_names);
1753
1754 /* Pour le domaine : */
1755 fprintf(DST, " If :\n" );
1756 poly = Constraints2Polyhedron(PV->Domain,200);
1757 Print_Domain(DST,poly,param_names);
1758 Domain_Free(poly);
1759 PV = PV->next;
1760 }
1761 return;
1762 } /* Param_Vertices_Print */
1763
1764 /*
1765 * Given a polyhedron 'Din' in combined data and parametre space, a context
1766 * polyhedron 'Cin' representing the constraints on the parameter space and
1767 * a working space size 'working_space', return a parametric polyhedron with
1768 * a list of distinct validity domains and a complete list of valid vertices
1769 * associated to each validity domain.
1770 */
Polyhedron2Param_Domain(Polyhedron * Din,Polyhedron * Cin,int working_space)1771 Param_Polyhedron *Polyhedron2Param_Domain(Polyhedron *Din,Polyhedron *Cin,int working_space) {
1772
1773 Param_Polyhedron *result;
1774 Param_Domain *D;
1775
1776 POL_ENSURE_FACETS(Din);
1777 POL_ENSURE_VERTICES(Din);
1778 POL_ENSURE_FACETS(Cin);
1779 POL_ENSURE_VERTICES(Cin);
1780
1781 if (emptyQ(Din) || emptyQ(Cin))
1782 return NULL;
1783
1784 #ifdef DEBUGPP
1785 fprintf(stderr,"Polyhedron2Param_Polyhedron algorithm starting at : %.2fs\n",
1786 (float)clock()/CLOCKS_PER_SEC);
1787 #endif
1788
1789 /* Find the m-faces, keeping the corresponding domains */
1790 /* in the linked list PDomains */
1791 result = Find_m_faces(&Din,Cin,1,working_space,NULL,NULL);
1792
1793 #ifdef DEBUGPP
1794 if(result) fprintf(stderr, "Number of vertices : %d\n",result->nbV);
1795 fprintf(stderr,"Vertices found at : %.2fs\n",(float)clock()/CLOCKS_PER_SEC);
1796 #endif
1797
1798 /* Processing of PVResult and PDomains */
1799 if(result && Cin->Dimension>0) /* at least 1 parameter */
1800 Compute_PDomains(result->D,result->nbV,working_space);
1801 if(result && CEqualities)
1802 for(D=result->D;D;D=D->next)
1803 D->Domain = Add_CEqualities(D->Domain);
1804 Polyhedron_Free(CEqualities);
1805
1806 #ifdef DEBUGPP
1807 fprintf(stderr, "domains found at : %.2fs\n", (float)clock()/CLOCKS_PER_SEC);
1808 #endif
1809
1810 return(result);
1811 } /* Polyhedon2Param_Domain */
1812
1813 /*
1814 *
1815 */
Polyhedron2Param_SimplifiedDomain(Polyhedron ** Din,Polyhedron * Cin,int working_space,Polyhedron ** CEq,Matrix ** CT)1816 Param_Polyhedron *Polyhedron2Param_SimplifiedDomain(Polyhedron **Din,Polyhedron *Cin,int working_space,Polyhedron **CEq,Matrix **CT) {
1817
1818 Param_Polyhedron *result;
1819
1820 assert(CEq != NULL);
1821 assert(CT != NULL);
1822
1823 POL_ENSURE_FACETS(*Din);
1824 POL_ENSURE_VERTICES(*Din);
1825 POL_ENSURE_FACETS(Cin);
1826 POL_ENSURE_VERTICES(Cin);
1827
1828 #ifdef DEBUGPP
1829 fprintf(stderr,"Polyhedron2Param_Polyhedron algorithm starting at : %.2fs\n",
1830 (float)clock()/CLOCKS_PER_SEC);
1831 #endif
1832
1833 /* Find the m-faces, keeping the corresponding domains */
1834 /* in the linked list PDomains */
1835 result = Find_m_faces(Din,Cin,1,working_space,CEq,CT);
1836
1837 #ifdef DEBUGPP
1838 if(result) fprintf(stderr, "Number of vertices : %d\n",result->nbV);
1839 fprintf(stderr,"Vertices found at : %.2fs\n",(float)clock()/CLOCKS_PER_SEC);
1840 #endif
1841
1842 /* Processing of PVResult and PDomains */
1843 if(result && Cin->Dimension>0) /* at least 1 parameter */
1844 Compute_PDomains(result->D,result->nbV,working_space);
1845
1846 /* Removed this, Vin100, March 01 */
1847 /* if(result && CEqualities )
1848 for(D=result->D;D;D=D->next)
1849 D->Domain = Add_CEqualities(D->Domain);
1850 */
1851
1852 #ifdef DEBUGPP
1853 fprintf(stderr, "domains found at : %.2fs\n", (float)clock()/CLOCKS_PER_SEC);
1854 #endif
1855
1856 return(result);
1857 } /* Polyhedron2Param_SimplifiedDomain */
1858
1859 /*
1860 * Free the memory allocated to a list of validity domain of a parametrized
1861 * polyhedron.
1862 */
Param_Domain_Free(Param_Domain * PD)1863 void Param_Domain_Free(Param_Domain *PD) {
1864
1865 Param_Domain *next_pd;
1866
1867 while(PD) {
1868 free(PD->F);
1869 Domain_Free(PD->Domain);
1870 next_pd = PD->next;
1871 free(PD);
1872 PD = next_pd;
1873 }
1874 return;
1875 } /* Param_Domain_Free */
1876
1877 /*
1878 * Free the memory allocated to a parametric polyhedron 'P'
1879 */
Param_Polyhedron_Free(Param_Polyhedron * P)1880 void Param_Polyhedron_Free(Param_Polyhedron *P) {
1881
1882 if (!P) return;
1883 Param_Vertices_Free(P->V);
1884 Param_Domain_Free(P->D);
1885 if (P->Constraints)
1886 Matrix_Free(P->Constraints);
1887 if (P->Rays)
1888 Matrix_Free(P->Rays);
1889 free(P);
1890 return;
1891 } /* Param_Polyhedron_Free */
1892
1893 /*
1894 * Scales the parametric polyhedron such that all vertices are integer.
1895 */
Param_Polyhedron_Scale_Integer(Param_Polyhedron * PP,Polyhedron ** P,Value * det,unsigned MaxRays)1896 void Param_Polyhedron_Scale_Integer(Param_Polyhedron *PP, Polyhedron **P,
1897 Value *det, unsigned MaxRays)
1898 {
1899 int i;
1900 int nb_param, nb_vars;
1901 Vector *denoms;
1902 Param_Vertices *V;
1903 Value global_var_lcm;
1904 Matrix *expansion;
1905
1906 value_set_si(*det, 1);
1907 if (!PP->nbV)
1908 return;
1909
1910 nb_param = PP->D->Domain->Dimension;
1911 nb_vars = PP->V->Vertex->NbRows;
1912
1913 /* Scan the vertices and make an orthogonal expansion of the variable
1914 space */
1915 /* a- prepare the array of common denominators */
1916 denoms = Vector_Alloc(nb_vars);
1917 value_init(global_var_lcm);
1918
1919 /* b- scan the vertices and compute the variables' global lcms */
1920 for (V = PP->V; V; V = V->next)
1921 for (i = 0; i < nb_vars; i++)
1922 value_lcm(denoms->p[i], denoms->p[i], V->Vertex->p[i][nb_param+1]);
1923
1924 value_set_si(global_var_lcm, 1);
1925 for (i = 0; i < nb_vars; i++) {
1926 value_multiply(*det, *det, denoms->p[i]);
1927 value_lcm(global_var_lcm, global_var_lcm, denoms->p[i]);
1928 }
1929
1930 /* scale vertices */
1931 for (V = PP->V; V; V = V->next)
1932 for (i = 0; i < nb_vars; i++) {
1933 Vector_Scale(V->Vertex->p[i], V->Vertex->p[i], denoms->p[i], nb_param+1);
1934 Vector_Normalize(V->Vertex->p[i], nb_param+2);
1935 }
1936
1937 /* the expansion can be actually writen as global_var_lcm.L^{-1} */
1938 /* this is equivalent to multiply the rows of P by denoms_det */
1939 for (i = 0; i < nb_vars; i++)
1940 value_division(denoms->p[i], global_var_lcm, denoms->p[i]);
1941
1942 /* OPT : we could use a vector instead of a diagonal matrix here (c- and d-).*/
1943 /* c- make the quick expansion matrix */
1944 expansion = Matrix_Alloc(nb_vars+nb_param+1, nb_vars+nb_param+1);
1945 for (i = 0; i < nb_vars; i++)
1946 value_assign(expansion->p[i][i], denoms->p[i]);
1947 for (i = nb_vars; i < nb_vars+nb_param+1; i++)
1948 value_assign(expansion->p[i][i], global_var_lcm);
1949
1950 /* d- apply the variable expansion to the polyhedron */
1951 if (P)
1952 *P = Polyhedron_Preimage(*P, expansion, MaxRays);
1953
1954 Matrix_Free(expansion);
1955 value_clear(global_var_lcm);
1956 Vector_Free(denoms);
1957 }
1958