evalue_convert.cc

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