1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glplib01.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,287 @@
1.4 +/* glplib01.c (bignum 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 +* Two routines below are intended to multiply and divide unsigned
1.33 +* integer numbers of arbitrary precision.
1.34 +*
1.35 +* The routines assume that an unsigned integer number is represented in
1.36 +* the positional numeral system with the base 2^16 = 65536, i.e. each
1.37 +* "digit" of the number is in the range [0, 65535] and represented as
1.38 +* a 16-bit value of the unsigned short type. In other words, a number x
1.39 +* has the following representation:
1.40 +*
1.41 +* n-1
1.42 +* x = sum d[j] * 65536^j,
1.43 +* j=0
1.44 +*
1.45 +* where n is the number of places (positions), and d[j] is j-th "digit"
1.46 +* of x, 0 <= d[j] <= 65535.
1.47 +***********************************************************************/
1.48 +
1.49 +/***********************************************************************
1.50 +* NAME
1.51 +*
1.52 +* bigmul - multiply unsigned integer numbers of arbitrary precision
1.53 +*
1.54 +* SYNOPSIS
1.55 +*
1.56 +* #include "glplib.h"
1.57 +* void bigmul(int n, int m, unsigned short x[], unsigned short y[]);
1.58 +*
1.59 +* DESCRIPTION
1.60 +*
1.61 +* The routine bigmul multiplies unsigned integer numbers of arbitrary
1.62 +* precision.
1.63 +*
1.64 +* n is the number of digits of multiplicand, n >= 1;
1.65 +*
1.66 +* m is the number of digits of multiplier, m >= 1;
1.67 +*
1.68 +* x is an array containing digits of the multiplicand in elements
1.69 +* x[m], x[m+1], ..., x[n+m-1]. Contents of x[0], x[1], ..., x[m-1] are
1.70 +* ignored on entry.
1.71 +*
1.72 +* y is an array containing digits of the multiplier in elements y[0],
1.73 +* y[1], ..., y[m-1].
1.74 +*
1.75 +* On exit digits of the product are stored in elements x[0], x[1], ...,
1.76 +* x[n+m-1]. The array y is not changed. */
1.77 +
1.78 +void bigmul(int n, int m, unsigned short x[], unsigned short y[])
1.79 +{ int i, j;
1.80 + unsigned int t;
1.81 + xassert(n >= 1);
1.82 + xassert(m >= 1);
1.83 + for (j = 0; j < m; j++) x[j] = 0;
1.84 + for (i = 0; i < n; i++)
1.85 + { if (x[i+m])
1.86 + { t = 0;
1.87 + for (j = 0; j < m; j++)
1.88 + { t += (unsigned int)x[i+m] * (unsigned int)y[j] +
1.89 + (unsigned int)x[i+j];
1.90 + x[i+j] = (unsigned short)t;
1.91 + t >>= 16;
1.92 + }
1.93 + x[i+m] = (unsigned short)t;
1.94 + }
1.95 + }
1.96 + return;
1.97 +}
1.98 +
1.99 +/***********************************************************************
1.100 +* NAME
1.101 +*
1.102 +* bigdiv - divide unsigned integer numbers of arbitrary precision
1.103 +*
1.104 +* SYNOPSIS
1.105 +*
1.106 +* #include "glplib.h"
1.107 +* void bigdiv(int n, int m, unsigned short x[], unsigned short y[]);
1.108 +*
1.109 +* DESCRIPTION
1.110 +*
1.111 +* The routine bigdiv divides one unsigned integer number of arbitrary
1.112 +* precision by another with the algorithm described in [1].
1.113 +*
1.114 +* n is the difference between the number of digits of dividend and the
1.115 +* number of digits of divisor, n >= 0.
1.116 +*
1.117 +* m is the number of digits of divisor, m >= 1.
1.118 +*
1.119 +* x is an array containing digits of the dividend in elements x[0],
1.120 +* x[1], ..., x[n+m-1].
1.121 +*
1.122 +* y is an array containing digits of the divisor in elements y[0],
1.123 +* y[1], ..., y[m-1]. The highest digit y[m-1] must be non-zero.
1.124 +*
1.125 +* On exit n+1 digits of the quotient are stored in elements x[m],
1.126 +* x[m+1], ..., x[n+m], and m digits of the remainder are stored in
1.127 +* elements x[0], x[1], ..., x[m-1]. The array y is changed but then
1.128 +* restored.
1.129 +*
1.130 +* REFERENCES
1.131 +*
1.132 +* 1. D. Knuth. The Art of Computer Programming. Vol. 2: Seminumerical
1.133 +* Algorithms. Stanford University, 1969. */
1.134 +
1.135 +void bigdiv(int n, int m, unsigned short x[], unsigned short y[])
1.136 +{ int i, j;
1.137 + unsigned int t;
1.138 + unsigned short d, q, r;
1.139 + xassert(n >= 0);
1.140 + xassert(m >= 1);
1.141 + xassert(y[m-1] != 0);
1.142 + /* special case when divisor has the only digit */
1.143 + if (m == 1)
1.144 + { d = 0;
1.145 + for (i = n; i >= 0; i--)
1.146 + { t = ((unsigned int)d << 16) + (unsigned int)x[i];
1.147 + x[i+1] = (unsigned short)(t / y[0]);
1.148 + d = (unsigned short)(t % y[0]);
1.149 + }
1.150 + x[0] = d;
1.151 + goto done;
1.152 + }
1.153 + /* multiply dividend and divisor by a normalizing coefficient in
1.154 + order to provide the condition y[m-1] >= base / 2 */
1.155 + d = (unsigned short)(0x10000 / ((unsigned int)y[m-1] + 1));
1.156 + if (d == 1)
1.157 + x[n+m] = 0;
1.158 + else
1.159 + { t = 0;
1.160 + for (i = 0; i < n+m; i++)
1.161 + { t += (unsigned int)x[i] * (unsigned int)d;
1.162 + x[i] = (unsigned short)t;
1.163 + t >>= 16;
1.164 + }
1.165 + x[n+m] = (unsigned short)t;
1.166 + t = 0;
1.167 + for (j = 0; j < m; j++)
1.168 + { t += (unsigned int)y[j] * (unsigned int)d;
1.169 + y[j] = (unsigned short)t;
1.170 + t >>= 16;
1.171 + }
1.172 + }
1.173 + /* main loop */
1.174 + for (i = n; i >= 0; i--)
1.175 + { /* estimate and correct the current digit of quotient */
1.176 + if (x[i+m] < y[m-1])
1.177 + { t = ((unsigned int)x[i+m] << 16) + (unsigned int)x[i+m-1];
1.178 + q = (unsigned short)(t / (unsigned int)y[m-1]);
1.179 + r = (unsigned short)(t % (unsigned int)y[m-1]);
1.180 + if (q == 0) goto putq; else goto test;
1.181 + }
1.182 + q = 0;
1.183 + r = x[i+m-1];
1.184 +decr: q--; /* if q = 0 then q-- = 0xFFFF */
1.185 + t = (unsigned int)r + (unsigned int)y[m-1];
1.186 + r = (unsigned short)t;
1.187 + if (t > 0xFFFF) goto msub;
1.188 +test: t = (unsigned int)y[m-2] * (unsigned int)q;
1.189 + if ((unsigned short)(t >> 16) > r) goto decr;
1.190 + if ((unsigned short)(t >> 16) < r) goto msub;
1.191 + if ((unsigned short)t > x[i+m-2]) goto decr;
1.192 +msub: /* now subtract divisor multiplied by the current digit of
1.193 + quotient from the current dividend */
1.194 + if (q == 0) goto putq;
1.195 + t = 0;
1.196 + for (j = 0; j < m; j++)
1.197 + { t += (unsigned int)y[j] * (unsigned int)q;
1.198 + if (x[i+j] < (unsigned short)t) t += 0x10000;
1.199 + x[i+j] -= (unsigned short)t;
1.200 + t >>= 16;
1.201 + }
1.202 + if (x[i+m] >= (unsigned short)t) goto putq;
1.203 + /* perform correcting addition, because the current digit of
1.204 + quotient is greater by one than its correct value */
1.205 + q--;
1.206 + t = 0;
1.207 + for (j = 0; j < m; j++)
1.208 + { t += (unsigned int)x[i+j] + (unsigned int)y[j];
1.209 + x[i+j] = (unsigned short)t;
1.210 + t >>= 16;
1.211 + }
1.212 +putq: /* store the current digit of quotient */
1.213 + x[i+m] = q;
1.214 + }
1.215 + /* divide divisor and remainder by the normalizing coefficient in
1.216 + order to restore their original values */
1.217 + if (d > 1)
1.218 + { t = 0;
1.219 + for (i = m-1; i >= 0; i--)
1.220 + { t = (t << 16) + (unsigned int)x[i];
1.221 + x[i] = (unsigned short)(t / (unsigned int)d);
1.222 + t %= (unsigned int)d;
1.223 + }
1.224 + t = 0;
1.225 + for (j = m-1; j >= 0; j--)
1.226 + { t = (t << 16) + (unsigned int)y[j];
1.227 + y[j] = (unsigned short)(t / (unsigned int)d);
1.228 + t %= (unsigned int)d;
1.229 + }
1.230 + }
1.231 +done: return;
1.232 +}
1.233 +
1.234 +/**********************************************************************/
1.235 +
1.236 +#if 0
1.237 +#include <assert.h>
1.238 +#include <stdio.h>
1.239 +#include <stdlib.h>
1.240 +#include "glprng.h"
1.241 +
1.242 +#define N_MAX 7
1.243 +/* maximal number of digits in multiplicand */
1.244 +
1.245 +#define M_MAX 5
1.246 +/* maximal number of digits in multiplier */
1.247 +
1.248 +#define N_TEST 1000000
1.249 +/* number of tests */
1.250 +
1.251 +int main(void)
1.252 +{ RNG *rand;
1.253 + int d, j, n, m, test;
1.254 + unsigned short x[N_MAX], y[M_MAX], z[N_MAX+M_MAX];
1.255 + rand = rng_create_rand();
1.256 + for (test = 1; test <= N_TEST; test++)
1.257 + { /* x[0,...,n-1] := multiplicand */
1.258 + n = 1 + rng_unif_rand(rand, N_MAX-1);
1.259 + assert(1 <= n && n <= N_MAX);
1.260 + for (j = 0; j < n; j++)
1.261 + { d = rng_unif_rand(rand, 65536);
1.262 + assert(0 <= d && d <= 65535);
1.263 + x[j] = (unsigned short)d;
1.264 + }
1.265 + /* y[0,...,m-1] := multiplier */
1.266 + m = 1 + rng_unif_rand(rand, M_MAX-1);
1.267 + assert(1 <= m && m <= M_MAX);
1.268 + for (j = 0; j < m; j++)
1.269 + { d = rng_unif_rand(rand, 65536);
1.270 + assert(0 <= d && d <= 65535);
1.271 + y[j] = (unsigned short)d;
1.272 + }
1.273 + if (y[m-1] == 0) y[m-1] = 1;
1.274 + /* z[0,...,n+m-1] := x * y */
1.275 + for (j = 0; j < n; j++) z[m+j] = x[j];
1.276 + bigmul(n, m, z, y);
1.277 + /* z[0,...,m-1] := z mod y, z[m,...,n+m-1] := z div y */
1.278 + bigdiv(n, m, z, y);
1.279 + /* z mod y must be 0 */
1.280 + for (j = 0; j < m; j++) assert(z[j] == 0);
1.281 + /* z div y must be x */
1.282 + for (j = 0; j < n; j++) assert(z[m+j] == x[j]);
1.283 + }
1.284 + fprintf(stderr, "%d tests successfully passed\n", N_TEST);
1.285 + rng_delete_rand(rand);
1.286 + return 0;
1.287 +}
1.288 +#endif
1.289 +
1.290 +/* eof */