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 *zz2matrix(const mat_ZZ& mat)
  74 {
  75     Matrix *M = Matrix_Alloc(mat.NumRows(), mat.NumCols());
  76     assert(M);
  77
  78     for (int i = 0; i < mat.NumRows(); ++i)
  79         zz2values(mat[i], M->p[i]);
  80     return M;
  81 }
  82
  83 Matrix *rays(Polyhedron *C)
  84 {
  85     unsigned dim = C->NbRays - 1; /* don't count zero vertex */
  86     assert(C->NbRays - 1 == C->Dimension);
  87
  88     Matrix *M = Matrix_Alloc(dim+1, dim+1);
  89     assert(M);
  90
  91     int i, c;
  92     for (i = 0, c = 0; i <= dim && c < dim; ++i)
  93         if (value_zero_p(C->Ray[i][dim+1])) {
  94             Vector_Copy(C->Ray[i] + 1, M->p[c], dim);
  95             value_set_si(M->p[c++][dim], 0);
  96         }
  97     assert(c == dim);
  98     value_set_si(M->p[dim][dim], 1);
  99
 100     return M;
 101 }
 102
 103 Matrix *rays2(Polyhedron *C)
 104 {
 105     unsigned dim = C->NbRays - 1; /* don't count zero vertex */
 106     assert(C->NbRays - 1 == C->Dimension);
 107
 108     Matrix *M = Matrix_Alloc(dim, dim);
 109     assert(M);
 110
 111     int i, c;
 112     for (i = 0, c = 0; i <= dim && c < dim; ++i)
 113         if (value_zero_p(C->Ray[i][dim+1]))
 114             Vector_Copy(C->Ray[i] + 1, M->p[c++], dim);
 115     assert(c == dim);
 116
 117     return M;
 118 }
 119
 120 void rays(Polyhedron *C, mat_ZZ& rays)
 121 {
 122     unsigned dim = C->NbRays - 1; /* don't count zero vertex */
 123     assert(C->NbRays - 1 == C->Dimension);
 124     rays.SetDims(dim, dim);
 125
 126     int i, j;
 127     for (i = 0, j = 0; i < C->NbRays; ++i) {
 128         if (value_notzero_p(C->Ray[i][dim+1]))
 129             continue;
 130         values2zz(C->Ray[i]+1, rays[j], dim);
 131         ++j;
 132     }
 133 }
 134
 135 void randomvector(Polyhedron *P, vec_ZZ& lambda, int nvar, int n_try)
 136 {
 137     Value tmp;
 138     int max = 5 * 16;
 139     unsigned int dim = P->Dimension;
 140     value_init(tmp);
 141
 142     for (int i = 0; i < P->NbRays; ++i) {
 143         for (int j = 1; j <= dim; ++j) {
 144             value_absolute(tmp, P->Ray[i][j]);
 145             int t = VALUE_TO_LONG(tmp) * 16;
 146             if (t > max)
 147                 max = t;
 148         }
 149     }
 150     for (int i = 0; i < P->NbConstraints; ++i) {
 151         for (int j = 1; j <= dim; ++j) {
 152             value_absolute(tmp, P->Constraint[i][j]);
 153             int t = VALUE_TO_LONG(tmp) * 16;
 154             if (t > max)
 155                 max = t;
 156         }
 157     }
 158     value_clear(tmp);
 159
 160     max += max << n_try;
 161
 162     lambda.SetLength(nvar);
 163     for (int k = 0; k < nvar; ++k) {
 164         int r = random_int(max*dim)+2;
 165         int v = (2*(r%2)-1) * (max/2*dim + (r >> 1));
 166         lambda[k] = v;
 167     }
 168 }