lemon-project-template-glpk

annotate deps/glpk/src/glplib01.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 /* glplib01.c (bignum 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 * Two routines below are intended to multiply and divide unsigned
alpar@9 30 * integer numbers of arbitrary precision.
alpar@9 31 *
alpar@9 32 * The routines assume that an unsigned integer number is represented in
alpar@9 33 * the positional numeral system with the base 2^16 = 65536, i.e. each
alpar@9 34 * "digit" of the number is in the range [0, 65535] and represented as
alpar@9 35 * a 16-bit value of the unsigned short type. In other words, a number x
alpar@9 36 * has the following representation:
alpar@9 37 *
alpar@9 38 * n-1
alpar@9 39 * x = sum d[j] * 65536^j,
alpar@9 40 * j=0
alpar@9 41 *
alpar@9 42 * where n is the number of places (positions), and d[j] is j-th "digit"
alpar@9 43 * of x, 0 <= d[j] <= 65535.
alpar@9 44 ***********************************************************************/
alpar@9 45
alpar@9 46 /***********************************************************************
alpar@9 47 * NAME
alpar@9 48 *
alpar@9 49 * bigmul - multiply unsigned integer numbers of arbitrary precision
alpar@9 50 *
alpar@9 51 * SYNOPSIS
alpar@9 52 *
alpar@9 53 * #include "glplib.h"
alpar@9 54 * void bigmul(int n, int m, unsigned short x[], unsigned short y[]);
alpar@9 55 *
alpar@9 56 * DESCRIPTION
alpar@9 57 *
alpar@9 58 * The routine bigmul multiplies unsigned integer numbers of arbitrary
alpar@9 59 * precision.
alpar@9 60 *
alpar@9 61 * n is the number of digits of multiplicand, n >= 1;
alpar@9 62 *
alpar@9 63 * m is the number of digits of multiplier, m >= 1;
alpar@9 64 *
alpar@9 65 * x is an array containing digits of the multiplicand in elements
alpar@9 66 * x[m], x[m+1], ..., x[n+m-1]. Contents of x[0], x[1], ..., x[m-1] are
alpar@9 67 * ignored on entry.
alpar@9 68 *
alpar@9 69 * y is an array containing digits of the multiplier in elements y[0],
alpar@9 70 * y[1], ..., y[m-1].
alpar@9 71 *
alpar@9 72 * On exit digits of the product are stored in elements x[0], x[1], ...,
alpar@9 73 * x[n+m-1]. The array y is not changed. */
alpar@9 74
alpar@9 75 void bigmul(int n, int m, unsigned short x[], unsigned short y[])
alpar@9 76 { int i, j;
alpar@9 77 unsigned int t;
alpar@9 78 xassert(n >= 1);
alpar@9 79 xassert(m >= 1);
alpar@9 80 for (j = 0; j < m; j++) x[j] = 0;
alpar@9 81 for (i = 0; i < n; i++)
alpar@9 82 { if (x[i+m])
alpar@9 83 { t = 0;
alpar@9 84 for (j = 0; j < m; j++)
alpar@9 85 { t += (unsigned int)x[i+m] * (unsigned int)y[j] +
alpar@9 86 (unsigned int)x[i+j];
alpar@9 87 x[i+j] = (unsigned short)t;
alpar@9 88 t >>= 16;
alpar@9 89 }
alpar@9 90 x[i+m] = (unsigned short)t;
alpar@9 91 }
alpar@9 92 }
alpar@9 93 return;
alpar@9 94 }
alpar@9 95
alpar@9 96 /***********************************************************************
alpar@9 97 * NAME
alpar@9 98 *
alpar@9 99 * bigdiv - divide unsigned integer numbers of arbitrary precision
alpar@9 100 *
alpar@9 101 * SYNOPSIS
alpar@9 102 *
alpar@9 103 * #include "glplib.h"
alpar@9 104 * void bigdiv(int n, int m, unsigned short x[], unsigned short y[]);
alpar@9 105 *
alpar@9 106 * DESCRIPTION
alpar@9 107 *
alpar@9 108 * The routine bigdiv divides one unsigned integer number of arbitrary
alpar@9 109 * precision by another with the algorithm described in [1].
alpar@9 110 *
alpar@9 111 * n is the difference between the number of digits of dividend and the
alpar@9 112 * number of digits of divisor, n >= 0.
alpar@9 113 *
alpar@9 114 * m is the number of digits of divisor, m >= 1.
alpar@9 115 *
alpar@9 116 * x is an array containing digits of the dividend in elements x[0],
alpar@9 117 * x[1], ..., x[n+m-1].
alpar@9 118 *
alpar@9 119 * y is an array containing digits of the divisor in elements y[0],
alpar@9 120 * y[1], ..., y[m-1]. The highest digit y[m-1] must be non-zero.
alpar@9 121 *
alpar@9 122 * On exit n+1 digits of the quotient are stored in elements x[m],
alpar@9 123 * x[m+1], ..., x[n+m], and m digits of the remainder are stored in
alpar@9 124 * elements x[0], x[1], ..., x[m-1]. The array y is changed but then
alpar@9 125 * restored.
alpar@9 126 *
alpar@9 127 * REFERENCES
alpar@9 128 *
alpar@9 129 * 1. D. Knuth. The Art of Computer Programming. Vol. 2: Seminumerical
alpar@9 130 * Algorithms. Stanford University, 1969. */
alpar@9 131
alpar@9 132 void bigdiv(int n, int m, unsigned short x[], unsigned short y[])
alpar@9 133 { int i, j;
alpar@9 134 unsigned int t;
alpar@9 135 unsigned short d, q, r;
alpar@9 136 xassert(n >= 0);
alpar@9 137 xassert(m >= 1);
alpar@9 138 xassert(y[m-1] != 0);
alpar@9 139 /* special case when divisor has the only digit */
alpar@9 140 if (m == 1)
alpar@9 141 { d = 0;
alpar@9 142 for (i = n; i >= 0; i--)
alpar@9 143 { t = ((unsigned int)d << 16) + (unsigned int)x[i];
alpar@9 144 x[i+1] = (unsigned short)(t / y[0]);
alpar@9 145 d = (unsigned short)(t % y[0]);
alpar@9 146 }
alpar@9 147 x[0] = d;
alpar@9 148 goto done;
alpar@9 149 }
alpar@9 150 /* multiply dividend and divisor by a normalizing coefficient in
alpar@9 151 order to provide the condition y[m-1] >= base / 2 */
alpar@9 152 d = (unsigned short)(0x10000 / ((unsigned int)y[m-1] + 1));
alpar@9 153 if (d == 1)
alpar@9 154 x[n+m] = 0;
alpar@9 155 else
alpar@9 156 { t = 0;
alpar@9 157 for (i = 0; i < n+m; i++)
alpar@9 158 { t += (unsigned int)x[i] * (unsigned int)d;
alpar@9 159 x[i] = (unsigned short)t;
alpar@9 160 t >>= 16;
alpar@9 161 }
alpar@9 162 x[n+m] = (unsigned short)t;
alpar@9 163 t = 0;
alpar@9 164 for (j = 0; j < m; j++)
alpar@9 165 { t += (unsigned int)y[j] * (unsigned int)d;
alpar@9 166 y[j] = (unsigned short)t;
alpar@9 167 t >>= 16;
alpar@9 168 }
alpar@9 169 }
alpar@9 170 /* main loop */
alpar@9 171 for (i = n; i >= 0; i--)
alpar@9 172 { /* estimate and correct the current digit of quotient */
alpar@9 173 if (x[i+m] < y[m-1])
alpar@9 174 { t = ((unsigned int)x[i+m] << 16) + (unsigned int)x[i+m-1];
alpar@9 175 q = (unsigned short)(t / (unsigned int)y[m-1]);
alpar@9 176 r = (unsigned short)(t % (unsigned int)y[m-1]);
alpar@9 177 if (q == 0) goto putq; else goto test;
alpar@9 178 }
alpar@9 179 q = 0;
alpar@9 180 r = x[i+m-1];
alpar@9 181 decr: q--; /* if q = 0 then q-- = 0xFFFF */
alpar@9 182 t = (unsigned int)r + (unsigned int)y[m-1];
alpar@9 183 r = (unsigned short)t;
alpar@9 184 if (t > 0xFFFF) goto msub;
alpar@9 185 test: t = (unsigned int)y[m-2] * (unsigned int)q;
alpar@9 186 if ((unsigned short)(t >> 16) > r) goto decr;
alpar@9 187 if ((unsigned short)(t >> 16) < r) goto msub;
alpar@9 188 if ((unsigned short)t > x[i+m-2]) goto decr;
alpar@9 189 msub: /* now subtract divisor multiplied by the current digit of
alpar@9 190 quotient from the current dividend */
alpar@9 191 if (q == 0) goto putq;
alpar@9 192 t = 0;
alpar@9 193 for (j = 0; j < m; j++)
alpar@9 194 { t += (unsigned int)y[j] * (unsigned int)q;
alpar@9 195 if (x[i+j] < (unsigned short)t) t += 0x10000;
alpar@9 196 x[i+j] -= (unsigned short)t;
alpar@9 197 t >>= 16;
alpar@9 198 }
alpar@9 199 if (x[i+m] >= (unsigned short)t) goto putq;
alpar@9 200 /* perform correcting addition, because the current digit of
alpar@9 201 quotient is greater by one than its correct value */
alpar@9 202 q--;
alpar@9 203 t = 0;
alpar@9 204 for (j = 0; j < m; j++)
alpar@9 205 { t += (unsigned int)x[i+j] + (unsigned int)y[j];
alpar@9 206 x[i+j] = (unsigned short)t;
alpar@9 207 t >>= 16;
alpar@9 208 }
alpar@9 209 putq: /* store the current digit of quotient */
alpar@9 210 x[i+m] = q;
alpar@9 211 }
alpar@9 212 /* divide divisor and remainder by the normalizing coefficient in
alpar@9 213 order to restore their original values */
alpar@9 214 if (d > 1)
alpar@9 215 { t = 0;
alpar@9 216 for (i = m-1; i >= 0; i--)
alpar@9 217 { t = (t << 16) + (unsigned int)x[i];
alpar@9 218 x[i] = (unsigned short)(t / (unsigned int)d);
alpar@9 219 t %= (unsigned int)d;
alpar@9 220 }
alpar@9 221 t = 0;
alpar@9 222 for (j = m-1; j >= 0; j--)
alpar@9 223 { t = (t << 16) + (unsigned int)y[j];
alpar@9 224 y[j] = (unsigned short)(t / (unsigned int)d);
alpar@9 225 t %= (unsigned int)d;
alpar@9 226 }
alpar@9 227 }
alpar@9 228 done: return;
alpar@9 229 }
alpar@9 230
alpar@9 231 /**********************************************************************/
alpar@9 232
alpar@9 233 #if 0
alpar@9 234 #include <assert.h>
alpar@9 235 #include <stdio.h>
alpar@9 236 #include <stdlib.h>
alpar@9 237 #include "glprng.h"
alpar@9 238
alpar@9 239 #define N_MAX 7
alpar@9 240 /* maximal number of digits in multiplicand */
alpar@9 241
alpar@9 242 #define M_MAX 5
alpar@9 243 /* maximal number of digits in multiplier */
alpar@9 244
alpar@9 245 #define N_TEST 1000000
alpar@9 246 /* number of tests */
alpar@9 247
alpar@9 248 int main(void)
alpar@9 249 { RNG *rand;
alpar@9 250 int d, j, n, m, test;
alpar@9 251 unsigned short x[N_MAX], y[M_MAX], z[N_MAX+M_MAX];
alpar@9 252 rand = rng_create_rand();
alpar@9 253 for (test = 1; test <= N_TEST; test++)
alpar@9 254 { /* x[0,...,n-1] := multiplicand */
alpar@9 255 n = 1 + rng_unif_rand(rand, N_MAX-1);
alpar@9 256 assert(1 <= n && n <= N_MAX);
alpar@9 257 for (j = 0; j < n; j++)
alpar@9 258 { d = rng_unif_rand(rand, 65536);
alpar@9 259 assert(0 <= d && d <= 65535);
alpar@9 260 x[j] = (unsigned short)d;
alpar@9 261 }
alpar@9 262 /* y[0,...,m-1] := multiplier */
alpar@9 263 m = 1 + rng_unif_rand(rand, M_MAX-1);
alpar@9 264 assert(1 <= m && m <= M_MAX);
alpar@9 265 for (j = 0; j < m; j++)
alpar@9 266 { d = rng_unif_rand(rand, 65536);
alpar@9 267 assert(0 <= d && d <= 65535);
alpar@9 268 y[j] = (unsigned short)d;
alpar@9 269 }
alpar@9 270 if (y[m-1] == 0) y[m-1] = 1;
alpar@9 271 /* z[0,...,n+m-1] := x * y */
alpar@9 272 for (j = 0; j < n; j++) z[m+j] = x[j];
alpar@9 273 bigmul(n, m, z, y);
alpar@9 274 /* z[0,...,m-1] := z mod y, z[m,...,n+m-1] := z div y */
alpar@9 275 bigdiv(n, m, z, y);
alpar@9 276 /* z mod y must be 0 */
alpar@9 277 for (j = 0; j < m; j++) assert(z[j] == 0);
alpar@9 278 /* z div y must be x */
alpar@9 279 for (j = 0; j < n; j++) assert(z[m+j] == x[j]);
alpar@9 280 }
alpar@9 281 fprintf(stderr, "%d tests successfully passed\n", N_TEST);
alpar@9 282 rng_delete_rand(rand);
alpar@9 283 return 0;
alpar@9 284 }
alpar@9 285 #endif
alpar@9 286
alpar@9 287 /* eof */