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