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
 343     CS = check_poly_context_scan(NULL, &P, P->Dimension, options);
 344
 345     check_poly_init(P, options);
 346
 347     if (!(CS && emptyQ2(CS))) {
 348         check_EP_set_scan(data, P, MaxRays);
 349         ok = check_poly(CS, &data->cp, nparam, 0, data->cp.z+1+nvar, options);
 350         check_EP_clear_scan(data);
 351     }
 352
 353     if (!options->print_all)
 354         printf("\n");
 355
 356     if (CS) {
 357         Domain_Free(CS);
 358         Domain_Free(P);
 359     }
 360
 361     return ok;
 362 }
 363
 364 /*
 365  * Project on final dim dimensions
 366  */
 367 Polyhedron *DomainProject(Polyhedron *D, unsigned dim, unsigned MaxRays)
 368 {
 369     Polyhedron *P;
 370     Polyhedron *R;
 371
 372     R = Polyhedron_Project(D, dim);
 373     for (P = D->next; P; P = P->next) {
 374         Polyhedron *R2 = Polyhedron_Project(P, dim);
 375         Polyhedron *R3 = DomainUnion(R, R2, MaxRays);
 376         Polyhedron_Free(R2);
 377         Domain_Free(R);
 378         R = R3;
 379     }
 380     return R;
 381 }
 382
 383 static Polyhedron *evalue_parameter_domain(const evalue *e, unsigned nparam,
 384                                              unsigned MaxRays)
 385 {
 386     int i;
 387     Polyhedron *U = NULL;
 388
 389     if (EVALUE_IS_ZERO(*e))
 390         return Universe_Polyhedron(0);
 391
 392     assert(value_zero_p(e->d));
 393     assert(e->x.p->type == partition);
 394     assert(e->x.p->size >= 2);
 395
 396     for (i = 0; i < e->x.p->size/2; ++i) {
 397         Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
 398         Polyhedron *P = DomainProject(D, nparam, MaxRays);
 399         if (!U) {
 400             U = P;
 401         } else {
 402             Polyhedron *T = U;
 403             U = DomainUnion(U, P, MaxRays);
 404             Domain_Free(P);
 405             Domain_Free(T);
 406         }
 407     }
 408     return U;
 409 }
 410
 411 int check_EP(struct check_EP_data *data, unsigned nvar, unsigned nparam,
 412              struct verify_options *options)
 413 {
 414     Vector *p;
 415     int ok = 1;
 416     Polyhedron *D, *P;
 417
 418     p = Vector_Alloc(nvar+nparam+2);
 419     value_set_si(p->p[nvar+nparam+1], 1);
 420
 421     data->cp.z = p->p;
 422
 423     D = evalue_parameter_domain(data->EP, nparam, options->barvinok->MaxRays);
 424
 425     for (P = D; P; P = P->next) {
 426         ok = check_EP_on_poly(P, data, nvar, nparam, options);
 427         if (!ok && !options->continue_on_error)
 428             break;
 429     }
 430
 431     Domain_Free(D);
 432     Vector_Free(p);
 433
 434     return ok;
 435 }
 436
 437 __isl_give isl_set *verify_context_set_bounds(__isl_take isl_set *set,
 438         const struct verify_options *options)
 439 {
 440         int i;
 441         unsigned nparam;
 442         isl_point *pt, *pt2;
 443         isl_set *box;
 444
 445         nparam = isl_set_dim(set, isl_dim_param);
 446
 447         if (options->r > 0) {
 448                 pt = isl_set_sample_point(isl_set_copy(set));
 449                 pt2 = isl_point_copy(pt);
 450
 451                 for (i = 0; i < nparam; ++i) {
 452                         pt = isl_point_add_ui(pt, isl_dim_param, i, options->r);
 453                         pt2 = isl_point_sub_ui(pt2, isl_dim_param, i, options->r);
 454                 }
 455         } else {
 456                 isl_ctx *ctx;
 457                 isl_val *v;
 458
 459                 ctx = isl_set_get_ctx(set);
 460                 pt = isl_point_zero(isl_set_get_space(set));
 461                 v = isl_val_int_from_si(ctx, options->m);
 462                 for (i = 0; i < nparam; ++i)
 463                         pt = isl_point_set_coordinate_val(pt, isl_dim_param, i,
 464                                                             isl_val_copy(v));
 465                 isl_val_free(v);
 466
 467                 pt2 = isl_point_zero(isl_set_get_space(set));
 468                 v = isl_val_int_from_si(ctx, options->M);
 469                 for (i = 0; i < nparam; ++i)
 470                         pt2 = isl_point_set_coordinate_val(pt2, isl_dim_param,
 471                                                             i, isl_val_copy(v));
 472                 isl_val_free(v);
 473         }
 474
 475         box = isl_set_box_from_points(pt, pt2);
 476
 477         return isl_set_intersect(set, box);
 478 }
 479
 480 int verify_point_data_init(struct verify_point_data *vpd,
 481         __isl_keep isl_set *context)
 482 {
 483         isl_val *v;
 484         int i;
 485
 486         context = isl_set_copy(context);
 487         context = isl_set_move_dims(context, isl_dim_set, 0, isl_dim_param, 0,
 488                                         isl_set_dim(context, isl_dim_param));
 489         v = isl_pw_qpolynomial_max(isl_set_card(context));
 490         vpd->n = isl_val_cmp_si(v, 200) < 0 ? isl_val_get_num_si(v) : 200;
 491         isl_val_free(v);
 492
 493         if (!vpd->options->print_all) {
 494                 vpd->s = vpd->n < 80 ? 1 : 1 + vpd->n/80;
 495                 for (i = 0; i < vpd->n; i += vpd->s)
 496                         printf(".");
 497                 printf("\r");
 498                 fflush(stdout);
 499         }
 500
 501         vpd->error = !v ? -1 : 0;
 502
 503         return vpd->error;
 504 }
 505
 506 void verify_point_data_fini(struct verify_point_data *vpd)
 507 {
 508         if (!vpd->options->print_all)
 509                 printf("\n");
 510
 511         if (vpd->error)
 512                 fprintf(stderr, "Check failed !\n");
 513 }