evalue_convert.cc

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