conversion.cc

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