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
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #undef POLY_DEBUG
00038 #undef POLY_RR_DEBUG
00039 #undef POLY_CH_DEBUG
00040
00041 #include <stdio.h>
00042 #include <stdlib.h>
00043 #include <string.h>
00044 #include <assert.h>
00045 #include <polylib/polylib.h>
00046
00047 #ifdef MAC_OS
00048 #define abs __abs
00049 #endif
00050
00051
00052 #define WSIZE (8*sizeof(int))
00053
00054 #define bexchange(a, b, l)\
00055 {\
00056 char *t = (char *)malloc(l*sizeof(char));\
00057 memcpy((t), (char *)(a), (int)(l));\
00058 memcpy((char *)(a), (char *)(b), (int)(l));\
00059 memcpy((char *)(b), (t), (int)(l));\
00060 free(t); \
00061 }
00062
00063 #define exchange(a, b, t)\
00064 { (t)=(a); (a)=(b); (b)=(t); }
00065
00066
00067
00068
00069
00070 void errormsg1(const char *f, const char *msgname, const char *msg);
00071
00072 int Pol_status;
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 typedef struct {
00083 unsigned int NbRows;
00084 unsigned int NbColumns;
00085 int **p;
00086 int *p_init;
00087 } SatMatrix;
00088
00089
00090
00091
00092 static SatMatrix *SMAlloc(int rows,int cols) {
00093
00094 int **q, *p, i;
00095 SatMatrix *result;
00096
00097 result = (SatMatrix *) malloc (sizeof(SatMatrix));
00098 if(!result) {
00099 errormsg1("SMAlloc", "outofmem", "out of memory space");
00100 return 0;
00101 }
00102 result->NbRows = rows;
00103 result->NbColumns = cols;
00104 if(rows == 0 || cols == 0) {
00105 result->p = NULL;
00106 return result;
00107 }
00108 result->p = q = (int **)malloc(rows * sizeof(int *));
00109 if(!result->p) {
00110 errormsg1("SMAlloc", "outofmem", "out of memory space");
00111 return 0;
00112 }
00113 result->p_init = p = (int *)malloc (rows * cols * sizeof (int));
00114 if(!result->p_init) {
00115 errormsg1("SMAlloc", "outofmem", "out of memory space");
00116 return 0;
00117 }
00118 for (i=0; i<rows; i++) {
00119 *q++ = p;
00120 p += cols;
00121 }
00122 return result;
00123 }
00124
00125
00126
00127
00128 static void SMFree (SatMatrix **matrix) {
00129 SatMatrix *SM = *matrix;
00130
00131 if (SM) {
00132 if (SM->p) {
00133 free ((char *) SM->p_init);
00134 free ((char *) SM->p);
00135 }
00136 free ((char *) SM);
00137 *matrix = NULL;
00138 }
00139 }
00140
00141
00142
00143
00144
00145 static void SMPrint (SatMatrix *matrix) {
00146
00147 int *p;
00148 int i, j;
00149 unsigned NbRows, NbColumns;
00150
00151 fprintf(stderr,"%d %d\n",NbRows=matrix->NbRows, NbColumns=matrix->NbColumns);
00152 for (i=0;i<NbRows;i++) {
00153 p = *(matrix->p+i);
00154 for (j=0;j<NbColumns;j++)
00155 fprintf(stderr, " %10X ", *p++);
00156 fprintf(stderr, "\n");
00157 }
00158 }
00159
00160
00161
00162
00163 static void SatVector_OR(int *p1,int *p2,int *p3,unsigned length) {
00164
00165 int *cp1, *cp2, *cp3;
00166 int i;
00167
00168 cp1=p1;
00169 cp2=p2;
00170 cp3=p3;
00171 for (i=0;i<length;i++) {
00172 *cp3 = *cp1 | *cp2;
00173 cp3++;
00174 cp1++;
00175 cp2++;
00176 }
00177 }
00178
00179
00180
00181
00182 #define SMVector_Copy(p1, p2, length) \
00183 memcpy((char *)(p2), (char *)(p1), (int)((length)*sizeof(int)))
00184
00185
00186
00187
00188 #define SMVector_Init(p1, length) \
00189 memset((char *)(p1), 0, (int)((length)*sizeof(int)))
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 static void Combine(Value *p1, Value *p2, Value *p3, int pos, unsigned length) {
00202
00203 Value a1, a2, gcd;
00204 Value abs_a1,abs_a2,neg_a1;
00205
00206
00207 value_init(a1); value_init(a2); value_init(gcd);
00208 value_init(abs_a1); value_init(abs_a2); value_init(neg_a1);
00209
00210
00211 value_assign(a1,p1[pos]);
00212
00213
00214 value_assign(a2,p2[pos]);
00215
00216
00217 value_absolute(abs_a1,a1);
00218
00219
00220 value_absolute(abs_a2,a2);
00221
00222
00223 value_gcd(gcd, abs_a1, abs_a2);
00224
00225
00226 value_divexact(a1, a1, gcd);
00227
00228
00229 value_divexact(a2, a2, gcd);
00230
00231
00232 value_oppose(neg_a1,a1);
00233
00234 Vector_Combine(p1+1,p2+1,p3+1,a2,neg_a1,length);
00235 Vector_Normalize(p3+1,length);
00236
00237
00238 value_clear(a1); value_clear(a2); value_clear(gcd);
00239 value_clear(abs_a1); value_clear(abs_a2); value_clear(neg_a1);
00240
00241 return;
00242 }
00243
00244
00245
00246
00247
00248
00249 static SatMatrix *TransformSat(Matrix *Mat, Matrix *Ray, SatMatrix *Sat) {
00250
00251 int i, j, sat_nbcolumns;
00252 unsigned jx1, jx2, bx1, bx2;
00253 SatMatrix *result;
00254
00255 if (Mat->NbRows != 0)
00256 sat_nbcolumns = (Mat->NbRows-1) /(sizeof(int)*8) + 1;
00257 else
00258 sat_nbcolumns = 0;
00259
00260 result = SMAlloc(Ray->NbRows, sat_nbcolumns);
00261 SMVector_Init(result->p_init, Ray->NbRows * sat_nbcolumns);
00262
00263 for(i=0,jx1=0,bx1=MSB; i<Ray->NbRows; i++) {
00264 for(j=0,jx2=0,bx2=MSB; j<Mat->NbRows; j++) {
00265 if (Sat->p[j][jx1] & bx1)
00266 result->p[i][jx2] |= bx2;
00267 NEXT(jx2,bx2);
00268 }
00269 NEXT(jx1, bx1);
00270 }
00271 return result;
00272 }
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288 static void RaySort(Matrix *Ray,SatMatrix *Sat,int NbBid,int NbRay,int *equal_bound,int *sup_bound,unsigned RowSize1, unsigned RowSize2, unsigned bx, unsigned jx) {
00289
00290 int inf_bound;
00291 Value **uni_eq, **uni_sup, **uni_inf;
00292 int **inc_eq, **inc_sup, **inc_inf;
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 *sup_bound = *equal_bound = NbBid;
00304 uni_sup = uni_eq = Ray->p+NbBid;
00305 inc_sup = inc_eq = Sat->p+NbBid;
00306 inf_bound = NbRay;
00307 uni_inf = Ray->p+NbRay;
00308 inc_inf = Sat->p+NbRay;
00309
00310 while (inf_bound>*sup_bound) {
00311 if (value_zero_p(**uni_sup)) {
00312 if (inc_eq != inc_sup) {
00313 Vector_Exchange(*uni_eq,*uni_sup,RowSize1);
00314 bexchange(*inc_eq,*inc_sup,RowSize2);
00315 }
00316 (*equal_bound)++; uni_eq++; inc_eq++;
00317 (*sup_bound)++; uni_sup++; inc_sup++;
00318 }
00319 else {
00320 *((*inc_sup)+jx)|=bx;
00321
00322
00323 if (value_neg_p(**uni_sup)) {
00324 inf_bound--; uni_inf--; inc_inf--;
00325 if (inc_inf != inc_sup) {
00326 Vector_Exchange(*uni_inf,*uni_sup,RowSize1);
00327 bexchange(*inc_inf,*inc_sup,RowSize2);
00328 }
00329 }
00330 else {
00331 (*sup_bound)++; uni_sup++; inc_sup++;
00332 }
00333 }
00334 }
00335 }
00336
00337 static void SatMatrix_Extend(SatMatrix *Sat, Matrix* Mat, unsigned rows)
00338 {
00339 int i;
00340 unsigned cols;
00341 cols = (Mat->NbRows - 1)/(sizeof(int)*8) + 1;
00342
00343 Sat->p = (int **)realloc(Sat->p, rows * sizeof(int *));
00344 if(!Sat->p) {
00345 errormsg1("SatMatrix_Extend", "outofmem", "out of memory space");
00346 return;
00347 }
00348 Sat->p_init = (int *)realloc(Sat->p_init, rows * cols * sizeof (int));
00349 if(!Sat->p_init) {
00350 errormsg1("SatMatrix_Extend", "outofmem", "out of memory space");
00351 return;
00352 }
00353 for (i = 0; i < rows; ++i)
00354 Sat->p[i] = Sat->p_init + (i * cols);
00355 Sat->NbRows = rows;
00356 }
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369 static int Chernikova (Matrix *Mat,Matrix *Ray,SatMatrix *Sat, unsigned NbBid, unsigned NbMaxRays, unsigned FirstConstraint,unsigned dual) {
00370
00371 unsigned NbRay, Dimension, NbConstraints, RowSize1, RowSize2, sat_nbcolumns;
00372 int sup_bound, equal_bound, index_non_zero, bound;
00373 int i, j, k, l, redundant, rayonly, nbcommonconstraints;
00374 int *Temp, aux;
00375 int *ip1, *ip2;
00376 unsigned bx, m, jx;
00377 Value *p1, *p2, *p3;
00378
00379 #ifdef POLY_CH_DEBUG
00380 fprintf(stderr, "[Chernikova: Input]\nRay = ");
00381 Matrix_Print(stderr,0,Ray);
00382 fprintf(stderr, "\nConstraints = ");
00383 Matrix_Print(stderr,0,Mat);
00384 fprintf(stderr, "\nSat = ");
00385 SMPrint(Sat);
00386 #endif
00387
00388 NbConstraints=Mat->NbRows;
00389 NbRay = Ray->NbRows;
00390 Dimension = Mat->NbColumns-1;
00391 sat_nbcolumns=Sat->NbColumns;
00392
00393 RowSize1=(Dimension+1);
00394 RowSize2=sat_nbcolumns * sizeof(int);
00395
00396 Temp=(int *)malloc(RowSize2);
00397 if(!Temp) {
00398 errormsg1("Chernikova", "outofmem", "out of memory space");
00399 return 0;
00400 }
00401 CATCH(any_exception_error) {
00402
00403
00404
00405
00406
00407 free(Temp);
00408 RETHROW();
00409 }
00410 TRY {
00411 jx = FirstConstraint/WSIZE;
00412 bx = MSB; bx >>= FirstConstraint%WSIZE;
00413 for (k=FirstConstraint; k<NbConstraints; k++) {
00414
00415
00416
00417
00418
00419
00420 index_non_zero = NbRay;
00421 for (i=0; i<NbRay; i++) {
00422 p1 = Ray->p[i]+1;
00423 p2 = Mat->p[k]+1;
00424 p3 = Ray->p[i];
00425
00426
00427 value_multiply(*p3,*p1,*p2);
00428 p1++; p2++;
00429 for (j=1; j<Dimension; j++) {
00430
00431
00432 value_addmul(*p3, *p1, *p2);
00433 p1++; p2++;
00434 }
00435 if (value_notzero_p(*p3) && (i<index_non_zero))
00436 index_non_zero=i;
00437 }
00438
00439 #ifdef POLY_CH_DEBUG
00440 fprintf(stderr, "[Chernikova: A]\nRay = ");
00441 Matrix_Print(stderr,0,Ray);
00442 fprintf(stderr, "\nConstraints = ");
00443 Matrix_Print(stderr,0,Mat);
00444 fprintf(stderr, "\nSat = ");
00445 SMPrint (Sat);
00446 #endif
00447
00448
00449 if (index_non_zero<NbBid) {
00450
00451
00452 NbBid--;
00453 if (NbBid!=index_non_zero)
00454 Vector_Exchange(Ray->p[index_non_zero],Ray->p[NbBid],RowSize1);
00455
00456 #ifdef POLY_CH_DEBUG
00457 fprintf(stderr,"************\n");
00458 for(i=0;i<RowSize1;i++) {
00459 value_print(stderr,P_VALUE_FMT,Ray->p[index_non_zero][i]);
00460 }
00461 fprintf(stderr,"\n******\n");
00462 for(i=0;i<RowSize1;i++) {
00463 value_print(stderr,P_VALUE_FMT,Ray->p[NbBid][i]);
00464 }
00465 fprintf(stderr,"\n*******\n");
00466 #endif
00467
00468
00469 for (i=0; i<NbBid; i++)
00470 if (value_notzero_p(Ray->p[i][0]))
00471 Combine(Ray->p[i],Ray->p[NbBid],Ray->p[i],0,Dimension);
00472
00473
00474
00475
00476 if (value_neg_p(Ray->p[NbBid][0])) {
00477 p1=Ray->p[NbBid];
00478 for (j=0;j<Dimension+1; j++) {
00479
00480
00481 value_oppose(*p1,*p1);
00482 p1++;
00483 }
00484 }
00485
00486 #ifdef POLY_CH_DEBUG
00487 fprintf(stderr, "[Chernikova: B]\nRay = ");
00488 Ray->NbRows=NbRay;
00489 Matrix_Print(stderr,0,Ray);
00490 fprintf(stderr, "\nConstraints = ");
00491 Matrix_Print(stderr,0,Mat);
00492 fprintf(stderr, "\nSat = ");
00493 SMPrint(Sat);
00494 #endif
00495
00496
00497 for (i=NbBid+1; i<NbRay; i++)
00498 if (value_notzero_p(Ray->p[i][0]))
00499 Combine(Ray->p[i],Ray->p[NbBid],Ray->p[i],0,Dimension);
00500
00501
00502 if (value_notzero_p(Mat->p[k][0])) {
00503 for (j=0;j<sat_nbcolumns;j++) {
00504 Sat->p[NbBid][j] = 0;
00505 }
00506
00507 Sat->p[NbBid][jx] |= bx;
00508 }
00509 else {
00510 if (--NbRay != NbBid) {
00511 Vector_Copy(Ray->p[NbRay],Ray->p[NbBid],Dimension+1);
00512 SMVector_Copy(Sat->p[NbRay],Sat->p[NbBid],sat_nbcolumns);
00513 }
00514 }
00515
00516 #ifdef POLY_CH_DEBUG
00517 fprintf(stderr, "[Chernikova: C]\nRay = ");
00518 Ray->NbRows=NbRay;
00519 Matrix_Print(stderr,0,Ray);
00520 fprintf(stderr, "\nConstraints = ");
00521 Matrix_Print(stderr,0,Mat);
00522 fprintf(stderr, "\nSat = ");
00523 SMPrint (Sat);
00524 #endif
00525
00526 }
00527 else {
00528 RaySort(Ray, Sat, NbBid, NbRay, &equal_bound, &sup_bound,
00529 RowSize1, RowSize2,bx,jx);
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542 #ifdef POLY_CH_DEBUG
00543 fprintf(stderr, "[Chernikova: D]\nRay = ");
00544 Ray->NbRows=NbRay;
00545 Matrix_Print(stderr,0,Ray);
00546 fprintf(stderr, "\nConstraints = ");
00547 Matrix_Print(stderr,0,Mat);
00548 fprintf(stderr, "\nSat = ");
00549 SMPrint (Sat);
00550 #endif
00551
00552
00553 bound=NbRay;
00554 for (i=equal_bound; i<sup_bound; i++)
00555 for(j=sup_bound; j<bound; j++) {
00556
00557
00558
00559
00560
00561
00562 nbcommonconstraints = 0;
00563 for (l=0; l<jx; l++) {
00564 aux = Temp[l] = Sat->p[i][l] | Sat->p[j][l];
00565 for (m=MSB; m!=0; m>>=1)
00566 if (!(aux&m))
00567 nbcommonconstraints++;
00568 }
00569 aux = Temp[jx] = Sat->p[i][jx] | Sat->p[j][jx];
00570 for (m=MSB; m!=bx; m>>=1)
00571 if (!(aux&m))
00572 nbcommonconstraints++;
00573 rayonly = (value_zero_p(Ray->p[i][Dimension]) &&
00574 value_zero_p(Ray->p[j][Dimension]) &&
00575 (dual == 0));
00576 if(rayonly)
00577 nbcommonconstraints++;
00578
00579
00580
00581
00582
00583 if (nbcommonconstraints+NbBid>=Dimension-2) {
00584
00585 redundant=0;
00586 for (m=NbBid; m<bound; m++)
00587 if ((m!=i)&&(m!=j)) {
00588
00589
00590
00591
00592
00593 if (rayonly && value_notzero_p(Ray->p[m][Dimension]))
00594 continue;
00595
00596
00597
00598
00599 ip1 = Temp;
00600 ip2 = Sat->p[m];
00601 for (l=0; l<=jx; l++,ip2++,ip1++)
00602 if (*ip2 & ~*ip1)
00603 break;
00604 if (l>jx) {
00605 redundant=1;
00606 break;
00607 }
00608 }
00609
00610 #ifdef POLY_CH_DEBUG
00611 fprintf(stderr, "[Chernikova: E]\nRay = ");
00612 Ray->NbRows=NbRay;
00613 Matrix_Print(stderr,0,Ray);
00614 fprintf(stderr, "\nConstraints = ");
00615 Matrix_Print(stderr,0,Mat);
00616 fprintf(stderr, "\nSat = ");
00617 SMPrint (Sat);
00618 #endif
00619
00620
00621
00622
00623
00624 if (!redundant) {
00625 if (NbRay==NbMaxRays) {
00626 NbMaxRays *= 2;
00627 Ray->NbRows = NbRay;
00628 Matrix_Extend(Ray, NbMaxRays);
00629 SatMatrix_Extend(Sat, Mat, NbMaxRays);
00630 }
00631
00632
00633 Combine(Ray->p[j],Ray->p[i],Ray->p[NbRay],0,Dimension);
00634 SatVector_OR(Sat->p[j],Sat->p[i],Sat->p[NbRay],sat_nbcolumns);
00635 Sat->p[NbRay][jx] &= ~bx;
00636 NbRay++;
00637 }
00638 }
00639 }
00640
00641 #ifdef POLY_CH_DEBUG
00642 fprintf(stderr,
00643 "[Chernikova: F]\n"
00644 "sup_bound=%d\n"
00645 "equal_bound=%d\n"
00646 "bound=%d\n"
00647 "NbRay=%d\n"
00648 "Dimension = %d\n"
00649 "Ray = ",sup_bound,equal_bound,bound,NbRay,Dimension);
00650 #endif
00651 #ifdef POLY_CH_DEBUG
00652 Ray->NbRows=NbRay;
00653 fprintf(stderr, "[Chernikova: F]:\nRay = ");
00654 Matrix_Print(stderr,0,Ray);
00655 #endif
00656
00657
00658
00659
00660 j = (value_notzero_p(Mat->p[k][0])) ?
00661 sup_bound : equal_bound;
00662
00663 i = NbRay;
00664 #ifdef POLY_CH_DEBUG
00665 fprintf(stderr, "i = %d\nj = %d \n", i, j);
00666 #endif
00667 while ((j<bound)&&(i>bound)) {
00668 i--;
00669 Vector_Copy(Ray->p[i],Ray->p[j],Dimension+1);
00670 SMVector_Copy(Sat->p[i],Sat->p[j],sat_nbcolumns);
00671 j++;
00672 }
00673
00674 #ifdef POLY_CH_DEBUG
00675 fprintf(stderr, "i = %d\nj = %d \n", i, j);
00676 fprintf(stderr,
00677 "[Chernikova: F]\n"
00678 "sup_bound=%d\n"
00679 "equal_bound=%d\n"
00680 "bound=%d\n"
00681 "NbRay=%d\n"
00682 "Dimension = %d\n"
00683 "Ray = ",sup_bound,equal_bound,bound,NbRay, Dimension);
00684 #endif
00685 #ifdef POLY_CH_DEBUG
00686 Ray->NbRows=NbRay;
00687 fprintf(stderr, "[Chernikova: G]\nRay = ");
00688 Matrix_Print(stderr,0,Ray);
00689 #endif
00690 if (j==bound)
00691 NbRay=i;
00692 else
00693 NbRay=j;
00694 }
00695 NEXT(jx,bx);
00696 }
00697 Ray->NbRows=NbRay;
00698 Sat->NbRows=NbRay;
00699
00700 }
00701
00702 UNCATCH(any_exception_error);
00703 free(Temp);
00704
00705 #ifdef POLY_CH_DEBUG
00706 fprintf(stderr, "[Chernikova: Output]\nRay = ");
00707 Matrix_Print(stderr,0,Ray);
00708 fprintf(stderr, "\nConstraints = ");
00709 Matrix_Print(stderr,0,Mat);
00710 fprintf(stderr, "\nSat = ");
00711 SMPrint (Sat);
00712 #endif
00713
00714 return 0;
00715 }
00716
00717 static int Gauss4(Value **p, int NbEq, int NbRows, int Dimension)
00718 {
00719 int i, j, k, pivot, Rank;
00720 int *column_index = NULL;
00721 Value gcd;
00722
00723 value_init(gcd);
00724 column_index=(int *)malloc(Dimension * sizeof(int));
00725 if(!column_index) {
00726 errormsg1("Gauss","outofmem","out of memory space");
00727 value_clear(gcd);
00728 return 0;
00729 }
00730 Rank=0;
00731
00732 CATCH(any_exception_error) {
00733 if (column_index)
00734 free(column_index);
00735 value_clear(gcd);
00736 RETHROW();
00737 }
00738 TRY {
00739
00740 for (j=1; j<=Dimension; j++) {
00741 for (i=Rank; i<NbEq; i++)
00742
00743
00744 if (value_notzero_p(p[i][j]))
00745 break;
00746 if (i!=NbEq) {
00747 if (i!=Rank)
00748 Vector_Exchange(p[Rank]+1,p[i]+1,Dimension);
00749
00750
00751
00752 Vector_Gcd(p[Rank]+1,Dimension,&gcd);
00753
00754 if (value_cmp_si(gcd, 2) >= 0)
00755 Vector_AntiScale(p[Rank]+1, p[Rank]+1, gcd, Dimension);
00756
00757 if (value_neg_p(p[Rank][j]))
00758 Vector_Oppose(p[Rank]+1, p[Rank]+1, Dimension);
00759
00760 pivot=i;
00761 for (i=pivot+1; i<NbEq; i++) {
00762
00763
00764 if (value_notzero_p(p[i][j]))
00765 Combine(p[i],p[Rank],p[i],j,Dimension);
00766 }
00767
00768
00769
00770
00771
00772 column_index[Rank]=j;
00773 Rank++;
00774 }
00775 }
00776
00777
00778 for (k=Rank-1; k>=0; k--) {
00779 j = column_index[k];
00780
00781
00782 for (i=0; i<k; i++) {
00783
00784
00785 if (value_notzero_p(p[i][j]))
00786 Combine(p[i],p[k],p[i],j,Dimension);
00787 }
00788
00789
00790 for (i=NbEq;i<NbRows;i++) {
00791
00792
00793 if (value_notzero_p(p[i][j]))
00794 Combine(p[i],p[k],p[i],j,Dimension);
00795 }
00796 }
00797 }
00798
00799 UNCATCH(any_exception_error);
00800 free(column_index), column_index = NULL;
00801
00802 value_clear(gcd);
00803 return Rank;
00804 }
00805
00806
00807
00808
00809
00810
00811
00812 int Gauss(Matrix *Mat, int NbEq, int Dimension)
00813 {
00814 int Rank;
00815
00816 #ifdef POLY_DEBUG
00817 fprintf(stderr, "[Gauss : Input]\nRay =");
00818 Matrix_Print(stderr,0,Mat);
00819 #endif
00820
00821 Rank = Gauss4(Mat->p, NbEq, Mat->NbRows, Dimension);
00822
00823 #ifdef POLY_DEBUG
00824 fprintf(stderr, "[Gauss : Output]\nRay =");
00825 Matrix_Print(stderr,0,Mat);
00826 #endif
00827
00828 return Rank;
00829 }
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842 static Polyhedron *Remove_Redundants(Matrix *Mat,Matrix *Ray,SatMatrix *Sat,unsigned *Filter) {
00843
00844 int i, j, k;
00845 unsigned Dimension, sat_nbcolumns, NbRay, NbConstraints, RowSize2,
00846 *Trace = NULL, *bx = NULL, *jx = NULL, Dim_RaySpace, b;
00847 unsigned NbBid, NbUni, NbEq, NbIneq;
00848 unsigned NbBid2, NbUni2, NbEq2, NbIneq2;
00849 int Redundant;
00850 int aux, *temp2 = NULL;
00851 Polyhedron *Pol = NULL;
00852 Vector *temp1 = NULL;
00853 unsigned Status;
00854
00855 Dimension = Mat->NbColumns-1;
00856 NbRay = Ray->NbRows;
00857 sat_nbcolumns = Sat->NbColumns;
00858 NbConstraints = Mat->NbRows;
00859 RowSize2=sat_nbcolumns * sizeof(int);
00860
00861 temp1 = Vector_Alloc(Dimension+1);
00862 if (!temp1) {
00863 errormsg1("Remove_Redundants", "outofmem", "out of memory space");
00864 return 0;
00865 }
00866
00867 if (Filter) {
00868 temp2 = (int *)calloc(sat_nbcolumns, sizeof(int));
00869 if (!temp2)
00870 goto oom;
00871 }
00872
00873
00874
00875 bx = (unsigned *)malloc(NbConstraints * sizeof(unsigned));
00876 if (!bx)
00877 goto oom;
00878 jx = (unsigned *)malloc(NbConstraints * sizeof(unsigned));
00879 if (!jx)
00880 goto oom;
00881 CATCH(any_exception_error) {
00882
00883 Vector_Free(temp1);
00884 if (temp2) free(temp2);
00885 if (bx) free(bx);
00886 if (jx) free(jx);
00887 if (Trace) free(Trace);
00888 if (Pol) Polyhedron_Free(Pol);
00889
00890 RETHROW();
00891 }
00892 TRY {
00893
00894
00895
00896
00897
00898
00899 i = 0;
00900 b = MSB;
00901 for (j=0; j<NbConstraints; j++) {
00902 jx[j] = i;
00903 bx[j] = b;
00904 NEXT(i,b);
00905 }
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915 aux = 0;
00916 for (i=0; i<NbRay; i++) {
00917
00918
00919 value_set_si(Ray->p[i][0],0);
00920
00921
00922 if (value_notzero_p(Ray->p[i][Dimension]))
00923 aux++;
00924 }
00925
00926
00927 if (!aux)
00928 goto empty;
00929
00930 #ifdef POLY_RR_DEBUG
00931 fprintf(stderr, "[Remove_redundants : Init]\nConstraints =");
00932 Matrix_Print(stderr,0,Mat);
00933 fprintf(stderr, "\nRays =");
00934 Matrix_Print(stderr,0,Ray);
00935 #endif
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946 NbEq=0;
00947
00948 #ifdef POLY_RR_DEBUG
00949 fprintf (stderr, " j = ");
00950 #endif
00951
00952 for (j=0; j<NbConstraints; j++) {
00953
00954 #ifdef POLY_RR_DEBUG
00955 fprintf (stderr, " %i ", j);
00956 fflush (stderr);
00957 #endif
00958
00959
00960 if (Filter && value_zero_p(Mat->p[j][0]))
00961 temp2[jx[j]] |= bx[j];
00962
00963 value_set_si(Mat->p[j][0],0);
00964
00965
00966 i = First_Non_Zero(Mat->p[j]+1, Dimension-1);
00967
00968 #ifdef POLY_RR_DEBUG
00969 fprintf(stderr, "[Remove_redundants : IntoStep1]\nConstraints =");
00970 Matrix_Print(stderr,0,Mat);
00971 fprintf (stderr, " j = %i \n", j);
00972 #endif
00973
00974
00975
00976
00977
00978 if (i == -1) {
00979 for (i=0; i<NbRay; i++)
00980 if (!(Sat->p[i][jx[j]]&bx[j])) {
00981
00982
00983 value_increment(Mat->p[j][0],Mat->p[j][0]);
00984 }
00985
00986
00987
00988 if ((value_cmp_si(Mat->p[j][0], NbRay) == 0) &&
00989 (value_notzero_p(Mat->p[j][Dimension])))
00990 goto empty;
00991
00992
00993 NbConstraints--;
00994 if (j==NbConstraints) continue;
00995 Vector_Exchange(Mat->p[j], Mat->p[NbConstraints], temp1->Size);
00996 exchange(jx[j], jx[NbConstraints], aux);
00997 exchange(bx[j], bx[NbConstraints], aux);
00998 j--; continue;
00999 }
01000
01001
01002
01003 for (i=0; i<NbRay; i++)
01004 if (!(Sat->p[i][jx[j]]&bx[j])) {
01005
01006
01007 value_increment(Mat->p[j][0],Mat->p[j][0]);
01008
01009
01010 value_increment (Ray->p[i][0],Ray->p[i][0]);
01011 }
01012
01013
01014 if (value_cmp_si(Mat->p[j][0], NbRay) == 0)
01015 NbEq++;
01016 }
01017 Mat->NbRows = NbConstraints;
01018
01019 NbBid=0;
01020 for (i=0; i<NbRay; i++) {
01021
01022
01023 if (value_zero_p(Ray->p[i][Dimension]))
01024
01025
01026 value_increment(Ray->p[i][0],Ray->p[i][0]);
01027
01028
01029
01030
01031
01032
01033 if (value_cmp_si(Ray->p[i][0], NbConstraints+1) == 0)
01034 NbBid++;
01035 }
01036
01037 #ifdef POLY_RR_DEBUG
01038 fprintf(stderr, "[Remove_redundants : Step1]\nConstraints =");
01039 Matrix_Print(stderr,0,Mat);
01040 fprintf(stderr, "\nRay =");
01041 Matrix_Print(stderr,0,Ray);
01042 #endif
01043
01044
01045
01046
01047
01048
01049
01050
01051 for (i=0; i<NbEq; i++) {
01052
01053
01054 if (value_cmp_si(Mat->p[i][0], NbRay) != 0) {
01055
01056
01057 for (k=i+1; value_cmp_si(Mat->p[k][0], NbRay) != 0 && k < NbConstraints; k++)
01058 ;
01059 if (k==NbConstraints) break;
01060
01061
01062
01063 Vector_Copy(Mat->p[k], temp1->p, temp1->Size);
01064 aux = jx[k];
01065 j = bx[k];
01066 for (;k>i;k--) {
01067 Vector_Copy(Mat->p[k-1], Mat->p[k], temp1->Size);
01068 jx[k] = jx[k-1];
01069 bx[k] = bx[k-1];
01070 }
01071 Vector_Copy(temp1->p, Mat->p[i], temp1->Size);
01072 jx[i] = aux;
01073 bx[i] = j;
01074 }
01075 }
01076
01077
01078 if (Filter) {
01079 Value mone;
01080 value_init(mone);
01081 value_set_si(mone, -1);
01082
01083
01084
01085 for (i = 0; i < NbEq; i++) {
01086 int pos = First_Non_Zero(Mat->p[i]+1, Dimension);
01087 if (pos == -1)
01088 continue;
01089 if (value_neg_p(Mat->p[i][1+pos]))
01090 Vector_Scale(Mat->p[i]+1, Mat->p[i]+1, mone, Dimension);
01091 }
01092 value_clear(mone);
01093 for (i=0; i<NbEq; i++) {
01094
01095
01096 Redundant = 0;
01097 for (j=i+1; j<NbEq; j++) {
01098
01099 if (!(temp2[jx[j]] & bx[j]))
01100 continue;
01101
01102 if (Vector_Equal(Mat->p[i]+1, Mat->p[j]+1, Dimension)) {
01103 Redundant=1;
01104 break;
01105 }
01106 }
01107
01108
01109 if (!Redundant) Filter[jx[i]] |= bx[i];
01110 }
01111 }
01112
01113 #ifdef POLY_RR_DEBUG
01114 fprintf(stderr, "[Remove_redundants : Step2]\nConstraints =");
01115 Matrix_Print(stderr,0,Mat);
01116 fprintf(stderr, "\nRay =");
01117 Matrix_Print(stderr,0,Ray);
01118 #endif
01119
01120
01121
01122
01123
01124
01125
01126
01127 NbEq2 = Gauss(Mat,NbEq,Dimension);
01128
01129
01130
01131
01132 if (NbEq2 >= Dimension)
01133 goto empty;
01134
01135 #ifdef POLY_RR_DEBUG
01136 fprintf(stderr, "[Remove_redundants : Step3]\nConstraints =");
01137 Matrix_Print(stderr,0,Mat);
01138 fprintf(stderr, "\nRay =");
01139 Matrix_Print(stderr,0,Ray);
01140 #endif
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150 for (i=0, k=NbRay; i<NbBid && k>i; i++) {
01151
01152 if (value_cmp_si(Ray->p[i][0], NbConstraints+1) != 0) {
01153
01154
01155 while (--k > i && value_cmp_si(Ray->p[k][0], NbConstraints+1) != 0)
01156 ;
01157
01158
01159
01160 Vector_Exchange(Ray->p[i], Ray->p[k], temp1->Size);
01161 bexchange(Sat->p[i], Sat->p[k], RowSize2);
01162 }
01163 }
01164
01165 #ifdef POLY_RR_DEBUG
01166 fprintf(stderr, "[Remove_redundants : Step4]\nConstraints =");
01167 Matrix_Print(stderr,0,Mat);
01168 fprintf(stderr, "\nRay =");
01169 Matrix_Print(stderr,0,Ray);
01170 #endif
01171
01172
01173
01174
01175
01176
01177
01178
01179 NbBid2 = Gauss(Ray, NbBid, Dimension);
01180
01181 #ifdef POLY_RR_DEBUG
01182 fprintf(stderr, "[Remove_redundants : After Gauss]\nRay =");
01183 Matrix_Print(stderr,0,Ray);
01184 #endif
01185
01186
01187
01188 if (NbBid2>=Dimension) {
01189 errormsg1("RemoveRedundants", "rmrdt", "dimension error");
01190 goto empty;
01191 }
01192
01193
01194 Dim_RaySpace = Dimension-1-NbEq2-NbBid2;
01195
01196 #ifdef POLY_RR_DEBUG
01197 fprintf(stderr, "[Remove_redundants : Step5]\nConstraints =");
01198 Matrix_Print(stderr,0,Mat);
01199 fprintf(stderr, "\nRay =");
01200 Matrix_Print(stderr,0,Ray);
01201 #endif
01202
01203
01204
01205
01206
01207
01208
01209
01210 NbIneq=0;
01211 for (j=0; j<NbConstraints; j++) {
01212
01213
01214 i = First_Non_Zero(Mat->p[j]+1, Dimension-1);
01215
01216
01217
01218 if (i == -1) {
01219
01220
01221 if ((value_cmp_si(Mat->p[j][0], NbRay) == 0) &&
01222 (value_notzero_p(Mat->p[j][Dimension])))
01223 goto empty;
01224
01225
01226
01227 value_set_si(Mat->p[j][0],2);
01228 continue;
01229 }
01230
01231 Status = VALUE_TO_INT(Mat->p[j][0]);
01232
01233 if (Status == 0)
01234 value_set_si(Mat->p[j][0], 2);
01235 else if (Status < Dim_RaySpace)
01236 value_set_si(Mat->p[j][0], 2);
01237 else if (Status == NbRay)
01238 value_set_si(Mat->p[j][0], 0);
01239 else {
01240 NbIneq++;
01241 value_set_si(Mat->p[j][0], 1);
01242 }
01243 }
01244
01245 #ifdef POLY_RR_DEBUG
01246 fprintf(stderr, "[Remove_redundants : Step6]\nConstraints =");
01247 Matrix_Print(stderr,0,Mat);
01248 fprintf(stderr, "\nRay =");
01249 Matrix_Print(stderr,0,Ray);
01250 #endif
01251
01252
01253
01254
01255
01256
01257 NbUni=0;
01258 for (j=0; j<NbRay; j++) {
01259 Status = VALUE_TO_INT(Ray->p[j][0]);
01260
01261 if (Status < Dim_RaySpace)
01262 value_set_si(Ray->p[j][0], 2);
01263 else if (Status == NbConstraints+1)
01264 value_set_si(Ray->p[j][0], 0);
01265 else {
01266 NbUni++;
01267 value_set_si(Ray->p[j][0], 1);
01268 }
01269 }
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279 Pol = Polyhedron_Alloc(Dimension-1, NbIneq+NbEq2+1, NbUni+NbBid2);
01280 if (!Pol) {
01281 UNCATCH(any_exception_error);
01282 goto oom;
01283 }
01284 Pol->NbBid = NbBid2;
01285 Pol->NbEq = NbEq2;
01286
01287
01288 if (NbBid2) Vector_Copy(Ray->p[0], Pol->Ray[0], (Dimension+1)*NbBid2);
01289 if (NbEq2) Vector_Copy(Mat->p[0], Pol->Constraint[0], (Dimension+1)*NbEq2);
01290
01291 #ifdef POLY_RR_DEBUG
01292 fprintf(stderr, "[Remove_redundants : Step7]\nConstraints =");
01293 Matrix_Print(stderr,0,Mat);
01294 fprintf(stderr, "\nRay =");
01295 Matrix_Print(stderr,0,Ray);
01296 #endif
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311 Trace=(unsigned *)malloc(sat_nbcolumns * sizeof(unsigned));
01312 if(!Trace) {
01313 UNCATCH(any_exception_error);
01314 goto oom;
01315 }
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338 NbIneq2 = 0;
01339 for (j=NbEq; j<NbConstraints; j++) {
01340
01341
01342 if (value_one_p (Mat->p[j][0])) {
01343 for (k=0; k<sat_nbcolumns; k++) Trace[k]=0;
01344
01345
01346
01347 for (i=NbBid; i<NbRay; i++)
01348
01349
01350 if (value_one_p(Ray->p[i][0])) {
01351 if (!(Sat->p[i][jx[j]]&bx[j]))
01352 for (k=0; k<sat_nbcolumns; k++) Trace[k] |= Sat->p[i][k];
01353 }
01354
01355
01356
01357
01358 Redundant=0;
01359 for (i=NbEq; i<NbConstraints; i++) {
01360
01361
01362 if (value_one_p(Mat->p[i][0]) && (i!=j) && !(Trace[jx[i]] & bx[i])) {
01363 Redundant=1;
01364 break;
01365 }
01366 }
01367 if (Redundant) {
01368 value_set_si(Mat->p[j][0],2);
01369 }
01370 else {
01371 Vector_Copy(Mat->p[j], Pol->Constraint[NbEq2+NbIneq2], Dimension+1);
01372 if (Filter) Filter[jx[j]] |= bx[j];
01373 NbIneq2++;
01374 }
01375 }
01376 }
01377 free(Trace), Trace = NULL;
01378
01379 #ifdef POLY_RR_DEBUG
01380 fprintf(stderr, "[Remove_redundants : Step8]\nConstraints =");
01381 Matrix_Print(stderr,0,Mat);
01382 fprintf(stderr, "\nRay =");
01383 Matrix_Print(stderr,0,Ray);
01384 #endif
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395 Trace=(unsigned *)malloc(NbRay * sizeof(unsigned));
01396 if(!Trace) {
01397 UNCATCH(any_exception_error);
01398 goto oom;
01399 }
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418 NbUni2 = 0;
01419
01420
01421 aux = 0;
01422 for (i=NbBid; i<NbRay; i++) {
01423
01424
01425 if (value_one_p (Ray->p[i][0])) {
01426
01427
01428 if (value_notzero_p (Ray->p[i][Dimension]))
01429 for (k=NbBid; k<NbRay; k++) Trace[k]=0;
01430 else
01431
01432
01433
01434
01435 for (k=NbBid; k<NbRay; k++)
01436
01437
01438 Trace[k] = (value_notzero_p (Ray->p[k][Dimension]));
01439
01440
01441
01442 for (j=NbEq; j<NbConstraints; j++)
01443
01444
01445 if (value_one_p (Mat->p[j][0])) {
01446 if (!(Sat->p[i][jx[j]]&bx[j]))
01447 for (k=NbBid; k<NbRay; k++) Trace[k] |= Sat->p[k][jx[j]]&bx[j];
01448 }
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458 Redundant = 0;
01459 for (j=NbBid; j<NbRay; j++) {
01460
01461
01462 if (value_one_p (Ray->p[j][0]) && (i!=j) && !Trace[j]) {
01463 Redundant=1;
01464 break;
01465 }
01466 }
01467 if (Redundant)
01468 value_set_si(Ray->p[i][0],2);
01469 else {
01470 Vector_Copy(Ray->p[i], Pol->Ray[NbBid2+NbUni2], Dimension+1);
01471 NbUni2++;
01472
01473
01474 if (value_zero_p (Ray->p[i][Dimension]))
01475 aux++;
01476 }
01477 }
01478 }
01479
01480
01481 if (aux>=Dim_RaySpace) {
01482 Vector_Set(Pol->Constraint[NbEq2+NbIneq2],0,Dimension+1);
01483 value_set_si(Pol->Constraint[NbEq2+NbIneq2][0],1);
01484 value_set_si(Pol->Constraint[NbEq2+NbIneq2][Dimension],1);
01485 NbIneq2++;
01486 }
01487 }
01488
01489 UNCATCH(any_exception_error);
01490
01491 #ifdef POLY_RR_DEBUG
01492 fprintf(stderr, "[Remove_redundants : Step9]\nConstraints =");
01493 Matrix_Print(stderr,0,Mat);
01494 fprintf(stderr, "\nRay =");
01495 Matrix_Print(stderr,0,Ray);
01496 #endif
01497
01498 free(Trace);
01499 free(bx);
01500 free(jx);
01501 if (temp2)
01502 free(temp2);
01503
01504 Pol->NbConstraints = NbEq2 + NbIneq2;
01505 Pol->NbRays = NbBid2 + NbUni2;
01506
01507 Vector_Free(temp1);
01508 F_SET(Pol,
01509 POL_VALID | POL_INEQUALITIES | POL_FACETS | POL_POINTS | POL_VERTICES);
01510 return Pol;
01511
01512 oom:
01513 errormsg1("Remove_Redundants", "outofmem", "out of memory space");
01514
01515 Vector_Free(temp1);
01516 if (temp2)
01517 free(temp2);
01518 if (bx)
01519 free(bx);
01520 if (jx)
01521 free(jx);
01522 return NULL;
01523
01524 empty:
01525 Vector_Free(temp1);
01526 if (temp2)
01527 free(temp2);
01528 if (bx)
01529 free(bx);
01530 if (jx)
01531 free(jx);
01532 UNCATCH(any_exception_error);
01533 return Empty_Polyhedron(Dimension-1);
01534 }
01535
01536
01537
01538
01539 Polyhedron* Polyhedron_Alloc(unsigned Dimension,unsigned NbConstraints,unsigned NbRays) {
01540
01541 Polyhedron *Pol;
01542 unsigned NbRows,NbColumns;
01543 int i,j;
01544 Value *p, **q;
01545
01546 Pol=(Polyhedron *)malloc(sizeof(Polyhedron));
01547 if(!Pol) {
01548 errormsg1("Polyhedron_Alloc", "outofmem", "out of memory space");
01549 return 0;
01550 }
01551
01552 Pol->next = (Polyhedron *)0;
01553 Pol->Dimension = Dimension;
01554 Pol->NbConstraints = NbConstraints;
01555 Pol->NbRays = NbRays;
01556 Pol->NbEq = 0;
01557 Pol->NbBid = 0;
01558 Pol->flags = 0;
01559 NbRows = NbConstraints + NbRays;
01560 NbColumns = Dimension + 2;
01561
01562 q = (Value **)malloc(NbRows * sizeof(Value *));
01563 if(!q) {
01564 errormsg1("Polyhedron_Alloc", "outofmem", "out of memory space");
01565 return 0;
01566 }
01567 p = value_alloc(NbRows*NbColumns, &Pol->p_Init_size);
01568 if(!p) {
01569 free(q);
01570 errormsg1("Polyhedron_Alloc", "outofmem", "out of memory space");
01571 return 0;
01572 }
01573 Pol->Constraint = q;
01574 Pol->Ray = q + NbConstraints;
01575 Pol->p_Init = p;
01576 for (i=0;i<NbRows;i++) {
01577 *q++ = p;
01578 p += NbColumns;
01579 }
01580 return Pol;
01581 }
01582
01583
01584
01585
01586 void Polyhedron_Free(Polyhedron *Pol)
01587 {
01588 if(!Pol)
01589 return;
01590 value_free(Pol->p_Init, Pol->p_Init_size);
01591 free(Pol->Constraint);
01592 free(Pol);
01593 return;
01594 }
01595
01596
01597
01598
01599 void Domain_Free(Polyhedron *Pol)
01600 {
01601 Polyhedron *Next;
01602
01603 for (; Pol; Pol = Next) {
01604 Next = Pol->next;
01605 Polyhedron_Free(Pol);
01606 }
01607 return;
01608 }
01609
01610
01611
01612
01613 void Polyhedron_Print(FILE *Dst, const char *Format, const Polyhedron *Pol)
01614 {
01615 unsigned Dimension, NbConstraints, NbRays;
01616 int i, j;
01617 Value *p;
01618
01619 if (!Pol) {
01620 fprintf(Dst, "<null polyhedron>\n");
01621 return;
01622 }
01623
01624 Dimension = Pol->Dimension + 2;
01625 NbConstraints = Pol->NbConstraints;
01626 NbRays = Pol->NbRays;
01627 fprintf(Dst, "POLYHEDRON Dimension:%d\n", Pol->Dimension);
01628 fprintf(Dst," Constraints:%d Equations:%d Rays:%d Lines:%d\n",
01629 Pol->NbConstraints, Pol->NbEq, Pol->NbRays, Pol->NbBid);
01630 fprintf(Dst,"Constraints %d %d\n", NbConstraints, Dimension);
01631
01632 for (i=0;i<NbConstraints;i++) {
01633 p=Pol->Constraint[i];
01634
01635
01636 if (value_notzero_p (*p))
01637 fprintf(Dst,"Inequality: [");
01638 else
01639 fprintf(Dst,"Equality: [");
01640 p++;
01641 for (j=1;j<Dimension;j++) {
01642 value_print(Dst,Format,*p++);
01643 }
01644 (void)fprintf(Dst," ]\n");
01645 }
01646
01647 (void)fprintf(Dst, "Rays %d %d\n", NbRays, Dimension);
01648 for (i=0;i<NbRays;i++) {
01649 p=Pol->Ray[i];
01650
01651
01652 if (value_notzero_p (*p)) {
01653 p++;
01654
01655
01656 if (value_notzero_p (p[Dimension-2]))
01657 fprintf(Dst, "Vertex: [");
01658 else
01659 fprintf(Dst, "Ray: [");
01660 }
01661 else {
01662 p++;
01663 fprintf(Dst, "Line: [");
01664 }
01665 for (j=1; j < Dimension-1; j++) {
01666 value_print(Dst,Format,*p++);
01667 }
01668
01669
01670 if (value_notzero_p (*p)) {
01671 fprintf( Dst, " ]/" );
01672 value_print(Dst,VALUE_FMT,*p);
01673 fprintf( Dst, "\n" );
01674 }
01675 else
01676 fprintf(Dst, " ]\n");
01677 }
01678 if (Pol->next) {
01679 fprintf(Dst, "UNION ");
01680 Polyhedron_Print(Dst,Format,Pol->next);
01681 }
01682 }
01683
01684
01685
01686
01687 void PolyPrint (Polyhedron *Pol) {
01688 Polyhedron_Print(stderr,"%4d",Pol);
01689 }
01690
01691
01692
01693
01694
01695
01696
01697
01698 Polyhedron *Empty_Polyhedron(unsigned Dimension) {
01699
01700 Polyhedron *Pol;
01701 int i;
01702
01703 Pol = Polyhedron_Alloc(Dimension, Dimension+1, 0);
01704 if (!Pol) {
01705 errormsg1("Empty_Polyhedron", "outofmem", "out of memory space");
01706 return 0;
01707 }
01708 Vector_Set(Pol->Constraint[0],0,(Dimension+1)*(Dimension+2));
01709 for (i=0; i<=Dimension; i++) {
01710
01711
01712 value_set_si(Pol->Constraint[i][i+1],1);
01713 }
01714 Pol->NbEq = Dimension+1;
01715 Pol->NbBid = 0;
01716 F_SET(Pol,
01717 POL_VALID | POL_INEQUALITIES | POL_FACETS | POL_POINTS | POL_VERTICES);
01718 return Pol;
01719 }
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731 Polyhedron *Universe_Polyhedron(unsigned Dimension) {
01732
01733 Polyhedron *Pol;
01734 int i;
01735
01736 Pol = Polyhedron_Alloc(Dimension,1,Dimension+1);
01737 if (!Pol) {
01738 errormsg1("Universe_Polyhedron", "outofmem", "out of memory space");
01739 return 0;
01740 }
01741 Vector_Set(Pol->Constraint[0],0,(Dimension+2));
01742
01743
01744 value_set_si(Pol->Constraint[0][0],1);
01745
01746
01747 value_set_si(Pol->Constraint[0][Dimension+1],1);
01748 Vector_Set(Pol->Ray[0],0,(Dimension+1)*(Dimension+2));
01749 for (i=0;i<=Dimension;i++) {
01750
01751
01752 value_set_si(Pol->Ray[i][i+1],1);
01753 }
01754
01755
01756 value_set_si(Pol->Ray[Dimension][0],1);
01757 Pol->NbEq = 0;
01758 Pol->NbBid = Dimension;
01759 F_SET(Pol,
01760 POL_VALID | POL_INEQUALITIES | POL_FACETS | POL_POINTS | POL_VERTICES);
01761 return Pol;
01762 }
01763
01764
01765
01766
01767
01768
01769 static void SortConstraints(Matrix *Constraints, unsigned NbEq)
01770 {
01771 int i, j, k;
01772 for (i = NbEq; i < Constraints->NbRows; ++i) {
01773 int max = i;
01774 for (k = i+1; k < Constraints->NbRows; ++k) {
01775 for (j = 1; j < Constraints->NbColumns-1; ++j) {
01776 if (value_eq(Constraints->p[k][j],
01777 Constraints->p[max][j]))
01778 continue;
01779 if (value_abs_lt(Constraints->p[k][j],
01780 Constraints->p[max][j]))
01781 break;
01782 if (value_abs_eq(Constraints->p[k][j],
01783 Constraints->p[max][j]) &&
01784 value_pos_p(Constraints->p[max][j]))
01785 break;
01786 max = k;
01787 break;
01788 }
01789
01790
01791
01792 if (j == Constraints->NbColumns-1) {
01793 if (value_lt(Constraints->p[k][j], Constraints->p[max][j]))
01794 Vector_Exchange(Constraints->p[k],
01795 Constraints->p[max],
01796 Constraints->NbColumns);
01797 Constraints->NbRows--;
01798 if (k < Constraints->NbRows)
01799 Vector_Exchange(Constraints->p[k],
01800 Constraints->p[Constraints->NbRows],
01801 Constraints->NbColumns);
01802 k--;
01803 }
01804 }
01805 if (max != i)
01806 Vector_Exchange(Constraints->p[max], Constraints->p[i],
01807 Constraints->NbColumns);
01808 }
01809 }
01810
01811
01812
01813
01814
01815
01816
01817
01818 static int ImplicitEqualities(Matrix *Constraints, unsigned NbEq)
01819 {
01820 int row, nrow, k;
01821 int found = 0;
01822 Value tmp;
01823 for (row = NbEq; row < Constraints->NbRows; ++row) {
01824 int d = First_Non_Zero(Constraints->p[row]+1, Constraints->NbColumns-2);
01825 if (d == -1) {
01826
01827 if (value_neg_p(Constraints->p[row][Constraints->NbColumns-1])) {
01828 value_set_si(Constraints->p[row][0], 0);
01829 found = 1;
01830 }
01831 break;
01832 }
01833 if (value_neg_p(Constraints->p[row][1+d]))
01834 continue;
01835 for (nrow = row+1; nrow < Constraints->NbRows; ++nrow) {
01836 if (value_zero_p(Constraints->p[nrow][1+d]))
01837 break;
01838 if (value_pos_p(Constraints->p[nrow][1+d]))
01839 continue;
01840 for (k = d; k < Constraints->NbColumns-1; ++k) {
01841 if (value_abs_ne(Constraints->p[row][1+k],
01842 Constraints->p[nrow][1+k]))
01843 break;
01844 if (value_zero_p(Constraints->p[row][1+k]))
01845 continue;
01846 if (value_eq(Constraints->p[row][1+k], Constraints->p[nrow][1+k]))
01847 break;
01848 }
01849 if (k == Constraints->NbColumns-1) {
01850 value_set_si(Constraints->p[row][0], 0);
01851 value_set_si(Constraints->p[nrow][0], 0);
01852 found = 1;
01853 break;
01854 }
01855 if (k != Constraints->NbColumns-2)
01856 continue;
01857
01858
01859
01860 value_init(tmp);
01861 value_addto(tmp, Constraints->p[row][1+k],
01862 Constraints->p[nrow][1+k]);
01863 if (value_sign(tmp) < 0) {
01864 Vector_Set(Constraints->p[row], 0, Constraints->NbColumns-1);
01865 Vector_Set(Constraints->p[nrow], 0, Constraints->NbColumns-1);
01866 value_set_si(Constraints->p[row][1+k], 1);
01867 value_set_si(Constraints->p[nrow][1+k], 1);
01868 found = 1;
01869 }
01870 value_clear(tmp);
01871 if (found)
01872 break;
01873 }
01874 }
01875 return found;
01876 }
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889 Polyhedron *Constraints2Polyhedron(Matrix *Constraints,unsigned NbMaxRays) {
01890
01891 Polyhedron *Pol = NULL;
01892 Matrix *Ray = NULL;
01893 SatMatrix *Sat = NULL;
01894 unsigned Dimension, nbcolumns;
01895 int i;
01896
01897 Dimension = Constraints->NbColumns - 1;
01898 if (Dimension < 1) {
01899 errormsg1("Constraints2Polyhedron","invalidpoly","invalid polyhedron dimension");
01900 return 0;
01901 }
01902
01903
01904
01905 if (Constraints->NbRows==0) {
01906 Pol = Universe_Polyhedron(Dimension-1);
01907 return Pol;
01908 }
01909
01910 if (POL_ISSET(NbMaxRays, POL_NO_DUAL)) {
01911 unsigned NbEq;
01912 unsigned Rank;
01913 Value tmp;
01914 if (POL_ISSET(NbMaxRays, POL_INTEGER))
01915 value_init(tmp);
01916 do {
01917 NbEq = 0;
01918
01919 for (i = 0; i < Constraints->NbRows; ++i)
01920 if (value_zero_p(Constraints->p[i][0])) {
01921 if (POL_ISSET(NbMaxRays, POL_INTEGER) &&
01922 ConstraintSimplify(Constraints->p[i],
01923 Constraints->p[i], Dimension+1, &tmp)) {
01924 value_clear(tmp);
01925 return Empty_Polyhedron(Dimension-1);
01926 }
01927
01928 if (First_Non_Zero(Constraints->p[i]+1, Dimension-1) == -1 &&
01929 value_notzero_p(Constraints->p[i][Dimension])) {
01930 if (POL_ISSET(NbMaxRays, POL_INTEGER))
01931 value_clear(tmp);
01932 return Empty_Polyhedron(Dimension-1);
01933 }
01934 if (i != NbEq)
01935 ExchangeRows(Constraints, i, NbEq);
01936 ++NbEq;
01937 }
01938 Rank = Gauss(Constraints, NbEq, Dimension);
01939 if (POL_ISSET(NbMaxRays, POL_INTEGER))
01940 for (i = NbEq; i < Constraints->NbRows; ++i)
01941 ConstraintSimplify(Constraints->p[i],
01942 Constraints->p[i], Dimension+1, &tmp);
01943 SortConstraints(Constraints, NbEq);
01944 } while (ImplicitEqualities(Constraints, NbEq));
01945 if (POL_ISSET(NbMaxRays, POL_INTEGER))
01946 value_clear(tmp);
01947 Pol = Polyhedron_Alloc(Dimension-1, Constraints->NbRows - (NbEq-Rank), 0);
01948 if (Rank > 0)
01949 Vector_Copy(Constraints->p[0], Pol->Constraint[0],
01950 Rank * Constraints->NbColumns);
01951 if (Constraints->NbRows > NbEq)
01952 Vector_Copy(Constraints->p[NbEq], Pol->Constraint[Rank],
01953 (Constraints->NbRows - NbEq) * Constraints->NbColumns);
01954 Pol->NbEq = Rank;
01955
01956 Pol->Ray = 0;
01957 F_SET(Pol, POL_VALID | POL_INEQUALITIES);
01958 return Pol;
01959 }
01960
01961 if (Dimension > NbMaxRays)
01962 NbMaxRays = Dimension;
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974 Ray = Matrix_Alloc(NbMaxRays, Dimension+1);
01975 if(!Ray) {
01976 errormsg1("Constraints2Polyhedron","outofmem","out of memory space");
01977 return 0;
01978 }
01979 Vector_Set(Ray->p_Init,0, NbMaxRays * (Dimension+1));
01980 for (i=0; i<Dimension; i++) {
01981
01982
01983 value_set_si(Ray->p[i][i+1],1);
01984 }
01985
01986
01987 value_set_si(Ray->p[Dimension-1][0],1);
01988 Ray->NbRows = Dimension;
01989
01990
01991 nbcolumns = (Constraints->NbRows - 1)/(sizeof(int)*8) + 1;
01992 Sat = SMAlloc(NbMaxRays, nbcolumns);
01993 SMVector_Init(Sat->p_init,Dimension*nbcolumns);
01994 Sat->NbRows = Dimension;
01995
01996 CATCH(any_exception_error) {
01997
01998
01999 if (Sat) SMFree(&Sat);
02000 if (Ray) Matrix_Free(Ray);
02001 if (Pol) Polyhedron_Free(Pol);
02002 RETHROW();
02003 }
02004 TRY {
02005
02006
02007 Chernikova(Constraints,Ray,Sat,Dimension-1,NbMaxRays,0,0);
02008
02009 #ifdef POLY_DEBUG
02010 fprintf(stderr, "[constraints2polyhedron]\nConstraints = ");
02011 Matrix_Print(stderr,0,Constraints);
02012 fprintf(stderr, "\nRay = ");
02013 Matrix_Print(stderr,0,Ray);
02014 fprintf(stderr, "\nSat = ");
02015 SMPrint(Sat);
02016 #endif
02017
02018
02019 Pol = Remove_Redundants(Constraints,Ray,Sat,0);
02020 }
02021
02022 UNCATCH(any_exception_error);
02023
02024 #ifdef POLY_DEBUG
02025 fprintf(stderr, "\nPol = ");
02026 Polyhedron_Print(stderr,"%4d",Pol);
02027 #endif
02028
02029 SMFree(&Sat), Sat = NULL;
02030 Matrix_Free(Ray), Ray = NULL;
02031 return Pol;
02032 }
02033
02034 #undef POLY_DEBUG
02035
02036
02037
02038
02039 Matrix *Polyhedron2Constraints(Polyhedron *Pol) {
02040
02041 Matrix *Mat;
02042 unsigned NbConstraints,Dimension;
02043
02044 POL_ENSURE_INEQUALITIES(Pol);
02045
02046 NbConstraints = Pol->NbConstraints;
02047 Dimension = Pol->Dimension+2;
02048 Mat = Matrix_Alloc(NbConstraints,Dimension);
02049 if(!Mat) {
02050 errormsg1("Polyhedron2Constraints", "outofmem", "out of memory space");
02051 return 0;
02052 }
02053 Vector_Copy(Pol->Constraint[0],Mat->p_Init,NbConstraints * Dimension);
02054 return Mat;
02055 }
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065
02066 Polyhedron *Rays2Polyhedron(Matrix *Ray,unsigned NbMaxConstrs) {
02067
02068 Polyhedron *Pol = NULL;
02069 Matrix *Mat = NULL;
02070 SatMatrix *Sat = NULL, *SatTranspose = NULL;
02071 unsigned Dimension, nbcolumns;
02072 int i;
02073
02074 Dimension = Ray->NbColumns-1;
02075 Sat = NULL;
02076 SatTranspose = NULL;
02077 Mat = NULL;
02078
02079 if (Ray->NbRows==0) {
02080
02081
02082 Pol = Empty_Polyhedron(Dimension-1);
02083 return(Pol);
02084 }
02085
02086
02087 if (POL_ISSET(NbMaxConstrs, POL_NO_DUAL))
02088 NbMaxConstrs = 0;
02089
02090 if (Dimension > NbMaxConstrs)
02091 NbMaxConstrs = Dimension;
02092
02093
02094 Mat = Matrix_Alloc(NbMaxConstrs,Dimension+1);
02095 if(!Mat) {
02096 errormsg1("Rays2Polyhedron","outofmem","out of memory space");
02097 return 0;
02098 }
02099
02100
02101 Vector_Set(Mat->p_Init,0,NbMaxConstrs * (Dimension+1));
02102 for (i=0; i<Dimension; i++) {
02103
02104
02105 value_set_si(Mat->p[i][i+1],1);
02106 }
02107
02108
02109
02110 Mat->NbRows = Dimension;
02111 nbcolumns = (Ray->NbRows -1)/(sizeof(int)*8) + 1;
02112 SatTranspose = SMAlloc(NbMaxConstrs,nbcolumns);
02113 SMVector_Init(SatTranspose->p[0],Dimension * nbcolumns);
02114 SatTranspose->NbRows = Dimension;
02115
02116 #ifdef POLY_DEBUG
02117 fprintf(stderr, "[ray2polyhedron: Before]\nRay = ");
02118 Matrix_Print(stderr,0,Ray);
02119 fprintf(stderr, "\nConstraints = ");
02120 Matrix_Print(stderr,0,Mat);
02121 fprintf(stderr, "\nSatTranspose = ");
02122 SMPrint (SatTranspose);
02123 #endif
02124
02125 CATCH(any_exception_error) {
02126
02127
02128
02129
02130 if (SatTranspose) SMFree(&SatTranspose);
02131 if (Sat) SMFree(&Sat);
02132 if (Mat) Matrix_Free(Mat);
02133 if (Pol) Polyhedron_Free(Pol);
02134 RETHROW();
02135 }
02136 TRY {
02137
02138
02139 Chernikova(Ray,Mat,SatTranspose,Dimension,NbMaxConstrs,0,1);
02140
02141 #ifdef POLY_DEBUG
02142 fprintf(stderr, "[ray2polyhedron: After]\nRay = ");
02143 Matrix_Print(stderr,0,Ray);
02144 fprintf(stderr, "\nConstraints = ");
02145 Matrix_Print(stderr,0,Mat);
02146 fprintf(stderr, "\nSatTranspose = ");
02147 SMPrint (SatTranspose);
02148 #endif
02149
02150
02151
02152 Sat = TransformSat(Mat,Ray,SatTranspose);
02153
02154 #ifdef POLY_DEBUG
02155 fprintf(stderr, "\nSat =");
02156 SMPrint(Sat);
02157 #endif
02158
02159 SMFree(&SatTranspose), SatTranspose = NULL;
02160
02161
02162 Pol = Remove_Redundants(Mat,Ray,Sat,0);
02163 }
02164
02165 UNCATCH(any_exception_error);
02166
02167 #ifdef POLY_DEBUG
02168 fprintf(stderr, "\nPol = ");
02169 Polyhedron_Print(stderr,"%4d",Pol);
02170 #endif
02171
02172 SMFree(&Sat);
02173 Matrix_Free(Mat);
02174 return Pol;
02175 }
02176
02177
02178
02179
02180
02181 void Polyhedron_Compute_Dual(Polyhedron *P)
02182 {
02183 if (!F_ISSET(P, POL_VALID))
02184 return;
02185
02186 if (F_ISSET(P, POL_FACETS | POL_VERTICES))
02187 return;
02188
02189 if (F_ISSET(P, POL_INEQUALITIES)) {
02190 Matrix M;
02191 Polyhedron tmp, *Q;
02192
02193 M.NbRows = P->NbConstraints;
02194 M.NbColumns = P->Dimension+2;
02195 M.p_Init = P->p_Init;
02196 M.p = P->Constraint;
02197 Q = Constraints2Polyhedron(&M, 0);
02198
02199
02200 tmp = *Q;
02201 *Q = *P;
02202 *P = tmp;
02203
02204 P->next = Q->next;
02205 Polyhedron_Free(Q);
02206 return;
02207 }
02208
02209
02210 assert(0);
02211 }
02212
02213
02214
02215
02216
02217
02218
02219
02220 static SatMatrix *BuildSat(Matrix *Mat,Matrix *Ray,unsigned NbConstraints,unsigned NbMaxRays) {
02221
02222 SatMatrix *Sat = NULL;
02223 int i, j, k, jx;
02224 Value *p1, *p2, *p3;
02225 unsigned Dimension, NbRay, bx, nbcolumns;
02226
02227 CATCH(any_exception_error) {
02228 if (Sat)
02229 SMFree(&Sat);
02230 RETHROW();
02231 }
02232 TRY {
02233 NbRay = Ray->NbRows;
02234 Dimension = Mat->NbColumns-1;
02235
02236
02237 nbcolumns = (Mat->NbRows - 1)/(sizeof(int)*8) + 1;
02238 Sat = SMAlloc(NbMaxRays,nbcolumns);
02239 Sat->NbRows = NbRay;
02240 SMVector_Init(Sat->p_init, nbcolumns * NbRay);
02241 jx=0; bx=MSB;
02242 for (k=0; k<NbConstraints; k++) {
02243 for (i=0; i<NbRay; i++) {
02244
02245
02246
02247 p1 = Ray->p[i]+1;
02248 p2 = Mat->p[k]+1;
02249 p3 = Ray->p[i];
02250 value_set_si(*p3,0);
02251 for (j=0; j<Dimension; j++) {
02252 value_addmul(*p3, *p1, *p2);
02253 p1++; p2++;
02254 }
02255 }
02256 for (j=0; j<NbRay; j++) {
02257
02258
02259
02260 if (value_notzero_p(Ray->p[j][0]))
02261 Sat->p[j][jx]|=bx;
02262 }
02263 NEXT(jx, bx);
02264 }
02265 }
02266
02267 UNCATCH(any_exception_error);
02268 return Sat;
02269 }
02270
02271
02272
02273
02274
02275
02276 Polyhedron *AddConstraints(Value *Con,unsigned NbConstraints,Polyhedron *Pol,unsigned NbMaxRays) {
02277
02278 Polyhedron *NewPol = NULL;
02279 Matrix *Mat = NULL, *Ray = NULL;
02280 SatMatrix *Sat = NULL;
02281 unsigned NbRay, NbCon, Dimension;
02282
02283 if (NbConstraints == 0)
02284 return Polyhedron_Copy(Pol);
02285
02286 POL_ENSURE_INEQUALITIES(Pol);
02287 Dimension = Pol->Dimension + 2;
02288
02289 if (POL_ISSET(NbMaxRays, POL_NO_DUAL)) {
02290 NbCon = Pol->NbConstraints + NbConstraints;
02291
02292 Mat = Matrix_Alloc(NbCon, Dimension);
02293 if (!Mat) {
02294 errormsg1("AddConstraints", "outofmem", "out of memory space");
02295 return NULL;
02296 }
02297 Vector_Copy(Pol->Constraint[0], Mat->p[0], Pol->NbConstraints * Dimension);
02298 Vector_Copy(Con, Mat->p[Pol->NbConstraints], NbConstraints * Dimension);
02299 NewPol = Constraints2Polyhedron(Mat, NbMaxRays);
02300 Matrix_Free(Mat);
02301 return NewPol;
02302 }
02303
02304 POL_ENSURE_VERTICES(Pol);
02305
02306 CATCH(any_exception_error) {
02307 if (NewPol) Polyhedron_Free(NewPol);
02308 if (Mat) Matrix_Free(Mat);
02309 if (Ray) Matrix_Free(Ray);
02310 if (Sat) SMFree(&Sat);
02311 RETHROW();
02312 }
02313 TRY {
02314 NbRay = Pol->NbRays;
02315 NbCon = Pol->NbConstraints + NbConstraints;
02316
02317 if (NbRay > NbMaxRays)
02318 NbMaxRays = NbRay;
02319
02320 Mat = Matrix_Alloc(NbCon, Dimension);
02321 if(!Mat) {
02322 errormsg1("AddConstraints", "outofmem", "out of memory space");
02323 UNCATCH(any_exception_error);
02324 return 0;
02325 }
02326
02327
02328 Vector_Copy(Pol->Constraint[0], Mat->p[0], Pol->NbConstraints * Dimension);
02329
02330
02331 Vector_Copy(Con, Mat->p[Pol->NbConstraints], NbConstraints * Dimension);
02332
02333
02334 Ray = Matrix_Alloc(NbMaxRays, Dimension);
02335 if(!Ray) {
02336 errormsg1("AddConstraints", "outofmem", "out of memory space");
02337 UNCATCH(any_exception_error);
02338 return 0;
02339 }
02340 Ray->NbRows = NbRay;
02341
02342
02343 if (NbRay)
02344 Vector_Copy(Pol->Ray[0], Ray->p[0], NbRay * Dimension);
02345
02346
02347
02348 Sat = BuildSat(Mat, Ray, Pol->NbConstraints, NbMaxRays);
02349
02350
02351 Pol_status = Chernikova(Mat, Ray, Sat, Pol->NbBid, NbMaxRays, Pol->NbConstraints,0);
02352
02353
02354 NewPol = Remove_Redundants(Mat, Ray, Sat, 0);
02355
02356 }
02357
02358 UNCATCH(any_exception_error);
02359 SMFree(&Sat);
02360 Matrix_Free(Ray);
02361 Matrix_Free(Mat);
02362 return NewPol;
02363 }
02364
02365
02366
02367
02368
02369
02370
02371 int PolyhedronIncludes(Polyhedron *Pol1,Polyhedron *Pol2) {
02372
02373 int Dimension = Pol1->Dimension + 1;
02374 int i, j, k;
02375 Value *p1, *p2, p3;
02376
02377 POL_ENSURE_FACETS(Pol1);
02378 POL_ENSURE_VERTICES(Pol1);
02379 POL_ENSURE_FACETS(Pol2);
02380 POL_ENSURE_VERTICES(Pol2);
02381
02382 value_init(p3);
02383 for (k=0; k<Pol1->NbConstraints; k++) {
02384 for (i=0;i<Pol2->NbRays;i++) {
02385
02386
02387 p1 = Pol2->Ray[i]+1;
02388 p2 = Pol1->Constraint[k]+1;
02389 value_set_si(p3,0);
02390 for(j=0;j<Dimension;j++) {
02391 value_addmul(p3, *p1,*p2);
02392 p1++; p2++;
02393 }
02394
02395
02396
02397 if(value_neg_p(p3) ||
02398 (value_notzero_p(p3)
02399 && (value_zero_p(Pol1->Constraint[k][0]) || (value_zero_p(Pol2->Ray[i][0])) ) )) {
02400 value_clear(p3);
02401 return 0;
02402 }
02403 }
02404 }
02405 value_clear(p3);
02406 return 1;
02407 }
02408
02409
02410
02411
02412
02413
02414
02415 Polyhedron *AddPolyToDomain(Polyhedron *Pol,Polyhedron *PolDomain) {
02416
02417 Polyhedron *p, *pnext, *p_domain_end = (Polyhedron *) 0;
02418 int Redundant;
02419
02420 if (!Pol)
02421 return PolDomain;
02422 if (!PolDomain)
02423 return Pol;
02424
02425 POL_ENSURE_FACETS(Pol);
02426 POL_ENSURE_VERTICES(Pol);
02427
02428
02429 if (emptyQ(Pol)) {
02430 Polyhedron_Free(Pol);
02431 return PolDomain;
02432 }
02433
02434 POL_ENSURE_FACETS(PolDomain);
02435 POL_ENSURE_VERTICES(PolDomain);
02436
02437
02438 if (emptyQ(PolDomain)) {
02439 Polyhedron_Free(PolDomain);
02440 return Pol;
02441 }
02442
02443
02444 Redundant = 0;
02445 for (p=PolDomain,PolDomain=(Polyhedron *)0; p; p=pnext) {
02446
02447
02448 if (PolyhedronIncludes(Pol, p))
02449 {
02450
02451 pnext = p->next;
02452 Polyhedron_Free( p );
02453 continue;
02454 }
02455
02456
02457 if (!PolDomain) PolDomain = p; else p_domain_end->next = p;
02458 p_domain_end = p;
02459
02460
02461 if (PolyhedronIncludes(p,Pol)) {
02462 Redundant = 1;
02463 break;
02464 }
02465 pnext = p->next;
02466 }
02467 if (!Redundant) {
02468
02469
02470
02471 if (!PolDomain) PolDomain = Pol; else p_domain_end->next = Pol;
02472 }
02473 else {
02474
02475
02476 Polyhedron_Free(Pol);
02477 }
02478 return PolDomain;
02479 }
02480
02481
02482
02483
02484
02485
02486
02487
02488
02489
02490
02491 Polyhedron *SubConstraint(Value *Con,Polyhedron *Pol,unsigned NbMaxRays,int Pass) {
02492
02493 Polyhedron *NewPol = NULL;
02494 Matrix *Mat = NULL, *Ray = NULL;
02495 SatMatrix *Sat = NULL;
02496 unsigned NbRay, NbCon, NbEle1, Dimension;
02497 int i;
02498
02499 POL_ENSURE_FACETS(Pol);
02500 POL_ENSURE_VERTICES(Pol);
02501
02502 CATCH(any_exception_error) {
02503 if (NewPol) Polyhedron_Free(NewPol);
02504 if (Mat) Matrix_Free(Mat);
02505 if (Ray) Matrix_Free(Ray);
02506 if (Sat) SMFree(&Sat);
02507 RETHROW();
02508 }
02509 TRY {
02510
02511
02512 Dimension = Pol->Dimension+1;
02513 for (i=1; i<Dimension; i++)
02514 if (value_notzero_p(Con[i])) break;
02515 if (i==Dimension) {
02516 UNCATCH(any_exception_error);
02517 return (Polyhedron *) 0;
02518 }
02519
02520 NbRay = Pol->NbRays;
02521 NbCon = Pol->NbConstraints;
02522 Dimension = Pol->Dimension+2;
02523 NbEle1 = NbCon * Dimension;
02524
02525
02526 if (POL_ISSET(NbMaxRays, POL_NO_DUAL))
02527 NbMaxRays = 0;
02528
02529 if (NbRay > NbMaxRays)
02530 NbMaxRays = NbRay;
02531
02532 Mat = Matrix_Alloc(NbCon + 1, Dimension);
02533 if(!Mat) {
02534 errormsg1("SubConstraint", "outofmem", "out of memory space");
02535 UNCATCH(any_exception_error);
02536 return 0;
02537 }
02538
02539
02540 Vector_Copy(Pol->Constraint[0], Mat->p[0], NbEle1);
02541
02542
02543 value_set_si(Mat->p[NbCon][0],1);
02544 if (!(Pass&1))
02545 for(i=1; i<Dimension; i++)
02546 value_oppose(Mat->p[NbCon][i],Con[i]);
02547 else
02548 for(i=1; i<Dimension; i++)
02549 value_assign(Mat->p[NbCon][i],Con[i]);
02550 if (!(Pass&2))
02551 value_decrement(Mat->p[NbCon][Dimension-1],Mat->p[NbCon][Dimension-1]);
02552
02553
02554 Ray = Matrix_Alloc(NbMaxRays, Dimension);
02555 if(!Ray) {
02556 errormsg1("SubConstraint", "outofmem", "out of memory space");
02557 UNCATCH(any_exception_error);
02558 return 0;
02559 }
02560
02561
02562 Ray->NbRows = NbRay;
02563 if (NbRay)
02564 Vector_Copy(Pol->Ray[0], Ray->p[0], NbRay * Dimension);
02565
02566
02567
02568 Sat = BuildSat(Mat, Ray, NbCon, NbMaxRays);
02569
02570
02571 Pol_status = Chernikova(Mat, Ray, Sat, Pol->NbBid, NbMaxRays, NbCon,0);
02572
02573
02574 NewPol = Remove_Redundants(Mat, Ray, Sat, 0);
02575
02576 }
02577
02578 UNCATCH(any_exception_error);
02579
02580 SMFree(&Sat);
02581 Matrix_Free(Ray);
02582 Matrix_Free(Mat);
02583 return NewPol;
02584 }
02585
02586
02587
02588
02589
02590 Polyhedron *DomainIntersection(Polyhedron *Pol1,Polyhedron *Pol2,unsigned NbMaxRays) {
02591
02592 Polyhedron *p1, *p2, *p3, *d;
02593
02594 if (!Pol1 || !Pol2) return (Polyhedron*) 0;
02595 if (Pol1->Dimension != Pol2->Dimension) {
02596 errormsg1( "DomainIntersection", "diffdim",
02597 "operation on different dimensions");
02598 return (Polyhedron*) 0;
02599 }
02600
02601
02602
02603
02604 d = (Polyhedron *)0;
02605 for (p1=Pol1; p1; p1=p1->next) {
02606 for (p2=Pol2; p2; p2=p2->next) {
02607 p3 = AddConstraints(p2->Constraint[0],
02608 p2->NbConstraints, p1, NbMaxRays);
02609 d = AddPolyToDomain(p3,d);
02610 }
02611 }
02612 if (!d)
02613 return Empty_Polyhedron(Pol1->Dimension);
02614 else
02615 return d;
02616
02617 }
02618
02619
02620
02621
02622 Matrix *Polyhedron2Rays(Polyhedron *Pol) {
02623
02624 Matrix *Ray;
02625 unsigned NbRays, Dimension;
02626
02627 POL_ENSURE_POINTS(Pol);
02628
02629 NbRays = Pol->NbRays;
02630 Dimension = Pol->Dimension+2;
02631 Ray = Matrix_Alloc(NbRays, Dimension);
02632 if(!Ray) {
02633 errormsg1("Polyhedron2Rays", "outofmem", "out of memory space");
02634 return 0;
02635 }
02636 Vector_Copy(Pol->Ray[0], Ray->p_Init, NbRays*Dimension);
02637 return Ray;
02638 }
02639
02640
02641
02642
02643
02644 Polyhedron *AddRays(Value *AddedRays,unsigned NbAddedRays,Polyhedron *Pol,unsigned NbMaxConstrs) {
02645
02646 Polyhedron *NewPol = NULL;
02647 Matrix *Mat = NULL, *Ray = NULL;
02648 SatMatrix *Sat = NULL, *SatTranspose = NULL;
02649 unsigned NbCon, NbRay,NbEle1, Dimension;
02650
02651 POL_ENSURE_FACETS(Pol);
02652 POL_ENSURE_VERTICES(Pol);
02653
02654 CATCH(any_exception_error) {
02655 if (NewPol) Polyhedron_Free(NewPol);
02656 if (Mat) Matrix_Free(Mat);
02657 if (Ray) Matrix_Free(Ray);
02658 if (Sat) SMFree(&Sat);
02659 if (SatTranspose) SMFree(&SatTranspose);
02660 RETHROW();
02661 }
02662 TRY {
02663
02664 NbCon = Pol->NbConstraints;
02665 NbRay = Pol->NbRays;
02666 Dimension = Pol->Dimension + 2;
02667 NbEle1 = NbRay * Dimension;
02668
02669 Ray = Matrix_Alloc(NbAddedRays + NbRay, Dimension);
02670 if(!Ray) {
02671 errormsg1("AddRays", "outofmem", "out of memory space");
02672 UNCATCH(any_exception_error);
02673 return 0;
02674 }
02675
02676
02677 if (NbRay)
02678 Vector_Copy(Pol->Ray[0], Ray->p_Init, NbEle1);
02679
02680
02681 Vector_Copy(AddedRays, Ray->p_Init+NbEle1, NbAddedRays * Dimension);
02682
02683
02684 if (POL_ISSET(NbMaxConstrs, POL_NO_DUAL))
02685 NbMaxConstrs = 0;
02686
02687
02688 if (NbMaxConstrs < NbCon)
02689 NbMaxConstrs = NbCon;
02690
02691
02692 Mat = Matrix_Alloc(NbMaxConstrs, Dimension);
02693 if(!Mat) {
02694 errormsg1("AddRays", "outofmem", "out of memory space");
02695 UNCATCH(any_exception_error);
02696 return 0;
02697 }
02698 Mat->NbRows = NbCon;
02699
02700
02701 Vector_Copy(Pol->Constraint[0], Mat->p_Init, NbCon*Dimension);
02702
02703
02704
02705
02706 SatTranspose = BuildSat(Ray, Mat, NbRay, NbMaxConstrs);
02707
02708
02709 Pol_status = Chernikova(Ray, Mat, SatTranspose, Pol->NbEq, NbMaxConstrs, NbRay,1);
02710
02711
02712
02713 Sat = TransformSat(Mat, Ray, SatTranspose);
02714 SMFree(&SatTranspose), SatTranspose = NULL;
02715
02716
02717 NewPol = Remove_Redundants(Mat, Ray, Sat, 0);
02718
02719 SMFree(&Sat), Sat = NULL;
02720 Matrix_Free(Mat), Mat = NULL;
02721 Matrix_Free(Ray), Ray = NULL;
02722 }
02723
02724 UNCATCH(any_exception_error);
02725 return NewPol;
02726 }
02727
02728
02729
02730
02731
02732
02733 Polyhedron *DomainAddRays(Polyhedron *Pol,Matrix *Ray,unsigned NbMaxConstrs) {
02734
02735 Polyhedron *PolA, *PolEndA, *p1, *p2, *p3;
02736 int Redundant;
02737
02738 if (!Pol) return (Polyhedron*) 0;
02739 if (!Ray || Ray->NbRows == 0)
02740 return Domain_Copy(Pol);
02741 if (Pol->Dimension != Ray->NbColumns-2) {
02742 errormsg1(
02743 "DomainAddRays", "diffdim", "operation on different dimensions");
02744 return (Polyhedron*) 0;
02745 }
02746
02747
02748 PolA = PolEndA = (Polyhedron *)0;
02749 for (p1=Pol; p1; p1=p1->next) {
02750 p3 = AddRays(Ray->p[0], Ray->NbRows, p1, NbMaxConstrs);
02751
02752
02753 Redundant = 0;
02754 for (p2=PolA; p2; p2=p2->next) {
02755 if (PolyhedronIncludes(p2, p3)) {
02756 Redundant = 1;
02757 break;
02758 }
02759 }
02760
02761
02762 if (Redundant)
02763 Polyhedron_Free(p3);
02764 else {
02765 if (!PolEndA)
02766 PolEndA = PolA = p3;
02767 else {
02768 PolEndA->next = p3;
02769 PolEndA = PolEndA->next;
02770 }
02771 }
02772 }
02773 return PolA;
02774 }
02775
02776
02777
02778
02779 Polyhedron *Polyhedron_Copy(Polyhedron *Pol) {
02780
02781 Polyhedron *Pol1;
02782
02783 if (!Pol) return (Polyhedron *)0;
02784
02785
02786 Pol1 = Polyhedron_Alloc(Pol->Dimension, Pol->NbConstraints, Pol->NbRays);
02787 if (!Pol1) {
02788 errormsg1("Polyhedron_Copy", "outofmem", "out of memory space");
02789 return 0;
02790 }
02791 if( Pol->NbConstraints )
02792 Vector_Copy(Pol->Constraint[0], Pol1->Constraint[0],
02793 Pol->NbConstraints*(Pol->Dimension+2));
02794 if( Pol->NbRays )
02795 Vector_Copy(Pol->Ray[0], Pol1->Ray[0],
02796 Pol->NbRays*(Pol->Dimension+2));
02797 Pol1->NbBid = Pol->NbBid;
02798 Pol1->NbEq = Pol->NbEq;
02799 Pol1->flags = Pol->flags;
02800 return Pol1;
02801 }
02802
02803
02804
02805
02806 Polyhedron *Domain_Copy(Polyhedron *Pol) {
02807
02808 Polyhedron *Pol1;
02809
02810 if (!Pol) return (Polyhedron *) 0;
02811 Pol1 = Polyhedron_Copy(Pol);
02812 if (Pol->next) Pol1->next = Domain_Copy(Pol->next);
02813 return Pol1;
02814 }
02815
02816
02817
02818
02819
02820
02821
02822
02823
02824
02825
02826
02827
02828
02829
02830
02831
02832
02833
02834
02835
02836
02837
02838
02839
02840
02841
02842
02843
02844 static void addToFilter(int k, unsigned *Filter, SatMatrix *Sat,Value *tmpR,Value *tmpC,int NbRays,int NbConstraints) {
02845
02846 int kj, i,j, jx;
02847 unsigned kb, bx;
02848
02849
02850 kj = k/WSIZE; kb = MSB; kb >>= k%WSIZE;
02851 Filter[kj]|=kb;
02852 value_set_si(tmpC[k],-1);
02853
02854
02855 for(i=0; i<NbRays; i++)
02856 if (value_posz_p(tmpR[i])) {
02857 if (Sat->p[i][kj]&kb)
02858 value_decrement(tmpR[i],tmpR[i]);
02859 else {
02860
02861
02862 value_set_si(tmpR[i],-1);
02863
02864
02865 jx=0; bx=MSB;
02866 for(j=0; j<NbConstraints; j++) {
02867 if (value_posz_p(tmpC[j]) && (Sat->p[i][jx]&bx) )
02868 value_decrement(tmpC[j],tmpC[j]);
02869 NEXT(jx,bx);
02870 }
02871 }
02872 }
02873 }
02874
02875
02876
02877
02878
02879
02880
02881
02882
02883
02884
02885 static void FindSimple(Polyhedron *P1,Polyhedron *P2,unsigned *Filter,unsigned NbMaxRays) {
02886
02887 Matrix *Mat = NULL;
02888 SatMatrix *Sat = NULL;
02889 int i, j, k, jx, found;
02890 Value *p1, *p2, p3;
02891 unsigned Dimension, NbRays, NbConstraints, bx, nc;
02892 Value NbConstraintsLeft, tmp;
02893 Value *tmpC = NULL, *tmpR = NULL;
02894 Polyhedron *Pol = NULL, *Pol2 = NULL;
02895
02896
02897 value_init(p3); value_init(NbConstraintsLeft);
02898 value_init(tmp);
02899
02900 CATCH(any_exception_error) {
02901 if (tmpC) free(tmpC);
02902 if (tmpR) free(tmpR);
02903 if (Mat) Matrix_Free(Mat);
02904 if (Sat) SMFree(&Sat);
02905 if (Pol2 && Pol2!=P2) Polyhedron_Free(Pol2);
02906 if (Pol && Pol!=Pol2 && Pol!=P2) Polyhedron_Free(Pol);
02907
02908
02909 value_clear(p3); value_clear(NbConstraintsLeft);
02910 value_clear(tmp);
02911 RETHROW();
02912 }
02913 TRY {
02914
02915 Dimension = P1->Dimension+2;
02916 Mat = Matrix_Alloc(P1->NbConstraints, Dimension);
02917 if(!Mat) {
02918 errormsg1("FindSimple", "outofmem", "out of memory space");
02919 UNCATCH(any_exception_error);
02920
02921
02922 value_clear(p3); value_clear(NbConstraintsLeft); value_clear(tmp);
02923 return;
02924 }
02925
02926
02927 jx = 0; bx = MSB; Mat->NbRows=0;
02928 value_set_si(NbConstraintsLeft,0);
02929 for (k=0; k<P1->NbConstraints; k++) {
02930 if (Filter[jx]&bx) {
02931 Vector_Copy(P1->Constraint[k], Mat->p[Mat->NbRows], Dimension);
02932 Mat->NbRows++;
02933 }
02934 else
02935 value_increment(NbConstraintsLeft,NbConstraintsLeft);
02936 NEXT(jx,bx);
02937 }
02938 Pol2 = P2;
02939
02940 for (;;) {
02941 if (Mat->NbRows==0)
02942 Pol = Polyhedron_Copy(Pol2);
02943 else {
02944 Pol = AddConstraints(Mat->p_Init, Mat->NbRows, Pol2, NbMaxRays);
02945 if (Pol2 != P2) Polyhedron_Free(Pol2), Pol2 = NULL;
02946 }
02947 if (emptyQ(Pol)) {
02948 Matrix_Free(Mat), Mat = NULL;
02949 Polyhedron_Free(Pol), Pol = NULL;
02950 UNCATCH(any_exception_error);
02951
02952
02953 value_clear(p3); value_clear(NbConstraintsLeft); value_clear(tmp);
02954 return;
02955 }
02956 Mat->NbRows = 0;
02957 Pol2 = Pol;
02958
02959
02960 Dimension = Pol->Dimension+1;
02961 NbRays = Pol->NbRays;
02962 NbConstraints = P1->NbConstraints;
02963 tmpR = (Value *)malloc(NbRays*sizeof(Value));
02964 if(!tmpR) {
02965 errormsg1("FindSimple", "outofmem", "out of memory space");
02966 UNCATCH(any_exception_error);
02967
02968
02969 value_clear(p3); value_clear(NbConstraintsLeft); value_clear(tmp);
02970 return;
02971 }
02972 for(i=0;i<NbRays;i++)
02973 value_init(tmpR[i]);
02974 tmpC = (Value *)malloc(NbConstraints*sizeof(Value));
02975 if(!tmpC) {
02976 errormsg1("FindSimple", "outofmem", "out of memory space");
02977 UNCATCH(any_exception_error);
02978
02979
02980 value_clear(p3); value_clear(NbConstraintsLeft);
02981 for(i=0;i<NbRays;i++)
02982 value_clear(tmpR[i]);
02983 free(tmpR);
02984 return;
02985 }
02986 for(i=0;i<NbConstraints;i++)
02987 value_init(tmpC[i]);
02988 Vector_Set(tmpR,0,NbRays);
02989 Vector_Set(tmpC,0,NbConstraints);
02990
02991
02992 nc = (NbConstraints - 1) / (sizeof(int)*8) + 1;
02993 Sat = SMAlloc(NbRays, nc);
02994 Sat->NbRows = NbRays;
02995 SMVector_Init(Sat->p_init, nc*NbRays);
02996
02997 jx=0; bx=MSB;
02998 for (k=0; k<NbConstraints; k++) {
02999 if (Filter[jx]&bx)
03000 value_set_si(tmpC[k],-1);
03001 else
03002 for (i=0; i<NbRays; i++) {
03003 p1 = Pol->Ray[i]+1;
03004 p2 = P1->Constraint[k]+1;
03005 value_set_si(p3,0);
03006 for (j=0; j<Dimension; j++) {
03007 value_addmul(p3, *p1, *p2);
03008 p1++; p2++;
03009 }
03010 if(value_zero_p(p3) ||
03011 (value_pos_p(p3) && value_notzero_p(P1->Constraint[k][0]))) {
03012 Sat->p[i][jx]|=bx;
03013 value_increment(tmpR[i],tmpR[i]);
03014 value_increment(tmpC[k],tmpC[k]);
03015 }
03016 }
03017 NEXT(jx, bx);
03018 }
03019
03020 do {
03021 found = 0;
03022 for(i=0; i<NbRays; i++)
03023 if(value_posz_p(tmpR[i])) {
03024 value_add_int(tmp,tmpR[i],1);
03025 if(value_eq(tmp,NbConstraintsLeft)) {
03026
03027
03028 jx = 0; bx = MSB;
03029 for(k=0; k<NbConstraints; k++) {
03030 if(value_posz_p(tmpC[k]) && ((Sat->p[i][jx]&bx)==0)) {
03031 addToFilter(k, Filter, Sat, tmpR, tmpC,
03032 NbRays, NbConstraints);
03033 Vector_Copy(P1->Constraint[k],
03034 Mat->p[Mat->NbRows],Dimension+1);
03035 Mat->NbRows++;
03036 value_decrement(NbConstraintsLeft,NbConstraintsLeft);
03037 found=1;
03038 break;
03039 }
03040 NEXT(jx,bx);
03041 }
03042 break;
03043 }
03044 }
03045 }
03046 while (found);
03047
03048 if (!Mat->NbRows) {
03049
03050 Value cmax;
03051 value_init(cmax);
03052
03053 #ifndef LINEAR_VALUE_IS_CHARS
03054 value_set_si(cmax,(NbRays * NbConstraints+1));
03055 #else
03056 value_set_si(cmax,1);
03057 #endif
03058
03059 j = -1;
03060 for(k=0; k<NbConstraints; k++)
03061 if (value_posz_p(tmpC[k])) {
03062 if (value_gt(cmax,tmpC[k])) {
03063 value_assign(cmax,tmpC[k]);
03064 j = k;
03065 }
03066 }
03067 value_clear(cmax);
03068 if (j<0) {
03069 errormsg1("DomSimplify","logerror","logic error");
03070 }
03071 else {
03072 addToFilter(j, Filter, Sat, tmpR, tmpC, NbRays, NbConstraints);
03073 Vector_Copy(P1->Constraint[j],Mat->p[Mat->NbRows],Dimension+1);
03074 Mat->NbRows++;
03075 value_decrement(NbConstraintsLeft,NbConstraintsLeft);
03076 }
03077 }
03078 SMFree(&Sat), Sat = NULL;
03079 free(tmpC), tmpC = NULL;
03080 free(tmpR), tmpR = NULL;
03081 }
03082 }
03083
03084
03085 value_clear(p3); value_clear(NbConstraintsLeft);
03086 value_clear(tmp);
03087 for(i=0;i<NbRays;i++)
03088 value_clear(tmpR[i]);
03089 for(i=0;i<NbRays;i++)
03090 value_clear(tmpC[i]);
03091
03092 UNCATCH(any_exception_error);
03093 }
03094
03095
03096
03097
03098
03099
03100
03101
03102
03103 static int SimplifyConstraints(Polyhedron *Pol1,Polyhedron *Pol2,unsigned *Filter,unsigned NbMaxRays) {
03104
03105 Polyhedron *Pol = NULL;
03106 Matrix *Mat = NULL, *Ray = NULL;
03107 SatMatrix *Sat = NULL;
03108 unsigned NbRay, NbCon, NbCon1, NbCon2, NbEle1, Dimension, notempty;
03109
03110 CATCH(any_exception_error) {
03111 if (Pol) Polyhedron_Free(Pol);
03112 if (Mat) Matrix_Free(Mat);
03113 if (Ray) Matrix_Free(Ray);
03114 if (Sat) SMFree(&Sat);
03115 RETHROW();
03116 }
03117 TRY {
03118
03119 NbRay = Pol1->NbRays;
03120 NbCon1 = Pol1->NbConstraints;
03121 NbCon2 = Pol2->NbConstraints;
03122 NbCon = NbCon1 + NbCon2;
03123 Dimension = Pol1->Dimension+2;
03124 NbEle1 = NbCon1*Dimension;
03125
03126
03127 if (POL_ISSET(NbMaxRays, POL_NO_DUAL))
03128 NbMaxRays = 0;
03129
03130 if (NbRay > NbMaxRays)
03131 NbMaxRays = NbRay;
03132
03133
03134 Mat = Matrix_Alloc(NbCon, Dimension);
03135 if(!Mat) {
03136 errormsg1("SimplifyConstraints", "outofmem", "out of memory space");
03137 UNCATCH(any_exception_error);
03138 return 0;
03139 }
03140
03141
03142 Vector_Copy(Pol1->Constraint[0], Mat->p_Init, NbEle1);
03143
03144
03145 Vector_Copy(Pol2->Constraint[0], Mat->p_Init+NbEle1, NbCon2*Dimension);
03146
03147
03148 Ray = Matrix_Alloc(NbMaxRays, Dimension);
03149 if(!Ray) {
03150 errormsg1("SimplifyConstraints", "outofmem", "out of memory space");
03151 UNCATCH(any_exception_error);
03152 return 0;
03153 }
03154 Ray->NbRows = NbRay;
03155
03156
03157 Vector_Copy(Pol1->Ray[0], Ray->p_Init, NbRay*Dimension);
03158
03159
03160
03161 Sat = BuildSat(Mat, Ray, NbCon1, NbMaxRays);
03162
03163
03164 Pol_status = Chernikova(Mat, Ray, Sat, Pol1->NbBid, NbMaxRays, NbCon1,0);
03165
03166
03167 Pol = Remove_Redundants(Mat, Ray, Sat, Filter);
03168 notempty = 1;
03169 if (Filter && emptyQ(Pol)) {
03170 notempty = 0;
03171 FindSimple(Pol1, Pol2, Filter, NbMaxRays);
03172 }
03173
03174
03175 Polyhedron_Free(Pol), Pol = NULL;
03176 SMFree(&Sat), Sat = NULL;
03177 Matrix_Free(Ray), Ray = NULL;
03178 Matrix_Free(Mat), Mat = NULL;
03179
03180 }
03181
03182 UNCATCH(any_exception_error);
03183 return notempty;
03184 }
03185
03186
03187
03188
03189
03190 static void SimplifyEqualities(Polyhedron *Pol1, Polyhedron *Pol2, unsigned *Filter) {
03191
03192 int i,j;
03193 unsigned ix, bx, NbEqn, NbEqn1, NbEqn2, NbEle2, Dimension;
03194 Matrix *Mat;
03195
03196 NbEqn1 = Pol1->NbEq;
03197 NbEqn2 = Pol2->NbEq;
03198 NbEqn = NbEqn1 + NbEqn2;
03199 Dimension = Pol1->Dimension+2;
03200 NbEle2 = NbEqn2*Dimension;
03201
03202 Mat = Matrix_Alloc(NbEqn, Dimension);
03203 if (!Mat) {
03204 errormsg1("SimplifyEqualities", "outofmem", "out of memory space");
03205 Pol_status = 1;
03206 return;
03207 }
03208
03209
03210 Vector_Copy(Pol2->Constraint[0], Mat->p_Init, NbEle2);
03211
03212
03213 Vector_Copy(Pol1->Constraint[0], Mat->p_Init+NbEle2, NbEqn1*Dimension);
03214
03215 Gauss(Mat, NbEqn2, Dimension-1);
03216
03217 ix = 0;
03218 bx = MSB;
03219 for (i=NbEqn2; i<NbEqn; i++) {
03220 for (j=1; j<Dimension; j++) {
03221 if (value_notzero_p(Mat->p[i][j])) {
03222
03223
03224
03225 Filter[ix] |= bx;
03226 break;
03227 }
03228 }
03229 NEXT(ix,bx);
03230 }
03231 Matrix_Free(Mat);
03232 return;
03233 }
03234
03235
03236
03237
03238
03239
03240
03241
03242
03243 Polyhedron *DomainSimplify(Polyhedron *Pol1, Polyhedron *Pol2, unsigned NbMaxRays) {
03244
03245 Polyhedron *p1, *p2, *p3, *d;
03246 unsigned k, jx, bx, nbentries, NbConstraints, Dimension, NbCon, empty;
03247 unsigned *Filter;
03248 Matrix *Constraints;
03249
03250
03251 if (!Pol1 || !Pol2) return Pol1;
03252 if (Pol1->Dimension != Pol2->Dimension) {
03253 errormsg1("DomSimplify","diffdim","operation on different dimensions");
03254 Pol_status = 1;
03255 return 0;
03256 }
03257 POL_ENSURE_VERTICES(Pol1);
03258 POL_ENSURE_VERTICES(Pol2);
03259 if (emptyQ(Pol1)||emptyQ(Pol2))
03260 return Empty_Polyhedron(Pol1->Dimension);
03261
03262
03263
03264 NbCon = 0;
03265 for (p2=Pol2; p2; p2=p2->next)
03266 if (p2->NbConstraints > NbCon)
03267 NbCon = p2->NbConstraints;
03268
03269 Dimension = Pol1->Dimension+2;
03270 d = (Polyhedron *)0;
03271 for (p1=Pol1; p1; p1=p1->next) {
03272
03273
03274
03275
03276
03277
03278 NbConstraints = p1->NbConstraints;
03279 nbentries = (NbConstraints + NbCon - 1) / (sizeof(int)*8) + 1;
03280
03281
03282 Filter = (unsigned *)malloc(nbentries * sizeof(int));
03283 if (!Filter) {
03284 errormsg1("DomSimplify", "outofmem", "out of memory space\n");
03285 Pol_status = 1;
03286 return 0;
03287 }
03288
03289
03290 SMVector_Init(Filter, nbentries);
03291
03292
03293 empty = 1;
03294 for (p2=Pol2; p2; p2=p2->next) {
03295
03296
03297
03298
03299
03300
03301 SimplifyEqualities(p1, p2, Filter);
03302 if (SimplifyConstraints(p1, p2, Filter, NbMaxRays))
03303 empty=0;
03304
03305
03306 }
03307
03308 if (!empty) {
03309
03310
03311 Constraints = Matrix_Alloc(NbConstraints, Dimension);
03312 if (!Constraints) {
03313 errormsg1("DomSimplify", "outofmem", "out of memory space\n");
03314 Pol_status = 1;
03315 return 0;
03316 }
03317 Constraints->NbRows = 0;
03318 for (k=0, jx=0, bx=MSB; k<NbConstraints; k++) {
03319
03320
03321
03322 if (Filter[jx]&bx) {
03323 Vector_Copy(p1->Constraint[k],
03324 Constraints->p[Constraints->NbRows],
03325 Dimension);
03326 Constraints->NbRows++;
03327 }
03328 NEXT(jx,bx);
03329 }
03330
03331
03332
03333 p3 = Constraints2Polyhedron(Constraints,NbMaxRays);
03334 Matrix_Free(Constraints);
03335
03336
03337 d = AddPolyToDomain (p3, d);
03338 p3 = NULL;
03339 }
03340 free(Filter);
03341 }
03342 if (!d)
03343 return Empty_Polyhedron(Pol1->Dimension);
03344 else return d;
03345
03346 }
03347
03348
03349
03350
03351 Polyhedron *Stras_DomainSimplify(Polyhedron *Pol1,Polyhedron *Pol2,unsigned NbMaxRays) {
03352
03353 Polyhedron *p1, *p2, *p3 = NULL, *d = NULL;
03354 unsigned k, jx, bx, nbentries, NbConstraints, Dimension, NbCon, empty;
03355 unsigned *Filter = NULL;
03356 Matrix *Constraints = NULL;
03357
03358 CATCH(any_exception_error) {
03359 if (Constraints) Matrix_Free(Constraints);
03360 if (Filter) free(Filter);
03361 if (d) Polyhedron_Free(d);
03362 if (p2) Polyhedron_Free(p3);
03363 RETHROW();
03364 }
03365 TRY {
03366 if (!Pol1 || !Pol2) {
03367 UNCATCH(any_exception_error);
03368 return Pol1;
03369 }
03370 if (Pol1->Dimension != Pol2->Dimension) {
03371 errormsg1("DomainSimplify","diffdim","operation on different dimensions");
03372 UNCATCH(any_exception_error);
03373 return 0;
03374 }
03375 POL_ENSURE_VERTICES(Pol1);
03376 POL_ENSURE_VERTICES(Pol2);
03377 if (emptyQ(Pol1)||emptyQ(Pol2)) {
03378 UNCATCH(any_exception_error);
03379 return Empty_Polyhedron(Pol1->Dimension);
03380 }
03381
03382
03383
03384 NbCon = 0;
03385 for (p2=Pol2; p2; p2=p2->next)
03386 if (p2->NbConstraints > NbCon)
03387 NbCon = p2->NbConstraints;
03388
03389 Dimension = Pol1->Dimension+2;
03390 d = (Polyhedron *)0;
03391 for (p1=Pol1; p1; p1=p1->next) {
03392
03393
03394
03395
03396
03397
03398 NbConstraints = p1->NbConstraints;
03399 nbentries = (NbConstraints + NbCon - 1)/(sizeof(int)*8) + 1;
03400
03401
03402 Filter = (unsigned *)malloc(nbentries * sizeof(int));
03403 if(!Filter) {
03404 errormsg1("DomainSimplify", "outofmem", "out of memory space");
03405 UNCATCH(any_exception_error);
03406 return 0;
03407 }
03408
03409
03410 SMVector_Init(Filter, nbentries);
03411
03412
03413 empty = 1;
03414 for (p2=Pol2; p2; p2=p2->next) {
03415
03416
03417
03418
03419
03420
03421 if (SimplifyConstraints(p1, p2, Filter, NbMaxRays))
03422 empty=0;
03423 }
03424
03425 if (!empty) {
03426
03427
03428 Constraints = Matrix_Alloc(NbConstraints,Dimension);
03429 if(!Constraints) {
03430 errormsg1("DomainSimplify", "outofmem", "out of memory space");
03431 UNCATCH(any_exception_error);
03432 return 0;
03433 }
03434 Constraints->NbRows = 0;
03435 for (k=0, jx=0, bx=MSB; k<NbConstraints; k++) {
03436
03437
03438
03439 if (Filter[jx]&bx) {
03440 Vector_Copy(p1->Constraint[k],
03441 Constraints->p[Constraints->NbRows],
03442 Dimension);
03443 Constraints->NbRows++;
03444 }
03445 NEXT(jx,bx);
03446 }
03447
03448
03449
03450 p3 = Constraints2Polyhedron(Constraints,NbMaxRays);
03451 Matrix_Free(Constraints), Constraints = NULL;
03452
03453
03454 d = AddPolyToDomain (p3, d);
03455 p3 = NULL;
03456 }
03457 free(Filter), Filter = NULL;
03458 }
03459 }
03460
03461 UNCATCH(any_exception_error);
03462 if (!d)
03463 return Empty_Polyhedron(Pol1->Dimension);
03464 else
03465 return d;
03466 }
03467
03468
03469
03470
03471
03472 Polyhedron *DomainUnion(Polyhedron *Pol1,Polyhedron *Pol2,unsigned NbMaxRays) {
03473
03474 Polyhedron *PolA, *PolEndA, *PolB, *PolEndB, *p1, *p2;
03475 int Redundant;
03476
03477 if (!Pol1 || !Pol2) return (Polyhedron *) 0;
03478 if (Pol1->Dimension != Pol2->Dimension) {
03479 errormsg1("DomainUnion","diffdim","operation on different dimensions");
03480 return (Polyhedron*) 0;
03481 }
03482
03483
03484
03485
03486
03487
03488
03489 PolA = PolEndA = (Polyhedron *)0;
03490 for (p1=Pol1; p1; p1=p1->next) {
03491
03492
03493 Redundant = 0;
03494 for (p2=Pol2; p2; p2=p2->next) {
03495 if (PolyhedronIncludes(p2, p1) ) {
03496 Redundant = 1;
03497
03498
03499 break;
03500
03501 }
03502 }
03503 if (!Redundant) {
03504
03505
03506 if (!PolEndA)
03507 PolEndA = PolA = Polyhedron_Copy(p1);
03508 else {
03509 PolEndA->next = Polyhedron_Copy(p1);
03510 PolEndA = PolEndA->next;
03511 }
03512
03513 }
03514 }
03515
03516
03517 PolB = PolEndB = (Polyhedron *)0;
03518 for (p2=Pol2; p2; p2=p2->next) {
03519
03520
03521 Redundant = 0;
03522 for (p1=PolA; p1; p1=p1->next) {
03523 if (PolyhedronIncludes(p1, p2)) {
03524 Redundant = 1;
03525 break;
03526 }
03527 }
03528 if (!Redundant) {
03529
03530
03531 if (!PolEndB)
03532 PolEndB = PolB = Polyhedron_Copy(p2);
03533 else {
03534 PolEndB->next = Polyhedron_Copy(p2);
03535 PolEndB = PolEndB->next;
03536 }
03537
03538
03539 }
03540 }
03541
03542 if (!PolA) return PolB;
03543 PolEndA->next = PolB;
03544 return PolA;
03545 }
03546
03547
03548
03549
03550
03551
03552
03553 Polyhedron *DomainConvex(Polyhedron *Pol,unsigned NbMaxConstrs) {
03554
03555 Polyhedron *p, *q, *NewPol = NULL;
03556
03557 CATCH(any_exception_error) {
03558 if (NewPol) Polyhedron_Free(NewPol);
03559 RETHROW();
03560 }
03561 TRY {
03562
03563 if (!Pol) {
03564 UNCATCH(any_exception_error);
03565 return (Polyhedron*) 0;
03566 }
03567
03568 NewPol = Polyhedron_Copy(Pol);
03569 for (p=Pol->next; p; p=p->next) {
03570 q = AddRays(p->Ray[0], p->NbRays, NewPol, NbMaxConstrs);
03571 Polyhedron_Free(NewPol);
03572 NewPol = q;
03573 }
03574 }
03575
03576 UNCATCH(any_exception_error);
03577
03578 return NewPol;
03579 }
03580
03581
03582
03583
03584
03585 Polyhedron *DomainDifference(Polyhedron *Pol1,Polyhedron *Pol2,unsigned NbMaxRays) {
03586
03587 Polyhedron *p1, *p2, *p3, *d;
03588 int i;
03589
03590 if (!Pol1 || !Pol2) return (Polyhedron*) 0;
03591 if (Pol1->Dimension != Pol2->Dimension) {
03592 errormsg1("DomainDifference",
03593 "diffdim", "operation on different dimensions");
03594 return (Polyhedron*) 0;
03595 }
03596 POL_ENSURE_FACETS(Pol1);
03597 POL_ENSURE_VERTICES(Pol1);
03598 POL_ENSURE_FACETS(Pol2);
03599 POL_ENSURE_VERTICES(Pol2);
03600 if (emptyQ(Pol1) || emptyQ(Pol2))
03601 return (Domain_Copy(Pol1));
03602 d = (Polyhedron *)0;
03603 for (p2=Pol2; p2; p2=p2->next) {
03604 for (p1=Pol1; p1; p1=p1->next) {
03605 for (i=0; i<p2->NbConstraints; i++) {
03606
03607
03608
03609 p3 = SubConstraint(p2->Constraint[i], p1, NbMaxRays,0);
03610
03611
03612 d = AddPolyToDomain (p3, d);
03613
03614
03615
03616
03617
03618 if( value_notzero_p(p2->Constraint[i][0]) )
03619 continue;
03620 p3 = SubConstraint(p2->Constraint[i], p1, NbMaxRays,1);
03621
03622
03623 d = AddPolyToDomain (p3, d);
03624 }
03625 }
03626 if (p2 != Pol2)
03627 Domain_Free(Pol1);
03628 Pol1 = d;
03629 d = (Polyhedron *)0;
03630 }
03631 if (!Pol1)
03632 return Empty_Polyhedron(Pol2->Dimension);
03633 else
03634 return Pol1;
03635 }
03636
03637
03638
03639
03640
03641
03642 Polyhedron *align_context(Polyhedron *Pol,int align_dimension,int NbMaxRays) {
03643
03644 int i, j, k;
03645 Polyhedron *p = NULL, **next, *result = NULL;
03646 unsigned dim;
03647
03648 CATCH(any_exception_error) {
03649 if (result) Polyhedron_Free(result);
03650 RETHROW();
03651 }
03652 TRY {
03653
03654 if (!Pol) return Pol;
03655 dim = Pol->Dimension;
03656 if (align_dimension < Pol->Dimension) {
03657 errormsg1("align_context", "diffdim", "context dimension exceeds data");
03658 UNCATCH(any_exception_error);
03659 return NULL;
03660 }
03661 if (align_dimension == Pol->Dimension) {
03662 UNCATCH(any_exception_error);
03663 return Domain_Copy(Pol);
03664 }
03665
03666
03667 k = align_dimension - Pol->Dimension;
03668 next = &result;
03669
03670
03671 for (; Pol; Pol=Pol->next) {
03672 int have_cons = !F_ISSET(Pol, POL_VALID) || F_ISSET(Pol, POL_INEQUALITIES);
03673 int have_rays = !F_ISSET(Pol, POL_VALID) || F_ISSET(Pol, POL_POINTS);
03674 unsigned NbCons = have_cons ? Pol->NbConstraints : 0;
03675 unsigned NbRays = have_rays ? Pol->NbRays + k : 0;
03676
03677 if (Pol->Dimension != dim) {
03678 Domain_Free(result);
03679 errormsg1("align_context", "diffdim", "context not of uniform dimension");
03680 UNCATCH(any_exception_error);
03681 return NULL;
03682 }
03683
03684 p = Polyhedron_Alloc(align_dimension, NbCons, NbRays);
03685 if (have_cons) {
03686 for (i = 0; i < NbCons; ++i) {
03687 value_assign(p->Constraint[i][0], Pol->Constraint[i][0]);
03688 Vector_Copy(Pol->Constraint[i]+1, p->Constraint[i]+k+1, Pol->Dimension+1);
03689 }
03690 p->NbEq = Pol->NbEq;
03691 }
03692
03693 if (have_rays) {
03694 for (i = 0; i < k; ++i)
03695 value_set_si(p->Ray[i][1+i], 1);
03696 for (i = 0; i < Pol->NbRays; ++i) {
03697 value_assign(p->Ray[k+i][0], Pol->Ray[i][0]);
03698 Vector_Copy(Pol->Ray[i]+1, p->Ray[i+k]+k+1, Pol->Dimension+1);
03699 }
03700 p->NbBid = Pol->NbBid + k;
03701 }
03702 p->flags = Pol->flags;
03703
03704 *next = p;
03705 next = &p->next;
03706 }
03707 }
03708
03709 UNCATCH(any_exception_error);
03710 return result;
03711 }
03712
03713
03714
03715
03716
03717
03718
03719
03720 Polyhedron *Polyhedron_Scan(Polyhedron *D, Polyhedron *C,unsigned NbMaxRays) {
03721
03722 int i, j, dim ;
03723 Matrix *Mat;
03724 Polyhedron *C1, *C2, *D1, *D2;
03725 Polyhedron *res, *last, *tmp;
03726
03727 dim = D->Dimension - C->Dimension;
03728 res = last = (Polyhedron *) 0;
03729 if (dim==0) return (Polyhedron *)0;
03730
03731 assert(!D->next);
03732
03733 POL_ENSURE_FACETS(D);
03734 POL_ENSURE_VERTICES(D);
03735 POL_ENSURE_FACETS(C);
03736 POL_ENSURE_VERTICES(C);
03737
03738
03739 Mat = Matrix_Alloc(D->Dimension, D->Dimension+2);
03740 if(!Mat) {
03741 errormsg1("Polyhedron_Scan", "outofmem", "out of memory space");
03742 return 0;
03743 }
03744 C1 = align_context(C,D->Dimension,NbMaxRays);
03745 if(!C1) {
03746 return 0;
03747 }
03748
03749 D2 = DomainIntersection( C1, D, NbMaxRays);
03750
03751 for (i=0; i<dim; i++)
03752 {
03753 Vector_Set(Mat->p_Init,0,D2->Dimension*(D2->Dimension + 2));
03754 for (j=i+1; j<dim; j++) {
03755 value_set_si(Mat->p[j-i-1][j+1],1);
03756 }
03757 Mat->NbRows = dim-i-1;
03758 D1 = Mat->NbRows ? DomainAddRays(D2, Mat, NbMaxRays) : D2;
03759 tmp = DomainSimplify(D1, C1, NbMaxRays);
03760 if (!last) res = last = tmp;
03761 else { last->next = tmp; last = tmp; }
03762 C2 = DomainIntersection(C1, D1, NbMaxRays);
03763 Domain_Free(C1);
03764 C1 = C2;
03765 if (Mat->NbRows) Domain_Free(D1);
03766 }
03767 Domain_Free(D2);
03768 Domain_Free(C1);
03769 Matrix_Free(Mat);
03770 return res;
03771 }
03772
03773
03774
03775
03776
03777
03778
03779
03780
03781
03782 int lower_upper_bounds(int pos,Polyhedron *P,Value *context,Value *LBp,Value *UBp) {
03783
03784 Value LB, UB;
03785 int flag, i;
03786 Value n, n1, d, tmp;
03787
03788 POL_ENSURE_FACETS(P);
03789 POL_ENSURE_VERTICES(P);
03790
03791
03792 value_init(LB); value_init(UB); value_init(tmp);
03793 value_init(n); value_init(n1); value_init(d);
03794
03795 value_set_si(LB,0);
03796 value_set_si(UB,0);
03797
03798
03799 flag = LB_INFINITY | UB_INFINITY;
03800 for (i=0; i<P->NbConstraints; i++) {
03801 value_assign(d,P->Constraint[i][pos]);
03802 Inner_Product(&context[1],&(P->Constraint[i][1]),P->Dimension+1,&n);
03803 if (value_zero_p(d)) {
03804
03805 if (value_neg_p(n) ||
03806 (value_zero_p(P->Constraint[i][0]) && value_notzero_p(n)))
03807 goto empty_loop;
03808 continue;
03809 }
03810 value_oppose(n,n);
03811
03812
03813
03814
03815
03816
03817
03818
03819
03820
03821
03822
03823
03824 if(value_zero_p(P->Constraint[i][0])) {
03825 value_modulus(tmp,n,d);
03826
03827
03828 if (value_notzero_p(tmp))
03829 goto empty_loop;
03830
03831 value_division(n1,n,d);
03832
03833
03834 if((flag&LB_INFINITY) || value_gt(n1,LB))
03835 value_assign(LB,n1);
03836 if((flag&UB_INFINITY) || value_lt(n1,UB))
03837 value_assign(UB,n1);
03838 flag = 0;
03839 }
03840
03841 if (value_pos_p(d)) {
03842 value_modulus(tmp,n,d);
03843
03844
03845 if (value_pos_p(n) && value_notzero_p(tmp)) {
03846 value_division(n1,n,d);
03847 value_add_int(n1,n1,1);
03848 }
03849 else
03850 value_division(n1,n,d);
03851 if (flag&LB_INFINITY) {
03852 value_assign(LB,n1);
03853 flag^=LB_INFINITY;
03854 }
03855 else if(value_gt(n1,LB))
03856 value_assign(LB,n1);
03857 }
03858
03859 if (value_neg_p(d)) {
03860 value_modulus(tmp,n,d);
03861
03862
03863 if (value_pos_p(n) && value_notzero_p(tmp)) {
03864 value_division(n1,n,d);
03865 value_sub_int(n1,n1,1);
03866 }
03867 else
03868 value_division(n1,n,d);
03869
03870 if (flag&UB_INFINITY) {
03871 value_assign(UB,n1);
03872 flag^=UB_INFINITY;
03873 }
03874 else if (value_lt(n1,UB))
03875 value_assign(UB, n1);
03876 }
03877 }
03878 if ((flag & LB_INFINITY)==0) value_assign(*LBp,LB);
03879 if ((flag & UB_INFINITY)==0) value_assign(*UBp,UB);
03880
03881 if (0) {
03882 empty_loop:
03883 flag = 0;
03884 value_set_si(*LBp, 1);
03885 value_set_si(*UBp, 0);
03886 }
03887
03888
03889 value_clear(LB); value_clear(UB); value_clear(tmp);
03890 value_clear(n); value_clear(n1); value_clear(d);
03891 return flag;
03892 }
03893
03894
03895
03896
03897 static void Rays_Mult(Value **A, Matrix *B, Value **C, unsigned NbRays)
03898 {
03899 int i, j, k;
03900 unsigned Dimension1, Dimension2;
03901 Value Sum, tmp;
03902
03903 value_init(Sum); value_init(tmp);
03904
03905 CATCH(any_exception_error) {
03906 value_clear(Sum); value_clear(tmp);
03907 RETHROW();
03908 }
03909 TRY {
03910 Dimension1 = B->NbRows;
03911 Dimension2 = B->NbColumns;
03912
03913 for (i=0; i<NbRays; i++) {
03914 value_assign(C[i][0],A[i][0]);
03915 for (j=0; j<Dimension2; j++) {
03916 value_set_si(Sum,0);
03917 for (k=0; k<Dimension1; k++) {
03918
03919
03920 value_addmul(Sum, A[i][k+1], B->p[k][j]);
03921 }
03922 value_assign(C[i][j+1],Sum);
03923 }
03924 Vector_Gcd(C[i]+1, Dimension2, &tmp);
03925 if (value_notone_p(tmp))
03926 Vector_AntiScale(C[i]+1, C[i]+1, tmp, Dimension2);
03927 }
03928 }
03929 UNCATCH(any_exception_error);
03930 value_clear(Sum); value_clear(tmp);
03931 }
03932
03933
03934
03935
03936 static void Rays_Mult_Transpose(Value **A, Matrix *B, Value **C,
03937 unsigned NbRays)
03938 {
03939 int i, j, k;
03940 unsigned Dimension1, Dimension2;
03941 Value Sum, tmp;
03942
03943 value_init(Sum); value_init(tmp);
03944
03945 CATCH(any_exception_error) {
03946 value_clear(Sum); value_clear(tmp);
03947 RETHROW();
03948 }
03949 TRY {
03950 Dimension1 = B->NbColumns;
03951 Dimension2 = B->NbRows;
03952
03953 for (i=0; i<NbRays; i++) {
03954 value_assign(C[i][0],A[i][0]);
03955 for (j=0; j<Dimension2; j++) {
03956 value_set_si(Sum,0);
03957 for (k=0; k<Dimension1; k++) {
03958
03959
03960 value_addmul(Sum, A[i][k+1], B->p[j][k]);
03961 }
03962 value_assign(C[i][j+1],Sum);
03963 }
03964 Vector_Gcd(C[i]+1, Dimension2, &tmp);
03965 if (value_notone_p(tmp))
03966 Vector_AntiScale(C[i]+1, C[i]+1, tmp, Dimension2);
03967 }
03968 }
03969 UNCATCH(any_exception_error);
03970 value_clear(Sum); value_clear(tmp);
03971 }
03972
03973
03974
03975
03976
03977
03978
03979 Polyhedron *Polyhedron_Preimage(Polyhedron *Pol,Matrix *Func,unsigned NbMaxRays) {
03980
03981 Matrix *Constraints = NULL;
03982 Polyhedron *NewPol = NULL;
03983 unsigned NbConstraints, Dimension1, Dimension2;
03984
03985 POL_ENSURE_INEQUALITIES(Pol);
03986
03987 CATCH(any_exception_error) {
03988 if (Constraints) Matrix_Free(Constraints);
03989 if (NewPol) Polyhedron_Free(NewPol);
03990 RETHROW();
03991 }
03992 TRY {
03993
03994 NbConstraints = Pol->NbConstraints;
03995 Dimension1 = Pol->Dimension+1;
03996 Dimension2 = Func->NbColumns;
03997 if (Dimension1!=(Func->NbRows)) {
03998 errormsg1("Polyhedron_Preimage", "dimincomp", "incompatable dimensions");
03999 UNCATCH(any_exception_error);
04000 return Empty_Polyhedron(Dimension2-1);
04001 }
04002
04003
04004
04005
04006
04007
04008
04009
04010
04011 Constraints = Matrix_Alloc(NbConstraints, Dimension2+1);
04012 if (!Constraints) {
04013 errormsg1("Polyhedron_Preimage", "outofmem", "out of memory space\n");
04014 Pol_status = 1;
04015 UNCATCH(any_exception_error);
04016 return 0;
04017 }
04018
04019
04020
04021 Rays_Mult(Pol->Constraint, Func, Constraints->p, NbConstraints);
04022 NewPol = Constraints2Polyhedron(Constraints, NbMaxRays);
04023 Matrix_Free(Constraints), Constraints = NULL;
04024
04025 }
04026
04027 UNCATCH(any_exception_error);
04028
04029 return NewPol;
04030 }
04031
04032
04033
04034
04035
04036
04037
04038 Polyhedron *DomainPreimage(Polyhedron *Pol,Matrix *Func,unsigned NbMaxRays) {
04039
04040 Polyhedron *p, *q, *d = NULL;
04041
04042 CATCH(any_exception_error) {
04043 if (d) Polyhedron_Free(d);
04044 RETHROW();
04045 }
04046 TRY {
04047 if (!Pol || !Func) {
04048 UNCATCH(any_exception_error);
04049 return (Polyhedron *) 0;
04050 }
04051 d = (Polyhedron *) 0;
04052 for (p=Pol; p; p=p->next) {
04053 q = Polyhedron_Preimage(p, Func, NbMaxRays);
04054 d = AddPolyToDomain (q, d);
04055 }
04056 }
04057 UNCATCH(any_exception_error);
04058 return d;
04059 }
04060
04061
04062
04063
04064
04065
04066 Polyhedron *Polyhedron_Image(Polyhedron *Pol, Matrix *Func,unsigned NbMaxConstrs) {
04067
04068 Matrix *Rays = NULL;
04069 Polyhedron *NewPol = NULL;
04070 unsigned NbRays, Dimension1, Dimension2;
04071
04072 POL_ENSURE_FACETS(Pol);
04073 POL_ENSURE_VERTICES(Pol);
04074
04075 CATCH(any_exception_error) {
04076 if (Rays) Matrix_Free(Rays);
04077 if (NewPol) Polyhedron_Free(NewPol);
04078 RETHROW();
04079 }
04080 TRY {
04081
04082 NbRays = Pol->NbRays;
04083 Dimension1 = Pol->Dimension+1;
04084 Dimension2 = Func->NbRows;
04085 if (Dimension1!=Func->NbColumns) {
04086 errormsg1("Polyhedron_Image", "dimincomp", "incompatible dimensions");
04087 UNCATCH(any_exception_error);
04088 return Empty_Polyhedron(Dimension2-1);
04089 }
04090
04091
04092
04093
04094
04095
04096
04097
04098
04099
04100 if (Dimension1 == Dimension2) {
04101 Matrix *M, *M2;
04102 int ok;
04103 M = Matrix_Copy(Func);
04104 M2 = Matrix_Alloc(Dimension2, Dimension1);
04105 if (!M2) {
04106 errormsg1("Polyhedron_Image", "outofmem", "out of memory space\n");
04107 UNCATCH(any_exception_error);
04108 return 0;
04109 }
04110
04111 ok = Matrix_Inverse(M, M2);
04112 Matrix_Free(M);
04113 if (ok) {
04114 NewPol = Polyhedron_Alloc(Pol->Dimension, Pol->NbConstraints,
04115 Pol->NbRays);
04116 if (!NewPol) {
04117 errormsg1("Polyhedron_Image", "outofmem",
04118 "out of memory space\n");
04119 UNCATCH(any_exception_error);
04120 return 0;
04121 }
04122 Rays_Mult_Transpose(Pol->Ray, Func, NewPol->Ray, NbRays);
04123 Rays_Mult(Pol->Constraint, M2, NewPol->Constraint,
04124 Pol->NbConstraints);
04125 NewPol->NbEq = Pol->NbEq;
04126 NewPol->NbBid = Pol->NbBid;
04127 if (NewPol->NbEq)
04128 Gauss4(NewPol->Constraint, NewPol->NbEq, NewPol->NbConstraints,
04129 NewPol->Dimension+1);
04130 if (NewPol->NbBid)
04131 Gauss4(NewPol->Ray, NewPol->NbBid, NewPol->NbRays,
04132 NewPol->Dimension+1);
04133 }
04134 Matrix_Free(M2);
04135 }
04136
04137 if (!NewPol) {
04138
04139 Rays = Matrix_Alloc(NbRays, Dimension2+1);
04140 if (!Rays) {
04141 errormsg1("Polyhedron_Image", "outofmem", "out of memory space\n");
04142 UNCATCH(any_exception_error);
04143 return 0;
04144 }
04145
04146
04147
04148 Rays_Mult_Transpose(Pol->Ray, Func, Rays->p, NbRays);
04149 NewPol = Rays2Polyhedron(Rays, NbMaxConstrs);
04150 Matrix_Free(Rays), Rays = NULL;
04151 }
04152
04153 }
04154
04155 UNCATCH(any_exception_error);
04156 return NewPol;
04157 }
04158
04159
04160
04161
04162
04163
04164 Polyhedron *DomainImage(Polyhedron *Pol,Matrix *Func,unsigned NbMaxConstrs) {
04165
04166 Polyhedron *p, *q, *d = NULL;
04167
04168 CATCH(any_exception_error) {
04169 if (d) Polyhedron_Free(d);
04170 RETHROW();
04171 }
04172 TRY {
04173
04174 if (!Pol || !Func) {
04175 UNCATCH(any_exception_error);
04176 return (Polyhedron *) 0;
04177 }
04178 d = (Polyhedron *) 0;
04179 for (p=Pol; p; p=p->next) {
04180 q = Polyhedron_Image(p, Func, NbMaxConstrs);
04181 d = AddPolyToDomain (q, d);
04182 }
04183 }
04184
04185 UNCATCH(any_exception_error);
04186
04187 return d;
04188 }
04189
04190
04191
04192
04193
04194
04195
04196
04197
04198
04199
04200 Interval *DomainCost(Polyhedron *Pol,Value *Cost) {
04201
04202 int i, j, NbRay, Dim;
04203 Value *p1, *p2, p3, d, status;
04204 Value tmp1, tmp2, tmp3;
04205 Value **Ray;
04206 Interval *I = NULL;
04207
04208 value_init(p3); value_init(d); value_init(status);
04209 value_init(tmp1); value_init(tmp2); value_init(tmp3);
04210
04211 POL_ENSURE_FACETS(Pol);
04212 POL_ENSURE_VERTICES(Pol);
04213
04214 CATCH(any_exception_error) {
04215 if (I) free(I);
04216 RETHROW();
04217 value_clear(p3); value_clear(d); value_clear(status);
04218 value_clear(tmp1); value_clear(tmp2); value_clear(tmp3);
04219 }
04220 TRY {
04221
04222 Ray = Pol->Ray;
04223 NbRay = Pol->NbRays;
04224 Dim = Pol->Dimension+1;
04225 I = (Interval *) malloc (sizeof(Interval));
04226 if (!I) {
04227 errormsg1("DomainCost", "outofmem", "out of memory space\n");
04228 UNCATCH(any_exception_error);
04229 value_clear(p3); value_clear(d); value_clear(status);
04230 value_clear(tmp1); value_clear(tmp2); value_clear(tmp3);
04231 return 0;
04232 }
04233
04234
04235
04236
04237
04238
04239
04240
04241 value_set_si(I->MaxN,-1);
04242 value_set_si(I->MaxD,0);
04243 I->MaxI = -1;
04244 value_set_si(I->MinN,1);
04245 value_set_si(I->MinD,0);
04246 I->MinI = -1;
04247
04248
04249 for (i=0; i<NbRay; i++) {
04250 p1 = Ray[i];
04251 value_assign(status, *p1);
04252 p1++;
04253 p2 = Cost;
04254
04255
04256 value_multiply(p3,*p1,*p2);
04257 p1++; p2++;
04258 for (j=1; j<Dim; j++) {
04259 value_multiply(tmp1,*p1,*p2);
04260
04261
04262 value_addto(p3,p3,tmp1);
04263 p1++; p2++;
04264 }
04265
04266
04267 p1--;
04268 value_assign(d,*p1);
04269 value_multiply(tmp1,p3,I->MaxD);
04270 value_multiply(tmp2,I->MaxN,d);
04271 value_set_si(tmp3,1);
04272
04273
04274 if (I->MaxI==-1 ||
04275 value_gt(tmp1,tmp2) ||
04276 (value_eq(tmp1,tmp2) &&
04277 value_eq(d,tmp3) && value_ne(I->MaxD,tmp3))) {
04278 value_assign(I->MaxN,p3);
04279 value_assign(I->MaxD,d);
04280 I->MaxI = i;
04281 }
04282 value_multiply(tmp1,p3,I->MinD);
04283 value_multiply(tmp2,I->MinN,d);
04284 value_set_si(tmp3,1);
04285
04286
04287 if (I->MinI==-1 ||
04288 value_lt(tmp1,tmp2) ||
04289 (value_eq(tmp1,tmp2) &&
04290 value_eq(d,tmp3) && value_ne(I->MinD,tmp3))) {
04291 value_assign(I->MinN, p3);
04292 value_assign(I->MinD, d);
04293 I->MinI = i;
04294 }
04295 value_multiply(tmp1,p3,I->MaxD);
04296 value_set_si(tmp2,0);
04297
04298
04299 if (value_zero_p(status)) {
04300 if (value_lt(tmp1,tmp2)) {
04301 value_oppose(I->MaxN,p3);
04302 value_set_si(I->MaxD,0);
04303 I->MaxI = i;
04304 }
04305 value_multiply(tmp1,p3,I->MinD);
04306 value_set_si(tmp2,0);
04307
04308 if (value_gt(tmp1,tmp2)) {
04309 value_oppose(I->MinN,p3);
04310 value_set_si(I->MinD,0);
04311 I->MinI = i;
04312 }
04313 }
04314 }
04315 }
04316
04317 UNCATCH(any_exception_error);
04318 value_clear(p3); value_clear(d); value_clear(status);
04319 value_clear(tmp1); value_clear(tmp2); value_clear(tmp3);
04320 return I;
04321 }
04322
04323
04324
04325
04326
04327
04328 Polyhedron *DomainAddConstraints(Polyhedron *Pol,Matrix *Mat,unsigned NbMaxRays) {
04329
04330 Polyhedron *PolA, *PolEndA, *p1, *p2, *p3;
04331 int Redundant;
04332
04333 if (!Pol) return (Polyhedron*) 0;
04334 if (!Mat) return Pol;
04335 if (Pol->Dimension != Mat->NbColumns-2) {
04336 errormsg1("DomainAddConstraints",
04337 "diffdim", "operation on different dimensions");
04338 return (Polyhedron*) 0;
04339 }
04340
04341
04342 PolA = PolEndA = (Polyhedron *)0;
04343 for (p1=Pol; p1; p1=p1->next) {
04344 p3 = AddConstraints(Mat->p_Init, Mat->NbRows, p1, NbMaxRays);
04345
04346
04347 Redundant = 0;
04348 for (p2=PolA; p2; p2=p2->next) {
04349 if (PolyhedronIncludes(p2, p3)) {
04350 Redundant = 1;
04351 break;
04352 }
04353 }
04354
04355
04356 if (Redundant)
04357 Polyhedron_Free(p3);
04358 else {
04359 if (!PolEndA)
04360 PolEndA = PolA = p3;
04361 else {
04362 PolEndA->next = p3;
04363 PolEndA = PolEndA->next;
04364 }
04365 }
04366 }
04367 return PolA;
04368 }
04369
04370
04371
04372
04373
04374
04375
04376
04377
04378
04379
04380
04381 Polyhedron *Disjoint_Domain( Polyhedron *P, int flag, unsigned NbMaxRays )
04382 {
04383 Polyhedron *lP, *tmp, *Result, *lR, *prec, *reste;
04384 Polyhedron *p1, *p2, *p3, *Pol1, *dx, *d1, *d2, *pi, *newpi;
04385 int i;
04386
04387 if( flag!=0 && flag!=1 )
04388 {
04389 errormsg1("Disjoint_Domain",
04390 "invalidarg", "flag should be equal to 0 or 1");
04391 return (Polyhedron*) 0;
04392 }
04393 if(!P) return (Polyhedron*) 0;
04394 if(!P->next) return Polyhedron_Copy(P);
04395
04396 Result = (Polyhedron *)0;
04397
04398 for(lP=P;lP;lP=lP->next)
04399 {
04400 reste = Polyhedron_Copy(lP);
04401 prec = (Polyhedron *)0;
04402
04403 lR=Result;
04404 while( lR && reste )
04405 {
04406
04407 dx = (Polyhedron *)0;
04408 for( p1=reste; p1; p1=p1->next )
04409 {
04410 p3 = AddConstraints(lR->Constraint[0], lR->NbConstraints, p1,
04411 NbMaxRays);
04412 dx = AddPolyToDomain(p3,dx);
04413 }
04414
04415
04416 if(!dx)
04417 { prec = lR;
04418 lR=lR->next;
04419 continue;
04420 }
04421 if (emptyQ(dx)) {
04422 Domain_Free(dx);
04423 prec = lR;
04424 lR=lR->next;
04425 continue;
04426 }
04427
04428
04429
04430
04431
04432
04433
04434
04435 d1 = (Polyhedron *)0;
04436 for (p1=reste; p1; p1=p1->next)
04437 {
04438 pi = p1;
04439 for (i=0; i<P->NbConstraints && pi ; i++)
04440 {
04441
04442
04443
04444 p3 = SubConstraint(P->Constraint[i], pi, NbMaxRays,2*flag);
04445
04446 d1 = AddPolyToDomain (p3, d1);
04447
04448
04449
04450
04451 if( value_zero_p(P->Constraint[i][0]) )
04452 {
04453 p3 = SubConstraint(P->Constraint[i], pi, NbMaxRays,1+2*flag);
04454
04455 d1 = AddPolyToDomain (p3, d1);
04456
04457
04458 newpi = AddConstraints( P->Constraint[i], 1, pi, NbMaxRays);
04459 }
04460 else
04461 {
04462
04463 newpi = SubConstraint(P->Constraint[i], pi, NbMaxRays,3);
04464 }
04465 if( newpi && emptyQ( newpi ) )
04466 {
04467 Domain_Free( newpi );
04468 newpi = (Polyhedron *)0;
04469 }
04470 if( pi != p1 )
04471 Domain_Free( pi );
04472 pi = newpi;
04473 }
04474 if( pi != p1 )
04475 Domain_Free( pi );
04476 }
04477
04478
04479 Pol1 = Polyhedron_Copy( lR );
04480 for (p2=reste; p2; p2=p2->next)
04481 {
04482 d2 = (Polyhedron *)0;
04483 for (p1=Pol1; p1; p1=p1->next)
04484 {
04485 pi = p1;
04486 for (i=0; i<p2->NbConstraints && pi ; i++)
04487 {
04488
04489
04490
04491 p3 = SubConstraint(p2->Constraint[i], pi, NbMaxRays,2*flag);
04492
04493 d2 = AddPolyToDomain (p3, d2);
04494
04495
04496
04497
04498 if( value_zero_p(p2->Constraint[i][0]) )
04499 {
04500 p3 = SubConstraint(p2->Constraint[i], pi, NbMaxRays,1+2*flag);
04501
04502 d2 = AddPolyToDomain (p3, d2);
04503
04504
04505 newpi = AddConstraints( p2->Constraint[i], 1, pi, NbMaxRays);
04506 }
04507 else
04508 {
04509
04510 newpi = SubConstraint(p2->Constraint[i], pi, NbMaxRays,3);
04511 }
04512 if( newpi && emptyQ( newpi ) )
04513 {
04514 Domain_Free( newpi );
04515 newpi = (Polyhedron *)0;
04516 }
04517 if( pi != p1 )
04518 Domain_Free( pi );
04519 pi = newpi;
04520 }
04521 if( pi && pi!=p1 )
04522 Domain_Free( pi );
04523 }
04524 if( Pol1 )
04525 Domain_Free( Pol1 );
04526 Pol1 = d2;
04527 }
04528
04529
04530
04531 if( d1 && emptyQ(d1) )
04532 {
04533 Domain_Free( d1 );
04534 d1 = NULL;
04535 }
04536 if( d2 && emptyQ(d2) )
04537 {
04538 Domain_Free( d2 );
04539 d2 = NULL;
04540 }
04541
04542
04543 Domain_Free( reste );
04544 reste = d1;
04545
04546
04547 if( d2 )
04548 {
04549 for( tmp=d2 ; tmp->next ; tmp=tmp->next )
04550 ;
04551 tmp->next = Result;
04552 Result = d2;
04553 if( !prec )
04554 prec = tmp;
04555 }
04556
04557
04558 for( tmp=dx ; tmp->next ; tmp=tmp->next )
04559 ;
04560 tmp->next = Result;
04561 Result = dx;
04562 if( !prec )
04563 prec = tmp;
04564
04565
04566 if( !prec )
04567 errormsg1( "Disjoint_Domain","internalerror","internal error");
04568 prec->next = lR->next;
04569 Polyhedron_Free( lR );
04570 lR = prec->next;
04571 }
04572
04573
04574 if(reste)
04575 {
04576 if(emptyQ(reste))
04577 {
04578 Domain_Free( reste );
04579 reste = NULL;
04580 }
04581 else
04582 {
04583 Polyhedron *tnext;
04584 for( tmp=reste ; tmp ; tmp=tnext )
04585 {
04586 tnext = tmp->next;
04587 tmp->next = NULL;
04588 Result = AddPolyToDomain(tmp, Result);
04589 }
04590 }
04591 }
04592 }
04593
04594 return( Result );
04595 }
04596
04597
04598
04599
04600 void Polyhedron_PrintConstraints(FILE *Dst, const char *Format, Polyhedron *Pol)
04601 {
04602 int i,j;
04603
04604 fprintf( Dst, "%d %d\n", Pol->NbConstraints, Pol->Dimension+2 );
04605 for( i=0 ; i<Pol->NbConstraints ; i++ )
04606 {
04607 for( j=0 ; j<Pol->Dimension+2 ; j++ )
04608 value_print( Dst, Format, Pol->Constraint[i][j] );
04609 fprintf( Dst, "\n" );
04610 }
04611
04612 }
04613
04614
04615 void Domain_PrintConstraints(FILE *Dst, const char *Format, Polyhedron *Pol)
04616 {
04617 Polyhedron *Q;
04618 for (Q = Pol; Q; Q = Q->next)
04619 Polyhedron_PrintConstraints(Dst, Format, Q);
04620 }
04621
04622 static Polyhedron *p_simplify_constraints(Polyhedron *P, Vector *row,
04623 Value *g, unsigned MaxRays)
04624 {
04625 Polyhedron *T, *R = P;
04626 int len = P->Dimension+2;
04627 int r;
04628
04629
04630
04631
04632
04633 for (r = 0; r < R->NbConstraints; ++r) {
04634 if (ConstraintSimplify(R->Constraint[r], row->p, len, g)) {
04635 T = R;
04636 if (value_zero_p(R->Constraint[r][0])) {
04637 R = Empty_Polyhedron(R->Dimension);
04638 r = R->NbConstraints;
04639 } else if (POL_ISSET(MaxRays, POL_NO_DUAL)) {
04640 R = Polyhedron_Copy(R);
04641 F_CLR(R, POL_FACETS | POL_VERTICES | POL_POINTS);
04642 Vector_Copy(row->p+1, R->Constraint[r]+1, R->Dimension+1);
04643 } else {
04644 R = AddConstraints(row->p, 1, R, MaxRays);
04645 r = -1;
04646 }
04647 if (T != P)
04648 Polyhedron_Free(T);
04649 }
04650 }
04651 if (R != P)
04652 Polyhedron_Free(P);
04653 return R;
04654 }
04655
04656
04657
04658
04659
04660
04661
04662 Polyhedron *DomainConstraintSimplify(Polyhedron *P, unsigned MaxRays)
04663 {
04664 Polyhedron **prev;
04665 int len = P->Dimension+2;
04666 Vector *row = Vector_Alloc(len);
04667 Value g;
04668 Polyhedron *R = P, *N;
04669 value_set_si(row->p[0], 1);
04670 value_init(g);
04671
04672 for (prev = &R; P; P = N) {
04673 Polyhedron *T;
04674 N = P->next;
04675 T = p_simplify_constraints(P, row, &g, MaxRays);
04676
04677 if (emptyQ(T) && prev != &R) {
04678 Polyhedron_Free(T);
04679 *prev = NULL;
04680 continue;
04681 }
04682
04683 if (T != P)
04684 T->next = N;
04685 *prev = T;
04686 prev = &T->next;
04687 }
04688
04689 if (R->next && emptyQ(R)) {
04690 N = R->next;
04691 Polyhedron_Free(R);
04692 R = N;
04693 }
04694
04695 value_clear(g);
04696 Vector_Free(row);
04697 return R;
04698 }