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 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 */ |
---|