verify.c

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