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