verify.c

   1 #include <assert.h>
   2 #include <stdlib.h>
   3 #include <barvinok/options.h>
   4 #include <barvinok/util.h>
   5 #include "verify.h"
   6
   7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
   8
   9 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
  10 #define RANGE 50
  11
  12 /* SRANGE : small range for evalutations */
  13 #define SRANGE 15
  14
  15 /* if dimension >= BIDDIM, use SRANGE */
  16 #define BIGDIM 5
  17
  18 /* VSRANGE : very small range for evalutations */
  19 #define VSRANGE 5
  20
  21 /* if dimension >= VBIDDIM, use VSRANGE */
  22 #define VBIGDIM 8
  23
  24 static int set_m(void *opt, long val)
  25 {
  26         struct verify_options *options = (struct verify_options *)opt;
  27         options->m = val;
  28         options->verify = 1;
  29 }
  30
  31 static int set_M(void *opt, long val)
  32 {
  33         struct verify_options *options = (struct verify_options *)opt;
  34         options->M = val;
  35         options->verify = 1;
  36 }
  37
  38 static int set_r(void *opt, long val)
  39 {
  40         struct verify_options *options = (struct verify_options *)opt;
  41         options->r = val;
  42         options->m = -val;
  43         options->M = val;
  44         options->verify = 1;
  45 }
  46
  47 struct isl_arg verify_options_arg[] = {
  48 ISL_ARG_CHILD(struct verify_options, barvinok, NULL, barvinok_options_arg, NULL)
  49 ISL_ARG_BOOL(struct verify_options, verify, 'T', "verify", 0, NULL)
  50 ISL_ARG_BOOL(struct verify_options, exact, 'E', "exact", 0, NULL)
  51 ISL_ARG_BOOL(struct verify_options, print_all, 'A', "print-all", 0, NULL)
  52 ISL_ARG_BOOL(struct verify_options, continue_on_error, 'C',
  53         "continue-on-error", 0, NULL)
  54 ISL_ARG_USER_LONG(struct verify_options, m, 'm', "min", set_m, INT_MAX, NULL)
  55 ISL_ARG_USER_LONG(struct verify_options, M, 'M', "max", set_M, INT_MIN, NULL)
  56 ISL_ARG_USER_LONG(struct verify_options, r, 'r', "range", set_r, -1, NULL)
  57 ISL_ARG_END
  58 };
  59
  60 void verify_options_set_range(struct verify_options *options, int dim)
  61 {
  62     int r;
  63
  64     if (dim >= VBIGDIM)
  65         r = VSRANGE;
  66     else if (dim >= BIGDIM)
  67         r = SRANGE;
  68     else
  69         r = RANGE;
  70     /* If the user didn't set m or M, then we try to adjust the window
  71      * to the context in check_poly_context_scan.
  72      */
  73     if (options->m == INT_MAX && options->M == INT_MIN)
  74         options->r = r;
  75     if (options->M == INT_MIN)
  76         options->M = r;
  77     if (options->m == INT_MAX)
  78         options->m = -r;
  79
  80     if (options->verify && options->m > options->M) {
  81         fprintf(stderr,"Nothing to do: min > max !\n");
  82         exit(0);
  83     }
  84 }
  85
  86 static Polyhedron *project_on(Polyhedron *P, int i)
  87 {
  88     unsigned dim = P->Dimension;
  89     Matrix *T;
  90     Polyhedron *I;
  91
  92     if (dim == 1)
  93         return Polyhedron_Copy(P);
  94
  95     T = Matrix_Alloc(2, dim+1);
  96     value_set_si(T->p[0][i], 1);
  97     value_set_si(T->p[1][dim], 1);
  98     I = Polyhedron_Image(P, T, P->NbConstraints);
  99     Matrix_Free(T);
 100     return I;
 101 }
 102
 103 static void set_bounds(Polyhedron *C, Value **rows, int i, unsigned nparam,
 104                         const struct verify_options *options)
 105 {
 106     Value min;
 107     Value max;
 108
 109     value_init(min);
 110     value_init(max);
 111     value_set_si(min, options->m);
 112     value_set_si(max, options->M);
 113
 114     if (options->r > 0) {
 115         Polyhedron *I = project_on(C, i);
 116         line_minmax(I, &min, &max);
 117         if (value_cmp_si(min, options->M) >= 0)
 118             value_add_int(max, min, options->r);
 119         else if (value_cmp_si(max, options->m) <= 0)
 120             value_sub_int(min, max, options->r);
 121         else {
 122             value_set_si(min, options->m);
 123             value_set_si(max, options->M);
 124         }
 125     }
 126
 127     value_set_si(rows[0][0], 1);
 128     value_set_si(rows[0][1+i], 1);
 129     value_oppose(rows[0][1+nparam], min);
 130     value_set_si(rows[1][0], 1);
 131     value_set_si(rows[1][1+i], -1);
 132     value_assign(rows[1][1+nparam], max);
 133
 134     value_clear(min);
 135     value_clear(max);
 136 }
 137
 138 Polyhedron *check_poly_context_scan(Polyhedron *P, Polyhedron **C,
 139                                     unsigned nparam,
 140                                     const struct verify_options *options)
 141 {
 142     int i;
 143     Matrix *MM;
 144     Polyhedron *CC, *CC2, *CS, *U;
 145     unsigned MaxRays = options->barvinok->MaxRays;
 146
 147     if (nparam <= 0)
 148         return NULL;
 149
 150     if (!P)
 151         CC = *C;
 152     else {
 153         CC = Polyhedron_Project(P, nparam);
 154         if (*C) {
 155             CC2 = DomainIntersection(*C, CC, MaxRays);
 156             Domain_Free(CC);
 157             CC = CC2;
 158         }
 159     }
 160
 161     /* Intersect context with range */
 162     MM = Matrix_Alloc(2*nparam, nparam+2);
 163     for (i = 0; i < nparam; ++i)
 164         set_bounds(CC, &MM->p[2*i], i, nparam, options);
 165     CC2 = AddConstraints(MM->p[0], 2*nparam, CC, options->barvinok->MaxRays);
 166     if (CC != *C)
 167         Domain_Free(CC);
 168     CC = CC2;
 169     U = Universe_Polyhedron(0);
 170     CS = Polyhedron_Scan(CC, U, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 171     Polyhedron_Free(U);
 172     *C = CC;
 173     Matrix_Free(MM);
 174     return CS;
 175 }
 176
 177 void check_poly_init(Polyhedron *C, struct verify_options *options)
 178 {
 179     int d, i;
 180
 181     if (options->print_all)
 182         return;
 183     if (C->Dimension <= 0)
 184         return;
 185
 186     d = options->M - options->m;
 187     if (d > 80)
 188         options->st = 1+d/80;
 189     else
 190         options->st = 1;
 191     for (i = options->m; i <= options->M; i += options->st)
 192         printf(".");
 193     printf( "\r" );
 194     fflush(stdout);
 195 }
 196
 197 static void print_rational(FILE *out, Value n, Value d)
 198 {
 199     value_print(out, VALUE_FMT, n);
 200     if (value_notone_p(d)) {
 201         fprintf(out, "/");
 202         value_print(out, VALUE_FMT, d);
 203     }
 204 }
 205
 206 void check_poly_print(int ok, int nparam, Value *z,
 207                       Value want_n, Value want_d,
 208                       Value got_n, Value got_d,
 209                       const char *op, const char *check,
 210                       const char *long_op,
 211                       const struct verify_options *options)
 212 {
 213     int k;
 214
 215     if (options->print_all) {
 216         printf("%s(", op);
 217         value_print(stdout, VALUE_FMT, z[0]);
 218         for (k = 1; k < nparam; ++k) {
 219             printf(", ");
 220             value_print(stdout, VALUE_FMT, z[k]);
 221         }
 222         printf(") = ");
 223         print_rational(stdout, got_n, got_d);
 224         printf(", %s = ", check);
 225         print_rational(stdout, want_n, want_d);
 226         printf(". ");
 227     }
 228
 229     if (!ok) {
 230         printf("\n");
 231         fflush(stdout);
 232         fprintf(stderr, "Error !\n");
 233         fprintf(stderr, "%s(", op);
 234         value_print(stderr, VALUE_FMT, z[0]);
 235         for (k = 1; k < nparam; ++k) {
 236             fprintf(stderr,", ");
 237             value_print(stderr, VALUE_FMT, z[k]);
 238         }
 239         fprintf(stderr, ") should be ");
 240         print_rational(stderr, want_n, want_d);
 241         fprintf(stderr, ", while %s gives ", long_op);
 242         print_rational(stderr, got_n, got_d);
 243         fprintf(stderr, ".\n");
 244     } else if (options->print_all)
 245         printf("OK.\n");
 246 }
 247
 248 /****************************************************/
 249 /* function check_poly :                            */
 250 /* scans the parameter space from Min to Max (all   */
 251 /* directions). Computes the number of points in    */
 252 /* the polytope using both methods, and compare them*/
 253 /* returns 1 on success                             */
 254 /****************************************************/
 255
 256 int check_poly(Polyhedron *CS, const struct check_poly_data *data,
 257                int nparam, int pos, Value *z,
 258                const struct verify_options *options)
 259 {
 260     if (pos == nparam) {
 261         if (!data->check(data, nparam, z, options) && !options->continue_on_error)
 262             return 0;
 263     } else {
 264         Value LB, UB;
 265         int ok;
 266         value_init(LB);
 267         value_init(UB);
 268         ok = !(lower_upper_bounds(1+pos, CS, z-1, &LB, &UB));
 269         assert(ok);
 270         for (; value_le(LB, UB); value_increment(LB, LB)) {
 271             if (!options->print_all) {
 272                 int k = VALUE_TO_INT(LB);
 273                 if (!pos && !(k % options->st)) {
 274                     printf("o");
 275                     fflush(stdout);
 276                 }
 277             }
 278
 279             value_assign(z[pos], LB);
 280             if (!check_poly(CS->next, data, nparam, pos+1, z, options)) {
 281                 value_clear(LB);
 282                 value_clear(UB);
 283                 return 0;
 284             }
 285         }
 286         value_set_si(z[pos], 0);
 287         value_clear(LB);
 288         value_clear(UB);
 289     }
 290     return 1;
 291 } /* check_poly */
 292
 293 void check_EP_set_scan(struct check_EP_data *data, Polyhedron *C,
 294                         unsigned MaxRays)
 295 {
 296     const evalue *EP = data->EP;
 297     int i;
 298     int n_S = 0;
 299
 300     for (i = 0; i < EP->x.p->size/2; ++i) {
 301         Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 302         for (; A; A = A->next)
 303             ++n_S;
 304     }
 305
 306     data->n_S = n_S;
 307     data->S = ALLOCN(Polyhedron *, n_S);
 308     n_S = 0;
 309     for (i = 0; i < EP->x.p->size/2; ++i) {
 310         Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 311         for (; A; A = A->next) {
 312             Polyhedron *next = A->next;
 313             A->next = NULL;
 314             data->S[n_S++] = Polyhedron_Scan(A, C,
 315                                     MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 316             A->next = next;
 317         }
 318     }
 319 }
 320
 321 void check_EP_clear_scan(struct check_EP_data *data)
 322 {
 323     int i;
 324
 325     for (i = 0; i < data->n_S; ++i)
 326         Domain_Free(data->S[i]);
 327     free(data->S);
 328 }
 329
 330 static int check_EP_on_poly(Polyhedron *P,
 331                             struct check_EP_data *data,
 332                             unsigned nvar, unsigned nparam,
 333                             struct verify_options *options)
 334 {
 335     Polyhedron *CS;
 336     unsigned MaxRays = options->barvinok->MaxRays;
 337     int ok = 1;
 338     const evalue *EP = data->EP;
 339
 340     CS = check_poly_context_scan(NULL, &P, P->Dimension, options);
 341
 342     check_poly_init(P, options);
 343
 344     if (!(CS && emptyQ2(CS))) {
 345         check_EP_set_scan(data, P, MaxRays);
 346         ok = check_poly(CS, &data->cp, nparam, 0, data->cp.z+1+nvar, options);
 347         check_EP_clear_scan(data);
 348     }
 349
 350     if (!options->print_all)
 351         printf("\n");
 352
 353     if (CS) {
 354         Domain_Free(CS);
 355         Domain_Free(P);
 356     }
 357
 358     return ok;
 359 }
 360
 361 /*
 362  * Project on final dim dimensions
 363  */
 364 Polyhedron *DomainProject(Polyhedron *D, unsigned dim, unsigned MaxRays)
 365 {
 366     Polyhedron *P;
 367     Polyhedron *R;
 368
 369     R = Polyhedron_Project(D, dim);
 370     for (P = D->next; P; P = P->next) {
 371         Polyhedron *R2 = Polyhedron_Project(P, dim);
 372         Polyhedron *R3 = DomainUnion(R, R2, MaxRays);
 373         Polyhedron_Free(R2);
 374         Domain_Free(R);
 375         R = R3;
 376     }
 377     return R;
 378 }
 379
 380 static Polyhedron *evalue_parameter_domain(const evalue *e, unsigned nparam,
 381                                              unsigned MaxRays)
 382 {
 383     int i;
 384     Polyhedron *U = NULL;
 385
 386     if (EVALUE_IS_ZERO(*e))
 387         return Universe_Polyhedron(0);
 388
 389     assert(value_zero_p(e->d));
 390     assert(e->x.p->type == partition);
 391     assert(e->x.p->size >= 2);
 392
 393     for (i = 0; i < e->x.p->size/2; ++i) {
 394         Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
 395         Polyhedron *P = DomainProject(D, nparam, MaxRays);
 396         if (!U) {
 397             U = P;
 398         } else {
 399             Polyhedron *T = U;
 400             U = DomainUnion(U, P, MaxRays);
 401             Domain_Free(P);
 402             Domain_Free(T);
 403         }
 404     }
 405     return U;
 406 }
 407
 408 int check_EP(struct check_EP_data *data, unsigned nvar, unsigned nparam,
 409              struct verify_options *options)
 410 {
 411     Vector *p;
 412     int i;
 413     int ok = 1;
 414     Polyhedron *D, *P;
 415
 416     p = Vector_Alloc(nvar+nparam+2);
 417     value_set_si(p->p[nvar+nparam+1], 1);
 418
 419     data->cp.z = p->p;
 420
 421     D = evalue_parameter_domain(data->EP, nparam, options->barvinok->MaxRays);
 422
 423     for (P = D; P; P = P->next) {
 424         ok = check_EP_on_poly(P, data, nvar, nparam, options);
 425         if (!ok && !options->continue_on_error)
 426             break;
 427     }
 428
 429     Domain_Free(D);
 430     Vector_Free(p);
 431
 432     return ok;
 433 }
 434
 435 __isl_give isl_set *verify_context_set_bounds(__isl_take isl_set *set,
 436         const struct verify_options *options)
 437 {
 438         int i;
 439         unsigned nparam;
 440         isl_point *pt, *pt2;
 441         isl_set *box;
 442
 443         nparam = isl_set_dim(set, isl_dim_param);
 444
 445         if (options->r > 0) {
 446                 pt = isl_set_sample_point(isl_set_copy(set));
 447                 pt2 = isl_point_copy(pt);
 448
 449                 for (i = 0; i < nparam; ++i) {
 450                         pt = isl_point_add_ui(pt, isl_dim_param, i, options->r);
 451                         pt2 = isl_point_sub_ui(pt2, isl_dim_param, i, options->r);
 452                 }
 453         } else {
 454                 isl_int v;
 455
 456                 isl_int_init(v);
 457                 pt = isl_point_zero(isl_set_get_dim(set));
 458                 isl_int_set_si(v, options->m);
 459                 for (i = 0; i < nparam; ++i)
 460                         pt = isl_point_set_coordinate(pt, isl_dim_param, i, v);
 461
 462                 pt2 = isl_point_zero(isl_set_get_dim(set));
 463                 isl_int_set_si(v, options->M);
 464                 for (i = 0; i < nparam; ++i)
 465                         pt2 = isl_point_set_coordinate(pt2, isl_dim_param, i, v);
 466
 467                 isl_int_clear(v);
 468         }
 469
 470         box = isl_set_box_from_points(pt, pt2);
 471
 472         return isl_set_intersect(set, box);
 473 }
 474
 475 int verify_point_data_init(struct verify_point_data *vpd,
 476         __isl_keep isl_set *context)
 477 {
 478         isl_int v;
 479         int i;
 480         int r;
 481
 482         isl_int_init(v);
 483         r = isl_set_count(context, &v);
 484         vpd->n = isl_int_cmp_si(v, 200) < 0 ? isl_int_get_si(v) : 200;
 485         isl_int_clear(v);
 486
 487         if (!vpd->options->print_all) {
 488                 vpd->s = vpd->n < 80 ? 1 : 1 + vpd->n/80;
 489                 for (i = 0; i < vpd->n; i += vpd->s)
 490                         printf(".");
 491                 printf("\r");
 492                 fflush(stdout);
 493         }
 494
 495         vpd->error = r < 0 ? -1 : 0;
 496
 497         return r;
 498 }
 499
 500 void verify_point_data_fini(struct verify_point_data *vpd)
 501 {
 502         if (!vpd->options->print_all)
 503                 printf("\n");
 504
 505         if (vpd->error)
 506                 fprintf(stderr, "Check failed !\n");
 507 }