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