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