COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glplib02.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 14 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 8.5 KB
Line 
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
42glp_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
62glp_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
84glp_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
106glp_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
125int 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
149glp_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
181glp_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
242double 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
257char *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
292static 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
314int 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 */
Note: See TracBrowser for help on using the repository browser.