conversion.cc

   1 #include <assert.h>
   2 #include <gmp.h>
   3 #include <NTL/mat_ZZ.h>
   4 #include <barvinok/util.h>
   5 #include "conversion.h"
   6
   7 #define SIZE(p) (((long *) (p))[1])
   8 #define DATA(p) ((mp_limb_t *) (((long *) (p)) + 2))
   9
  10 void value2zz(Value v, ZZ& z)
  11 {
  12     int sa = v[0]._mp_size;
  13     int abs_sa = sa < 0 ? -sa : sa;
  14
  15     _ntl_gsetlength(&z.rep, abs_sa);
  16     mp_limb_t * adata = DATA(z.rep);
  17     for (int i = 0; i < abs_sa; ++i)
  18         adata[i] = v[0]._mp_d[i];
  19     SIZE(z.rep) = sa;
  20 }
  21
  22 void zz2value(const ZZ& z, Value& v)
  23 {
  24     if (!z.rep) {
  25         value_set_si(v, 0);
  26         return;
  27     }
  28
  29     int sa = SIZE(z.rep);
  30     int abs_sa = sa < 0 ? -sa : sa;
  31
  32     mp_limb_t * adata = DATA(z.rep);
  33     _mpz_realloc(v, abs_sa);
  34     for (int i = 0; i < abs_sa; ++i)
  35         v[0]._mp_d[i] = adata[i];
  36     v[0]._mp_size = sa;
  37 }
  38
  39 void values2zz(Value *p, vec_ZZ& v, int len)
  40 {
  41     v.SetLength(len);
  42
  43     for (int i = 0; i < len; ++i) {
  44         value2zz(p[i], v[i]);
  45     }
  46 }
  47
  48 /*
  49  */
  50 void zz2values(const vec_ZZ& v, Value *p)
  51 {
  52     for (int i = 0; i < v.length(); ++i)
  53         zz2value(v[i], p[i]);
  54 }
  55
  56 /*
  57  * We just ignore the last column and row
  58  * If the final element is not equal to one
  59  * then the result will actually be a multiple of the input
  60  */
  61 void matrix2zz(Matrix *M, mat_ZZ& m, unsigned nr, unsigned nc)
  62 {
  63     m.SetDims(nr, nc);
  64
  65     for (int i = 0; i < nr; ++i) {
  66 //      assert(value_one_p(M->p[i][M->NbColumns - 1]));
  67         for (int j = 0; j < nc; ++j) {
  68             value2zz(M->p[i][j], m[i][j]);
  69         }
  70     }
  71 }
  72
  73 Matrix *rays2matrix(const mat_ZZ& rays)
  74 {
  75     unsigned dim = rays.NumRows();
  76
  77     Matrix *M = Matrix_Alloc(dim+1, dim+1);
  78     assert(M);
  79
  80     for (int i = 0; i < dim; ++i)
  81         zz2values(rays[i], M->p[i]);
  82     value_set_si(M->p[dim][dim], 1);
  83     return M;
  84 }
  85
  86 Matrix *rays2matrix2(const mat_ZZ& rays)
  87 {
  88     unsigned dim = rays.NumRows();
  89
  90     Matrix *M = Matrix_Alloc(dim, dim);
  91     assert(M);
  92
  93     for (int i = 0; i < dim; ++i)
  94         zz2values(rays[i], M->p[i]);
  95     return M;
  96 }
  97
  98 Matrix *rays(Polyhedron *C)
  99 {
 100     unsigned dim = C->NbRays - 1; /* don't count zero vertex */
 101     assert(C->NbRays - 1 == C->Dimension);
 102
 103     Matrix *M = Matrix_Alloc(dim+1, dim+1);
 104     assert(M);
 105
 106     int i, c;
 107     for (i = 0, c = 0; i <= dim && c < dim; ++i)
 108         if (value_zero_p(C->Ray[i][dim+1])) {
 109             Vector_Copy(C->Ray[i] + 1, M->p[c], dim);
 110             value_set_si(M->p[c++][dim], 0);
 111         }
 112     assert(c == dim);
 113     value_set_si(M->p[dim][dim], 1);
 114
 115     return M;
 116 }
 117
 118 Matrix *rays2(Polyhedron *C)
 119 {
 120     unsigned dim = C->NbRays - 1; /* don't count zero vertex */
 121     assert(C->NbRays - 1 == C->Dimension);
 122
 123     Matrix *M = Matrix_Alloc(dim, dim);
 124     assert(M);
 125
 126     int i, c;
 127     for (i = 0, c = 0; i <= dim && c < dim; ++i)
 128         if (value_zero_p(C->Ray[i][dim+1]))
 129             Vector_Copy(C->Ray[i] + 1, M->p[c++], dim);
 130     assert(c == dim);
 131
 132     return M;
 133 }
 134
 135 void rays(Polyhedron *C, mat_ZZ& rays)
 136 {
 137     unsigned dim = C->NbRays - 1; /* don't count zero vertex */
 138     assert(C->NbRays - 1 == C->Dimension);
 139     rays.SetDims(dim, dim);
 140
 141     int i, j;
 142     for (i = 0, j = 0; i < C->NbRays; ++i) {
 143         if (value_notzero_p(C->Ray[i][dim+1]))
 144             continue;
 145         values2zz(C->Ray[i]+1, rays[j], dim);
 146         ++j;
 147     }
 148 }
 149
 150 void randomvector(Polyhedron *P, vec_ZZ& lambda, int nvar)
 151 {
 152     Value tmp;
 153     int max = 10 * 16;
 154     unsigned int dim = P->Dimension;
 155     value_init(tmp);
 156
 157     for (int i = 0; i < P->NbRays; ++i) {
 158         for (int j = 1; j <= dim; ++j) {
 159             value_absolute(tmp, P->Ray[i][j]);
 160             int t = VALUE_TO_LONG(tmp) * 16;
 161             if (t > max)
 162                 max = t;
 163         }
 164     }
 165     for (int i = 0; i < P->NbConstraints; ++i) {
 166         for (int j = 1; j <= dim; ++j) {
 167             value_absolute(tmp, P->Constraint[i][j]);
 168             int t = VALUE_TO_LONG(tmp) * 16;
 169             if (t > max)
 170                 max = t;
 171         }
 172     }
 173     value_clear(tmp);
 174
 175     lambda.SetLength(nvar);
 176     for (int k = 0; k < nvar; ++k) {
 177         int r = random_int(max*dim)+2;
 178         int v = (2*(r%2)-1) * (max/2*dim + (r >> 1));
 179         lambda[k] = v;
 180     }
 181 }