summate.cc

   1 #include <assert.h>
   2 #include <iostream>
   3 #include <barvinok/barvinok.h>
   4 #include <barvinok/options.h>
   5 #include <barvinok/util.h>
   6 #include "argp.h"
   7 #include "progname.h"
   8 #include "evalue_convert.h"
   9 #include "evalue_read.h"
  10 #include "verify.h"
  11
  12 using std::cout;
  13 using std::cerr;
  14 using std::endl;
  15
  16 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  17
  18 #define OPT_VARS            (BV_OPT_LAST+1)
  19
  20 struct argp_option argp_options[] = {
  21     { "variables",          OPT_VARS,   "list", 0,
  22         "comma separated list of variables over which to sum" },
  23     { "verbose",            'v',        0,      0, },
  24     { 0 }
  25 };
  26
  27 struct options {
  28     struct convert_options   convert;
  29     struct verify_options    verify;
  30     char* var_list;
  31     int verbose;
  32 };
  33
  34 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  35 {
  36     struct options *options = (struct options*) state->input;
  37
  38     switch (key) {
  39     case ARGP_KEY_INIT:
  40         state->child_inputs[0] = &options->convert;
  41         state->child_inputs[1] = &options->verify;
  42         state->child_inputs[2] = options->verify.barvinok;
  43         options->var_list = NULL;
  44         options->verbose = 0;
  45         break;
  46     case 'v':
  47         options->verbose = 1;
  48         break;
  49     case OPT_VARS:
  50         options->var_list = strdup(arg);
  51         break;
  52     default:
  53         return ARGP_ERR_UNKNOWN;
  54     }
  55     return 0;
  56 }
  57
  58 static int check_poly_sum(const struct check_poly_data *data,
  59                           int nparam, Value *z,
  60                           const struct verify_options *options);
  61
  62 struct check_poly_sum_data : public check_poly_data {
  63     Polyhedron          **S;
  64     evalue               *EP;
  65     evalue               *sum;
  66
  67     check_poly_sum_data(Value *z, evalue *EP, evalue *sum) :
  68                         EP(EP), sum(sum) {
  69         this->z = z;
  70         this->check = check_poly_sum;
  71     }
  72 };
  73
  74 static void sum(Polyhedron *S, int pos, const check_poly_sum_data *data,
  75                 evalue *s, const struct verify_options *options)
  76 {
  77     if (!S) {
  78         evalue *e = evalue_eval(data->EP, data->z+1);
  79         eadd(e, s);
  80         evalue_free(e);
  81     } else {
  82         Value LB, UB;
  83         int ok;
  84         value_init(LB);
  85         value_init(UB);
  86         ok = !(lower_upper_bounds(1+pos, S, data->z, &LB, &UB));
  87         assert(ok);
  88         for (; value_le(LB, UB); value_increment(LB, LB)) {
  89             value_assign(data->z[1+pos], LB);
  90             sum(S->next, pos+1, data, s, options);
  91         }
  92         value_set_si(data->z[1+pos], 0);
  93         value_clear(LB);
  94         value_clear(UB);
  95     }
  96 }
  97
  98 static evalue *sum(const check_poly_sum_data *data,
  99                           const struct verify_options *options)
 100 {
 101     evalue *s = evalue_zero();
 102     for (int i = 0; i < data->EP->x.p->size/2; ++i)
 103         if (!emptyQ2(data->S[i]))
 104             sum(data->S[i], 0, data, s, options);
 105     return s;
 106 }
 107
 108 static int check_poly_sum(const struct check_poly_data *data,
 109                           int nparam, Value *z,
 110                           const struct verify_options *options)
 111 {
 112     const check_poly_sum_data *sum_data;
 113     sum_data = static_cast<const check_poly_sum_data *>(data);
 114     evalue *e, *s;
 115     int k;
 116     int ok;
 117
 118     e = evalue_eval(sum_data->sum, z);
 119     if (options->print_all) {
 120         printf("sum(");
 121         value_print(stdout, VALUE_FMT, z[0]);
 122         for (k = 1; k < nparam; ++k) {
 123             printf(", ");
 124             value_print(stdout, VALUE_FMT, z[k]);
 125         }
 126         printf(") = ");
 127         value_print(stdout, VALUE_FMT, e->x.n);
 128         if (value_notone_p(e->d)) {
 129             printf("/");
 130             value_print(stdout, VALUE_FMT, e->d);
 131         }
 132     }
 133
 134     s = sum(sum_data, options);
 135
 136     if (options->print_all) {
 137         printf(", sum(EP) = ");
 138         value_print(stdout, VALUE_FMT, s->x.n);
 139         if (value_notone_p(s->d)) {
 140             printf("/");
 141             value_print(stdout, VALUE_FMT, s->d);
 142         }
 143         printf(". ");
 144     }
 145
 146     ok = eequal(e, s);
 147
 148     if (!ok) {
 149         printf("\n");
 150         fflush(stdout);
 151         fprintf(stderr,"Error !\n");
 152         fprintf(stderr,"sum( ");
 153         value_print(stderr, VALUE_FMT, z[0]);
 154         for (k = 1; k < nparam; ++k) {
 155             fprintf(stderr, ", ");
 156             value_print(stderr, VALUE_FMT, z[k]);
 157         }
 158         fprintf(stderr," ) should be ");
 159         value_print(stderr, VALUE_FMT, s->x.n);
 160         if (value_notone_p(s->d)) {
 161             fprintf(stderr, "/");
 162             value_print(stderr, VALUE_FMT, s->d);
 163         }
 164         fprintf(stderr,", while summation gives ");
 165         value_print(stderr, VALUE_FMT, e->x.n);
 166         if (value_notone_p(e->d)) {
 167             fprintf(stderr, "/");
 168             value_print(stderr, VALUE_FMT, e->d);
 169         }
 170         fprintf(stderr, ".\n");
 171     } else if (options->print_all)
 172         printf("OK.\n");
 173
 174     evalue_free(s);
 175     evalue_free(e);
 176
 177     return ok;
 178 }
 179
 180 static int verify(Polyhedron *P, evalue *sum, evalue *EP,
 181                   unsigned nvar, unsigned nparam, Vector *p,
 182                   struct verify_options *options)
 183 {
 184     Polyhedron *CS;
 185     unsigned MaxRays = options->barvinok->MaxRays;
 186     int error = 0;
 187
 188     CS = check_poly_context_scan(NULL, &P, P->Dimension, options);
 189
 190     check_poly_init(P, options);
 191
 192     if (!(CS && emptyQ2(CS))) {
 193         check_poly_sum_data data(p->p, EP, sum);
 194         data.S = ALLOCN(Polyhedron *, EP->x.p->size/2);
 195         for (int i = 0; i < EP->x.p->size/2; ++i) {
 196             Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 197             data.S[i] = Polyhedron_Scan(A, P, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 198         }
 199         error = !check_poly(CS, &data, nparam, 0, p->p+1+nvar, options);
 200         for (int i = 0; i < EP->x.p->size/2; ++i)
 201             Domain_Free(data.S[i]);
 202         free(data.S);
 203     }
 204
 205     if (!options->print_all)
 206         printf("\n");
 207
 208     if (CS) {
 209         Domain_Free(CS);
 210         Domain_Free(P);
 211     }
 212
 213     return error;
 214 }
 215
 216 /*
 217  * Project on final dim dimensions
 218  */
 219 Polyhedron *DomainProject(Polyhedron *D, unsigned dim, unsigned MaxRays)
 220 {
 221     Polyhedron *P;
 222     Polyhedron *R;
 223
 224     R = Polyhedron_Project(D, dim);
 225     for (P = D->next; P; P = P->next) {
 226         Polyhedron *R2 = Polyhedron_Project(P, dim);
 227         Polyhedron *R3 = DomainUnion(R, R2, MaxRays);
 228         Polyhedron_Free(R2);
 229         Domain_Free(R);
 230         R = R3;
 231     }
 232     return R;
 233 }
 234
 235 static int verify(evalue *EP, evalue *sum, unsigned nvar, unsigned nparam,
 236                   struct verify_options *options)
 237 {
 238     Vector *p;
 239
 240     p = Vector_Alloc(nvar+nparam+2);
 241     value_set_si(p->p[nvar+nparam+1], 1);
 242
 243     assert(value_zero_p(EP->d));
 244     assert(EP->x.p->type == partition);
 245
 246     Polyhedron *EP_D = EVALUE_DOMAIN(EP->x.p->arr[0]);
 247     Polyhedron *D = Polyhedron_Project(EP_D, nparam);
 248
 249     for (int i = 1; i < EP->x.p->size/2; ++i) {
 250         Polyhedron *D2 = D;
 251         EP_D = DomainProject(EVALUE_DOMAIN(EP->x.p->arr[2*i]), nparam,
 252                              options->barvinok->MaxRays);
 253         D = DomainUnion(EP_D, D, options->barvinok->MaxRays);
 254         Domain_Free(D2);
 255     }
 256
 257     int error = 0;
 258
 259     for (Polyhedron *P = D; P; P = P->next) {
 260         error = verify(P, sum, EP, nvar, nparam, p, options);
 261         if (error && !options->continue_on_error)
 262             break;
 263     }
 264
 265     Domain_Free(D);
 266     Vector_Free(p);
 267
 268     return error;
 269 }
 270
 271 int main(int argc, char **argv)
 272 {
 273     evalue *EP;
 274     char **all_vars = NULL;
 275     unsigned nvar;
 276     unsigned nparam;
 277     struct options options;
 278     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 279     static struct argp_child argp_children[] = {
 280         { &convert_argp,        0,      "input conversion",     1 },
 281         { &verify_argp,         0,      "verification",         2 },
 282         { &barvinok_argp,       0,      "barvinok options",     3 },
 283         { 0 }
 284     };
 285     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 286     int result = 0;
 287
 288     options.verify.barvinok = bv_options;
 289     set_program_name(argv[0]);
 290     argp_parse(&argp, argc, argv, 0, 0, &options);
 291
 292     EP = evalue_read_from_file(stdin, options.var_list, &all_vars,
 293                                &nvar, &nparam, bv_options->MaxRays);
 294     assert(EP);
 295
 296     if (options.verify.verify)
 297         verify_options_set_range(&options.verify, nvar+nparam);
 298
 299     evalue_convert(EP, &options.convert, options.verbose, nparam, all_vars);
 300
 301     if (EVALUE_IS_ZERO(*EP))
 302         print_evalue(stdout, EP, all_vars);
 303     else {
 304         evalue *sum = barvinok_summate(EP, nvar, bv_options);
 305         assert(sum);
 306         if (options.verify.verify)
 307             result = verify(EP, sum, nvar, nparam, &options.verify);
 308         else
 309             print_evalue(stdout, sum, all_vars+nvar);
 310         evalue_free(sum);
 311     }
 312
 313     evalue_free(EP);
 314
 315     if (options.var_list)
 316         free(options.var_list);
 317     Free_ParamNames(all_vars, nvar+nparam);
 318     barvinok_options_free(bv_options);
 319     return result;
 320 }