bound.c

   1 #include <assert.h>
   2 #include <isl/stream.h>
   3 #include <barvinok/options.h>
   4 #include <barvinok/util.h>
   5 #include <barvinok/barvinok.h>
   6 #include "verify.h"
   7
   8 struct options {
   9         struct verify_options    *verify;
  10         long split;
  11         int lower;
  12         long iterate;
  13 };
  14
  15 ISL_ARGS_START(struct options, options_args)
  16 ISL_ARG_CHILD(struct options, verify, NULL,
  17         &verify_options_args, "verification")
  18 ISL_ARG_LONG(struct options, split, 0, "split", 0, NULL)
  19 ISL_ARG_OPT_LONG(struct options, iterate, 0, "iterate", 0, -1,
  20         "exact result by iterating over domain (of specified maximal size)")
  21 ISL_ARG_BOOL(struct options, lower, 0, "lower", 0,
  22         "compute lower bound instead of upper bound")
  23 ISL_ARGS_END
  24
  25 ISL_ARG_DEF(options, struct options, options_args)
  26
  27 struct verify_point_bound {
  28         struct verify_point_data vpd;
  29         isl_pw_qpolynomial *pwqp;
  30         isl_pw_qpolynomial_fold *pwf;
  31         enum isl_fold type;
  32 };
  33
  34 static int verify_point(__isl_take isl_point *pnt, void *user)
  35 {
  36         int i;
  37         unsigned nparam;
  38         struct verify_point_bound *vpb = (struct verify_point_bound *) user;
  39         isl_int v, n, d, b, t;
  40         isl_pw_qpolynomial *pwqp;
  41         isl_qpolynomial *bound;
  42         isl_qpolynomial *opt;
  43         const char *minmax;
  44         int sign;
  45         int ok;
  46         int cst;
  47         FILE *out = vpb->vpd.options->print_all ? stdout : stderr;
  48
  49         vpb->vpd.n--;
  50
  51         if (vpb->type == isl_fold_max) {
  52                 minmax = "max";
  53                 sign = 1;
  54         } else {
  55                 minmax = "min";
  56                 sign = -1;
  57         }
  58
  59         isl_int_init(b);
  60         isl_int_init(t);
  61         isl_int_init(v);
  62         isl_int_init(n);
  63         isl_int_init(d);
  64
  65         pwqp = isl_pw_qpolynomial_copy(vpb->pwqp);
  66
  67         nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
  68         for (i = 0; i < nparam; ++i) {
  69                 isl_point_get_coordinate(pnt, isl_dim_param, i, &v);
  70                 pwqp = isl_pw_qpolynomial_fix_dim(pwqp, isl_dim_param, i, v);
  71         }
  72
  73         bound = isl_pw_qpolynomial_fold_eval(isl_pw_qpolynomial_fold_copy(vpb->pwf),
  74                                                 isl_point_copy(pnt));
  75
  76         if (sign > 0)
  77                 opt = isl_pw_qpolynomial_max(pwqp);
  78         else
  79                 opt = isl_pw_qpolynomial_min(pwqp);
  80
  81         cst = isl_qpolynomial_is_cst(opt, &n, &d);
  82         if (cst != 1)
  83                 goto error;
  84         if (sign > 0)
  85                 isl_int_fdiv_q(v, n, d);
  86         else
  87                 isl_int_cdiv_q(v, n, d);
  88
  89         cst = isl_qpolynomial_is_cst(bound, &n, &d);
  90         if (cst != 1)
  91                 goto error;
  92         if (sign > 0)
  93                 isl_int_fdiv_q(b, n, d);
  94         else
  95                 isl_int_cdiv_q(b, n, d);
  96
  97         if (sign > 0)
  98                 ok = value_ge(b, v);
  99         else
 100                 ok = value_le(b, v);
 101
 102         if (vpb->vpd.options->print_all || !ok) {
 103                 fprintf(out, "%s(", minmax);
 104                 for (i = 0; i < nparam; ++i) {
 105                         if (i)
 106                                 fprintf(out, ", ");
 107                         isl_point_get_coordinate(pnt, isl_dim_param, i, &t);
 108                         isl_int_print(out, t, 0);
 109                 }
 110                 fprintf(out, ") = ");
 111                 isl_int_print(out, n, 0);
 112                 if (!isl_int_is_one(d)) {
 113                         fprintf(out, "/");
 114                         isl_int_print(out, d, 0);
 115                         fprintf(out, " (");
 116                         isl_int_print(out, b, 0);
 117                         fprintf(out, ")");
 118                 }
 119                 fprintf(out, ", %s(EP) = ", minmax);
 120                 isl_int_print(out, v, 0);
 121                 if (ok)
 122                         fprintf(out, ". OK\n");
 123                 else
 124                         fprintf(out, ". NOT OK\n");
 125         } else if ((vpb->vpd.n % vpb->vpd.s) == 0) {
 126                 printf("o");
 127                 fflush(stdout);
 128         }
 129
 130         if (0) {
 131 error:
 132                 ok = 0;
 133         }
 134
 135         isl_qpolynomial_free(bound);
 136         isl_qpolynomial_free(opt);
 137         isl_point_free(pnt);
 138
 139         isl_int_clear(t);
 140         isl_int_clear(d);
 141         isl_int_clear(n);
 142         isl_int_clear(v);
 143         isl_int_clear(b);
 144
 145         if (!ok)
 146                 vpb->vpd.error = 1;
 147
 148         if (vpb->vpd.options->continue_on_error)
 149                 ok = 1;
 150
 151         return (vpb->vpd.n >= 1 && ok) ? 0 : -1;
 152 }
 153
 154 static int verify(__isl_keep isl_pw_qpolynomial_fold *pwf,
 155         __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type,
 156         struct verify_options *options)
 157 {
 158         struct verify_point_bound vpb = { { options } };
 159         isl_set *context;
 160         int r;
 161
 162         vpb.pwf = pwf;
 163         vpb.type = type;
 164         vpb.pwqp = pwqp;
 165         context = isl_pw_qpolynomial_fold_domain(
 166                                         isl_pw_qpolynomial_fold_copy(vpb.pwf));
 167         context = verify_context_set_bounds(context, options);
 168
 169         r = verify_point_data_init(&vpb.vpd, context);
 170
 171         if (r == 0)
 172                 isl_set_foreach_point(context, verify_point, &vpb);
 173         if (vpb.vpd.error)
 174                 r = -1;
 175
 176         isl_set_free(context);
 177         isl_pw_qpolynomial_free(pwqp);
 178
 179         verify_point_data_fini(&vpb.vpd);
 180
 181         return r;
 182 }
 183
 184 static __isl_give isl_pw_qpolynomial_fold *iterate(
 185         __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type)
 186 {
 187         isl_space *dim = isl_pw_qpolynomial_get_space(pwqp);
 188         isl_set *set;
 189         isl_qpolynomial *qp;
 190         isl_qpolynomial_fold *fold;
 191         unsigned nvar;
 192
 193         assert(isl_space_dim(dim, isl_dim_param) == 0);
 194         nvar = isl_space_dim(dim, isl_dim_set);
 195
 196         if (type == isl_fold_min)
 197                 qp = isl_pw_qpolynomial_min(pwqp);
 198         else
 199                 qp = isl_pw_qpolynomial_max(pwqp);
 200
 201         qp = isl_qpolynomial_drop_dims(qp, isl_dim_set, 0, nvar);
 202         fold = isl_qpolynomial_fold_alloc(type, qp);
 203         dim = isl_space_drop_dims(dim, isl_dim_set, 0, nvar);
 204         set = isl_set_universe(dim);
 205
 206         return isl_pw_qpolynomial_fold_alloc(type, set, fold);
 207 }
 208
 209 struct bv_split_data {
 210         int size;
 211         __isl_give isl_pw_qpolynomial *pwqp_less;
 212         __isl_give isl_pw_qpolynomial *pwqp_more;
 213 };
 214
 215 static int split_on_size(__isl_take isl_set *set,
 216         __isl_take isl_qpolynomial *qp, void *user)
 217 {
 218         struct bv_split_data *data = (struct bv_split_data *)user;
 219         int bounded;
 220         isl_set *set_np;
 221         isl_pw_qpolynomial *pwqp;
 222         int nparam;
 223
 224         nparam = isl_set_dim(set, isl_dim_param);
 225         set_np = isl_set_move_dims(isl_set_copy(set), isl_dim_set, 0,
 226                                         isl_dim_param, 0, nparam);
 227         bounded = isl_set_is_bounded(set_np);
 228         assert(bounded >= 0);
 229         if (bounded) {
 230                 isl_pw_qpolynomial *pwqp;
 231                 isl_qpolynomial *cst;
 232                 isl_int m;
 233                 int is_cst;
 234
 235                 pwqp = isl_set_card(set_np);
 236                 cst = isl_pw_qpolynomial_max(pwqp);
 237                 isl_int_init(m);
 238                 is_cst = isl_qpolynomial_is_cst(cst, &m, NULL);
 239                 isl_qpolynomial_free(cst);
 240                 assert(is_cst);
 241                 bounded = isl_int_cmp_si(m, data->size) <= 0;
 242                 isl_int_clear(m);
 243         } else
 244                 isl_set_free(set_np);
 245
 246         pwqp = isl_pw_qpolynomial_alloc(set, qp);
 247         if (bounded)
 248                 data->pwqp_less = isl_pw_qpolynomial_add_disjoint(
 249                                                 data->pwqp_less, pwqp);
 250         else
 251                 data->pwqp_more = isl_pw_qpolynomial_add_disjoint(
 252                                                 data->pwqp_more, pwqp);
 253
 254         return 0;
 255 }
 256
 257 /*
 258  * Split (partition) pwpq into a partition with (sub)domains containing
 259  * size integer points or less and a partition with (sub)domains
 260  * containing more integer points.
 261  */
 262 static int split_on_domain_size(__isl_take isl_pw_qpolynomial *pwqp,
 263         __isl_give isl_pw_qpolynomial **pwqp_less,
 264         __isl_give isl_pw_qpolynomial **pwqp_more,
 265         int size)
 266 {
 267         isl_space *dim;
 268         struct bv_split_data data = { size };
 269         int r;
 270
 271         dim = isl_pw_qpolynomial_get_space(pwqp);
 272         data.pwqp_less = isl_pw_qpolynomial_zero(isl_space_copy(dim));
 273         data.pwqp_more = isl_pw_qpolynomial_zero(dim);
 274         r = isl_pw_qpolynomial_foreach_piece(pwqp, &split_on_size, &data);
 275         *pwqp_less = data.pwqp_less;
 276         *pwqp_more = data.pwqp_more;
 277         isl_pw_qpolynomial_free(pwqp);
 278
 279         return r;
 280 }
 281
 282 static __isl_give isl_pw_qpolynomial_fold *optimize(
 283         __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type,
 284         struct options *options)
 285 {
 286     isl_pw_qpolynomial_fold *pwf;
 287
 288     if (options->iterate > 0) {
 289         isl_pw_qpolynomial *pwqp_less, *pwqp_more;
 290         isl_pw_qpolynomial_fold *pwf_less, *pwf_more;
 291         split_on_domain_size(pwqp, &pwqp_less, &pwqp_more, options->iterate);
 292         pwf_less = iterate(pwqp_less, type);
 293         pwf_more = isl_pw_qpolynomial_bound(pwqp_more, type, NULL);
 294         pwf = isl_pw_qpolynomial_fold_fold(pwf_less, pwf_more);
 295     } else if (options->iterate)
 296         pwf = iterate(pwqp, type);
 297     else
 298         pwf = isl_pw_qpolynomial_bound(pwqp, type, NULL);
 299     return pwf;
 300 }
 301
 302 static int optimize_and_print(__isl_take isl_pw_qpolynomial *pwqp,
 303         struct options *options)
 304 {
 305     int print_solution = 1;
 306     int result = 0;
 307     isl_pw_qpolynomial_fold *pwf;
 308     enum isl_fold type = options->lower ? isl_fold_min : isl_fold_max;
 309
 310     if (options->verify->verify) {
 311         isl_space *dim = isl_pw_qpolynomial_get_space(pwqp);
 312         unsigned total = isl_space_dim(dim, isl_dim_all);
 313         isl_space_free(dim);
 314         verify_options_set_range(options->verify, total);
 315         if (!options->verify->barvinok->verbose)
 316             print_solution = 0;
 317     }
 318
 319     pwf = optimize(isl_pw_qpolynomial_copy(pwqp), type, options);
 320     assert(pwf);
 321     if (print_solution) {
 322         isl_ctx *ctx = isl_pw_qpolynomial_get_ctx(pwqp);
 323         isl_printer *p = isl_printer_to_file(ctx, stdout);
 324         p = isl_printer_print_pw_qpolynomial_fold(p, pwf);
 325         p = isl_printer_end_line(p);
 326         isl_printer_free(p);
 327     }
 328     if (options->verify->verify) {
 329         enum isl_fold type = options->lower ? isl_fold_min : isl_fold_max;
 330         result = verify(pwf, pwqp, type, options->verify);
 331     } else
 332         isl_pw_qpolynomial_free(pwqp);
 333     isl_pw_qpolynomial_fold_free(pwf);
 334
 335     return result;
 336 }
 337
 338 int main(int argc, char **argv)
 339 {
 340     struct options *options = options_new_with_defaults();
 341     isl_ctx *ctx;
 342     isl_pw_qpolynomial *pwqp;
 343     struct isl_stream *s;
 344     int result = 0;
 345
 346     argc = options_parse(options, argc, argv, ISL_ARG_ALL);
 347     ctx = isl_ctx_alloc_with_options(&options_args, options);
 348
 349     s = isl_stream_new_file(ctx, stdin);
 350     pwqp = isl_stream_read_pw_qpolynomial(s);
 351
 352     if (options->split)
 353         pwqp = isl_pw_qpolynomial_split_periods(pwqp, options->split);
 354
 355     result = optimize_and_print(pwqp, options);
 356
 357     isl_stream_free(s);
 358     isl_ctx_free(ctx);
 359     return result;
 360 }