Annotation of OpenXM_contrib/gmp/mpz/divexact.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpz_divexact -- finds quotient when known that quot * den == num && den != 0.
2:
1.1.1.2 ! maekawa 3: Copyright (C) 1991, 1993, 1994, 1995, 1996, 1997, 1998, 2000 Free Software
! 4: Foundation, Inc.
1.1 maekawa 5:
6: This file is part of the GNU MP Library.
7:
8: The GNU MP Library is free software; you can redistribute it and/or modify
1.1.1.2 ! maekawa 9: it under the terms of the GNU Lesser General Public License as published by
! 10: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 11: option) any later version.
12:
13: The GNU MP Library is distributed in the hope that it will be useful, but
14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! maekawa 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 16: License for more details.
17:
1.1.1.2 ! maekawa 18: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 19: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21: MA 02111-1307, USA. */
22:
23: /* Ken Weber (kweber@mat.ufrgs.br, kweber@mcs.kent.edu)
24:
25: Funding for this work has been partially provided by Conselho Nacional
26: de Desenvolvimento Cienti'fico e Tecnolo'gico (CNPq) do Brazil, Grant
27: 301314194-2, and was done while I was a visiting reseacher in the Instituto
28: de Matema'tica at Universidade Federal do Rio Grande do Sul (UFRGS).
29:
30: References:
1.1.1.2 ! maekawa 31: T. Jebelean, An algorithm for exact division, Journal of Symbolic
! 32: Computation, v. 15, 1993, pp. 169-180. */
1.1 maekawa 33:
34: #include "gmp.h"
35: #include "gmp-impl.h"
36: #include "longlong.h"
37:
38: void
39: #if __STDC__
40: mpz_divexact (mpz_ptr quot, mpz_srcptr num, mpz_srcptr den)
41: #else
42: mpz_divexact (quot, num, den)
43: mpz_ptr quot;
44: mpz_srcptr num;
45: mpz_srcptr den;
46: #endif
47: {
48: mp_ptr qp, tp;
49: mp_size_t qsize, tsize;
1.1.1.2 ! maekawa 50: mp_srcptr np, dp;
! 51: mp_size_t nsize, dsize;
1.1 maekawa 52: TMP_DECL (marker);
53:
1.1.1.2 ! maekawa 54: nsize = ABS (num->_mp_size);
! 55: dsize = ABS (den->_mp_size);
! 56:
! 57: qsize = nsize - dsize + 1;
! 58: if (quot->_mp_alloc < qsize)
! 59: _mpz_realloc (quot, qsize);
! 60:
! 61: np = num->_mp_d;
! 62: dp = den->_mp_d;
! 63: qp = quot->_mp_d;
1.1 maekawa 64:
65: if (nsize == 0)
66: {
1.1.1.2 ! maekawa 67: if (dsize == 0)
! 68: DIVIDE_BY_ZERO;
1.1 maekawa 69: quot->_mp_size = 0;
70: return;
71: }
72:
1.1.1.2 ! maekawa 73: if (dsize <= 1)
! 74: {
! 75: if (dsize == 1)
! 76: {
! 77: mpn_divmod_1 (qp, np, nsize, dp[0]);
! 78: qsize -= qp[qsize - 1] == 0;
! 79: quot->_mp_size = (num->_mp_size ^ den->_mp_size) >= 0 ? qsize : -qsize;
! 80: return;
! 81: }
! 82:
! 83: /* Generate divide-by-zero error since dsize == 0. */
! 84: DIVIDE_BY_ZERO;
! 85: }
1.1 maekawa 86:
87: TMP_MARK (marker);
88:
89: /* QUOT <-- NUM/2^r, T <-- DEN/2^r where = r number of twos in DEN. */
90: while (dp[0] == 0)
91: np += 1, nsize -= 1, dp += 1, dsize -= 1;
92: tsize = MIN (qsize, dsize);
1.1.1.2 ! maekawa 93: if ((dp[0] & 1) != 0)
1.1 maekawa 94: {
1.1.1.2 ! maekawa 95: if (quot == den) /* QUOT and DEN overlap. */
1.1 maekawa 96: {
1.1.1.2 ! maekawa 97: tp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);
1.1 maekawa 98: MPN_COPY (tp, dp, tsize);
99: }
100: else
101: tp = (mp_ptr) dp;
1.1.1.2 ! maekawa 102: if (qp != np)
! 103: MPN_COPY_INCR (qp, np, qsize);
1.1 maekawa 104: }
105: else
106: {
1.1.1.2 ! maekawa 107: unsigned int r;
! 108: tp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);
1.1 maekawa 109: count_trailing_zeros (r, dp[0]);
110: mpn_rshift (tp, dp, tsize, r);
111: if (dsize > tsize)
1.1.1.2 ! maekawa 112: tp[tsize - 1] |= dp[tsize] << (BITS_PER_MP_LIMB - r);
1.1 maekawa 113: mpn_rshift (qp, np, qsize, r);
114: if (nsize > qsize)
1.1.1.2 ! maekawa 115: qp[qsize - 1] |= np[qsize] << (BITS_PER_MP_LIMB - r);
1.1 maekawa 116: }
117:
118: /* Now QUOT <-- QUOT/T. */
119: mpn_bdivmod (qp, qp, qsize, tp, tsize, qsize * BITS_PER_MP_LIMB);
120: MPN_NORMALIZE (qp, qsize);
121:
1.1.1.2 ! maekawa 122: quot->_mp_size = (num->_mp_size ^ den->_mp_size) >= 0 ? qsize : -qsize;
1.1 maekawa 123:
124: TMP_FREE (marker);
125: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>