evalue_convert.cc

   1 #include <assert.h>
   2 #include <unistd.h>
   3 #include <sstream>
   4 #include <isl_dim.h>
   5 #include <barvinok/util.h>
   6 #include "conversion.h"
   7 #include "evalue_convert.h"
   8 #include "lattice_point.h"
   9 #include "config.h"
  10 #ifdef USE_FDSTREAM
  11 #include "fdstream.h"
  12 #endif
  13
  14 using std::cout;
  15 using std::cerr;
  16 using std::endl;
  17
  18 using std::string;
  19
  20 static struct argp_option argp_options[] = {
  21     { "convert",            'c', 0, 0, "convert fractionals to periodics" },
  22     { "combine",            'C', 0, 0 },
  23     { "floor",              'f', 0, 0, "convert fractionals to floorings" },
  24     { "list",               'l', 0, 0 },
  25     { "latex",              'L', 0, 0 },
  26     { "to-isl",             'I', 0, 0 },
  27     { "range-reduction",    'R',    0,      0 },
  28     0
  29 };
  30
  31 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  32 {
  33     struct convert_options *options = (struct convert_options *)state->input;
  34
  35     switch (key) {
  36     case ARGP_KEY_INIT:
  37         options->floor = 0;
  38         options->convert = 0;
  39         options->combine = 0;
  40         options->range = 0;
  41         options->list = 0;
  42         options->latex = 0;
  43         options->isl = 0;
  44         break;
  45     case ARGP_KEY_FINI:
  46         break;
  47     case 'f':
  48         options->floor = 1;
  49         break;
  50     case 'c':
  51         options->convert = 1;
  52         break;
  53     case 'C':
  54         options->combine = 1;
  55         break;
  56     case 'l':
  57         options->list = 1;
  58         break;
  59     case 'L':
  60         options->latex = 1;
  61         break;
  62     case 'R':
  63         options->range = 1;
  64         break;
  65     case 'I':
  66         options->isl = 1;
  67         break;
  68     default:
  69         return ARGP_ERR_UNKNOWN;
  70     }
  71     return 0;
  72 }
  73
  74 struct argp convert_argp = {
  75     argp_options, parse_opt, 0, 0
  76 };
  77
  78 static int type_offset(enode *p)
  79 {
  80    return p->type == fractional ? 1 :
  81           p->type == flooring ? 1 : 0;
  82 }
  83
  84 static Lattice *extract_lattice(evalue *e, int nparam)
  85 {
  86     int i;
  87     Lattice *L;
  88     Matrix *U;
  89     Vector *X;
  90     /* For some mysterious reason, SolveDiophantine expects an extra
  91      * [0 0 0 1] row in its input matrix.
  92      */
  93     Matrix *M = Matrix_Alloc(2, nparam+1+1);
  94     value_set_si(M->p[1][nparam+1], 1);
  95     evalue_extract_affine(e, M->p[0], M->p[0]+nparam+1, M->p[0]+nparam);
  96     /* ignore constant */
  97     value_set_si(M->p[0][nparam+1], 0);
  98     SolveDiophantine(M, &U, &X);
  99     Matrix_Free(M);
 100     Vector_Free(X);
 101     L = Matrix_Alloc(nparam+1, nparam+1);
 102     for (i = 0; i < nparam; ++i)
 103         Vector_Copy(U->p[i], L->p[i], nparam);
 104     value_set_si(L->p[nparam][nparam], 1);
 105     Matrix_Free(U);
 106     return L;
 107 }
 108
 109 /* Returns a lattice such that the quasi-polynomial e can be represented
 110  * by a list of polynomials, one for each point in the fundamental
 111  * parallelepiped of the lattice.
 112  * If e is a polynomial, then this function returns NULL.
 113  */
 114 static Lattice *extract_common_lattice(evalue *e, Lattice *L, int nparam)
 115 {
 116     int i, offset;
 117
 118     if (value_notzero_p(e->d))
 119         return L;
 120
 121     assert(e->x.p->type != partition);
 122
 123     if (e->x.p->type == fractional) {
 124         Lattice *L2 = extract_lattice(&e->x.p->arr[0], nparam);
 125         if (!L)
 126             L = L2;
 127         else {
 128             Lattice *L3 = LatticeIntersection(L, L2);
 129             Matrix_Free(L);
 130             Matrix_Free(L2);
 131             L = L3;
 132         }
 133     }
 134
 135     offset = type_offset(e->x.p);
 136     for (i = e->x.p->size-1; i >= offset; --i)
 137         L = extract_common_lattice(&e->x.p->arr[i], L, nparam);
 138     return L;
 139 }
 140
 141 /* Construct an evalue dst from src corresponding to the coset represented
 142  * by coset, a vector of size number of parameters plus one.
 143  * The final value in this vector should be 1.
 144  */
 145 static void evalue_coset(const evalue *src, const Vector *coset, evalue *dst)
 146 {
 147     if (value_notzero_p(src->d)) {
 148         value_assign(dst->d, src->d);
 149         value_init(dst->x.n);
 150         value_assign(dst->x.n, src->x.n);
 151         return;
 152     }
 153
 154     if (src->x.p->type == fractional) {
 155         evalue f;
 156         evalue t;
 157         Vector *c = Vector_Alloc(coset->Size);
 158         value_init(f.d);
 159         value_init(f.x.n);
 160         evalue_extract_affine(&src->x.p->arr[0], c->p, c->p+c->Size-1, &f.d);
 161         Inner_Product(coset->p, c->p, c->Size, &f.x.n);
 162         Vector_Free(c);
 163         mpz_fdiv_r(f.x.n, f.x.n, f.d);
 164
 165         evalue_set_si(dst, 0, 1);
 166         for (int i = src->x.p->size-1; i >= 1; --i) {
 167             emul(&f, dst);
 168             value_init(t.d);
 169             evalue_coset(&src->x.p->arr[i], coset, &t);
 170             eadd(&t, dst);
 171             free_evalue_refs(&t);
 172         }
 173         free_evalue_refs(&f);
 174         return;
 175     }
 176
 177     if (src->x.p->type == relation) {
 178         evalue *arg = evalue_eval(&src->x.p->arr[0], coset->p);
 179         if (value_zero_p(arg->x.n))
 180             evalue_coset(&src->x.p->arr[1], coset, dst);
 181         else if (src->x.p->size > 2)
 182             evalue_coset(&src->x.p->arr[2], coset, dst);
 183         else
 184             evalue_set_si(dst, 0, 1);
 185         evalue_free(arg);
 186         return;
 187     }
 188
 189     assert(src->x.p->type == polynomial);
 190     value_set_si(dst->d, 0);
 191     dst->x.p = new_enode(src->x.p->type, src->x.p->size, src->x.p->pos);
 192     for (int i = 0; i < src->x.p->size; ++i)
 193         evalue_coset(&src->x.p->arr[i], coset, &dst->x.p->arr[i]);
 194 }
 195
 196 #ifndef USE_FDSTREAM
 197 static void evalue_print_list_evalue(FILE *out, evalue *e, int nparam,
 198                                      const char **params)
 199 {
 200         cerr << "not supported" << endl;
 201 }
 202 #else
 203 static void evalue_print_list_evalue(FILE *out, evalue *e, int nparam,
 204                                      const char **params)
 205 {
 206     Lattice *L;
 207     L = extract_common_lattice(e, NULL, nparam);
 208     if (!L)
 209         print_evalue(out, e, params);
 210     else {
 211         fdostream os(dup(fileno(out)));
 212         Vector *coset = Vector_Alloc(nparam+1);
 213         value_set_si(coset->p[nparam], 1);
 214         mat_ZZ R;
 215         Matrix_Transposition(L);
 216         matrix2zz(L, R, nparam, nparam);
 217         fprintf(out, "Lattice:\n");
 218         os << R << endl;
 219         unsigned long det = to_ulong(abs(determinant(R)));
 220         mat_ZZ vertices;
 221         Matrix *points = Matrix_Alloc(det, nparam);
 222         lattice_points_fixed(coset->p, coset->p, L, L, points, det);
 223         matrix2zz(points, vertices, points->NbRows, points->NbColumns);
 224         Matrix_Free(points);
 225         Matrix_Free(L);
 226         for (int i = 0; i < vertices.NumRows(); ++i) {
 227             evalue t;
 228             os << vertices[i] << endl;
 229             zz2values(vertices[i], coset->p);
 230             value_init(t.d);
 231             evalue_coset(e, coset, &t);
 232             print_evalue(out, &t, params);
 233             free_evalue_refs(&t);
 234         }
 235         Vector_Free(coset);
 236     }
 237 }
 238 #endif
 239
 240 static void evalue_print_list(FILE *out, evalue *e, int nparam,
 241                                 const char **params)
 242 {
 243     int i;
 244     assert(value_zero_p(e->d));
 245     assert(e->x.p->type == partition);
 246
 247     for (i = 0; i < e->x.p->size/2; i++) {
 248         Print_Domain(out, EVALUE_DOMAIN(e->x.p->arr[2*i]), params);
 249         evalue_print_list_evalue(out, &e->x.p->arr[2*i+1], nparam, params);
 250     }
 251 }
 252
 253 static void print_domain_latex(std::ostream& o, Polyhedron *D, int nparam,
 254                                const char **params)
 255 {
 256     int fr = 1;
 257     for (int i = 0; i < D->NbConstraints; ++i) {
 258         if (First_Non_Zero(D->Constraint[i]+1, D->Dimension) == -1)
 259             continue;
 260         int fc = 1;
 261         if (!fr)
 262             o << " \\wedge\n";
 263         for (int j = 0; j < D->Dimension; ++j) {
 264             if (value_zero_p(D->Constraint[i][1+j]))
 265                 continue;
 266             o << " ";
 267             if (!fc && value_pos_p(D->Constraint[i][1+j]))
 268                 o << "+";
 269             if (value_mone_p(D->Constraint[i][1+j]))
 270                 o << "-";
 271             else if (!value_one_p(D->Constraint[i][1+j]))
 272                 o << VALUE_TO_INT(D->Constraint[i][1+j]);
 273             o << " " << params[j];
 274             fc = 0;
 275         }
 276         if (!fc && value_pos_p(D->Constraint[i][1+D->Dimension]))
 277             o << "+";
 278         if (value_notzero_p(D->Constraint[i][1+D->Dimension]))
 279             o << VALUE_TO_INT(D->Constraint[i][1+D->Dimension]);
 280         if (value_zero_p(D->Constraint[i][0]))
 281             o << " = 0";
 282         else
 283             o << " \\ge 0";
 284         fr = 0;
 285     }
 286 }
 287
 288 static void evalue_print_latex(std::ostream& o, const evalue *e,
 289                                int first, int nested,
 290                                const string& suffix1, const string& suffix2,
 291                                int nparam, const char **params);
 292
 293 static void evalue_print_poly_latex1(std::ostream& o, const evalue *e,
 294                                      int first, int nested, const string& base,
 295                                      const string& suffix1, const string& suffix2,
 296                                      int nparam, const char **params)
 297 {
 298     int offset = type_offset(e->x.p);
 299     for (int i = e->x.p->size-1; i >= offset; --i) {
 300         std::ostringstream strm;
 301         strm << suffix1;
 302         if (i-offset)
 303             strm << " " << base;
 304         if (i-offset > 1)
 305             strm << "^" << (i-offset);
 306         evalue_print_latex(o, &e->x.p->arr[i], first, nested,
 307                            strm.str(), suffix2, nparam, params);
 308         first = 0;
 309     }
 310 }
 311
 312 static void evalue_print_poly_latex2(std::ostream& o, const evalue *e,
 313                                      int first, int nested, const string& base,
 314                                      const string& suffix1, const string& suffix2,
 315                                      int nparam, const char **params)
 316 {
 317     int offset = type_offset(e->x.p);
 318     for (int i = e->x.p->size-1; i >= offset; --i) {
 319         std::ostringstream strm;
 320         strm << suffix2;
 321         if (i-offset)
 322             strm << " " << base;
 323         if (i-offset > 1)
 324             strm << "^" << (i-offset);
 325         evalue_print_latex(o, &e->x.p->arr[i], first, nested,
 326                            suffix1, strm.str(), nparam, params);
 327         first = 0;
 328     }
 329 }
 330
 331 static void evalue_print_latex(std::ostream& o, const evalue *e,
 332                                int first, int nested,
 333                                const string& suffix1, const string &suffix2,
 334                                int nparam, const char **params)
 335 {
 336     if (value_notzero_p(e->d)) {
 337         if (value_zero_p(e->x.n)) {
 338             if (first)
 339                 o << "0";
 340             return;
 341         }
 342         Value tmp;
 343         value_init(tmp);
 344         value_absolute(tmp, e->x.n);
 345         if (!first && value_pos_p(e->x.n))
 346             o << " + ";
 347         if (value_neg_p(e->x.n))
 348             o << " - ";
 349         if (value_one_p(e->d)) {
 350             if (!value_one_p(tmp) ||
 351                     (suffix1.length() == 0 && suffix2.length() == 0))
 352                 o << VALUE_TO_INT(tmp);
 353         } else {
 354             o << "\\frac{";
 355             if (value_one_p(tmp) && suffix1.length() != 0)
 356                 o << suffix1;
 357             else
 358                 o << VALUE_TO_INT(tmp);
 359             o << "}{"
 360               << VALUE_TO_INT(e->d) << "}";
 361         }
 362         if (!value_one_p(tmp)) {
 363             o << suffix1;
 364             o << " ";
 365         }
 366         value_clear(tmp);
 367         o << suffix2;
 368         if (!nested)
 369             o << endl;
 370         return;
 371     }
 372     switch (e->x.p->type) {
 373     case partition:
 374         o << "\\begin{cases}\n";
 375         for (int i = 0; i < e->x.p->size/2; ++i) {
 376             if (i)
 377                 o << "\\\\\n";
 378             evalue_print_latex(o, &e->x.p->arr[2*i+1], 1, 0,
 379                                 suffix1, suffix2, nparam, params);
 380             o << "& \\text{if $";
 381             print_domain_latex(o, EVALUE_DOMAIN(e->x.p->arr[2*i]), nparam, params);
 382             o << "$}\n";
 383         }
 384         o << "\\end{cases}\n";
 385         break;
 386     case polynomial:
 387         evalue_print_poly_latex1(o, e, first, nested, params[e->x.p->pos-1],
 388                                 suffix1, suffix2, nparam, params);
 389         break;
 390     case fractional: {
 391         std::ostringstream strm;
 392         strm << "\\fractional{";
 393         evalue_print_latex(strm, &e->x.p->arr[0], 1, 1, "", "", nparam, params);
 394         strm << "}";
 395         evalue_print_poly_latex2(o, e, first, nested,
 396                                 strm.str(), suffix1, suffix2, nparam, params);
 397         break;
 398     }
 399     default:
 400         assert(0);
 401     }
 402 }
 403
 404 #ifndef USE_FDSTREAM
 405 static void evalue_print_latex(FILE *out, const evalue *e, int nparam,
 406                                const char **params)
 407 {
 408         cerr << "not supported" << endl;
 409 }
 410 #else
 411 static void evalue_print_latex(FILE *out, const evalue *e, int nparam,
 412                                const char **params)
 413 {
 414     fdostream os(dup(fileno(out)));
 415     evalue_print_latex(os, e, 1, 0, "", "", nparam, params);
 416 }
 417 #endif
 418
 419 static void evalue_print_isl(FILE *out, const evalue *e, int nparam,
 420                                const char **params)
 421 {
 422         int i;
 423         isl_ctx *ctx = isl_ctx_alloc();
 424         isl_dim *dim = isl_dim_set_alloc(ctx, nparam, 0);
 425         isl_pw_qpolynomial *pwqp;
 426
 427         for (i = 0; i < nparam; ++i)
 428                 dim = isl_dim_set_name(dim, isl_dim_param, i, params[i]);
 429
 430         pwqp = isl_pw_qpolynomial_from_evalue(dim, e);
 431
 432         isl_pw_qpolynomial_print(pwqp, out, ISL_FORMAT_ISL);
 433
 434         isl_pw_qpolynomial_free(pwqp);
 435
 436         isl_ctx_free(ctx);
 437 }
 438
 439 int evalue_convert(evalue *EP, struct convert_options *options,
 440                    int verbose, unsigned nparam, const char **params)
 441 {
 442     int printed = 0;
 443     if (options->combine)
 444         evalue_combine(EP);
 445     if (options->range)
 446         evalue_range_reduction(EP);
 447     if (verbose) {
 448         print_evalue(stdout, EP, params);
 449         printed = 1;
 450     }
 451     if (options->floor) {
 452         fprintf(stderr, "WARNING: floor conversion not supported\n");
 453         evalue_frac2floor(EP);
 454         if (params)
 455             print_evalue(stdout, EP, params);
 456     } else if (options->list && params) {
 457         evalue_print_list(stdout, EP, nparam, params);
 458         printed = 1;
 459     } else if (options->latex && params) {
 460         evalue_print_latex(stdout, EP, nparam, params);
 461         printed = 1;
 462     } else if (options->isl && params) {
 463         evalue_print_isl(stdout, EP, nparam, params);
 464         printed = 1;
 465     } else if (options->convert) {
 466         evalue_mod2table(EP, nparam);
 467         if (verbose)
 468             print_evalue(stdout, EP, params);
 469     }
 470     return printed;
 471 }