Annotation of OpenXM_contrib/gmp/mpz/dmincl.c, Revision 1.1.1.1
1.1 maekawa 1: /* dmincl.c -- include file for tdiv_qr.c, tdiv_r.c.
2:
3: Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
4:
5: This file is part of the GNU MP Library.
6:
7: The GNU MP Library is free software; you can redistribute it and/or modify
8: it under the terms of the GNU Library General Public License as published by
9: the Free Software Foundation; either version 2 of the License, or (at your
10: option) any later version.
11:
12: The GNU MP Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15: License for more details.
16:
17: You should have received a copy of the GNU Library General Public License
18: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: /* If den == quot, den needs temporary storage.
23: If den == rem, den needs temporary storage.
24: If num == quot, num needs temporary storage.
25: If den has temporary storage, it can be normalized while being copied,
26: i.e no extra storage should be allocated. */
27:
28: /* This is the function body of mdiv, mpz_divmod, and mpz_mod.
29:
30: If COMPUTE_QUOTIENT is defined, the quotient is put in the MP_INT
31: object quot, otherwise that variable is not referenced at all.
32:
33: The remainder is always computed, and the result is put in the MP_INT
34: object rem. */
35:
36: {
37: mp_ptr np, dp;
38: mp_ptr qp, rp;
39: mp_size_t nsize = num->_mp_size;
40: mp_size_t dsize = den->_mp_size;
41: mp_size_t qsize, rsize;
42: mp_size_t sign_remainder = nsize;
43: #ifdef COMPUTE_QUOTIENT
44: mp_size_t sign_quotient = nsize ^ dsize;
45: #endif
46: unsigned normalization_steps;
47: mp_limb_t q_limb;
48: TMP_DECL (marker);
49:
50: nsize = ABS (nsize);
51: dsize = ABS (dsize);
52:
53: /* Ensure space is enough for quotient and remainder. */
54:
55: /* We need space for an extra limb in the remainder, because it's
56: up-shifted (normalized) below. */
57: rsize = nsize + 1;
58: if (rem->_mp_alloc < rsize)
59: _mpz_realloc (rem, rsize);
60:
61: qsize = rsize - dsize; /* qsize cannot be bigger than this. */
62: if (qsize <= 0)
63: {
64: if (num != rem)
65: {
66: rem->_mp_size = num->_mp_size;
67: MPN_COPY (rem->_mp_d, num->_mp_d, nsize);
68: }
69: #ifdef COMPUTE_QUOTIENT
70: /* This needs to follow the assignment to rem, in case the
71: numerator and quotient are the same. */
72: quot->_mp_size = 0;
73: #endif
74: return;
75: }
76:
77: #ifdef COMPUTE_QUOTIENT
78: if (quot->_mp_alloc < qsize)
79: _mpz_realloc (quot, qsize);
80: #endif
81:
82: /* Read pointers here, when reallocation is finished. */
83: np = num->_mp_d;
84: dp = den->_mp_d;
85: rp = rem->_mp_d;
86:
87: /* Optimize division by a single-limb divisor. */
88: if (dsize == 1)
89: {
90: mp_limb_t rlimb;
91: #ifdef COMPUTE_QUOTIENT
92: qp = quot->_mp_d;
93: rlimb = mpn_divmod_1 (qp, np, nsize, dp[0]);
94: qsize -= qp[qsize - 1] == 0;
95: quot->_mp_size = sign_quotient >= 0 ? qsize : -qsize;
96: #else
97: rlimb = mpn_mod_1 (np, nsize, dp[0]);
98: #endif
99: rp[0] = rlimb;
100: rsize = rlimb != 0;
101: rem->_mp_size = sign_remainder >= 0 ? rsize : -rsize;
102: return;
103: }
104:
105: TMP_MARK (marker);
106:
107: #ifdef COMPUTE_QUOTIENT
108: qp = quot->_mp_d;
109:
110: /* Make sure QP and NP point to different objects. Otherwise the
111: numerator would be gradually overwritten by the quotient limbs. */
112: if (qp == np)
113: {
114: /* Copy NP object to temporary space. */
115: np = (mp_ptr) TMP_ALLOC (nsize * BYTES_PER_MP_LIMB);
116: MPN_COPY (np, qp, nsize);
117: }
118:
119: #else
120: /* Put quotient at top of remainder. */
121: qp = rp + dsize;
122: #endif
123:
124: count_leading_zeros (normalization_steps, dp[dsize - 1]);
125:
126: /* Normalize the denominator, i.e. make its most significant bit set by
127: shifting it NORMALIZATION_STEPS bits to the left. Also shift the
128: numerator the same number of steps (to keep the quotient the same!). */
129: if (normalization_steps != 0)
130: {
131: mp_ptr tp;
132: mp_limb_t nlimb;
133:
134: /* Shift up the denominator setting the most significant bit of
135: the most significant word. Use temporary storage not to clobber
136: the original contents of the denominator. */
137: tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
138: mpn_lshift (tp, dp, dsize, normalization_steps);
139: dp = tp;
140:
141: /* Shift up the numerator, possibly introducing a new most
142: significant word. Move the shifted numerator in the remainder
143: meanwhile. */
144: nlimb = mpn_lshift (rp, np, nsize, normalization_steps);
145: if (nlimb != 0)
146: {
147: rp[nsize] = nlimb;
148: rsize = nsize + 1;
149: }
150: else
151: rsize = nsize;
152: }
153: else
154: {
155: /* The denominator is already normalized, as required. Copy it to
156: temporary space if it overlaps with the quotient or remainder. */
157: #ifdef COMPUTE_QUOTIENT
158: if (dp == rp || dp == qp)
159: #else
160: if (dp == rp)
161: #endif
162: {
163: mp_ptr tp;
164:
165: tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
166: MPN_COPY (tp, dp, dsize);
167: dp = tp;
168: }
169:
170: /* Move the numerator to the remainder. */
171: if (rp != np)
172: MPN_COPY (rp, np, nsize);
173:
174: rsize = nsize;
175: }
176:
177: q_limb = mpn_divmod (qp, rp, rsize, dp, dsize);
178:
179: #ifdef COMPUTE_QUOTIENT
180: qsize = rsize - dsize;
181: if (q_limb)
182: {
183: qp[qsize] = q_limb;
184: qsize += 1;
185: }
186:
187: quot->_mp_size = sign_quotient >= 0 ? qsize : -qsize;
188: #endif
189:
190: rsize = dsize;
191: MPN_NORMALIZE (rp, rsize);
192:
193: if (normalization_steps != 0 && rsize != 0)
194: {
195: mpn_rshift (rp, rp, rsize, normalization_steps);
196: rsize -= rp[rsize - 1] == 0;
197: }
198:
199: rem->_mp_size = sign_remainder >= 0 ? rsize : -rsize;
200: TMP_FREE (marker);
201: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dmincl.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz