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                                      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);
 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_frac2floor(EP);
 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 }