00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include <stdio.h>
00027 #include <stdlib.h>
00028 #include <string.h>
00029 #include <assert.h>
00030 #ifdef DEBUGPP
00031 #include <time.h>
00032 #endif
00033
00034 #include <polylib/polylib.h>
00035
00036 static void traite_m_face(Polyhedron *, unsigned int *, unsigned int *);
00037 static void scan_m_face(int,int,Polyhedron *,unsigned int *);
00038
00039
00040
00041
00042
00043
00044
00045 Polyhedron *PDomainIntersection(Polyhedron *Pol1,Polyhedron *Pol2,unsigned NbMaxRays) {
00046
00047 Polyhedron *p1, *p2, *p3, *d;
00048
00049 if (!Pol1 || !Pol2) return (Polyhedron*) 0;
00050 if((Pol1->Dimension != Pol2->Dimension) || (Pol1->NbEq != Pol2->NbEq)) {
00051 fprintf(stderr,
00052 "? PDomainIntersection: operation on different dimensions\n");
00053 return (Polyhedron*) 0;
00054 }
00055
00056 POL_ENSURE_FACETS(Pol1);
00057 POL_ENSURE_VERTICES(Pol1);
00058 POL_ENSURE_FACETS(Pol2);
00059 POL_ENSURE_VERTICES(Pol2);
00060
00061 d = (Polyhedron *)0;
00062 for (p1=Pol1; p1; p1=p1->next) {
00063 for (p2=Pol2; p2; p2=p2->next) {
00064 p3 = AddConstraints(p2->Constraint[0],
00065 p2->NbConstraints,p1,NbMaxRays);
00066 if (!p3) continue;
00067
00068
00069 if (p3->NbEq!=Pol1->NbEq)
00070 Polyhedron_Free(p3) ;
00071
00072
00073 else
00074 d = AddPolyToDomain(p3,d);
00075 }
00076 }
00077 return d;
00078 }
00079
00080
00081
00082
00083
00084
00085
00086 Polyhedron *PDomainDifference(Polyhedron *Pol1,Polyhedron *Pol2,unsigned NbMaxRays) {
00087
00088 Polyhedron *p1, *p2, *p3, *d;
00089 int i;
00090
00091 if (!Pol1 || !Pol2)
00092 return (Polyhedron*) 0;
00093 if((Pol1->Dimension != Pol2->Dimension) || (Pol1->NbEq != Pol2->NbEq)) {
00094 fprintf(stderr,
00095 "? PDomainDifference: operation on different dimensions\n");
00096 return (Polyhedron*) 0;
00097 }
00098
00099 POL_ENSURE_FACETS(Pol1);
00100 POL_ENSURE_VERTICES(Pol1);
00101 POL_ENSURE_FACETS(Pol2);
00102 POL_ENSURE_VERTICES(Pol2);
00103
00104 d = (Polyhedron *)0;
00105 for (p2=Pol2; p2; p2=p2->next) {
00106 for (p1=Pol1; p1; p1=p1->next) {
00107 for (i=0; i<p2->NbConstraints; i++) {
00108
00109
00110 p3 = SubConstraint(p2->Constraint[i],p1,NbMaxRays,2);
00111 if (!p3) continue;
00112
00113
00114
00115 if (emptyQ(p3) || p3->NbEq!=Pol1->NbEq)
00116 Polyhedron_Free(p3);
00117
00118
00119 else
00120 d = AddPolyToDomain(p3,d);
00121 }
00122 }
00123 if (p2 != Pol2)
00124 Domain_Free(Pol1);
00125 Pol1 = d;
00126 d = (Polyhedron *)0;
00127 }
00128 return Pol1;
00129 }
00130
00131
00132
00133
00134 static int TestRank(Matrix *Mat) {
00135
00136 int i,j,k;
00137 Value m1,m2,m3,gcd,tmp;
00138
00139
00140 value_init(m1); value_init(m2);
00141 value_init(m3); value_init(gcd); value_init(tmp);
00142
00143 for(k=0;k<Mat->NbColumns;++k) {
00144
00145
00146
00147 if(value_zero_p(Mat->p[k][k])) {
00148 for(j=k+1;j<Mat->NbRows;++j) {
00149
00150
00151 if(value_notzero_p(Mat->p[j][k])) {
00152
00153
00154 for(i=k;i<Mat->NbColumns;++i) {
00155 value_assign(tmp,Mat->p[j][i]);
00156 value_assign(Mat->p[j][i],Mat->p[k][i]);
00157 value_assign(Mat->p[k][i],tmp);
00158 }
00159 break;
00160 }
00161 }
00162
00163
00164
00165 if(j>=Mat->NbRows) {
00166
00167
00168 value_clear(m1); value_clear(m2);
00169 value_clear(m3); value_clear(gcd); value_clear(tmp);
00170 return 0;
00171 }
00172 }
00173
00174
00175 for(j=k+1;j<Mat->NbRows;++j) {
00176
00177
00178 value_gcd(gcd, Mat->p[j][k], Mat->p[k][k]);
00179 for(i=k+1;i<Mat->NbColumns;++i) {
00180
00181
00182 value_multiply(m1,Mat->p[j][i],Mat->p[k][k]);
00183 value_multiply(m2,Mat->p[j][k],Mat->p[k][i]);
00184 value_subtract(m3,m1,m2);
00185 value_division(Mat->p[j][i],m3,gcd);
00186 }
00187 }
00188 }
00189
00190
00191 value_clear(m1); value_clear(m2);
00192 value_clear(m3); value_clear(gcd); value_clear(tmp);
00193
00194
00195 return 1;
00196 }
00197
00198
00199
00200
00201
00202
00203
00204
00205 typedef struct {
00206 unsigned int NbRows;
00207 unsigned int NbColumns;
00208 unsigned int **p;
00209 unsigned int *p_init;
00210 } SatMatrix;
00211
00212 static SatMatrix *SMAlloc(int rows,int cols) {
00213
00214 unsigned int **q, *p;
00215 int i;
00216
00217 SatMatrix *result = (SatMatrix *)malloc(sizeof(SatMatrix));
00218 assert (result != NULL);
00219
00220 result->NbRows = rows;
00221 result->NbColumns = cols;
00222 result->p = q = (unsigned int **)malloc(rows * sizeof(unsigned int *));
00223 assert (result->p != NULL);
00224 result->p_init = p = (unsigned int *)malloc(rows * cols * sizeof(unsigned int));
00225 assert (result->p_init != NULL);
00226
00227 for (i=0;i<rows;i++) {
00228 *q++ = p;
00229 p += cols;
00230 }
00231
00232 return result;
00233 }
00234
00235 static void SMPrint (SatMatrix *matrix) {
00236
00237 unsigned int *p;
00238 int i, j;
00239 unsigned NbRows, NbColumns;
00240
00241 fprintf(stderr,"%d %d\n",NbRows=matrix->NbRows, NbColumns=matrix->NbColumns);
00242 for (i=0;i<NbRows;i++) {
00243 p = *(matrix->p+i);
00244 for (j=0;j<NbColumns;j++)
00245 fprintf(stderr, " %10X ", *p++);
00246 fprintf(stderr, "\n");
00247 }
00248 }
00249
00250
00251 static void SMFree (SatMatrix *matrix) {
00252
00253 free ((char *) matrix->p_init);
00254 free ((char *) matrix->p);
00255 free ((char *) matrix);
00256 return;
00257 }
00258
00259
00260
00261
00262
00263
00264
00265 static int m;
00266 static int m_dim;
00267 static int n;
00268 static int ws;
00269 static int nr;
00270
00271 static Polyhedron *CEqualities;
00272 static SatMatrix *Sat;
00273 static unsigned int *egalite;
00274 static Matrix *Xi, *Pi;
00275 static Matrix *PiTest;
00276 static Matrix *CTest;
00277 static Matrix *PiInv;
00278
00279 static Matrix *RaysDi;
00280
00281 static int KD;
00282
00283
00284 static int nbPV;
00285 static Param_Vertices *PV_Result;
00286 static Param_Domain *PDomains;
00287
00288 #ifdef DEBUGPP
00289 static int nbfaces;
00290 #endif
00291
00292
00293
00294
00295
00296
00297 static Polyhedron *Add_CEqualities(Polyhedron *D) {
00298
00299 Polyhedron *d,*r,*tmp;
00300
00301 if(!CEqualities)
00302 return D;
00303 else {
00304 if(!D || emptyQ(D)) {
00305 if(D)
00306 Domain_Free(D);
00307 return(Polyhedron_Copy(CEqualities));
00308 }
00309 r = AddConstraints(D->Constraint[0],D->NbConstraints,
00310 CEqualities,ws);
00311 tmp = r;
00312 for(d=D->next;d;d=d->next) {
00313 tmp->next = AddConstraints(d->Constraint[0],d->NbConstraints,
00314 CEqualities,ws);
00315 tmp = tmp->next;
00316 }
00317 Domain_Free(D);
00318 return(r);
00319 }
00320 }
00321
00322 #define INT_BITS (sizeof(unsigned) * 8)
00323
00324 unsigned int *int_array2bit_vector(unsigned int *array, int n)
00325 {
00326 int i, ix;
00327 unsigned bx;
00328 int words = (n+INT_BITS-1)/INT_BITS;
00329 unsigned int *bv = (unsigned int *)calloc(words, sizeof(unsigned));
00330 assert(bv);
00331 for (i = 0, ix = 0, bx = MSB; i < n; ++i) {
00332 if (array[i])
00333 bv[ix] |= bx;
00334 NEXT(ix, bx);
00335 }
00336 return bv;
00337 }
00338
00339
00340
00341
00342
00343
00344
00345
00346 static void traite_m_face(Polyhedron *D, unsigned int *mf,
00347 unsigned int *egalite)
00348 {
00349 Matrix *Si;
00350 Polyhedron *PDi;
00351 Param_Vertices *PV;
00352 int j,k,c,r;
00353 unsigned kx, bx;
00354
00355 #ifdef DEBUGPP
00356 ++nbfaces;
00357 #endif
00358
00359
00360 RaysDi->NbRows = 0;
00361 for(k=0,c=0,kx=0,bx=MSB;k<D->NbRays;++k) {
00362 if(mf[kx]&bx) {
00363 if(c<m+1) {
00364 int i;
00365
00366
00367
00368
00369
00370
00371 for(j=0;j<m+1;++j) {
00372 for(i=0;i<c;++i)
00373
00374
00375 value_assign(PiTest->p[j][i],Pi->p[j][i]);
00376
00377
00378 value_assign(PiTest->p[j][c],D->Ray[k][j+1+n]);
00379 }
00380 PiTest->NbColumns = c+1;
00381 r = TestRank(PiTest);
00382 if(r ) {
00383
00384
00385 for (j=0;j<n;j++)
00386 value_assign(Xi->p[j][c],D->Ray[k][j+1]);
00387 for (j=0;j<m;j++)
00388 value_assign(Pi->p[j][c],D->Ray[k][j+1+n]);
00389 value_assign(Xi->p[n][c],D->Ray[k][n+m+1]);
00390 value_assign(Pi->p[m][c],D->Ray[k][n+m+1]);
00391 c++;
00392 }
00393 }
00394
00395
00396 value_assign(RaysDi->p[RaysDi->NbRows][0],D->Ray[k][0]);
00397 Vector_Copy(&D->Ray[k][n+1],&RaysDi->p[RaysDi->NbRows][1],(m+1));
00398 ++RaysDi->NbRows;
00399 }
00400 NEXT(kx,bx);
00401 }
00402
00403 #ifdef DEBUGPP41
00404 fprintf(stderr, "\nRaysDi=\n");
00405 Matrix_Print(stderr,P_VALUE_FMT,RaysDi);
00406 if(c < m+1)
00407 fprintf(stderr, "Invalid ");
00408 fprintf(stderr, "Pi=\n");
00409 Matrix_Print(stderr,P_VALUE_FMT,Pi);
00410 #endif
00411
00412 #ifdef DEBUGPP4
00413 if(c < m+1)
00414 fprintf(stderr,"Eliminated because of no vertex\n");
00415 #endif
00416
00417 if(c < m+1)
00418 return;
00419
00420
00421
00422
00423
00424
00425
00426
00427 #ifdef DEBUGPP4
00428 fprintf(stderr,"Xi = ");
00429 Matrix_Print(stderr,P_VALUE_FMT,Xi);
00430 fprintf(stderr,"Pi = ");
00431 Matrix_Print(stderr,P_VALUE_FMT,Pi);
00432 #endif
00433
00434
00435
00436 if(!MatInverse(Pi,PiInv)) {
00437
00438 #ifdef DEBUGPP4
00439 fprintf(stderr, "Eliminated because of no inverse Pi\n");
00440 #endif
00441
00442 return;
00443 }
00444
00445 #ifdef DEBUGPP4
00446 fprintf(stderr,"FACE GENERATED!\n");
00447 fprintf(stderr,"PiInv = ");
00448 Matrix_Print(stderr,P_VALUE_FMT,PiInv);
00449 #endif
00450
00451
00452 Si = Matrix_Alloc(Xi->NbRows,PiInv->NbColumns);
00453 rat_prodmat(Si,Xi,PiInv);
00454
00455 #ifdef DEBUGPP4
00456 fprintf(stderr,"Si = ");
00457 Matrix_Print(stderr,P_VALUE_FMT,Si);
00458 #endif
00459
00460 Si->NbRows--;
00461
00462
00463 PV = (Param_Vertices *) malloc(sizeof(Param_Vertices));
00464 PV->next = PV_Result;
00465 PV->Vertex = Si;
00466 PV->Domain = NULL;
00467 PV->Facets = int_array2bit_vector(egalite, D->NbConstraints);
00468 PV_Result = PV;
00469 nbPV++;
00470
00471
00472 PDi = Rays2Polyhedron(RaysDi,ws);
00473
00474 #ifdef DEBUGPP3
00475 fprintf(stderr,"RaysDi = ");
00476 Matrix_Print(stderr,P_VALUE_FMT,RaysDi);
00477 fprintf(stderr,"PDi = ");
00478 Polyhedron_Print(stderr,P_VALUE_FMT,PDi);
00479 #endif
00480
00481 if(KD==0) {
00482
00483
00484 PDi = Add_CEqualities(PDi);
00485 PV->Domain = Polyhedron2Constraints(PDi);
00486 Polyhedron_Free(PDi);
00487 }
00488 else {
00489 Param_Domain *PD;
00490 PD = (Param_Domain *) malloc(sizeof(Param_Domain));
00491 PD->Domain = PDi;
00492 PD->F = NULL;
00493 PD->next = PDomains;
00494 PDomains = PD;
00495 }
00496 return;
00497 }
00498
00499
00500
00501
00502
00503
00504 int cntbit[256] = {
00505 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
00506 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
00507 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
00508 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
00509
00510 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
00511 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
00512 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
00513 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
00514
00515 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
00516 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
00517 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
00518 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
00519
00520 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
00521 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
00522 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
00523 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8 };
00524
00525 static int count_sat (unsigned int *mf) {
00526
00527 register unsigned int i, tmp, cnt=0;
00528
00529 for (i=0; i<nr; i++) {
00530 tmp = mf[i];
00531 cnt = cnt
00532 + cntbit[ tmp & 0xff ]
00533 + cntbit[ (tmp>>8) & 0xff ]
00534 + cntbit[ (tmp>>16) & 0xff ]
00535 + cntbit[ (tmp>>24) & 0xff ]
00536 ;
00537 }
00538 return cnt;
00539 }
00540
00541
00542 static int bit_vector_includes(unsigned int *bv, int len, unsigned int *part)
00543 {
00544 int j;
00545
00546 for (j = 0; j < len; j++) {
00547 #ifdef DEBUGPP4
00548 fprintf(stderr, "mf=%08X Sat=%08X &=%08X\n",
00549 part[j],bv[j], (part[j] & bv[j]));
00550 #endif
00551 if ((part[j] & bv[j]) != part[j])
00552 return 0;
00553 }
00554 return 1;
00555 }
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598 static void scan_m_face(int pos,int nb_un,Polyhedron *D,unsigned int *mf) {
00599
00600
00601
00602
00603
00604 unsigned int *new_mf;
00605
00606 #ifdef DEBUGPP4
00607 fprintf(stderr,"Start scan_m_face(pos=%d, nb_un=%d, n=%d, m=%d\n",
00608 pos,nb_un,n,m);
00609 fprintf(stderr,"mf = ");
00610 {
00611 int i;
00612 for(i=0;i<nr;i++)
00613 fprintf(stderr,"%08X", mf[i]);
00614 fprintf(stderr,"\nequality = [");
00615 for(i=0;i<D->NbConstraints;i++)
00616 fprintf(stderr," %1d",egalite[i]);
00617 fprintf(stderr,"]\n");
00618 }
00619 #endif
00620
00621 if(nb_un == 0) {
00622 int i;
00623
00624
00625
00626
00627
00628 for(i=0;i<pos-1;i++)
00629 {
00630 if(egalite[i])
00631 continue;
00632
00633
00634 #ifdef DEBUGPP4
00635 fprintf(stderr, "Sat[%d]\n", i);
00636 #endif
00637 if (bit_vector_includes(Sat->p[i], nr, mf)) {
00638 #ifdef DEBUGPP4
00639 fprintf(stderr, "Redundant with constraint %d\n", i);
00640 #endif
00641 return;
00642 }
00643 }
00644
00645
00646 for (i = pos; i < D->NbConstraints; ++i) {
00647 if (bit_vector_includes(Sat->p[i], nr, mf))
00648 egalite[i] = 1;
00649 }
00650
00651
00652 traite_m_face(D, mf, egalite);
00653 for (i = pos; i < D->NbConstraints; ++i)
00654 egalite[i] = 0;
00655 return;
00656 }
00657
00658
00659 if((pos+nb_un)>D->NbConstraints) return;
00660
00661
00662
00663 {
00664 int k;
00665 new_mf = (unsigned int *)malloc(nr*sizeof(unsigned int));
00666 for (k=0; k<nr; k++)
00667 new_mf[k] = mf[k] & Sat->p[pos][k];
00668 }
00669 #ifdef DEBUGPP4
00670 fprintf(stderr,"new_mf = ");
00671 {
00672 int i;
00673 for(i=0;i<nr;i++) {
00674 fprintf(stderr,"%08X", new_mf[i]);
00675 }
00676 fprintf(stderr,"\ncount(new_mf) = %d\n",count_sat(new_mf));
00677 }
00678 #endif
00679
00680 {
00681 int c;
00682 c = count_sat(new_mf);
00683
00684 if (c>m_dim )
00685 {
00686 int redundant = 0;
00687
00688 egalite[pos]=1;
00689
00690
00691
00692
00693
00694
00695
00696 if( c==count_sat(mf) ) {
00697 int i, c, j;
00698
00699 for (i = 0, c = 0; i < D->NbConstraints; ++i) {
00700 if (egalite[i] == 0 || egalite[i] == -1)
00701 continue;
00702 for (j = 0; j < D->Dimension+1; ++j)
00703 value_assign(CTest->p[j][c],
00704 D->Constraint[i][j+1]);
00705 ++c;
00706 }
00707 CTest->NbColumns = c;
00708 #ifdef DEBUGPP41
00709 Matrix_Print(stderr,P_VALUE_FMT,CTest);
00710 #endif
00711 redundant = !TestRank(CTest);
00712 }
00713
00714
00715 if( redundant )
00716 {
00717 egalite[pos]=-1;
00718
00719 scan_m_face(pos+1,nb_un,D,new_mf);
00720 }
00721 else
00722 {
00723 scan_m_face(pos+1,nb_un-1,D,new_mf);
00724 }
00725 }
00726 }
00727 free(new_mf);
00728 egalite[pos]=0;
00729 if ((pos+nb_un)>=D->NbConstraints) return;
00730 scan_m_face(pos+1,nb_un,D,mf);
00731 return;
00732 }
00733
00734
00735
00736
00737
00738
00739 static SatMatrix *Poly2Sat(Polyhedron *Pol,unsigned int **L) {
00740
00741 SatMatrix *Sat;
00742 int i, j, k, kx;
00743 unsigned int *Temp;
00744 Value *p1, *p2, p3,tmp;
00745 unsigned Dimension, NbRay, NbCon, bx;
00746
00747
00748 value_init(p3); value_init(tmp);
00749
00750 NbRay = Pol->NbRays;
00751 NbCon = Pol->NbConstraints;
00752 Dimension = Pol->Dimension+1;
00753
00754
00755 nr = (NbRay - 1)/(sizeof(int)*8) + 1;
00756 Sat = SMAlloc(NbCon,nr);
00757 Temp = (unsigned int *)malloc(nr*sizeof(unsigned int));
00758 memset(Sat->p_init,0,nr*NbCon*sizeof(int));
00759 memset(Temp,0,nr*sizeof(unsigned int));
00760 kx=0; bx=MSB;
00761 for (k=0; k<NbRay; k++) {
00762 for (i=0; i<NbCon; i++) {
00763 p1 = &Pol->Constraint[i][1];
00764 p2 = &Pol->Ray[k][1];
00765 value_set_si(p3,0);
00766 for (j=0;j<Dimension;j++) {
00767 value_multiply(tmp,*p1,*p2);
00768 value_addto(p3,p3,tmp);
00769 p1++; p2++;
00770 }
00771 if (value_zero_p(p3))
00772 Sat->p[i][kx]|=bx;
00773 }
00774 Temp[kx] |= bx;
00775 NEXT(kx, bx);
00776 }
00777
00778
00779
00780 *L = Temp;
00781
00782
00783 value_clear(p3); value_clear(tmp);
00784
00785 return Sat;
00786 }
00787
00788
00789
00790
00791
00792 Param_Polyhedron *GenParamPolyhedron(Polyhedron *Pol, Matrix *Rays) {
00793
00794 Param_Polyhedron *result;
00795 int nbRows, nbColumns;
00796 int i, size, rays;
00797
00798 nbRows=Pol->NbRays;
00799 nbColumns=Pol->Dimension+2;
00800
00801
00802 for(i=0, rays=0; i<nbRows; i++)
00803 if(value_notone_p(Pol->Ray[i][0]) ||
00804 value_zero_p(Pol->Ray[i][nbColumns-1]))
00805 ++rays;
00806
00807
00808 result=(Param_Polyhedron *)malloc(sizeof(Param_Polyhedron));
00809 result->nbV=nbRows-rays;
00810 result->V=NULL;
00811 result->Constraints = Polyhedron2Constraints(Pol);
00812 result->Rays = Rays;
00813
00814
00815 for(i=0;i<nbRows;i++) {
00816 Matrix *vertex;
00817 Param_Vertices *paramVertex;
00818 int j;
00819
00820 if (value_notone_p(Pol->Ray[i][0]) ||
00821 value_zero_p(Pol->Ray[i][nbColumns-1]))
00822 continue;
00823
00824 vertex=Matrix_Alloc(nbColumns-2,2);
00825 for(j=1;j<nbColumns-1;j++) {
00826 value_assign(vertex->p[j-1][0],Pol->Ray[i][j]);
00827 value_assign(vertex->p[j-1][1],Pol->Ray[i][nbColumns-1]);
00828 }
00829 paramVertex=(Param_Vertices *)malloc(sizeof(Param_Vertices));
00830 paramVertex->Vertex=vertex;
00831
00832
00833 paramVertex->Domain=Matrix_Alloc(1,2);
00834 value_set_si(paramVertex->Domain->p[0][0],1);
00835 value_set_si(paramVertex->Domain->p[0][1],1);
00836 paramVertex->Facets = NULL;
00837 paramVertex->next=result->V;
00838 result->V=paramVertex;
00839 }
00840
00841
00842 if (nbRows > 1)
00843 size=(nbRows-1)/(8*sizeof(int))+1;
00844 else
00845 size = 1;
00846 result->D=(Param_Domain *)malloc(sizeof(Param_Domain));
00847 result->D->next=NULL;
00848 result->D->Domain=Universe_Polyhedron(0);
00849 result->D->F=(unsigned int *)malloc(size*sizeof(int));
00850 memset(&result->D->F[0],0xFF,size*sizeof(int));
00851
00852 return result;
00853 }
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865 Matrix *PreElim_Columns(Polyhedron *E,int *p,int *ref,int m) {
00866
00867 int i,j,l;
00868 Matrix *T;
00869
00870
00871
00872
00873
00874 memset(p,0,sizeof(int)*(E->NbEq));
00875
00876 #ifdef DEBUGPP32
00877 fprintf(stderr,"\n\nPreElim_Columns starting\n");
00878 fprintf(stderr,"E =\n");
00879 Polyhedron_Print(stderr,P_VALUE_FMT,E);
00880 #endif
00881
00882 for(l=0;l<E->NbEq;++l) {
00883 if(value_notzero_p(E->Constraint[l][0])) {
00884 fprintf(stderr,"Internal error: Elim_Columns (polyparam.c), expecting equalities first in E.\n");
00885 free(p);
00886 return(NULL);
00887 }
00888 for(i=1;value_zero_p(E->Constraint[l][i]);++i) {
00889 if(i==E->Dimension+1) {
00890 fprintf(stderr,"Internal error: Elim_Columns (polyparam.c), expecting non-empty constraint in E.\n");
00891 free(p);
00892 return( NULL );
00893 }
00894 }
00895 p[l] = i;
00896
00897 #ifdef DEBUGPP32
00898 fprintf(stderr,"p[%d] = %d,",l,p[l]);
00899 #endif
00900 }
00901
00902
00903 for(i=0;i<E->Dimension+2-E->NbEq;++i) {
00904 ref[i]=i;
00905 for(j=0;j<E->NbEq;++j)
00906 if(p[j]<=ref[i])
00907 ref[i]++;
00908
00909 #ifdef DEBUGPP32
00910 fprintf(stderr,"ref[%d] = %d,",i,ref[i]);
00911 #endif
00912 }
00913
00914
00915 T = Matrix_Alloc(m+1-E->NbEq,m+1);
00916 for(i=0;i<T->NbColumns;i++)
00917 for(j=0;j<T->NbRows;j++)
00918 if(ref[E->Dimension-m+j+1] == E->Dimension-m+i+1)
00919 value_set_si(T->p[j][i],1);
00920 else
00921 value_set_si(T->p[j][i],0);
00922
00923 #ifdef DEBUGPP32
00924 fprintf(stderr,"\nTransMatrix =\n");
00925 Matrix_Print(stderr,P_VALUE_FMT,T);
00926 #endif
00927
00928 return(T);
00929
00930 }
00931
00932
00933
00934
00935
00936
00937
00938 Polyhedron *Elim_Columns(Polyhedron *A,Polyhedron *E,int *p,int *ref) {
00939
00940 int i,l,c;
00941 Matrix *M, *C;
00942 Polyhedron *R;
00943 Value tmp1,tmp2;
00944
00945
00946 value_init(tmp1); value_init(tmp2);
00947
00948 #ifdef DEBUGPP32
00949 fprintf(stderr,"\nElim_Columns starting\n");
00950 fprintf(stderr,"A =\n" );
00951 Polyhedron_Print(stderr,P_VALUE_FMT,A);
00952 #endif
00953
00954
00955 M = Polyhedron2Constraints(A);
00956 for(l=0;l<E->NbEq;++l) {
00957 for(c=0;c<M->NbRows;++c) {
00958 if(value_notzero_p(M->p[c][p[l]])) {
00959
00960
00961 for(i=1;i<M->NbColumns;++i) {
00962 value_multiply(tmp1,M->p[c][i],E->Constraint[l][p[l]]);
00963 value_multiply(tmp2,E->Constraint[l][i],M->p[c][p[l]]);
00964 value_subtract(M->p[c][i],tmp1,tmp2);
00965 }
00966 }
00967 }
00968 }
00969
00970 #ifdef DEBUGPP32
00971 fprintf(stderr,"\nElim_Columns after zeroing columns of A.\n");
00972 fprintf(stderr,"M =\n");
00973 Matrix_Print(stderr,P_VALUE_FMT,M);
00974 #endif
00975
00976
00977 C = Matrix_Alloc(M->NbRows,M->NbColumns - E->NbEq);
00978 for(l=0;l<C->NbRows;++l)
00979 for(c=0;c<C->NbColumns;++c) {
00980 value_assign(C->p[l][c],M->p[l][ref[c]]);
00981 }
00982
00983 #ifdef DEBUGPP32
00984 fprintf(stderr,"\nElim_Columns after eliminating columns of A.\n");
00985 fprintf(stderr,"C =\n");
00986 Matrix_Print(stderr,P_VALUE_FMT,C);
00987 #endif
00988
00989 R = Constraints2Polyhedron(C,ws);
00990 Matrix_Free(C);
00991 Matrix_Free(M);
00992 value_clear(tmp1); value_clear(tmp2);
00993 return(R);
00994 }
00995
00996
00997
00998
00999
01000
01001
01002 Param_Polyhedron *Find_m_faces(Polyhedron **Di,Polyhedron *C,int keep_dom,int working_space,Polyhedron **CEq,Matrix **CT) {
01003
01004 unsigned int *mf;
01005 int i, j, c1r, nbr, dimfaces;
01006 Polyhedron *D=*Di,
01007 *C1,
01008 *D1,
01009 *npP;
01010 Matrix *M;
01011 Matrix *Rays;
01012 Param_Polyhedron *res;
01013 int *p, *ref;
01014 Value tmp;
01015
01016 CEqualities = NULL;
01017
01018 if(CT) {
01019 *CEq = NULL;
01020 *CT = NULL;
01021 }
01022
01023 if(!D || !C)
01024 return (Param_Polyhedron *) 0;
01025
01026 ws = working_space;
01027 m = C->Dimension;
01028 n = D->Dimension - m;
01029 if(n<0) {
01030 fprintf(stderr,
01031 "Find_m_faces: ?%d parameters of a %d-polyhedron !\n",m,n);
01032 return (Param_Polyhedron *) 0;
01033 }
01034 if (m==0)
01035 return GenParamPolyhedron(D,Matrix_Alloc(0,2));
01036
01037
01038 C1 = align_context(C,D->Dimension,ws);
01039
01040 #ifdef DEBUGPP31
01041 fprintf(stderr,"m = %d\n",m);
01042 fprintf(stderr, "D = ");
01043 Polyhedron_Print(stderr,P_VALUE_FMT,D);
01044 fprintf(stderr,"C1 = ");
01045 Polyhedron_Print(stderr,P_VALUE_FMT,C1);
01046 #endif
01047
01048 D1 = DomainIntersection(D,C1,ws);
01049
01050 #ifdef DEBUGPP31
01051 fprintf(stderr,"D1 = ");
01052 Polyhedron_Print(stderr,P_VALUE_FMT,D1);
01053 #endif
01054
01055 Domain_Free(C1);
01056
01057 if (!D1)
01058 return(NULL);
01059 if (emptyQ(D1)) {
01060 Polyhedron_Free(D1);
01061 return(NULL);
01062 }
01063
01064
01065
01066 M = Matrix_Alloc(n, D1->Dimension+2);
01067 for (i=0; i<n; i++)
01068 value_set_si(M->p[i][i+1],1);
01069 C1 = DomainAddRays(D1,M,ws);
01070 Matrix_Free(M);
01071
01072 #ifdef DEBUGPP31
01073 fprintf(stderr,"True context C1 = ");
01074 Polyhedron_Print(stderr,P_VALUE_FMT,C1);
01075 #endif
01076
01077
01078
01079
01080
01081 for(c1r=0 ; c1r<C1->NbRays ; c1r++)
01082 if( value_notzero_p(C1->Ray[c1r][0]) &&
01083 value_notzero_p(C1->Ray[c1r][C1->Dimension+1]) )
01084 break;
01085 assert(c1r != C1->NbRays);
01086 value_init(tmp);
01087 M = Matrix_Alloc(D1->NbConstraints, n+2);
01088 for(i=0 ; i<D1->NbConstraints ; i++)
01089 {
01090 value_assign(M->p[i][0], D1->Constraint[i][0]);
01091 for(j=1 ; j<=n ; j++)
01092 value_multiply(M->p[i][j], D1->Constraint[i][j],
01093 C1->Ray[c1r][C1->Dimension+1]);
01094 value_multiply(M->p[i][j], D1->Constraint[i][D1->Dimension+1],
01095 C1->Ray[c1r][C1->Dimension+1]);
01096 Inner_Product (&D1->Constraint[i][j],
01097 &C1->Ray[c1r][n+1], m, &tmp);
01098 value_addto (M->p[i][j], tmp, M->p[i][j]);
01099 }
01100 value_clear(tmp);
01101 #ifdef DEBUGPP31
01102 fprintf(stderr,"Parameters instanciated = ");
01103 Matrix_Print(stderr,P_VALUE_FMT,M);
01104 #endif
01105
01106 npP = Constraints2Polyhedron(M,ws);
01107 Matrix_Free(M);
01108 nbr=0;
01109 for(i=0 ; i<npP->NbRays ; i++)
01110 if(value_zero_p(npP->Ray[i][n+1]))
01111 nbr++;
01112 Rays=Matrix_Alloc(nbr, n+2);
01113 for(i=0,j=0 ; j<nbr ; i++)
01114 if(value_zero_p(npP->Ray[i][n+1]))
01115 Vector_Copy(npP->Ray[i], Rays->p[j++], n+2);
01116 dimfaces = npP->NbBid;
01117 Polyhedron_Free(npP);
01118 #ifdef DEBUGPP31
01119 fprintf(stderr,"Rays = ");
01120 Matrix_Print(stderr,P_VALUE_FMT,Rays);
01121 fprintf(stderr,"dimfaces = %d\n", dimfaces);
01122 #endif
01123
01124
01125
01126 if(!CT) {
01127 if (C1->NbEq == 0) {
01128 Polyhedron_Free(C1);
01129 C1 = NULL;
01130 } else {
01131 Polyhedron *CEq1,
01132 *D2;
01133
01134
01135
01136
01137 M = Matrix_Alloc(C1->NbEq,m+2);
01138 for(j=0,i=0;i<C1->NbEq;++i,++j) {
01139 while(value_notzero_p(C1->Constraint[j][0]))
01140 ++j;
01141 value_assign(M->p[i][0],C1->Constraint[j][0]);
01142 Vector_Copy(&C1->Constraint[j][D->Dimension-m+1],&M->p[i][1],(m+1));
01143 }
01144 CEqualities = Constraints2Polyhedron(M,ws);
01145 Matrix_Free(M);
01146 CEq1 = align_context(CEqualities,D->Dimension,ws);
01147
01148
01149 D2 = DomainSimplify(D1,CEq1,ws);
01150 Polyhedron_Free(D1);
01151 Polyhedron_Free(C1);
01152 Polyhedron_Free(CEq1);
01153 D1 = D2;
01154 C1 = NULL;
01155 }
01156 }
01157 else {
01158 Polyhedron *CEq1,
01159 *D2;
01160
01161
01162
01163
01164 CEq1 = C1;
01165 M = Matrix_Alloc(C1->NbConstraints,m+2);
01166 for(i=0;i<C1->NbConstraints;++i) {
01167 value_assign(M->p[i][0],C1->Constraint[i][0]);
01168 Vector_Copy(&C1->Constraint[i][D->Dimension-m+1],&M->p[i][1],(m+1));
01169 }
01170 CEqualities = Constraints2Polyhedron( M, ws );
01171 Matrix_Free(M);
01172
01173 D2 = DomainSimplify(D1,CEq1,ws);
01174 Polyhedron_Free(D1);
01175 D1 = D2;
01176 C1 = Universe_Polyhedron(D2->Dimension);
01177
01178
01179
01180
01181 if( CEq1->NbEq )
01182 {
01183 m -= CEq1->NbEq;
01184 p = (int *)malloc(sizeof(int)*(CEq1->NbEq));
01185 }
01186 else
01187 p = NULL;
01188 ref = (int*) malloc(sizeof(int)*
01189 (CEq1->Dimension+2-CEq1->NbEq));
01190 *CT = PreElim_Columns(CEq1,p,ref,CEqualities->Dimension);
01191 D2 = Elim_Columns(D1,CEq1,p,ref);
01192 if (p)
01193 free(p);
01194 free(ref);
01195
01196 #ifdef DEBUGPP3
01197 fprintf(stderr,"D2\t Dim = %3d\tNbEq = %3d\tLines = %3d\n",
01198 D2->Dimension,D2->NbEq,D2->NbBid);
01199 C2 = Elim_Columns(C1,CEq1,p,ref);
01200 fprintf(stderr,"C2\t Dim = %3d\tNbEq = %3d\tLines = %3d\n",
01201 C2->Dimension,C2->NbEq,C2->NbBid);
01202 Polyhedron_Free(C2);
01203 #endif
01204
01205 Polyhedron_Free(D1);
01206 Polyhedron_Free(C1);
01207 D1 = D2;
01208 C1 = NULL;
01209 *CEq = CEqualities;
01210
01211 #ifdef DEBUGPP3
01212 fprintf(stderr,"Polyhedron CEq = ");
01213 Polyhedron_Print(stderr,P_VALUE_FMT,*CEq);
01214 fprintf(stderr,"Matrix CT = ");
01215 Matrix_Print(stderr,P_VALUE_FMT,*CT);
01216 #endif
01217
01218 Polyhedron_Free(CEq1);
01219 CEqualities = NULL;
01220
01221
01222 if(m==0) {
01223
01224 *Di = D1;
01225
01226 return GenParamPolyhedron(D1, Rays);
01227 }
01228 }
01229
01230 #ifdef DEBUGPP3
01231 fprintf(stderr,"Polyhedron D1 (D AND C) = ");
01232 Polyhedron_Print(stderr,P_VALUE_FMT, D1);
01233 fprintf(stderr,"Polyhedron CEqualities = ");
01234 if(CEqualities) Polyhedron_Print(stderr,P_VALUE_FMT, CEqualities);
01235 else fprintf(stderr,"NULL\n");
01236 #endif
01237
01238 KD = keep_dom;
01239 PDomains = NULL;
01240 PV_Result = NULL;
01241 nbPV = 0;
01242
01243 if (emptyQ(D1)) {
01244 Polyhedron_Free(D1);
01245 Matrix_Free(Rays);
01246 return NULL;
01247 }
01248
01249
01250
01251
01252 Sat = Poly2Sat(D1,&mf);
01253
01254 #ifdef DEBUGPP4
01255 fprintf(stderr,"Sat = ");
01256 SMPrint(Sat);
01257
01258 fprintf(stderr,"mf = ");
01259 for (i=0; i<nr; i++)
01260 fprintf(stderr,"%08X", mf[i]);
01261 fprintf(stderr, "\n");
01262 #endif
01263
01264
01265 egalite = (unsigned int *)malloc(sizeof(int)*D1->NbConstraints );
01266 memset(egalite,0, sizeof(int)*D1->NbConstraints);
01267
01268 for (i=0; i<D1->NbEq; i++)
01269 egalite[i] = 1;
01270
01271 Xi = Matrix_Alloc(n+1,m+1);
01272 Pi = Matrix_Alloc(m+1,m+1);
01273 PiTest = Matrix_Alloc(m+1,m+1);
01274 CTest = Matrix_Alloc(D->Dimension+1,D->NbConstraints);
01275 PiInv = Matrix_Alloc(m+1,m+2);
01276 RaysDi = Matrix_Alloc(D1->NbRays,m+2);
01277 m_dim = m;
01278
01279
01280
01281
01282 m_dim += dimfaces;
01283
01284
01285
01286
01287 #ifdef DEBUGPP3
01288 if (m_dim < D1->NbBid)
01289 fprintf(stderr, "m_dim (%d) < D1->NbBid (%d)\n", m_dim, D1->NbBid );
01290 #endif
01291 if (m_dim < D1->NbBid)
01292 m_dim = D1->NbBid;
01293
01294 #ifdef DEBUGPP
01295 nbfaces=0;
01296 #endif
01297 #ifdef DEBUGPP3
01298 fprintf(stderr, "m_dim = %d\n", m_dim );
01299 fprintf(stderr,
01300 "Target: find faces that saturate %d constraints and %d rays/lines\n",
01301 D1->Dimension - m_dim,m_dim+1);
01302 #endif
01303
01304
01305 scan_m_face(D1->NbEq,(D1->Dimension - m_dim - D1->NbEq),D1,mf);
01306
01307
01308
01309 #ifdef DEBUGPP
01310 fprintf( stderr, "Number of m-faces: %d\n", nbfaces );
01311 #endif
01312
01313 Matrix_Free(RaysDi);
01314 Matrix_Free(PiInv);
01315 Matrix_Free(PiTest);
01316 Matrix_Free(CTest);
01317 Matrix_Free(Pi);
01318 Matrix_Free(Xi);
01319 free(egalite);
01320 free(mf);
01321 SMFree(Sat);
01322
01323
01324
01325
01326
01327
01328 res = (Param_Polyhedron *) malloc (sizeof(Param_Polyhedron));
01329 res->nbV = nbPV;
01330 res->V = PV_Result;
01331 res->D = PDomains;
01332 res->Constraints = Polyhedron2Constraints(D1);
01333 res->Rays = Rays;
01334
01335 if(CT)
01336 *Di = D1;
01337 else
01338 Domain_Free(D1);
01339
01340 return(res);
01341 }
01342
01343
01344
01345
01346
01347 void Compute_PDomains(Param_Domain *PD,int nb_domains,int working_space) {
01348
01349 unsigned bx;
01350 int i, ix, nv;
01351 Polyhedron *dx, *d1, *d2;
01352 Param_Domain *p1, *p2, *p2prev, *PDNew;
01353
01354 if (nb_domains==0) {
01355
01356 #ifdef DEBUGPP5
01357 fprintf(stderr,"No domains\n");
01358 #endif
01359
01360 return;
01361 }
01362
01363
01364 if (!PD->next && PD->F)
01365 return;
01366
01367
01368 nv = (nb_domains - 1)/(8*sizeof(int)) + 1;
01369
01370 #ifdef DEBUGPP5
01371 fprintf(stderr,"nv = %d\n",nv);
01372 #endif
01373
01374 for(p1=PD,i=0,ix=0,bx=MSB;p1;p1=p1->next,i++) {
01375
01376
01377 p1->F = (unsigned *) malloc (nv * sizeof(unsigned));
01378
01379
01380 memset(p1->F,0,nv * sizeof(unsigned));
01381 p1->F[ix] |= bx;
01382 NEXT(ix, bx);
01383 }
01384
01385 #ifdef DEBUGPP5
01386 fprintf(stderr,"nb of vertices=%d\n",i);
01387 #endif
01388
01389
01390 ix = 0; bx=MSB;
01391 for (p1=PD;p1;p1=p1->next) {
01392 for (p2prev=p1,p2=p1->next;p2;p2prev=p2,p2=p2->next) {
01393
01394
01395 dx = PDomainIntersection(p1->Domain,p2->Domain,working_space);
01396
01397 if (!dx || emptyQ(dx)) {
01398
01399 #ifdef DEBUGPP5
01400 fprintf( stderr, "Empty dx (p1 inter p2). Continuing\n");
01401 #endif
01402 if(dx)
01403 Domain_Free(dx);
01404 continue;
01405 }
01406
01407 #ifdef DEBUGPP5
01408 fprintf(stderr,"Begin PDomainDifference\n");
01409 fprintf(stderr, "p1=");
01410 Polyhedron_Print(stderr,P_VALUE_FMT,p1->Domain);
01411 fprintf(stderr,"p2=");
01412 Polyhedron_Print(stderr,P_VALUE_FMT,p2->Domain);
01413 #endif
01414 d1 = PDomainDifference(p1->Domain,p2->Domain,working_space);
01415 d2 = PDomainDifference(p2->Domain,p1->Domain,working_space);
01416
01417 #ifdef DEBUGPP5
01418 fprintf(stderr,"p1\\p2=");
01419 Polyhedron_Print(stderr,P_VALUE_FMT,d1);
01420 fprintf(stderr,"p2\\p1=");
01421 Polyhedron_Print(stderr,P_VALUE_FMT,d2);
01422 fprintf(stderr,"END PDomainDifference\n\n");
01423 #endif
01424 if (!d1 || emptyQ(d1) || d1->NbEq!=0) {
01425
01426 #ifdef DEBUGPP5
01427 fprintf(stderr,"Empty d1\n");
01428 #endif
01429 if (d1)
01430 Domain_Free(d1);
01431 Domain_Free(dx);
01432
01433 if (!d2 || emptyQ(d2) || d2->NbEq!=0) {
01434
01435 #ifdef DEBUGPP5
01436 fprintf( stderr, "Empty d2 (deleting)\n");
01437 #endif
01438
01439 if (d2) Domain_Free(d2);
01440
01441
01442 for (i=0;i<nv;i++)
01443 p1->F[i] |= p2->F[i];
01444
01445
01446 p2prev->next = p2->next;
01447 Domain_Free(p2->Domain);
01448 free(p2->F);
01449 free(p2);
01450 p2 = p2prev;
01451 }
01452 else {
01453
01454 #ifdef DEBUGPP5
01455 fprintf( stderr, "p2 replaced by d2\n");
01456 #endif
01457
01458 for(i=0;i<nv;i++)
01459 p1->F[i] |= p2->F[i];
01460
01461
01462 Domain_Free( p2->Domain );
01463 p2->Domain = d2;
01464 }
01465 }
01466 else {
01467 if (!d2 || emptyQ(d2) || d2->NbEq!=0) {
01468
01469 #ifdef DEBUGPP5
01470 fprintf( stderr, "p1 replaced by d1\n");
01471 #endif
01472 if (d2) Domain_Free(d2);
01473
01474
01475 Domain_Free(dx);
01476
01477
01478 for(i=0;i<nv;i++)
01479 p2->F[i] |= p1->F[i];
01480
01481
01482 Domain_Free(p1->Domain);
01483 p1->Domain = d1;
01484 }
01485 else {
01486
01487 #ifdef DEBUGPP5
01488 fprintf(stderr,"Non-empty d1 and d2\nNew node created\n");
01489 #endif
01490
01491 PDNew = (Param_Domain *) malloc( sizeof(Param_Domain) );
01492 PDNew->F = (unsigned int *)malloc( nv*sizeof(int) );
01493 memset(PDNew->F,0,nv*sizeof(int));
01494 PDNew->Domain = dx;
01495
01496 for (i=0;i<nv;i++)
01497 PDNew->F[i] = p1->F[i] | p2->F[i];
01498
01499
01500 Domain_Free( p1->Domain );
01501 p1->Domain = d1;
01502
01503
01504 Domain_Free( p2->Domain );
01505 p2->Domain = d2;
01506
01507
01508 PDNew->next = p1->next;
01509 p1->next = PDNew;
01510 }
01511 }
01512 }
01513 if (p1->Domain->next) {
01514 Polyhedron *C = DomainConvex(p1->Domain, working_space);
01515 Domain_Free(p1->Domain);
01516 p1->Domain = C;
01517 }
01518 }
01519 }
01520
01521
01522
01523
01524
01525
01526
01527 Param_Polyhedron *Polyhedron2Param_Vertices(Polyhedron *Din,Polyhedron *Cin,int working_space) {
01528
01529 Param_Polyhedron *result;
01530
01531 POL_ENSURE_FACETS(Din);
01532 POL_ENSURE_VERTICES(Din);
01533 POL_ENSURE_FACETS(Cin);
01534 POL_ENSURE_VERTICES(Cin);
01535
01536 #ifdef DEBUGPP
01537 fprintf(stderr,"Polyhedron2Param_Vertices algorithm starting at : %.2fs\n",
01538 (float)clock()/CLOCKS_PER_SEC);
01539 #endif
01540
01541
01542 result = Find_m_faces(&Din,Cin,0,working_space,NULL,NULL);
01543
01544 #ifdef DEBUGPP
01545 fprintf(stderr, "nb of points : %d\n",result->nbV);
01546 #endif
01547
01548 #ifdef DEBUGPP
01549 fprintf(stderr, "end main loop : %.2fs\n", (float)clock()/CLOCKS_PER_SEC);
01550 #endif
01551
01552 return(result);
01553 }
01554
01555
01556
01557
01558 void Param_Vertices_Free(Param_Vertices *PV) {
01559
01560 Param_Vertices *next_pv;
01561
01562 while(PV) {
01563 next_pv = PV->next;
01564 if (PV->Vertex) Matrix_Free(PV->Vertex);
01565 if (PV->Domain) Matrix_Free(PV->Domain);
01566 if (PV->Facets) free(PV->Facets);
01567 free(PV);
01568 PV = next_pv;
01569 }
01570 }
01571
01572
01573
01574
01575 void Print_Vertex(FILE *DST, Matrix *V, const char **param_names)
01576 {
01577 int l, v;
01578 int first;
01579 Value gcd,tmp;
01580
01581 value_init(gcd); value_init(tmp);
01582
01583 fprintf(DST, "[" );
01584 for(l=0;l<V->NbRows;++l){
01585
01586
01587 first=1;
01588 fprintf(DST, " " );
01589 for(v=0;v < V->NbColumns-2;++v) {
01590 if(value_notzero_p(V->p[l][v])) {
01591 value_gcd(gcd, V->p[l][v], V->p[l][V->NbColumns-1]);
01592 value_divexact(tmp, V->p[l][v], gcd);
01593 if(value_posz_p(tmp)) {
01594 if(!first)
01595 fprintf(DST, "+");
01596 if(value_notone_p(tmp)) {
01597 value_print(DST,VALUE_FMT,tmp);
01598 }
01599 }
01600 else {
01601 if(value_mone_p(tmp))
01602 fprintf(DST, "-" );
01603 else {
01604 value_print(DST,VALUE_FMT,tmp);
01605 }
01606 }
01607 value_divexact(tmp, V->p[l][V->NbColumns-1], gcd);
01608 if(value_notone_p(tmp)) {
01609 fprintf(DST, "%s/", param_names[v]);
01610 value_print(DST,VALUE_FMT,tmp);
01611 }
01612 else
01613 fprintf(DST, "%s", param_names[v]);
01614 first=0;
01615 }
01616 }
01617
01618
01619 if(value_notzero_p(V->p[l][v]) || first) {
01620 if(value_posz_p(V->p[l][v]) && !first)
01621 fprintf(DST,"+");
01622 value_gcd(gcd, V->p[l][v], V->p[l][V->NbColumns-1]);
01623 value_divexact(tmp, V->p[l][v], gcd);
01624 value_print(DST,VALUE_FMT,tmp);
01625 value_divexact(tmp, V->p[l][V->NbColumns-1], gcd);
01626 if(value_notone_p(tmp)) {
01627 fprintf(DST,"/");
01628 value_print(DST,VALUE_FMT,tmp);
01629 fprintf(DST," ");
01630 }
01631 }
01632 if (l<V->NbRows-1)
01633 fprintf(DST, ", ");
01634 }
01635 fprintf(DST, " ]");
01636 value_clear(gcd); value_clear(tmp);
01637 return;
01638 }
01639
01640
01641
01642
01643
01644 Matrix *VertexCT(Matrix *V,Matrix *CT) {
01645
01646 Matrix *Vt;
01647 int i,j,k;
01648
01649 if(CT) {
01650
01651
01652 Vt = Matrix_Alloc(V->NbRows,CT->NbColumns+1);
01653 for(i=0;i<V->NbRows;++i) {
01654 value_assign(Vt->p[i][CT->NbColumns],V->p[i][V->NbColumns-1]);
01655 for(j=0;j<CT->NbColumns;j++) {
01656 for(k=0;k<CT->NbRows;k++)
01657 if(value_notzero_p(CT->p[k][j]))
01658 break;
01659 if(k<CT->NbRows)
01660 value_assign(Vt->p[i][j],V->p[i][k]);
01661 else
01662 value_set_si(Vt->p[i][j],0);
01663 }
01664 }
01665 return(Vt);
01666 }
01667 else
01668 return(NULL);
01669 }
01670
01671
01672
01673
01674 void Print_Domain(FILE *DST, Polyhedron *D, const char **pname)
01675 {
01676 int l, v;
01677 int first;
01678
01679 POL_ENSURE_FACETS(D);
01680 POL_ENSURE_VERTICES(D);
01681
01682 for(l=0;l<D->NbConstraints;++l) {
01683 fprintf(DST, " ");
01684 first = 1;
01685 for(v=1;v<=D->Dimension;++v) {
01686 if(value_notzero_p(D->Constraint[l][v])) {
01687 if(value_one_p(D->Constraint[l][v])) {
01688 if(first)
01689 fprintf(DST, "%s ", pname[v-1]);
01690 else
01691 fprintf(DST, "+ %s ", pname[v-1] );
01692 }
01693 else if(value_mone_p(D->Constraint[l][v]))
01694 fprintf(DST, "- %s ", pname[v-1] );
01695 else {
01696 if(value_pos_p(D->Constraint[l][v]) && !first )
01697 fprintf(DST, "+ " );
01698 value_print(DST,VALUE_FMT,D->Constraint[l][v]);
01699 fprintf(DST,"%s ",pname[v-1]);
01700 }
01701 first = 0;
01702 }
01703 }
01704 if(value_notzero_p(D->Constraint[l][v])) {
01705 if(value_pos_p(D->Constraint[l][v]) && !first)
01706 fprintf(DST,"+");
01707 fprintf(DST," ");
01708 value_print(DST,VALUE_FMT,D->Constraint[l][v]);
01709 }
01710 fprintf(DST,(value_notzero_p(D->Constraint[l][0])) ?" >= 0":" = 0");
01711 fprintf(DST, "\n" );
01712 }
01713 fprintf(DST, "\n");
01714
01715 if( D->next )
01716 {
01717 fprintf( DST, "UNION\n" );
01718 Print_Domain( DST, D->next, pname );
01719 }
01720 return;
01721 }
01722
01723
01724
01725
01726
01727 void Param_Vertices_Print(FILE *DST, Param_Vertices *PV, const char **param_names)
01728 {
01729 Polyhedron *poly;
01730
01731 while(PV) {
01732 fprintf(DST, "Vertex :\n" );
01733 Print_Vertex(DST,PV->Vertex,param_names);
01734
01735
01736 fprintf(DST, " If :\n" );
01737 poly = Constraints2Polyhedron(PV->Domain,200);
01738 Print_Domain(DST,poly,param_names);
01739 Domain_Free(poly);
01740 PV = PV->next;
01741 }
01742 return;
01743 }
01744
01745
01746
01747
01748
01749
01750
01751
01752 Param_Polyhedron *Polyhedron2Param_Domain(Polyhedron *Din,Polyhedron *Cin,int working_space) {
01753
01754 Param_Polyhedron *result;
01755 Param_Domain *D;
01756
01757 POL_ENSURE_FACETS(Din);
01758 POL_ENSURE_VERTICES(Din);
01759 POL_ENSURE_FACETS(Cin);
01760 POL_ENSURE_VERTICES(Cin);
01761
01762 if (emptyQ(Din) || emptyQ(Cin))
01763 return NULL;
01764
01765 #ifdef DEBUGPP
01766 fprintf(stderr,"Polyhedron2Param_Polyhedron algorithm starting at : %.2fs\n",
01767 (float)clock()/CLOCKS_PER_SEC);
01768 #endif
01769
01770
01771
01772 result = Find_m_faces(&Din,Cin,1,working_space,NULL,NULL);
01773
01774 #ifdef DEBUGPP
01775 if(result) fprintf(stderr, "Number of vertices : %d\n",result->nbV);
01776 fprintf(stderr,"Vertices found at : %.2fs\n",(float)clock()/CLOCKS_PER_SEC);
01777 #endif
01778
01779
01780 if(result && Cin->Dimension>0)
01781 Compute_PDomains(result->D,result->nbV,working_space);
01782 if(result && CEqualities)
01783 for(D=result->D;D;D=D->next)
01784 D->Domain = Add_CEqualities(D->Domain);
01785 Polyhedron_Free(CEqualities);
01786
01787 #ifdef DEBUGPP
01788 fprintf(stderr, "domains found at : %.2fs\n", (float)clock()/CLOCKS_PER_SEC);
01789 #endif
01790
01791 return(result);
01792 }
01793
01794
01795
01796
01797 Param_Polyhedron *Polyhedron2Param_SimplifiedDomain(Polyhedron **Din,Polyhedron *Cin,int working_space,Polyhedron **CEq,Matrix **CT) {
01798
01799 Param_Polyhedron *result;
01800
01801 assert(CEq != NULL);
01802 assert(CT != NULL);
01803
01804 POL_ENSURE_FACETS(*Din);
01805 POL_ENSURE_VERTICES(*Din);
01806 POL_ENSURE_FACETS(Cin);
01807 POL_ENSURE_VERTICES(Cin);
01808
01809 #ifdef DEBUGPP
01810 fprintf(stderr,"Polyhedron2Param_Polyhedron algorithm starting at : %.2fs\n",
01811 (float)clock()/CLOCKS_PER_SEC);
01812 #endif
01813
01814
01815
01816 result = Find_m_faces(Din,Cin,1,working_space,CEq,CT);
01817
01818 #ifdef DEBUGPP
01819 if(result) fprintf(stderr, "Number of vertices : %d\n",result->nbV);
01820 fprintf(stderr,"Vertices found at : %.2fs\n",(float)clock()/CLOCKS_PER_SEC);
01821 #endif
01822
01823
01824 if(result && Cin->Dimension>0)
01825 Compute_PDomains(result->D,result->nbV,working_space);
01826
01827
01828
01829
01830
01831
01832
01833 #ifdef DEBUGPP
01834 fprintf(stderr, "domains found at : %.2fs\n", (float)clock()/CLOCKS_PER_SEC);
01835 #endif
01836
01837 return(result);
01838 }
01839
01840
01841
01842
01843
01844 void Param_Domain_Free(Param_Domain *PD) {
01845
01846 Param_Domain *next_pd;
01847
01848 while(PD) {
01849 free(PD->F);
01850 Domain_Free(PD->Domain);
01851 next_pd = PD->next;
01852 free(PD);
01853 PD = next_pd;
01854 }
01855 return;
01856 }
01857
01858
01859
01860
01861 void Param_Polyhedron_Free(Param_Polyhedron *P) {
01862
01863 if (!P) return;
01864 Param_Vertices_Free(P->V);
01865 Param_Domain_Free(P->D);
01866 if (P->Constraints)
01867 Matrix_Free(P->Constraints);
01868 free(P);
01869 return;
01870 }
01871
01872
01873
01874
01875 void Param_Polyhedron_Scale_Integer(Param_Polyhedron *PP, Polyhedron **P,
01876 Value *det, unsigned MaxRays)
01877 {
01878 int i;
01879 int nb_param, nb_vars;
01880 Vector *denoms;
01881 Param_Vertices *V;
01882 Value global_var_lcm;
01883 Matrix *expansion;
01884
01885 value_set_si(*det, 1);
01886 if (!PP->nbV)
01887 return;
01888
01889 nb_param = PP->D->Domain->Dimension;
01890 nb_vars = PP->V->Vertex->NbRows;
01891
01892
01893
01894
01895 denoms = Vector_Alloc(nb_vars);
01896 value_init(global_var_lcm);
01897
01898
01899 for (V = PP->V; V; V = V->next)
01900 for (i = 0; i < nb_vars; i++)
01901 value_lcm(denoms->p[i], denoms->p[i], V->Vertex->p[i][nb_param+1]);
01902
01903 value_set_si(global_var_lcm, 1);
01904 for (i = 0; i < nb_vars; i++) {
01905 value_multiply(*det, *det, denoms->p[i]);
01906 value_lcm(global_var_lcm, global_var_lcm, denoms->p[i]);
01907 }
01908
01909
01910 for (V = PP->V; V; V = V->next)
01911 for (i = 0; i < nb_vars; i++) {
01912 Vector_Scale(V->Vertex->p[i], V->Vertex->p[i], denoms->p[i], nb_param+1);
01913 Vector_Normalize(V->Vertex->p[i], nb_param+2);
01914 }
01915
01916
01917
01918 for (i = 0; i < nb_vars; i++)
01919 value_division(denoms->p[i], global_var_lcm, denoms->p[i]);
01920
01921
01922
01923 expansion = Matrix_Alloc(nb_vars+nb_param+1, nb_vars+nb_param+1);
01924 for (i = 0; i < nb_vars; i++)
01925 value_assign(expansion->p[i][i], denoms->p[i]);
01926 for (i = nb_vars; i < nb_vars+nb_param+1; i++)
01927 value_assign(expansion->p[i][i], global_var_lcm);
01928
01929
01930 if (P)
01931 *P = Polyhedron_Preimage(*P, expansion, MaxRays);
01932
01933 Matrix_Free(expansion);
01934 value_clear(global_var_lcm);
01935 Vector_Free(denoms);
01936 }