evalue_convert.cc

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