1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 #include <iostream>
17 #include <assert.h>
18
19 #include "soplex/spxdefines.h"
|
(1) Event include_recursion: |
#include file "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scipoptsuite-8.0.1/soplex/src/soplex/ratrecon.h" includes itself: ratrecon.h -> ratrecon.hpp -> ratrecon.h |
|
(2) Event caretline: |
^ |
20 #include "soplex/ratrecon.h"
21 #include "soplex/rational.h"
22
23 namespace soplex
24 {
25
26 /** this reconstruction routine will set x equal to the mpq vector where each component is the best rational
27 * approximation of xnum / denom with where the GCD of denominators of x is at most Dbound; it will return true on
28 * success and false if more accuracy is required: specifically if componentwise rational reconstruction does not
29 * produce such a vector
30 */
31 static int Reconstruct(VectorRational& resvec, Integer* xnum, Integer denom, int dim,
32 const Rational& denomBoundSquared, const DIdxSet* indexSet = 0)
33 {
34 bool rval = true;
35 int done = 0;
36
37 /* denominator must be positive */
38 assert(denom > 0);
39 assert(denomBoundSquared > 0);
40
41 Integer temp = 0;
42 Integer td = 0;
43 Integer tn = 0;
44 Integer Dbound = 0;
45 Integer gcd = 1;
46
47 Dbound = numerator(denomBoundSquared) / denominator(
48 denomBoundSquared); /* this is the working bound on the denominator size */
49
50 Dbound = (Integer) sqrt(Dbound);
51
52 MSG_DEBUG(std::cout << "reconstructing " << dim << " dimensional vector with denominator bound " <<
53 Dbound << "\n");
54
55 /* if Dbound is below 2^24 increase it to this value, this avoids changing input vectors that have low denominator
56 * because they are floating point representable
57 */
58 if(Dbound < 16777216)
59 Dbound = 16777216;
60
61 /* The following represent a_i, the cont frac representation and p_i/q_i, the convergents */
62 Integer a0 = 0;
63 Integer ai = 0;
64
65 /* here we use p[2]=pk, p[1]=pk-1,p[0]=pk-2 and same for q */
66 Integer p[3];
67 Integer q[3];
68
69 for(int c = 0; (indexSet == 0 && c < dim) || (indexSet != 0 && c < indexSet->size()); c++)
70 {
71 int j = (indexSet == 0 ? c : indexSet->index(c));
72
73 assert(j >= 0);
74 assert(j < dim);
75
76 MSG_DEBUG(std::cout << " --> component " << j << " = " << &xnum[j] << " / denom\n");
77
78 /* if xnum =0 , then just leave x[j] as zero */
79 if(xnum[j] != 0)
80 {
81 /* setup n and d for computing a_i the cont. frac. rep */
82 tn = xnum[j];
83 td = denom;
84
85 /* divide tn and td by gcd */
86 SpxGcd(temp, tn, td);
87 tn = tn / temp;
88 td = td / temp;
89
90 if(td <= Dbound)
91 {
92 MSG_DEBUG(std::cout << "marker 1\n");
93
94 resvec[j] = Rational(tn, td);
95 }
96 else
97 {
98 MSG_DEBUG(std::cout << "marker 2\n");
99
100 temp = 1;
101
102 divide_qr(tn, td, a0, temp);
103
104 tn = td;
105 td = temp;
106
107 divide_qr(tn, td, ai, temp);
108 tn = td;
109 td = temp;
110
111 p[1] = a0;
112 p[2] = 1;
113 p[2] += a0 * ai;
114
115 q[1] = 1;
116 q[2] = ai;
117
118 done = 0;
119
120 /* if q is already big, skip loop */
121 if(q[2] > Dbound)
122 {
123 MSG_DEBUG(std::cout << "marker 3\n");
124 done = 1;
125 }
126
127 int cfcnt = 2;
128
129 while(!done && td != 0)
130 {
131 /* update everything: compute next ai, then update convergents */
132
133 /* update ai */
134 divide_qr(tn, td, ai, temp);
135 tn = td;
136 td = temp;
137
138 /* shift p,q */
139 q[0] = q[1];
140 q[1] = q[2];
141 p[0] = p[1];
142 p[1] = p[2];
143
144 /* compute next p,q */
145 p[2] = p[0];
146 p[2] += p[1] * ai;
147 q[2] = q[0];
148 q[2] += q[1] * ai;
149
150 if(q[2] > Dbound)
151 done = 1;
152
153 cfcnt++;
154
155 MSG_DEBUG(std::cout << " --> convergent denominator = " << &q[2] << "\n");
156 }
157
158 assert(q[1] != 0);
159
160 /* Assign the values */
161 if(q[1] >= 0)
162 resvec[j] = Rational(p[1], q[1]);
163 else
164 resvec[j] = Rational(-p[1], -q[1]);
165
166 SpxGcd(temp, gcd, denominator(resvec[j]));
167 gcd *= temp;
168
169 if(gcd > Dbound)
170 {
171 MSG_DEBUG(std::cout << "terminating with gcd " << &gcd << " exceeding Dbound " << &Dbound << "\n");
172 rval = false;
173 break;
174 }
175 }
176 }
177 }
178
179 return rval;
180 }
181
182 /** reconstruct a rational vector */
183 inline bool reconstructVector(VectorRational& input, const Rational& denomBoundSquared,
184 const DIdxSet* indexSet)
185 {
186 std::vector<Integer> xnum(input.dim()); /* numerator of input vector */
187 Integer denom = 1; /* common denominator of input vector */
188 int rval = true;
189 int dim;
190
191 dim = input.dim();
192
193 /* find common denominator */
194 if(indexSet == 0)
195 {
196 for(int i = 0; i < dim; i++)
197 SpxLcm(denom, denom, denominator(input[i]));
198
199 for(int i = 0; i < dim; i++)
200 {
201 xnum[i] = denom * Integer(numerator(input[i]));
202 xnum[i] = xnum[i] / Integer(denominator(input[i]));
203 }
204 }
205 else
206 {
207 for(int i = 0; i < indexSet->size(); i++)
208 {
209 assert(indexSet->index(i) >= 0);
210 assert(indexSet->index(i) < input.dim());
211 SpxLcm(denom, denom, denominator(input[indexSet->index(i)]));
212 }
213
214 for(int i = 0; i < indexSet->size(); i++)
215 {
216 int k = indexSet->index(i);
217 assert(k >= 0);
218 assert(k < input.dim());
219 xnum[k] = denom * Integer(numerator(input[k]));
220 xnum[k] = xnum[k] / Integer(denominator(input[k]));
221 }
222 }
223
224 MSG_DEBUG(std::cout << "LCM = " << mpz_get_str(0, 10, denom) << "\n");
225
226 /* reconstruct */
227 rval = Reconstruct(input, xnum.data(), denom, dim, denomBoundSquared, indexSet);
228
229 return rval;
230 }
231
232
233
234 /** reconstruct a rational solution */
235 /**@todo make this a method of class SoPlex */
236 inline bool reconstructSol(SolRational& solution)
237 {
238 #if 0
239 VectorRational buffer;
240
241 if(solution.hasPrimal())
242 {
243 buffer.reDim((solution._primal).dim());
244 solution.getPrimalSol(buffer);
245 reconstructVector(buffer);
246 solution._primal = buffer;
247
248 buffer.reDim((solution._slacks).dim());
249 solution.getSlacks(buffer);
250 reconstructVector(buffer);
251 solution._slacks = buffer;
252 }
253
254 if(solution.hasPrimalray())
255 {
256 buffer.reDim((solution._primalray).dim());
257 solution.getPrimalray(buffer);
258 reconstructVector(buffer);
259 solution._primalray = buffer;
260 }
261
262 if(solution.hasDual())
263 {
264 buffer.reDim((solution._dual).dim());
265 solution.getDualSol(buffer);
266 reconstructVector(buffer);
267 solution._dual = buffer;
268
269 buffer.reDim((solution._redcost).dim());
270 solution.getRedcost(buffer);
271 reconstructVector(buffer);
272 solution._redcost = buffer;
273 }
274
275 if(solution.hasDualfarkas())
276 {
277 buffer.reDim((solution._dualfarkas).dim());
278 solution.getDualfarkas(buffer);
279 reconstructVector(buffer);
280 solution._dualfarkas = buffer;
281 }
282
283 #endif
284 return true;
285 }
286 } // namespace soplex
287