evalue_convert.cc

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