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