1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glplib02.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,335 @@
1.4 +/* glplib02.c (64-bit arithmetic) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpenv.h"
1.29 +#include "glplib.h"
1.30 +
1.31 +/***********************************************************************
1.32 +* NAME
1.33 +*
1.34 +* xlset - expand integer to long integer
1.35 +*
1.36 +* SYNOPSIS
1.37 +*
1.38 +* #include "glplib.h"
1.39 +* glp_long xlset(int x);
1.40 +*
1.41 +* RETURNS
1.42 +*
1.43 +* The routine xlset returns x expanded to long integer. */
1.44 +
1.45 +glp_long xlset(int x)
1.46 +{ glp_long t;
1.47 + t.lo = x, t.hi = (x >= 0 ? 0 : -1);
1.48 + return t;
1.49 +}
1.50 +
1.51 +/***********************************************************************
1.52 +* NAME
1.53 +*
1.54 +* xlneg - negate long integer
1.55 +*
1.56 +* SYNOPSIS
1.57 +*
1.58 +* #include "glplib.h"
1.59 +* glp_long xlneg(glp_long x);
1.60 +*
1.61 +* RETURNS
1.62 +*
1.63 +* The routine xlneg returns the difference 0 - x. */
1.64 +
1.65 +glp_long xlneg(glp_long x)
1.66 +{ if (x.lo)
1.67 + x.lo = - x.lo, x.hi = ~x.hi;
1.68 + else
1.69 + x.hi = - x.hi;
1.70 + return x;
1.71 +}
1.72 +
1.73 +/***********************************************************************
1.74 +* NAME
1.75 +*
1.76 +* xladd - add long integers
1.77 +*
1.78 +* SYNOPSIS
1.79 +*
1.80 +* #include "glplib.h"
1.81 +* glp_long xladd(glp_long x, glp_long y);
1.82 +*
1.83 +* RETURNS
1.84 +*
1.85 +* The routine xladd returns the sum x + y. */
1.86 +
1.87 +glp_long xladd(glp_long x, glp_long y)
1.88 +{ if ((unsigned int)x.lo <= 0xFFFFFFFF - (unsigned int)y.lo)
1.89 + x.lo += y.lo, x.hi += y.hi;
1.90 + else
1.91 + x.lo += y.lo, x.hi += y.hi + 1;
1.92 + return x;
1.93 +}
1.94 +
1.95 +/***********************************************************************
1.96 +* NAME
1.97 +*
1.98 +* xlsub - subtract long integers
1.99 +*
1.100 +* SYNOPSIS
1.101 +*
1.102 +* #include "glplib.h"
1.103 +* glp_long xlsub(glp_long x, glp_long y);
1.104 +*
1.105 +* RETURNS
1.106 +*
1.107 +* The routine xlsub returns the difference x - y. */
1.108 +
1.109 +glp_long xlsub(glp_long x, glp_long y)
1.110 +{ return
1.111 + xladd(x, xlneg(y));
1.112 +}
1.113 +
1.114 +/***********************************************************************
1.115 +* NAME
1.116 +*
1.117 +* xlcmp - compare long integers
1.118 +*
1.119 +* SYNOPSIS
1.120 +*
1.121 +* #include "glplib.h"
1.122 +* int xlcmp(glp_long x, glp_long y);
1.123 +*
1.124 +* RETURNS
1.125 +*
1.126 +* The routine xlcmp returns the sign of the difference x - y. */
1.127 +
1.128 +int xlcmp(glp_long x, glp_long y)
1.129 +{ if (x.hi >= 0 && y.hi < 0) return +1;
1.130 + if (x.hi < 0 && y.hi >= 0) return -1;
1.131 + if ((unsigned int)x.hi < (unsigned int)y.hi) return -1;
1.132 + if ((unsigned int)x.hi > (unsigned int)y.hi) return +1;
1.133 + if ((unsigned int)x.lo < (unsigned int)y.lo) return -1;
1.134 + if ((unsigned int)x.lo > (unsigned int)y.lo) return +1;
1.135 + return 0;
1.136 +}
1.137 +
1.138 +/***********************************************************************
1.139 +* NAME
1.140 +*
1.141 +* xlmul - multiply long integers
1.142 +*
1.143 +* SYNOPSIS
1.144 +*
1.145 +* #include "glplib.h"
1.146 +* glp_long xlmul(glp_long x, glp_long y);
1.147 +*
1.148 +* RETURNS
1.149 +*
1.150 +* The routine xlmul returns the product x * y. */
1.151 +
1.152 +glp_long xlmul(glp_long x, glp_long y)
1.153 +{ unsigned short xx[8], yy[4];
1.154 + xx[4] = (unsigned short)x.lo;
1.155 + xx[5] = (unsigned short)(x.lo >> 16);
1.156 + xx[6] = (unsigned short)x.hi;
1.157 + xx[7] = (unsigned short)(x.hi >> 16);
1.158 + yy[0] = (unsigned short)y.lo;
1.159 + yy[1] = (unsigned short)(y.lo >> 16);
1.160 + yy[2] = (unsigned short)y.hi;
1.161 + yy[3] = (unsigned short)(y.hi >> 16);
1.162 + bigmul(4, 4, xx, yy);
1.163 + x.lo = (unsigned int)xx[0] | ((unsigned int)xx[1] << 16);
1.164 + x.hi = (unsigned int)xx[2] | ((unsigned int)xx[3] << 16);
1.165 + return x;
1.166 +}
1.167 +
1.168 +/***********************************************************************
1.169 +* NAME
1.170 +*
1.171 +* xldiv - divide long integers
1.172 +*
1.173 +* SYNOPSIS
1.174 +*
1.175 +* #include "glplib.h"
1.176 +* glp_ldiv xldiv(glp_long x, glp_long y);
1.177 +*
1.178 +* RETURNS
1.179 +*
1.180 +* The routine xldiv returns a structure of type glp_ldiv containing
1.181 +* members quot (the quotient) and rem (the remainder), both of type
1.182 +* glp_long. */
1.183 +
1.184 +glp_ldiv xldiv(glp_long x, glp_long y)
1.185 +{ glp_ldiv t;
1.186 + int m, sx, sy;
1.187 + unsigned short xx[8], yy[4];
1.188 + /* sx := sign(x) */
1.189 + sx = (x.hi < 0);
1.190 + /* sy := sign(y) */
1.191 + sy = (y.hi < 0);
1.192 + /* x := |x| */
1.193 + if (sx) x = xlneg(x);
1.194 + /* y := |y| */
1.195 + if (sy) y = xlneg(y);
1.196 + /* compute x div y and x mod y */
1.197 + xx[0] = (unsigned short)x.lo;
1.198 + xx[1] = (unsigned short)(x.lo >> 16);
1.199 + xx[2] = (unsigned short)x.hi;
1.200 + xx[3] = (unsigned short)(x.hi >> 16);
1.201 + yy[0] = (unsigned short)y.lo;
1.202 + yy[1] = (unsigned short)(y.lo >> 16);
1.203 + yy[2] = (unsigned short)y.hi;
1.204 + yy[3] = (unsigned short)(y.hi >> 16);
1.205 + if (yy[3])
1.206 + m = 4;
1.207 + else if (yy[2])
1.208 + m = 3;
1.209 + else if (yy[1])
1.210 + m = 2;
1.211 + else if (yy[0])
1.212 + m = 1;
1.213 + else
1.214 + xerror("xldiv: divide by zero\n");
1.215 + bigdiv(4 - m, m, xx, yy);
1.216 + /* remainder in x[0], x[1], ..., x[m-1] */
1.217 + t.rem.lo = (unsigned int)xx[0], t.rem.hi = 0;
1.218 + if (m >= 2) t.rem.lo |= (unsigned int)xx[1] << 16;
1.219 + if (m >= 3) t.rem.hi = (unsigned int)xx[2];
1.220 + if (m >= 4) t.rem.hi |= (unsigned int)xx[3] << 16;
1.221 + if (sx) t.rem = xlneg(t.rem);
1.222 + /* quotient in x[m], x[m+1], ..., x[4] */
1.223 + t.quot.lo = (unsigned int)xx[m], t.quot.hi = 0;
1.224 + if (m <= 3) t.quot.lo |= (unsigned int)xx[m+1] << 16;
1.225 + if (m <= 2) t.quot.hi = (unsigned int)xx[m+2];
1.226 + if (m <= 1) t.quot.hi |= (unsigned int)xx[m+3] << 16;
1.227 + if (sx ^ sy) t.quot = xlneg(t.quot);
1.228 + return t;
1.229 +}
1.230 +
1.231 +/***********************************************************************
1.232 +* NAME
1.233 +*
1.234 +* xltod - convert long integer to double
1.235 +*
1.236 +* SYNOPSIS
1.237 +*
1.238 +* #include "glplib.h"
1.239 +* double xltod(glp_long x);
1.240 +*
1.241 +* RETURNS
1.242 +*
1.243 +* The routine xltod returns x converted to double. */
1.244 +
1.245 +double xltod(glp_long x)
1.246 +{ double s, z;
1.247 + if (x.hi >= 0)
1.248 + s = +1.0;
1.249 + else
1.250 + s = -1.0, x = xlneg(x);
1.251 + if (x.hi >= 0)
1.252 + z = 4294967296.0 * (double)x.hi + (double)(unsigned int)x.lo;
1.253 + else
1.254 + { xassert(x.hi == 0x80000000 && x.lo == 0x00000000);
1.255 + z = 9223372036854775808.0; /* 2^63 */
1.256 + }
1.257 + return s * z;
1.258 +}
1.259 +
1.260 +char *xltoa(glp_long x, char *s)
1.261 +{ /* convert long integer to character string */
1.262 + static const char *d = "0123456789";
1.263 + glp_ldiv t;
1.264 + int neg, len;
1.265 + if (x.hi >= 0)
1.266 + neg = 0;
1.267 + else
1.268 + neg = 1, x = xlneg(x);
1.269 + if (x.hi >= 0)
1.270 + { len = 0;
1.271 + while (!(x.hi == 0 && x.lo == 0))
1.272 + { t = xldiv(x, xlset(10));
1.273 + xassert(0 <= t.rem.lo && t.rem.lo <= 9);
1.274 + s[len++] = d[t.rem.lo];
1.275 + x = t.quot;
1.276 + }
1.277 + if (len == 0) s[len++] = d[0];
1.278 + if (neg) s[len++] = '-';
1.279 + s[len] = '\0';
1.280 + strrev(s);
1.281 + }
1.282 + else
1.283 + strcpy(s, "-9223372036854775808"); /* -2^63 */
1.284 + return s;
1.285 +}
1.286 +
1.287 +/**********************************************************************/
1.288 +
1.289 +#if 0
1.290 +#include "glprng.h"
1.291 +
1.292 +#define N_TEST 1000000
1.293 +/* number of tests */
1.294 +
1.295 +static glp_long myrand(RNG *rand)
1.296 +{ glp_long x;
1.297 + int k;
1.298 + k = rng_unif_rand(rand, 4);
1.299 + xassert(0 <= k && k <= 3);
1.300 + x.lo = rng_unif_rand(rand, 65536);
1.301 + if (k == 1 || k == 3)
1.302 + { x.lo <<= 16;
1.303 + x.lo += rng_unif_rand(rand, 65536);
1.304 + }
1.305 + if (k <= 1)
1.306 + x.hi = 0;
1.307 + else
1.308 + x.hi = rng_unif_rand(rand, 65536);
1.309 + if (k == 3)
1.310 + { x.hi <<= 16;
1.311 + x.hi += rng_unif_rand(rand, 65536);
1.312 + }
1.313 + if (rng_unif_rand(rand, 2)) x = xlneg(x);
1.314 + return x;
1.315 +}
1.316 +
1.317 +int main(void)
1.318 +{ RNG *rand;
1.319 + glp_long x, y;
1.320 + glp_ldiv z;
1.321 + int test;
1.322 + rand = rng_create_rand();
1.323 + for (test = 1; test <= N_TEST; test++)
1.324 + { x = myrand(rand);
1.325 + y = myrand(rand);
1.326 + if (y.lo == 0 && y.hi == 0) y.lo = 1;
1.327 + /* z.quot := x div y, z.rem := x mod y */
1.328 + z = xldiv(x, y);
1.329 + /* x must be equal to y * z.quot + z.rem */
1.330 + xassert(xlcmp(x, xladd(xlmul(y, z.quot), z.rem)) == 0);
1.331 + }
1.332 + xprintf("%d tests successfully passed\n", N_TEST);
1.333 + rng_delete_rand(rand);
1.334 + return 0;
1.335 +}
1.336 +#endif
1.337 +
1.338 +/* eof */