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