alpar@9: /* glplib01.c (bignum arithmetic) */ alpar@9: alpar@9: /*********************************************************************** alpar@9: * This code is part of GLPK (GNU Linear Programming Kit). alpar@9: * alpar@9: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@9: * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, alpar@9: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@9: * E-mail: . alpar@9: * alpar@9: * GLPK is free software: you can redistribute it and/or modify it alpar@9: * under the terms of the GNU General Public License as published by alpar@9: * the Free Software Foundation, either version 3 of the License, or alpar@9: * (at your option) any later version. alpar@9: * alpar@9: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@9: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@9: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@9: * License for more details. alpar@9: * alpar@9: * You should have received a copy of the GNU General Public License alpar@9: * along with GLPK. If not, see . alpar@9: ***********************************************************************/ alpar@9: alpar@9: #include "glpenv.h" alpar@9: #include "glplib.h" alpar@9: alpar@9: /*********************************************************************** alpar@9: * Two routines below are intended to multiply and divide unsigned alpar@9: * integer numbers of arbitrary precision. alpar@9: * alpar@9: * The routines assume that an unsigned integer number is represented in alpar@9: * the positional numeral system with the base 2^16 = 65536, i.e. each alpar@9: * "digit" of the number is in the range [0, 65535] and represented as alpar@9: * a 16-bit value of the unsigned short type. In other words, a number x alpar@9: * has the following representation: alpar@9: * alpar@9: * n-1 alpar@9: * x = sum d[j] * 65536^j, alpar@9: * j=0 alpar@9: * alpar@9: * where n is the number of places (positions), and d[j] is j-th "digit" alpar@9: * of x, 0 <= d[j] <= 65535. alpar@9: ***********************************************************************/ alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * bigmul - multiply unsigned integer numbers of arbitrary precision alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glplib.h" alpar@9: * void bigmul(int n, int m, unsigned short x[], unsigned short y[]); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine bigmul multiplies unsigned integer numbers of arbitrary alpar@9: * precision. alpar@9: * alpar@9: * n is the number of digits of multiplicand, n >= 1; alpar@9: * alpar@9: * m is the number of digits of multiplier, m >= 1; alpar@9: * alpar@9: * x is an array containing digits of the multiplicand in elements alpar@9: * x[m], x[m+1], ..., x[n+m-1]. Contents of x[0], x[1], ..., x[m-1] are alpar@9: * ignored on entry. alpar@9: * alpar@9: * y is an array containing digits of the multiplier in elements y[0], alpar@9: * y[1], ..., y[m-1]. alpar@9: * alpar@9: * On exit digits of the product are stored in elements x[0], x[1], ..., alpar@9: * x[n+m-1]. The array y is not changed. */ alpar@9: alpar@9: void bigmul(int n, int m, unsigned short x[], unsigned short y[]) alpar@9: { int i, j; alpar@9: unsigned int t; alpar@9: xassert(n >= 1); alpar@9: xassert(m >= 1); alpar@9: for (j = 0; j < m; j++) x[j] = 0; alpar@9: for (i = 0; i < n; i++) alpar@9: { if (x[i+m]) alpar@9: { t = 0; alpar@9: for (j = 0; j < m; j++) alpar@9: { t += (unsigned int)x[i+m] * (unsigned int)y[j] + alpar@9: (unsigned int)x[i+j]; alpar@9: x[i+j] = (unsigned short)t; alpar@9: t >>= 16; alpar@9: } alpar@9: x[i+m] = (unsigned short)t; alpar@9: } alpar@9: } alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * bigdiv - divide unsigned integer numbers of arbitrary precision alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glplib.h" alpar@9: * void bigdiv(int n, int m, unsigned short x[], unsigned short y[]); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine bigdiv divides one unsigned integer number of arbitrary alpar@9: * precision by another with the algorithm described in [1]. alpar@9: * alpar@9: * n is the difference between the number of digits of dividend and the alpar@9: * number of digits of divisor, n >= 0. alpar@9: * alpar@9: * m is the number of digits of divisor, m >= 1. alpar@9: * alpar@9: * x is an array containing digits of the dividend in elements x[0], alpar@9: * x[1], ..., x[n+m-1]. alpar@9: * alpar@9: * y is an array containing digits of the divisor in elements y[0], alpar@9: * y[1], ..., y[m-1]. The highest digit y[m-1] must be non-zero. alpar@9: * alpar@9: * On exit n+1 digits of the quotient are stored in elements x[m], alpar@9: * x[m+1], ..., x[n+m], and m digits of the remainder are stored in alpar@9: * elements x[0], x[1], ..., x[m-1]. The array y is changed but then alpar@9: * restored. alpar@9: * alpar@9: * REFERENCES alpar@9: * alpar@9: * 1. D. Knuth. The Art of Computer Programming. Vol. 2: Seminumerical alpar@9: * Algorithms. Stanford University, 1969. */ alpar@9: alpar@9: void bigdiv(int n, int m, unsigned short x[], unsigned short y[]) alpar@9: { int i, j; alpar@9: unsigned int t; alpar@9: unsigned short d, q, r; alpar@9: xassert(n >= 0); alpar@9: xassert(m >= 1); alpar@9: xassert(y[m-1] != 0); alpar@9: /* special case when divisor has the only digit */ alpar@9: if (m == 1) alpar@9: { d = 0; alpar@9: for (i = n; i >= 0; i--) alpar@9: { t = ((unsigned int)d << 16) + (unsigned int)x[i]; alpar@9: x[i+1] = (unsigned short)(t / y[0]); alpar@9: d = (unsigned short)(t % y[0]); alpar@9: } alpar@9: x[0] = d; alpar@9: goto done; alpar@9: } alpar@9: /* multiply dividend and divisor by a normalizing coefficient in alpar@9: order to provide the condition y[m-1] >= base / 2 */ alpar@9: d = (unsigned short)(0x10000 / ((unsigned int)y[m-1] + 1)); alpar@9: if (d == 1) alpar@9: x[n+m] = 0; alpar@9: else alpar@9: { t = 0; alpar@9: for (i = 0; i < n+m; i++) alpar@9: { t += (unsigned int)x[i] * (unsigned int)d; alpar@9: x[i] = (unsigned short)t; alpar@9: t >>= 16; alpar@9: } alpar@9: x[n+m] = (unsigned short)t; alpar@9: t = 0; alpar@9: for (j = 0; j < m; j++) alpar@9: { t += (unsigned int)y[j] * (unsigned int)d; alpar@9: y[j] = (unsigned short)t; alpar@9: t >>= 16; alpar@9: } alpar@9: } alpar@9: /* main loop */ alpar@9: for (i = n; i >= 0; i--) alpar@9: { /* estimate and correct the current digit of quotient */ alpar@9: if (x[i+m] < y[m-1]) alpar@9: { t = ((unsigned int)x[i+m] << 16) + (unsigned int)x[i+m-1]; alpar@9: q = (unsigned short)(t / (unsigned int)y[m-1]); alpar@9: r = (unsigned short)(t % (unsigned int)y[m-1]); alpar@9: if (q == 0) goto putq; else goto test; alpar@9: } alpar@9: q = 0; alpar@9: r = x[i+m-1]; alpar@9: decr: q--; /* if q = 0 then q-- = 0xFFFF */ alpar@9: t = (unsigned int)r + (unsigned int)y[m-1]; alpar@9: r = (unsigned short)t; alpar@9: if (t > 0xFFFF) goto msub; alpar@9: test: t = (unsigned int)y[m-2] * (unsigned int)q; alpar@9: if ((unsigned short)(t >> 16) > r) goto decr; alpar@9: if ((unsigned short)(t >> 16) < r) goto msub; alpar@9: if ((unsigned short)t > x[i+m-2]) goto decr; alpar@9: msub: /* now subtract divisor multiplied by the current digit of alpar@9: quotient from the current dividend */ alpar@9: if (q == 0) goto putq; alpar@9: t = 0; alpar@9: for (j = 0; j < m; j++) alpar@9: { t += (unsigned int)y[j] * (unsigned int)q; alpar@9: if (x[i+j] < (unsigned short)t) t += 0x10000; alpar@9: x[i+j] -= (unsigned short)t; alpar@9: t >>= 16; alpar@9: } alpar@9: if (x[i+m] >= (unsigned short)t) goto putq; alpar@9: /* perform correcting addition, because the current digit of alpar@9: quotient is greater by one than its correct value */ alpar@9: q--; alpar@9: t = 0; alpar@9: for (j = 0; j < m; j++) alpar@9: { t += (unsigned int)x[i+j] + (unsigned int)y[j]; alpar@9: x[i+j] = (unsigned short)t; alpar@9: t >>= 16; alpar@9: } alpar@9: putq: /* store the current digit of quotient */ alpar@9: x[i+m] = q; alpar@9: } alpar@9: /* divide divisor and remainder by the normalizing coefficient in alpar@9: order to restore their original values */ alpar@9: if (d > 1) alpar@9: { t = 0; alpar@9: for (i = m-1; i >= 0; i--) alpar@9: { t = (t << 16) + (unsigned int)x[i]; alpar@9: x[i] = (unsigned short)(t / (unsigned int)d); alpar@9: t %= (unsigned int)d; alpar@9: } alpar@9: t = 0; alpar@9: for (j = m-1; j >= 0; j--) alpar@9: { t = (t << 16) + (unsigned int)y[j]; alpar@9: y[j] = (unsigned short)(t / (unsigned int)d); alpar@9: t %= (unsigned int)d; alpar@9: } alpar@9: } alpar@9: done: return; alpar@9: } alpar@9: alpar@9: /**********************************************************************/ alpar@9: alpar@9: #if 0 alpar@9: #include alpar@9: #include alpar@9: #include alpar@9: #include "glprng.h" alpar@9: alpar@9: #define N_MAX 7 alpar@9: /* maximal number of digits in multiplicand */ alpar@9: alpar@9: #define M_MAX 5 alpar@9: /* maximal number of digits in multiplier */ alpar@9: alpar@9: #define N_TEST 1000000 alpar@9: /* number of tests */ alpar@9: alpar@9: int main(void) alpar@9: { RNG *rand; alpar@9: int d, j, n, m, test; alpar@9: unsigned short x[N_MAX], y[M_MAX], z[N_MAX+M_MAX]; alpar@9: rand = rng_create_rand(); alpar@9: for (test = 1; test <= N_TEST; test++) alpar@9: { /* x[0,...,n-1] := multiplicand */ alpar@9: n = 1 + rng_unif_rand(rand, N_MAX-1); alpar@9: assert(1 <= n && n <= N_MAX); alpar@9: for (j = 0; j < n; j++) alpar@9: { d = rng_unif_rand(rand, 65536); alpar@9: assert(0 <= d && d <= 65535); alpar@9: x[j] = (unsigned short)d; alpar@9: } alpar@9: /* y[0,...,m-1] := multiplier */ alpar@9: m = 1 + rng_unif_rand(rand, M_MAX-1); alpar@9: assert(1 <= m && m <= M_MAX); alpar@9: for (j = 0; j < m; j++) alpar@9: { d = rng_unif_rand(rand, 65536); alpar@9: assert(0 <= d && d <= 65535); alpar@9: y[j] = (unsigned short)d; alpar@9: } alpar@9: if (y[m-1] == 0) y[m-1] = 1; alpar@9: /* z[0,...,n+m-1] := x * y */ alpar@9: for (j = 0; j < n; j++) z[m+j] = x[j]; alpar@9: bigmul(n, m, z, y); alpar@9: /* z[0,...,m-1] := z mod y, z[m,...,n+m-1] := z div y */ alpar@9: bigdiv(n, m, z, y); alpar@9: /* z mod y must be 0 */ alpar@9: for (j = 0; j < m; j++) assert(z[j] == 0); alpar@9: /* z div y must be x */ alpar@9: for (j = 0; j < n; j++) assert(z[m+j] == x[j]); alpar@9: } alpar@9: fprintf(stderr, "%d tests successfully passed\n", N_TEST); alpar@9: rng_delete_rand(rand); alpar@9: return 0; alpar@9: } alpar@9: #endif alpar@9: alpar@9: /* eof */