lemon-project-template-glpk

comparison deps/glpk/src/glplib02.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:e109bc17c145
1 /* glplib02.c (64-bit arithmetic) */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 *
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9 * E-mail: <mao@gnu.org>.
10 *
11 * GLPK is free software: you can redistribute it and/or modify it
12 * under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19 * License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23 ***********************************************************************/
24
25 #include "glpenv.h"
26 #include "glplib.h"
27
28 /***********************************************************************
29 * NAME
30 *
31 * xlset - expand integer to long integer
32 *
33 * SYNOPSIS
34 *
35 * #include "glplib.h"
36 * glp_long xlset(int x);
37 *
38 * RETURNS
39 *
40 * The routine xlset returns x expanded to long integer. */
41
42 glp_long xlset(int x)
43 { glp_long t;
44 t.lo = x, t.hi = (x >= 0 ? 0 : -1);
45 return t;
46 }
47
48 /***********************************************************************
49 * NAME
50 *
51 * xlneg - negate long integer
52 *
53 * SYNOPSIS
54 *
55 * #include "glplib.h"
56 * glp_long xlneg(glp_long x);
57 *
58 * RETURNS
59 *
60 * The routine xlneg returns the difference 0 - x. */
61
62 glp_long xlneg(glp_long x)
63 { if (x.lo)
64 x.lo = - x.lo, x.hi = ~x.hi;
65 else
66 x.hi = - x.hi;
67 return x;
68 }
69
70 /***********************************************************************
71 * NAME
72 *
73 * xladd - add long integers
74 *
75 * SYNOPSIS
76 *
77 * #include "glplib.h"
78 * glp_long xladd(glp_long x, glp_long y);
79 *
80 * RETURNS
81 *
82 * The routine xladd returns the sum x + y. */
83
84 glp_long xladd(glp_long x, glp_long y)
85 { if ((unsigned int)x.lo <= 0xFFFFFFFF - (unsigned int)y.lo)
86 x.lo += y.lo, x.hi += y.hi;
87 else
88 x.lo += y.lo, x.hi += y.hi + 1;
89 return x;
90 }
91
92 /***********************************************************************
93 * NAME
94 *
95 * xlsub - subtract long integers
96 *
97 * SYNOPSIS
98 *
99 * #include "glplib.h"
100 * glp_long xlsub(glp_long x, glp_long y);
101 *
102 * RETURNS
103 *
104 * The routine xlsub returns the difference x - y. */
105
106 glp_long xlsub(glp_long x, glp_long y)
107 { return
108 xladd(x, xlneg(y));
109 }
110
111 /***********************************************************************
112 * NAME
113 *
114 * xlcmp - compare long integers
115 *
116 * SYNOPSIS
117 *
118 * #include "glplib.h"
119 * int xlcmp(glp_long x, glp_long y);
120 *
121 * RETURNS
122 *
123 * The routine xlcmp returns the sign of the difference x - y. */
124
125 int xlcmp(glp_long x, glp_long y)
126 { if (x.hi >= 0 && y.hi < 0) return +1;
127 if (x.hi < 0 && y.hi >= 0) return -1;
128 if ((unsigned int)x.hi < (unsigned int)y.hi) return -1;
129 if ((unsigned int)x.hi > (unsigned int)y.hi) return +1;
130 if ((unsigned int)x.lo < (unsigned int)y.lo) return -1;
131 if ((unsigned int)x.lo > (unsigned int)y.lo) return +1;
132 return 0;
133 }
134
135 /***********************************************************************
136 * NAME
137 *
138 * xlmul - multiply long integers
139 *
140 * SYNOPSIS
141 *
142 * #include "glplib.h"
143 * glp_long xlmul(glp_long x, glp_long y);
144 *
145 * RETURNS
146 *
147 * The routine xlmul returns the product x * y. */
148
149 glp_long xlmul(glp_long x, glp_long y)
150 { unsigned short xx[8], yy[4];
151 xx[4] = (unsigned short)x.lo;
152 xx[5] = (unsigned short)(x.lo >> 16);
153 xx[6] = (unsigned short)x.hi;
154 xx[7] = (unsigned short)(x.hi >> 16);
155 yy[0] = (unsigned short)y.lo;
156 yy[1] = (unsigned short)(y.lo >> 16);
157 yy[2] = (unsigned short)y.hi;
158 yy[3] = (unsigned short)(y.hi >> 16);
159 bigmul(4, 4, xx, yy);
160 x.lo = (unsigned int)xx[0] | ((unsigned int)xx[1] << 16);
161 x.hi = (unsigned int)xx[2] | ((unsigned int)xx[3] << 16);
162 return x;
163 }
164
165 /***********************************************************************
166 * NAME
167 *
168 * xldiv - divide long integers
169 *
170 * SYNOPSIS
171 *
172 * #include "glplib.h"
173 * glp_ldiv xldiv(glp_long x, glp_long y);
174 *
175 * RETURNS
176 *
177 * The routine xldiv returns a structure of type glp_ldiv containing
178 * members quot (the quotient) and rem (the remainder), both of type
179 * glp_long. */
180
181 glp_ldiv xldiv(glp_long x, glp_long y)
182 { glp_ldiv t;
183 int m, sx, sy;
184 unsigned short xx[8], yy[4];
185 /* sx := sign(x) */
186 sx = (x.hi < 0);
187 /* sy := sign(y) */
188 sy = (y.hi < 0);
189 /* x := |x| */
190 if (sx) x = xlneg(x);
191 /* y := |y| */
192 if (sy) y = xlneg(y);
193 /* compute x div y and x mod y */
194 xx[0] = (unsigned short)x.lo;
195 xx[1] = (unsigned short)(x.lo >> 16);
196 xx[2] = (unsigned short)x.hi;
197 xx[3] = (unsigned short)(x.hi >> 16);
198 yy[0] = (unsigned short)y.lo;
199 yy[1] = (unsigned short)(y.lo >> 16);
200 yy[2] = (unsigned short)y.hi;
201 yy[3] = (unsigned short)(y.hi >> 16);
202 if (yy[3])
203 m = 4;
204 else if (yy[2])
205 m = 3;
206 else if (yy[1])
207 m = 2;
208 else if (yy[0])
209 m = 1;
210 else
211 xerror("xldiv: divide by zero\n");
212 bigdiv(4 - m, m, xx, yy);
213 /* remainder in x[0], x[1], ..., x[m-1] */
214 t.rem.lo = (unsigned int)xx[0], t.rem.hi = 0;
215 if (m >= 2) t.rem.lo |= (unsigned int)xx[1] << 16;
216 if (m >= 3) t.rem.hi = (unsigned int)xx[2];
217 if (m >= 4) t.rem.hi |= (unsigned int)xx[3] << 16;
218 if (sx) t.rem = xlneg(t.rem);
219 /* quotient in x[m], x[m+1], ..., x[4] */
220 t.quot.lo = (unsigned int)xx[m], t.quot.hi = 0;
221 if (m <= 3) t.quot.lo |= (unsigned int)xx[m+1] << 16;
222 if (m <= 2) t.quot.hi = (unsigned int)xx[m+2];
223 if (m <= 1) t.quot.hi |= (unsigned int)xx[m+3] << 16;
224 if (sx ^ sy) t.quot = xlneg(t.quot);
225 return t;
226 }
227
228 /***********************************************************************
229 * NAME
230 *
231 * xltod - convert long integer to double
232 *
233 * SYNOPSIS
234 *
235 * #include "glplib.h"
236 * double xltod(glp_long x);
237 *
238 * RETURNS
239 *
240 * The routine xltod returns x converted to double. */
241
242 double xltod(glp_long x)
243 { double s, z;
244 if (x.hi >= 0)
245 s = +1.0;
246 else
247 s = -1.0, x = xlneg(x);
248 if (x.hi >= 0)
249 z = 4294967296.0 * (double)x.hi + (double)(unsigned int)x.lo;
250 else
251 { xassert(x.hi == 0x80000000 && x.lo == 0x00000000);
252 z = 9223372036854775808.0; /* 2^63 */
253 }
254 return s * z;
255 }
256
257 char *xltoa(glp_long x, char *s)
258 { /* convert long integer to character string */
259 static const char *d = "0123456789";
260 glp_ldiv t;
261 int neg, len;
262 if (x.hi >= 0)
263 neg = 0;
264 else
265 neg = 1, x = xlneg(x);
266 if (x.hi >= 0)
267 { len = 0;
268 while (!(x.hi == 0 && x.lo == 0))
269 { t = xldiv(x, xlset(10));
270 xassert(0 <= t.rem.lo && t.rem.lo <= 9);
271 s[len++] = d[t.rem.lo];
272 x = t.quot;
273 }
274 if (len == 0) s[len++] = d[0];
275 if (neg) s[len++] = '-';
276 s[len] = '\0';
277 strrev(s);
278 }
279 else
280 strcpy(s, "-9223372036854775808"); /* -2^63 */
281 return s;
282 }
283
284 /**********************************************************************/
285
286 #if 0
287 #include "glprng.h"
288
289 #define N_TEST 1000000
290 /* number of tests */
291
292 static glp_long myrand(RNG *rand)
293 { glp_long x;
294 int k;
295 k = rng_unif_rand(rand, 4);
296 xassert(0 <= k && k <= 3);
297 x.lo = rng_unif_rand(rand, 65536);
298 if (k == 1 || k == 3)
299 { x.lo <<= 16;
300 x.lo += rng_unif_rand(rand, 65536);
301 }
302 if (k <= 1)
303 x.hi = 0;
304 else
305 x.hi = rng_unif_rand(rand, 65536);
306 if (k == 3)
307 { x.hi <<= 16;
308 x.hi += rng_unif_rand(rand, 65536);
309 }
310 if (rng_unif_rand(rand, 2)) x = xlneg(x);
311 return x;
312 }
313
314 int main(void)
315 { RNG *rand;
316 glp_long x, y;
317 glp_ldiv z;
318 int test;
319 rand = rng_create_rand();
320 for (test = 1; test <= N_TEST; test++)
321 { x = myrand(rand);
322 y = myrand(rand);
323 if (y.lo == 0 && y.hi == 0) y.lo = 1;
324 /* z.quot := x div y, z.rem := x mod y */
325 z = xldiv(x, y);
326 /* x must be equal to y * z.quot + z.rem */
327 xassert(xlcmp(x, xladd(xlmul(y, z.quot), z.rem)) == 0);
328 }
329 xprintf("%d tests successfully passed\n", N_TEST);
330 rng_delete_rand(rand);
331 return 0;
332 }
333 #endif
334
335 /* eof */