source/oldpolytest.c

   1 /*
   2     This file is part of PolyLib.
   3
   4     PolyLib is free software: you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation, either version 3 of the License, or
   7     (at your option) any later version.
   8
   9     PolyLib is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with PolyLib.  If not, see <http://www.gnu.org/licenses/>.
  16 */
  17
  18 /* polytest.c */
  19 #include <stdio.h>
  20 #include <polylib/polylib.h>
  21
  22 /* alpha.c
  23      COPYRIGHT
  24           Both this software and its documentation are
  25
  26               Copyright 1993 by IRISA /Universite de Rennes I - France,
  27               Copyright 1995,1996 by BYU, Provo, Utah
  28                          all rights reserved.
  29
  30 */
  31
  32 #include <stdio.h>
  33 #include <stdlib.h>
  34 #include <string.h>
  35
  36 /*---------------------------------------------------------------------*/
  37 /* int exist_points(pos,P,context)                                     */
  38 /*    pos : index position of current loop index (0..hdim-1)           */
  39 /*    P: loop domain                                                   */
  40 /*    context : context values for fixed indices                       */
  41 /* recursive procedure, recurs for each imbriquation                   */
  42 /* returns 1 if there exists any integer points in this polyhedron     */
  43 /* returns 0 if there are none                                         */
  44 /*---------------------------------------------------------------------*/
  45 static int exist_points(int pos,Polyhedron *Pol,Value *context) {
  46
  47   Value LB, UB, k,tmp;
  48
  49   value_init(LB); value_init(UB);
  50   value_init(k);  value_init(tmp);
  51   value_set_si(LB,0);
  52   value_set_si(UB,0);
  53
  54   /* Problem if UB or LB is INFINITY */
  55   if (lower_upper_bounds(pos,Pol,context,&LB,&UB) !=0) {
  56     errormsg1("exist_points", "infdom", "infinite domain");
  57     value_clear(LB);
  58     value_clear(UB);
  59     value_clear(k);
  60     value_clear(tmp);
  61     return -1;
  62   }
  63   value_set_si(context[pos],0);
  64   if(value_lt(UB,LB)) {
  65     value_clear(LB);
  66     value_clear(UB);
  67     value_clear(k);
  68     value_clear(tmp);
  69     return 0;
  70   }
  71   if (!Pol->next) {
  72     value_subtract(tmp,UB,LB);
  73     value_increment(tmp,tmp);
  74     value_clear(UB);
  75     value_clear(LB);
  76     value_clear(k);
  77     return (value_pos_p(tmp));
  78   }
  79
  80   for (value_assign(k,LB);value_le(k,UB);value_increment(k,k)) {
  81
  82     /* insert k in context */
  83     value_assign(context[pos],k);
  84     if (exist_points(pos+1,Pol->next,context) > 0 ) {
  85       value_clear(LB); value_clear(UB);
  86       value_clear(k); value_clear(tmp);
  87       return 1;
  88     }
  89   }
  90   /* Reset context */
  91   value_set_si(context[pos],0);
  92   value_clear(UB); value_clear(LB);
  93   value_clear(k); value_clear(tmp);
  94   return 0;
  95 }
  96
  97 /*--------------------------------------------------------------*/
  98 /* Test to see if there are any integral points in a polyhedron */
  99 /*--------------------------------------------------------------*/
 100 int Polyhedron_Not_Empty(Polyhedron *P,Polyhedron *C,int MAXRAYS) {
 101
 102   int res,i;
 103   Value *context;
 104   Polyhedron *L;
 105
 106   /* Create a context vector size dim+2 and set it to all zeros */
 107   context = (Value *) malloc((P->Dimension+2)*sizeof(Value));
 108
 109   /* Initialize array 'context' */
 110   for (i=0;i<(P->Dimension+2);i++)
 111     value_init(context[i]);
 112
 113   Vector_Set(context,0,(P->Dimension+2));
 114
 115   /* Set context[P->Dimension+1] = 1  (the constant) */
 116   value_set_si(context[P->Dimension+1],1);
 117
 118   L = Polyhedron_Scan(P,C,MAXRAYS);
 119   res = exist_points(1,L,context);
 120   Domain_Free(L);
 121
 122   /* Clear array 'context' */
 123   for (i=0;i<(P->Dimension+2);i++)
 124     value_clear(context[i]);
 125   free(context);
 126   return res;
 127 }
 128
 129 /* PolyhedronLTQ ( P1, P2 ) */
 130 /* P1 and P2 must be simple polyhedra */
 131 /* result =  0 --> not comparable */
 132 /* result = -1 --> P1 < P2        */
 133 /* result =  1 --> P1 > P2        */
 134 int PolyhedronLTQ (Polyhedron *Pol1,Polyhedron *Pol2,int NbMaxConstrs) {
 135
 136   int res, dim, i, j, k;
 137   Polyhedron *Q1, *Q2, *Q3, *Q;
 138   Matrix *Mat;
 139
 140 #define INDEX 1
 141
 142   if (Pol1->next || Pol2->next) {
 143     errormsg1("PolyhedronLTQ", "compoly", "Can only compare polyhedra");
 144     return 0;
 145   }
 146   if (Pol1->Dimension != Pol2->Dimension) {
 147     errormsg1("PolyhedronLTQ","diffdim","Polyhedra are not same dimension");
 148     return 0;
 149   }
 150   dim = Pol1->Dimension+2;
 151
 152 #ifdef DEBUG
 153   fprintf(stdout, "P1\n");
 154   Polyhedron_Print(stdout,P_VALUE_FMT,Pol1);
 155   fprintf(stdout, "P2\n");
 156   Polyhedron_Print(stdout,P_VALUE_FMT,Pol2);
 157 #endif
 158
 159   /* Create the Line to add */
 160   Mat = Matrix_Alloc(1,dim);
 161   Vector_Set(Mat->p_Init,0,dim);
 162   value_set_si(Mat->p[0][INDEX],1);  /* line in INDEX dimension */
 163
 164   Q1 = AddRays(Mat->p[0],1,Pol1,NbMaxConstrs);
 165   Q2 = AddRays(Mat->p[0],1,Pol2,NbMaxConstrs);
 166
 167 #ifdef DEBUG
 168   fprintf(stdout, "Q1\n");
 169   Polyhedron_Print(stdout,P_VALUE_FMT,Q1);
 170   fprintf(stdout, "Q2\n");
 171   Polyhedron_Print(stdout,P_VALUE_FMT,Q2);
 172 #endif
 173
 174   Matrix_Free(Mat);
 175   Q  = DomainIntersection(Q1,Q2,NbMaxConstrs);
 176
 177 #ifdef DEBUG
 178   fprintf(stdout, "Q\n");
 179   Polyhedron_Print(stdout,P_VALUE_FMT,Q);
 180 #endif
 181
 182   Domain_Free(Q1);
 183   Domain_Free(Q2);
 184
 185   if (emptyQ(Q)) res = 0;       /* not comparable */
 186   else {
 187     Q1 = DomainIntersection(Pol1,Q,NbMaxConstrs);
 188     Q2 = DomainIntersection(Pol2,Q,NbMaxConstrs);
 189
 190 #ifdef DEBUG
 191     fprintf(stdout, "Q1\n");
 192     Polyhedron_Print(stdout,P_VALUE_FMT,Q1);
 193     fprintf(stdout, "Q2\n");
 194     Polyhedron_Print(stdout,P_VALUE_FMT,Q2);
 195 #endif
 196
 197     /* Test if Q1 < Q2 */
 198     /* Build a constraint array for >= Q1 and <= Q2 */
 199     Mat = Matrix_Alloc(2,dim);
 200     Vector_Set(Mat->p_Init,0,2*dim);
 201
 202     /* Choose a contraint from Q1 */
 203     for (i=0; i<Q1->NbConstraints; i++) {
 204       if (value_zero_p(Q1->Constraint[i][0])) {
 205
 206         /* Equality */
 207         if (value_zero_p(Q1->Constraint[i][INDEX])) {
 208
 209           /* Ignore side constraint (they are in Q) */
 210           continue;
 211         }
 212         else if (value_neg_p(Q1->Constraint[i][INDEX])) {
 213
 214           /* copy -constraint to Mat */
 215           value_set_si(Mat->p[0][0],1);
 216           for (k=1; k<dim; k++)
 217             value_oppose(Mat->p[0][k],Q1->Constraint[i][k]);
 218         }
 219         else {
 220
 221           /* Copy constraint to Mat */
 222
 223           value_set_si(Mat->p[0][0],1);
 224           for (k=1; k<dim; k++)
 225             value_assign(Mat->p[0][k],Q1->Constraint[i][k]);
 226         }
 227       }
 228       else if(value_neg_p(Q1->Constraint[i][INDEX])) {
 229
 230         /* Upper bound -- make a lower bound from it */
 231         value_set_si(Mat->p[0][0],1);
 232         for (k=1; k<dim; k++)
 233           value_oppose(Mat->p[0][k],Q1->Constraint[i][k]);
 234       }
 235       else {
 236
 237         /* Lower or side bound -- ignore it */
 238         continue;
 239       }
 240
 241       /* Choose a constraint from Q2 */
 242       for (j=0; j<Q2->NbConstraints; j++) {
 243         if (value_zero_p(Q2->Constraint[j][0])) {   /* equality */
 244           if (value_zero_p(Q2->Constraint[j][INDEX])) {
 245
 246             /* Ignore side constraint (they are in Q) */
 247             continue;
 248           }
 249           else if (value_pos_p(Q2->Constraint[j][INDEX])) {
 250
 251             /* Copy -constraint to Mat */
 252             value_set_si(Mat->p[1][0],1);
 253             for (k=1; k<dim; k++)
 254               value_oppose(Mat->p[1][k],Q2->Constraint[j][k]);
 255           }
 256           else {
 257
 258             /* Copy constraint to Mat */
 259             value_set_si(Mat->p[1][0],1);
 260             for (k=1; k<dim; k++)
 261               value_assign(Mat->p[1][k],Q2->Constraint[j][k]);
 262           };
 263         }
 264         else if (value_pos_p(Q2->Constraint[j][INDEX])) {
 265
 266           /* Lower bound -- make an upper bound from it */
 267           value_set_si(Mat->p[1][0],1);
 268           for(k=1;k<dim;k++)
 269             value_oppose(Mat->p[1][k],Q2->Constraint[j][k]);
 270         }
 271         else {
 272
 273           /* Upper or side bound -- ignore it */
 274           continue;
 275         };
 276
 277 #ifdef DEBUG
 278         fprintf(stdout, "i=%d j=%d M=\n", i+1, j+1);
 279         Matrix_Print(stdout,P_VALUE_FMT,Mat);
 280 #endif
 281
 282         /* Add Mat to Q and see if anything is made */
 283         Q3 = AddConstraints(Mat->p[0],2,Q,NbMaxConstrs);
 284
 285 #ifdef DEBUG
 286         fprintf(stdout, "Q3\n");
 287         Polyhedron_Print(stdout,P_VALUE_FMT,Q3);
 288 #endif
 289
 290         if (!emptyQ(Q3)) {
 291           Domain_Free(Q3);
 292
 293 #ifdef DEBUG
 294           fprintf(stdout, "not empty\n");
 295 #endif
 296           res = -1;
 297           goto LTQdone;
 298         }
 299 #ifdef DEBUG
 300         fprintf(stdout,"empty\n");
 301 #endif
 302         Domain_Free(Q3);
 303       } /* end for j */
 304     } /* end for i */
 305     res = 1;
 306 LTQdone:
 307     Matrix_Free(Mat);
 308     Domain_Free(Q1);
 309     Domain_Free(Q2);
 310   }
 311   Domain_Free(Q);
 312
 313 #ifdef DEBUG
 314   fprintf(stdout, "res = %d\n", res);
 315 #endif
 316
 317   return res;
 318 } /* PolyhedronLTQ */
 319
 320 /* GaussSimplify --
 321    Given Mat1, a matrix of equalities, performs Gaussian elimination.
 322    Find a minimum basis, Returns the rank.
 323    Mat1 is context, Mat2 is reduced in context of Mat1
 324 */
 325 int GaussSimplify(Matrix *Mat1,Matrix *Mat2) {
 326
 327   int NbRows = Mat1->NbRows;
 328   int NbCols = Mat1->NbColumns;
 329   int *column_index;
 330   int i, j, k, n, pivot, Rank;
 331   Value gcd, tmp, *cp;
 332
 333   column_index=(int *)malloc(NbCols * sizeof(int));
 334   if (!column_index) {
 335     errormsg1("GaussSimplify", "outofmem", "out of memory space\n");
 336     Pol_status = 1;
 337     return 0;
 338   }
 339
 340   /* Initialize all the 'Value' variables */
 341   value_init(gcd); value_init(tmp);
 342
 343   Rank=0;
 344   for (j=0; j<NbCols; j++) {              /* for each column starting at */
 345     for (i=Rank; i<NbRows; i++)           /* diagonal, look down to find */
 346       if (value_notzero_p(Mat1->p[i][j])) /* the first non-zero entry    */
 347         break;
 348     if (i!=NbRows) {                      /* was one found ? */
 349       if (i!=Rank)                        /* was it found below the diagonal?*/
 350         Vector_Exchange(Mat1->p[Rank],Mat1->p[i],NbCols);
 351
 352       /* Normalize the pivot row */
 353       Vector_Gcd(Mat1->p[Rank],NbCols,&gcd);
 354
 355       /* If (gcd >= 2) */
 356       value_set_si(tmp,2);
 357       if (value_ge(gcd,tmp)) {
 358         cp = Mat1->p[Rank];
 359         for (k=0; k<NbCols; k++,cp++)
 360           value_division(*cp,*cp,gcd);
 361       }
 362       if (value_neg_p(Mat1->p[Rank][j])) {
 363         cp = Mat1->p[Rank];
 364         for (k=0; k<NbCols; k++,cp++)
 365           value_oppose(*cp,*cp);
 366       }
 367       /* End of normalize */
 368
 369       pivot=i;
 370       for (i=0;i<NbRows;i++)    /* Zero out the rest of the column */
 371         if (i!=Rank) {
 372           if (value_notzero_p(Mat1->p[i][j])) {
 373             Value a, a1, a2;
 374             value_init(a); value_init(a1); value_init(a2);
 375             value_absolute(a1,Mat1->p[i][j]);
 376             value_absolute(a2,Mat1->p[Rank][j]);
 377             value_gcd(a, a1, a2));
 378             value_divexact(a1, a1, a);
 379             value_divexact(a2, a2, a);
 380             value_oppose(a1,a1);
 381             Vector_Combine(Mat1->p[i],Mat1->p[Rank],Mat1->p[i],a2,
 382                            a1,NbCols);
 383             Vector_Normalize(Mat1->p[i],NbCols);
 384             value_clear(a); value_clear(a1); value_clear(a2);
 385           }
 386         }
 387       column_index[Rank]=j;
 388       Rank++;
 389     }
 390   } /* end of Gauss elimination */
 391
 392   if (Mat2) {  /* Mat2 is a transformation matrix  (i,j->f(i,j))....
 393                   can't scale it because can't scale both sides of -> */
 394     /* normalizes an affine transformation        */
 395     /* priority of forms                          */
 396     /*    1. i' -> i                (identity)    */
 397     /*    2. i' -> i + constant     (uniform)     */
 398     /*    3. i' -> constant         (broadcast)   */
 399     /*    4. i' -> j                (permutation) */
 400     /*    5. i' -> j + constant     (      )      */
 401     /*    6. i' -> i + j + constant (non-uniform) */
 402     for (k=0; k<Rank; k++) {
 403       j = column_index[k];
 404       for (i=0; i<(Mat2->NbRows-1);i++) {   /* all but the last row 0...0 1 */
 405         if ((i!=j) && value_notzero_p(Mat2->p[i][j])) {
 406
 407           /* Remove dependency of i' on j */
 408           Value a, a1, a2;
 409           value_init(a); value_init(a1); value_init(a2);
 410           value_absolute(a1,Mat2->p[i][j]);
 411           value_absolute(a2,Mat2->p[k][j]);
 412           value_gcd(a, a1, a2);
 413           value_divexact(a1, a1, a);
 414           value_divexact(a2, a2, a);
 415           value_oppose(a1,a1);
 416           if (value_one_p(a2)) {
 417             Vector_Combine(Mat2->p[i],Mat1->p[k],Mat2->p[i],a2,
 418                            a1,NbCols);
 419
 420             /* Vector_Normalize(Mat2->p[i],NbCols); -- can't do T        */
 421           } /* otherwise, can't do it without mult lhs prod (2i,3j->...) */
 422           value_clear(a); value_clear(a1); value_clear(a2);
 423         }
 424         else if ((i==j) && value_zero_p(Mat2->p[i][j])) {
 425
 426           /* 'i' does not depend on j */
 427           for (n=j+1; n < (NbCols-1); n++) {
 428             if (value_notzero_p(Mat2->p[i][n])) { /* i' depends on some n */
 429               value_set_si(tmp,1);
 430               Vector_Combine(Mat2->p[i],Mat1->p[k],Mat2->p[i],tmp,
 431                              tmp,NbCols);
 432               break;
 433             }  /* if 'i' depends on just a constant, then leave it alone.*/
 434           }
 435         }
 436       }
 437     }
 438
 439     /* Check last row of transformation Mat2 */
 440     for (j=0; j<(NbCols-1); j++)
 441       if (value_notzero_p(Mat2->p[Mat2->NbRows-1][j])) {
 442         errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
 443         break;
 444       }
 445
 446     if (value_notone_p(Mat2->p[Mat2->NbRows-1][NbCols-1])) {
 447       errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
 448     }
 449   }
 450   value_clear(gcd); value_clear(tmp);
 451   free(column_index);
 452   return Rank;
 453 } /* GaussSimplify */
 454
 455 char s[128];
 456
 457 int main() {
 458
 459   Matrix *a=NULL, *b=NULL, *c, *d, *e, *f;
 460   Polyhedron *A, *B, *C, *D, *last, *tmp;
 461   int i, nbPol, nbMat, func;
 462
 463   fgets(s, 128, stdin);
 464   nbPol = nbMat = 0;
 465   while ((*s=='#') ||
 466           ((sscanf(s, "D %d", &nbPol)<1) && (sscanf(s, "M %d", &nbMat)<1)) )
 467     fgets(s, 128, stdin);
 468
 469   for (i=0, A=last=(Polyhedron *)0; i<nbPol; i++) {
 470     a = Matrix_Read();
 471     tmp = Constraints2Polyhedron(a,600);
 472     Matrix_Free(a);
 473     if (!last) A = last = tmp;
 474     else {
 475       last->next = tmp;
 476       last = tmp;
 477     }
 478     }
 479
 480     if (nbMat) {
 481       a = Matrix_Read();
 482     }
 483     fgets(s,128,stdin);
 484     nbPol = nbMat = 0;
 485     while ( (*s=='#') ||
 486         ((sscanf(s, "D %d", &nbPol)<1) && (sscanf(s, "M %d", &nbMat)<1)) )
 487       fgets(s, 128, stdin);
 488
 489     for (i=0, B=last=(Polyhedron *)0; i<nbPol; i++) {
 490       b = Matrix_Read();
 491       tmp = Constraints2Polyhedron(b,200);
 492       Matrix_Free(b);
 493       if (!last) B = last = tmp;
 494       else {
 495         last->next = tmp;
 496         last = tmp;
 497       }
 498     }
 499
 500     if (nbMat) {
 501       b = Matrix_Read();
 502     }
 503     fgets(s, 128, stdin);
 504     while ((*s=='#') || (sscanf(s, "F %d", &func)<1))
 505       fgets(s, 128, stdin);
 506
 507     switch (func) {
 508     case 1:
 509       C = DomainUnion(A, B, 200);
 510       D = DomainConvex(C, 200);
 511       d = Polyhedron2Constraints(D);
 512       Matrix_Print(stdout,P_VALUE_FMT,d);
 513       Matrix_Free(d);
 514       Domain_Free(C);
 515       Domain_Free(D);
 516       break;
 517     case 2:
 518       D = DomainSimplify(A, B, 200);
 519       d = Polyhedron2Constraints(D);
 520       Matrix_Print(stdout,P_VALUE_FMT,d);
 521       Matrix_Free(d);
 522       Domain_Free(D);
 523       break;
 524     case 3:
 525       a = Polyhedron2Constraints(A);
 526       Matrix_Print(stdout,P_VALUE_FMT,a);
 527       b = Polyhedron2Constraints(B);
 528       Matrix_Print(stdout,P_VALUE_FMT,b);
 529       break;
 530     case 4:
 531       a = Polyhedron2Rays(A);
 532       Matrix_Print(stdout,P_VALUE_FMT,a);
 533       break;
 534     case 5:
 535
 536       /* a = ec , da = c , ed = 1 */
 537       right_hermite(a,&c,&d,&e);
 538       Matrix_Print(stdout,P_VALUE_FMT,c);
 539       Matrix_Print(stdout,P_VALUE_FMT,d);
 540       Matrix_Print(stdout,P_VALUE_FMT,e);
 541       f = Matrix_Alloc(e->NbRows,c->NbColumns);
 542       Matrix_Product(e,c,f);
 543       Matrix_Print(stdout,P_VALUE_FMT,f);
 544       Matrix_Free(f);
 545       f = Matrix_Alloc(d->NbRows,a->NbColumns);
 546       Matrix_Product(d,a,f);
 547       Matrix_Print(stdout,P_VALUE_FMT,f);
 548       Matrix_Free(f);
 549       f = Matrix_Alloc(e->NbRows, d->NbColumns);
 550       Matrix_Product(e,d,f);
 551       Matrix_Print(stdout,P_VALUE_FMT,f);
 552       break;
 553     case 6:
 554
 555       /* a = ce , ad = c , de = 1 */
 556       left_hermite(a,&c,&d,&e);
 557       Matrix_Print(stdout,P_VALUE_FMT,c);
 558       Matrix_Print(stdout,P_VALUE_FMT,d);
 559       Matrix_Print(stdout,P_VALUE_FMT,e);
 560       f = Matrix_Alloc(c->NbRows, e->NbColumns);
 561       Matrix_Product(c,e,f);
 562       Matrix_Print(stdout,P_VALUE_FMT,f);
 563       Matrix_Free(f);
 564       f = Matrix_Alloc(a->NbRows, d->NbColumns);
 565       Matrix_Product(a,d,f);
 566       Matrix_Print(stdout,P_VALUE_FMT,f);
 567       Matrix_Free(f);
 568       f = Matrix_Alloc(d->NbRows, e->NbColumns);
 569       Matrix_Product(d,e,f);
 570       Matrix_Print(stdout,P_VALUE_FMT,f);
 571       break;
 572     case 7:
 573
 574       /* Polyhedron_Print(stdout,"%5d", A); */
 575       /* Matrix_Print(stdout,"%4d", b);     */
 576
 577       C = Polyhedron_Image(A, b, 400);
 578       Polyhedron_Print(stdout,P_VALUE_FMT,C);
 579       break;
 580     case 8:
 581
 582       printf("%s\n",
 583              Polyhedron_Not_Empty(A,B,600) ? "Not Empty" : "Empty");
 584       break;
 585     case 9:
 586
 587       i = PolyhedronLTQ(A,B,600);
 588       printf("%s\n",
 589              i==-1 ? "A<B" : i==1 ? "A>B" : i==0 ? "A><B" : "error");
 590       i = PolyhedronLTQ(B,A,600);
 591       printf("%s\n",
 592              i==-1 ? "A<B" : i==1 ? "A>B" : i==0 ? "A><B" : "error");
 593       break;
 594     default:
 595       printf("? unknown function\n");
 596     }
 597
 598     Domain_Free(A);
 599     Domain_Free(B);
 600
 601     return 0;
 602 }
 603
 604