evalue_convert.cc

   1 #include "conversion.h"
   2 #include "evalue_convert.h"
   3 #include "lattice_point.h"
   4 #include "fdstream.h"
   5
   6 using std::cout;
   7 using std::cerr;
   8 using std::endl;
   9
  10 static struct argp_option argp_options[] = {
  11     { "convert",            'c', 0, 0, "convert fractionals to periodics" },
  12     { "combine",            'C', 0, 0 },
  13     { "floor",              'f', 0, 0, "convert fractionals to floorings" },
  14     { "list",               'l', 0, 0 },
  15     { "range-reduction",    'R',    0,      0 },
  16     0
  17 };
  18
  19 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  20 {
  21     struct convert_options *options = (struct convert_options *)state->input;
  22
  23     switch (key) {
  24     case ARGP_KEY_INIT:
  25         options->floor = 0;
  26         options->convert = 0;
  27         options->combine = 0;
  28         options->range = 0;
  29         options->list = 0;
  30         break;
  31     case ARGP_KEY_FINI:
  32         break;
  33     case 'f':
  34         options->floor = 1;
  35         break;
  36     case 'c':
  37         options->convert = 1;
  38         break;
  39     case 'C':
  40         options->combine = 1;
  41         break;
  42     case 'l':
  43         options->list = 1;
  44         break;
  45     case 'R':
  46         options->range = 1;
  47         break;
  48     default:
  49         return ARGP_ERR_UNKNOWN;
  50     }
  51     return 0;
  52 }
  53
  54 struct argp convert_argp = {
  55     argp_options, parse_opt, 0, 0
  56 };
  57
  58 static int type_offset(enode *p)
  59 {
  60    return p->type == fractional ? 1 :
  61           p->type == flooring ? 1 : 0;
  62 }
  63
  64 static Lattice *extract_lattice(evalue *e, int nparam)
  65 {
  66     int i;
  67     Lattice *L;
  68     Matrix *U;
  69     Vector *X;
  70     /* For some mysterious reason, SolveDiophantine expects an extra
  71      * [0 0 0 1] row in its input matrix.
  72      */
  73     Matrix *M = Matrix_Alloc(2, nparam+1+1);
  74     value_set_si(M->p[1][nparam+1], 1);
  75     evalue_extract_affine(e, M->p[0], M->p[0]+nparam+1, M->p[0]+nparam);
  76     /* ignore constant */
  77     value_set_si(M->p[0][nparam+1], 0);
  78     SolveDiophantine(M, &U, &X);
  79     Matrix_Free(M);
  80     Vector_Free(X);
  81     L = Matrix_Alloc(nparam+1, nparam+1);
  82     for (i = 0; i < nparam; ++i)
  83         Vector_Copy(U->p[i], L->p[i], nparam);
  84     value_set_si(L->p[nparam][nparam], 1);
  85     Matrix_Free(U);
  86     return L;
  87 }
  88
  89 /* Returns a lattice such that the quasi-polynomial e can be represented
  90  * by a list of polynomials, one for each point in the fundamental
  91  * parallelepiped of the lattice.
  92  * If e is a polynomial, then this function returns NULL.
  93  */
  94 static Lattice *extract_common_lattice(evalue *e, Lattice *L, int nparam)
  95 {
  96     int i, offset;
  97
  98     if (value_notzero_p(e->d))
  99         return L;
 100
 101     assert(e->x.p->type != partition);
 102
 103     if (e->x.p->type == fractional) {
 104         Lattice *L2 = extract_lattice(&e->x.p->arr[0], nparam);
 105         if (!L)
 106             L = L2;
 107         else {
 108             Lattice *L3 = LatticeIntersection(L, L2);
 109             Matrix_Free(L);
 110             Matrix_Free(L2);
 111             L = L3;
 112         }
 113     }
 114
 115     offset = type_offset(e->x.p);
 116     for (i = e->x.p->size-1; i >= offset; --i)
 117         L = extract_common_lattice(&e->x.p->arr[i], L, nparam);
 118     return L;
 119 }
 120
 121 /* Construct an evalue dst from src corresponding to the coset represented
 122  * by coset, a vector of size number of parameters plus one.
 123  * The final value in this vector should be 1.
 124  */
 125 static void evalue_coset(const evalue *src, const Vector *coset, evalue *dst)
 126 {
 127     if (value_notzero_p(src->d)) {
 128         value_assign(dst->d, src->d);
 129         value_init(dst->x.n);
 130         value_assign(dst->x.n, src->x.n);
 131         return;
 132     }
 133
 134     if (src->x.p->type == fractional) {
 135         evalue f;
 136         evalue t;
 137         Vector *c = Vector_Alloc(coset->Size);
 138         value_init(f.d);
 139         value_init(f.x.n);
 140         evalue_extract_affine(&src->x.p->arr[0], c->p, c->p+c->Size-1, &f.d);
 141         Inner_Product(coset->p, c->p, c->Size, &f.x.n);
 142         Vector_Free(c);
 143         mpz_fdiv_r(f.x.n, f.x.n, f.d);
 144
 145         evalue_set_si(dst, 0, 1);
 146         for (int i = src->x.p->size-1; i >= 1; --i) {
 147             emul(&f, dst);
 148             value_init(t.d);
 149             evalue_coset(&src->x.p->arr[i], coset, &t);
 150             eadd(&t, dst);
 151             free_evalue_refs(&t);
 152         }
 153         free_evalue_refs(&f);
 154         return;
 155     }
 156
 157     assert(src->x.p->type == polynomial);
 158     value_set_si(dst->d, 0);
 159     dst->x.p = new_enode(src->x.p->type, src->x.p->size, src->x.p->pos);
 160     for (int i = 0; i < src->x.p->size; ++i)
 161         evalue_coset(&src->x.p->arr[i], coset, &dst->x.p->arr[i]);
 162 }
 163
 164 static void evalue_print_list_evalue(FILE *out, evalue *e, int nparam,
 165                                      char **params)
 166 {
 167     Lattice *L;
 168     L = extract_common_lattice(e, NULL, nparam);
 169     if (!L)
 170         print_evalue(out, e, params);
 171     else {
 172         fdostream os(dup(fileno(out)));
 173         Vector *coset = Vector_Alloc(nparam+1);
 174         value_set_si(coset->p[nparam], 1);
 175         mat_ZZ RT;
 176         mat_ZZ R;
 177         matrix2zz(L, RT, nparam, nparam);
 178         R = transpose(RT);
 179         fprintf(out, "Lattice:\n");
 180         os << R << endl;
 181         mat_ZZ vertices;
 182         lattice_point(coset->p, R, vertices, to_ulong(abs(determinant(R))), NULL);
 183         Matrix_Free(L);
 184         for (int i = 0; i < vertices.NumRows(); ++i) {
 185             evalue t;
 186             os << vertices[i] << endl;
 187             zz2values(vertices[i], coset->p);
 188             value_init(t.d);
 189             evalue_coset(e, coset, &t);
 190             print_evalue(out, &t, params);
 191             free_evalue_refs(&t);
 192         }
 193         Vector_Free(coset);
 194     }
 195 }
 196
 197 static void evalue_print_list(FILE *out, evalue *e, int nparam, char **params)
 198 {
 199     int i;
 200     assert(value_zero_p(e->d));
 201     assert(e->x.p->type == partition);
 202
 203     for (i = 0; i < e->x.p->size/2; i++) {
 204         Print_Domain(out, EVALUE_DOMAIN(e->x.p->arr[2*i]), params);
 205         evalue_print_list_evalue(out, &e->x.p->arr[2*i+1], nparam, params);
 206     }
 207 }
 208
 209 int evalue_convert(evalue *EP, struct convert_options *options,
 210                    int verbose, unsigned nparam, char **params)
 211 {
 212     int printed = 0;
 213     if (options->combine)
 214         evalue_combine(EP);
 215     if (options->range)
 216         evalue_range_reduction(EP);
 217     if (verbose) {
 218         print_evalue(stdout, EP, params);
 219         printed = 1;
 220     }
 221     if (options->floor) {
 222         fprintf(stderr, "WARNING: floor conversion not supported\n");
 223         evalue_frac2floor2(EP, 0);
 224         if (params)
 225             print_evalue(stdout, EP, params);
 226     } else if (options->list && params) {
 227         evalue_print_list(stdout, EP, nparam, params);
 228         printed = 1;
 229     } else if (options->convert) {
 230         evalue_mod2table(EP, nparam);
 231         if (verbose)
 232             print_evalue(stdout, EP, params);
 233     }
 234     return printed;
 235 }