evalue_convert.cc

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