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