evalue_convert.cc

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