evalue_convert.cc

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