lemon-project-template-glpk

annotate deps/glpk/src/glplib02.c @ 9:33de93886c88

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