verify_lexsmaller.c

   1 #include <polylib/polylibgmp.h>
   2 #include <barvinok/util.h>
   3 #include "config.h"
   4
   5 #ifdef HAVE_GROWING_CHERNIKOVA
   6 #define MAXRAYS    POL_NO_DUAL
   7 #else
   8 #define MAXRAYS  600
   9 #endif
  10
  11 #include "config.h"
  12 #ifndef HAVE_GETOPT_H
  13 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
  14 #else
  15 #include <getopt.h>
  16 struct option options[] = {
  17     { "continue",  no_argument,  0,  'k' },
  18     { "max",  no_argument,  0,  'M' },
  19     { "live",  no_argument,  0,  'l' },
  20     { "verbose",  no_argument,  0,  'v' },
  21     { "version",   no_argument,  0,  'V' },
  22     { 0, 0, 0, 0 }
  23 };
  24 #endif
  25
  26 static int live = 0;
  27 static int print_max = 0;
  28 static int verbose = 0;
  29 static int keep_going = 0;
  30 static Value max;
  31
  32 #define LS_OK       1
  33 #define LS_P        2       /* continue searching P */
  34 #define LS_D        4       /* continue searching D */
  35
  36 static int check_lexsmaller(Polyhedron *SP, Polyhedron *SD, Enumeration *en,
  37                          int pos, int nvar, Value *zP, Value *zD, Value *zE,
  38                          Value *count)
  39 {
  40     int i;
  41     int ok;
  42     Value PLB, PUB, DLB, DUB, LB, UB, tmp, c;
  43
  44     if (!SP && !SD)
  45         return LS_OK | LS_P | LS_D;
  46
  47     value_init(PLB); value_init(PUB);
  48     value_init(DLB); value_init(DUB);
  49     value_init(LB); value_init(UB);
  50     value_init(tmp);
  51
  52     if (SP) {
  53         ok = !(lower_upper_bounds(1+pos, SP, zP, &PLB, &PUB));
  54         assert(ok);
  55     }
  56     if (SD) {
  57         ok = !(lower_upper_bounds(1+pos, SD, zD, &DLB, &DUB));
  58         assert(ok);
  59     }
  60     if (!SD || (SP && value_lt(PLB, DLB)))
  61         value_assign(LB, PLB);
  62     else
  63         value_assign(LB, DLB);
  64     if (!SD || (SP && value_gt(PUB, DUB)))
  65         value_assign(UB, PUB);
  66     else
  67         value_assign(UB, DUB);
  68
  69     if (SD && !SD->next)
  70         value_init(c);
  71
  72     ok = LS_OK | LS_P | LS_D;
  73
  74     for(value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
  75         int inP = SP && value_ge(tmp, PLB) && value_le(tmp, PUB);
  76         int inD = SD && value_ge(tmp, DLB) && value_le(tmp, DUB);
  77         if (!inP && !inD)
  78             continue;
  79
  80         if (inP)
  81             value_assign(zP[pos+1], tmp);
  82         if (inD)
  83             value_assign(zD[pos+1], tmp);
  84         if (inD && pos < nvar)
  85             value_assign(zE[pos], tmp);
  86
  87         if (inD && !SD->next) {
  88             Value *ctmp;
  89
  90             value_assign(c,*(ctmp=compute_poly(en, zE)));
  91             value_clear(*ctmp);
  92             free(ctmp);
  93
  94             if (verbose >= 2) {
  95                 printf("EP( ");
  96                 value_print(stdout, VALUE_FMT, zE[0]);
  97                 for (i = 1; i < nvar; ++i) {
  98                     printf(", ");
  99                     value_print(stdout, VALUE_FMT, zE[i]);
 100                 }
 101                 printf(" ) = ");
 102                 value_print(stdout, VALUE_FMT, c);
 103                 printf("; count = ");
 104                 value_print(stdout, VALUE_FMT, *count);
 105                 printf("\n");
 106             }
 107
 108             if (value_ne(c, *count)) {
 109                 printf("\n");
 110                 fflush(stdout);
 111                 fprintf(stderr,"Error !\n");
 112                 fprintf(stderr, "EP( ");
 113                 value_print(stderr, VALUE_FMT, zE[0]);
 114                 for (i = 1; i < nvar; ++i) {
 115                     fprintf(stderr, ", ");
 116                     value_print(stderr, VALUE_FMT, zE[i]);
 117                 }
 118                 fprintf(stderr, " ) = ");
 119                 value_print(stderr, VALUE_FMT, c);
 120                 fprintf(stderr, " but count = ");
 121                 value_print(stderr, VALUE_FMT, *count);
 122                 printf("\n");
 123                 ok = 0;
 124             }
 125
 126             if (live)
 127                 value_decrement(*count, *count);
 128
 129             ok &= ~LS_D;
 130         }
 131
 132         if (pos < nvar-1)
 133             ok &= check_lexsmaller(inP ? SP->next : NULL,
 134                                    inD ? SD->next : NULL,
 135                                    en, pos+1, nvar, zP, zD, zE, count);
 136         else {
 137             ok &= check_lexsmaller(NULL, inD ? SD->next : NULL,
 138                                    en, pos+1, nvar, zP, zD, zE, count)
 139                &  check_lexsmaller(inP ? SP->next : NULL, NULL,
 140                                    en, pos+1, nvar, zP, zD, zE, count);
 141             if (pos >= nvar && !(ok & LS_D))
 142                 break;
 143             if (pos >= nvar && !(ok & LS_P))
 144                 break;
 145         }
 146
 147         if (!ok && !keep_going)
 148             goto end;
 149
 150         if (inP && !SP->next) {
 151             value_increment(*count, *count);
 152             if (value_gt(*count, max))
 153                 value_assign(max, *count);
 154             ok &= ~LS_P;
 155         }
 156     }
 157     if (SP)
 158         value_set_si(zP[pos+1], 0);
 159     if (SD)
 160         value_set_si(zD[pos+1], 0);
 161
 162 end:
 163     if (SD && !SD->next)
 164         value_clear(c);
 165
 166     value_clear(PLB); value_clear(PUB);
 167     value_clear(DLB); value_clear(DUB);
 168     value_clear(LB); value_clear(UB);
 169     value_clear(tmp);
 170
 171     return ok;
 172 }
 173
 174 int main(int argc,char *argv[])
 175 {
 176     Matrix *M;
 177     Polyhedron *P, *D, *C;
 178     Polyhedron *SP, *SD;
 179     int nb_parms;
 180     char **param_name = NULL;
 181     evalue *EP;
 182     Enumeration *en;
 183     Vector *zP, *zD, *zE;
 184     Value count;
 185     int c, ind = 0;
 186     char s[128];
 187     unsigned dim;
 188
 189     while ((c = getopt_long(argc, argv, "klMvV", options, &ind)) != -1) {
 190         switch (c) {
 191         case 'k':
 192             keep_going = 1;
 193             break;
 194         case 'M':
 195             print_max = 1;
 196         case 'l':
 197             live = 1;
 198             break;
 199         case 'v':
 200             ++verbose;
 201             break;
 202         case 'V':
 203             printf(barvinok_version());
 204             exit(0);
 205             break;
 206         }
 207     }
 208
 209     M = Matrix_Read();
 210     P = Constraints2Polyhedron(M, MAXRAYS);
 211     assert(P != NULL);
 212     Matrix_Free(M);
 213     M = Matrix_Read();
 214     D = Constraints2Polyhedron(M, MAXRAYS);
 215     assert(D != NULL);
 216     Matrix_Free(M);
 217
 218     fgets(s, 128, stdin);
 219     while ((*s=='#') || (sscanf(s, "D %u", &dim)<1))
 220         fgets(s, 128, stdin);
 221
 222     M = Matrix_Read();
 223     C = Constraints2Polyhedron(M, MAXRAYS);
 224     assert(C != NULL);
 225     Matrix_Free(M);
 226
 227     nb_parms = D->Dimension;
 228     param_name = Read_ParamNames(stdin, nb_parms);
 229
 230     EP = barvinok_lexsmaller_ev(P, D, dim, C, MAXRAYS);
 231     if (live) {
 232         evalue mone;
 233         evalue *EC = barvinok_lexsmaller_ev(D, D, dim, C, MAXRAYS);
 234         if (verbose >= 2) {
 235             puts("EP");
 236             print_evalue(stdout, EP, param_name);
 237             puts("EC");
 238             print_evalue(stdout, EC, param_name);
 239         }
 240         value_init(mone.d);
 241         evalue_set_si(&mone, -1, 1);
 242         emul(&mone, EC);
 243         eadd(EC, EP);
 244         free_evalue_refs(&mone);
 245         free_evalue_refs(EC);
 246         free(EC);
 247         reduce_evalue(EP);
 248     }
 249     if (verbose >= 1) {
 250         puts("Enumeration");
 251         print_evalue(stdout, EP, param_name);
 252     }
 253     en = partition2enumeration(EP);
 254
 255     assert(C->Dimension == 0); /* for now */
 256
 257     /* S = scanning list of polyhedra */
 258     SP = Polyhedron_Scan(P, C, MAXRAYS);
 259     SD = Polyhedron_Scan(D, C, MAXRAYS);
 260
 261     zP = Vector_Alloc(1+P->Dimension+1);
 262     value_set_si(zP->p[1+P->Dimension], 1);
 263     zD = Vector_Alloc(1+D->Dimension+1);
 264     value_set_si(zD->p[1+D->Dimension], 1);
 265     zE = Vector_Alloc(dim+C->Dimension);
 266
 267     if (print_max)
 268         value_init(max);
 269     value_init(count);
 270     check_lexsmaller(SP, SD, en, 0, dim, zP->p, zD->p, zE->p, &count);
 271     value_clear(count);
 272
 273     if (print_max) {
 274         printf("max = ");
 275         value_print(stdout, VALUE_FMT, max);
 276         printf("\n");
 277         value_clear(max);
 278     }
 279
 280     Enumeration_Free(en);
 281     Free_ParamNames(param_name, nb_parms);
 282     Polyhedron_Free(P);
 283     Polyhedron_Free(D);
 284     Polyhedron_Free(C);
 285     Domain_Free(SP);
 286     Domain_Free(SD);
 287     Vector_Free(zP);
 288     Vector_Free(zD);
 289     Vector_Free(zE);
 290 }