24 #ifdef SOPLEX_WITH_GMP 27 for(
int i = 0; i < dim; i++)
33 for(
int i = 0; i < dim; i++)
39 static std::ostream&
operator<<(std::ostream& os,
const mpz_t* number)
41 os << mpz_get_str(0, 10, *number);
58 assert(mpz_sgn(denom) > 0);
59 assert(denomBoundSquared > 0);
78 mpz_set(Dbound, denom);
81 mpz_sqrt(Dbound, Dbound);
83 MSG_DEBUG(std::cout <<
"reconstructing " << dim <<
" dimensional vector with denominator bound " <<
84 mpz_get_str(0, 10, Dbound) <<
"\n");
89 if(mpz_cmp_ui(Dbound, 16777216) < 0)
90 mpz_set_ui(Dbound, 16777216);
104 for(
int c = 0; (indexSet == 0 && c < dim) || (indexSet != 0 && c < indexSet->size()); c++)
106 int j = (indexSet == 0 ? c : indexSet->index(c));
111 MSG_DEBUG(std::cout <<
" --> component " << j <<
" = " << &xnum[j] <<
" / denom\n");
114 if(mpz_sgn(xnum[j]) != 0)
117 mpz_set(tn, xnum[j]);
121 mpz_gcd(temp, tn, td);
122 mpz_divexact(tn, tn, temp);
123 mpz_divexact(td, td, temp);
125 if(mpz_cmp(td, Dbound) <= 0)
129 mpq_set_num(resvec[j].getMpqRef_w(), tn);
130 mpq_set_den(resvec[j].getMpqRef_w(), td);
138 mpz_fdiv_q(a0, tn, td);
139 mpz_fdiv_r(temp, tn, td);
142 mpz_fdiv_q(ai, tn, td);
143 mpz_fdiv_r(temp, tn, td);
149 mpz_addmul(p[2], a0, ai);
157 if(mpz_cmp(q[2], Dbound) > 0)
165 while(!done && mpz_cmp_ui(td, 0))
170 mpz_fdiv_q(ai, tn, td);
171 mpz_fdiv_r(temp, tn, td);
183 mpz_addmul(p[2], p[1], ai);
185 mpz_addmul(q[2], q[1], ai);
187 if(mpz_cmp(q[2], Dbound) > 0)
192 MSG_DEBUG(std::cout <<
" --> convergent denominator = " << &q[2] <<
"\n");
196 mpq_set_num(resvec[j].getMpqRef_w(), p[1]);
197 mpq_set_den(resvec[j].getMpqRef_w(), q[1]);
198 mpq_canonicalize(resvec[j].getMpqRef_w());
199 mpz_gcd(temp, gcd, mpq_denref(resvec[j].getMpqRef()));
200 mpz_mul(gcd, gcd, temp);
202 if(mpz_cmp(gcd, Dbound) > 0)
204 MSG_DEBUG(std::cout <<
"terminating with gcd " << &gcd <<
" exceeding Dbound " << &Dbound <<
"\n");
233 #ifdef SOPLEX_WITH_GMP 244 mpz_init_set_ui(denom, 1);
249 for(
int i = 0; i < dim; i++)
250 mpz_lcm(denom, denom, mpq_denref(input[i].getMpqRef()));
252 for(
int i = 0; i < dim; i++)
254 mpz_mul(xnum[i], denom, mpq_numref(input[i].getMpqRef()));
255 mpz_divexact(xnum[i], xnum[i], mpq_denref(input[i].getMpqRef()));
260 for(
int i = 0; i < indexSet->
size(); i++)
262 assert(indexSet->
index(i) >= 0);
263 assert(indexSet->
index(i) < input.
dim());
264 mpz_lcm(denom, denom, mpq_denref(input[indexSet->
index(i)].getMpqRef()));
267 for(
int i = 0; i < indexSet->
size(); i++)
269 int k = indexSet->
index(i);
271 assert(k < input.
dim());
272 mpz_mul(xnum[k], denom, mpq_numref(input[k].getMpqRef()));
273 mpz_divexact(xnum[k], xnum[k], mpq_denref(input[k].getMpqRef()));
277 MSG_DEBUG(std::cout <<
"LCM = " << mpz_get_str(0, 10, denom) <<
"\n");
280 rval =
Reconstruct(input, xnum, denom, dim, denomBoundSquared, indexSet);
301 if(solution.hasPrimal())
314 if(solution.hasPrimalray())
316 buffer.
reDim((solution._primalray).dim());
317 solution.getPrimalray(buffer);
319 solution._primalray = buffer;
322 if(solution.hasDual())
327 solution.
_dual = buffer;
329 buffer.
reDim((solution._redcost).dim());
330 solution.getRedcost(buffer);
332 solution._redcost = buffer;
335 if(solution.hasDualfarkas())
337 buffer.
reDim((solution._dualfarkas).dim());
338 solution.getDualfarkas(buffer);
340 solution._dualfarkas = buffer;
bool getDual(VectorBase< R > &vector) const
gets the dual solution vector; returns true on success
void reDim(int newdim, const bool setZero=true)
Resets DVectorBase's dimension to newdim.
Dense vector.Class VectorBase provides dense linear algebra vectors. It does not provide memory manag...
bool getSlacks(VectorBase< R > &vector) const
gets the vector of slack values; returns true on success
static void GMPv_init(mpz_t *vect, int dim)
static void GMPv_clear(mpz_t *vect, int dim)
Wrapper for GMP type mpq_class.We wrap mpq_class so that we can replace it by a double type if GMP is...
bool getPrimal(VectorBase< R > &vector) const
gets the primal solution vector; returns true on success
void spx_alloc(T &p, int n=1)
Allocate memory.
Rational reconstruction of solution vector.
bool reconstructVector(VectorRational &input, const Rational &denomBoundSquared, const DIdxSet *indexSet)
static int Reconstruct(VectorRational &resvec, mpz_t *xnum, mpz_t denom, int dim, const Rational &denomBoundSquared, const DIdxSet *indexSet=0)
bool reconstructSol(SolRational &solution)
std::ostream & operator<<(std::ostream &s, const VectorBase< R > &vec)
Output operator.
int size() const
returns the number of used indices.
Dynamic index set.Class DIdxSet provides dynamic IdxSet in the sense, that no restrictions are posed ...
const mpq_t & getMpqRef() const
provides read-only access to underlying mpq_t
Debugging, floating point type and parameter definitions.
int dim() const
Dimension of vector.
Everything should be within this namespace.
int index(int n) const
access n 'th index.
void spx_free(T &p)
Release memory.