COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glplib01.c @ 10:5545663ca997

subpack-glpk
Last change on this file since 10:5545663ca997 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import GLPK 4.47

File size: 9.4 KB
Line 
1/* glplib01.c (bignum arithmetic) */
2
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***********************************************************************/
24
25#include "glpenv.h"
26#include "glplib.h"
27
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***********************************************************************/
45
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. */
74
75void 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}
95
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. */
131
132void 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];
181decr:    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;
185test:    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;
189msub:    /* 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         }
209putq:    /* 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      }
228done: return;
229}
230
231/**********************************************************************/
232
233#if 0
234#include <assert.h>
235#include <stdio.h>
236#include <stdlib.h>
237#include "glprng.h"
238
239#define N_MAX 7
240/* maximal number of digits in multiplicand */
241
242#define M_MAX 5
243/* maximal number of digits in multiplier */
244
245#define N_TEST 1000000
246/* number of tests */
247
248int 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
286
287/* eof */
Note: See TracBrowser for help on using the repository browser.