evalue_convert.cc

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