verify.c

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