Annotation of OpenXM_contrib/gmp/mpn/generic/sb_divrem_mn.c, Revision 1.1.1.1
1.1 maekawa 1: /* mpn_sb_divrem_mn -- Divide natural numbers, producing both remainder and
2: quotient.
3:
4: THE FUNCTIONS IN THIS FILE ARE INTERNAL FUNCTIONS WITH MUTABLE
5: INTERFACES. IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.
6: IN FACT, IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A
7: FUTURE GNU MP RELEASE.
8:
9:
10: Copyright (C) 1993, 1994, 1995, 1996, 2000 Free Software Foundation, Inc.
11:
12: This file is part of the GNU MP Library.
13:
14: The GNU MP Library is free software; you can redistribute it and/or modify
15: it under the terms of the GNU Lesser General Public License as published by
16: the Free Software Foundation; either version 2.1 of the License, or (at your
17: option) any later version.
18:
19: The GNU MP Library is distributed in the hope that it will be useful, but
20: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
21: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
22: License for more details.
23:
24: You should have received a copy of the GNU Lesser General Public License
25: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
26: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
27: MA 02111-1307, USA. */
28:
29: #include "gmp.h"
30: #include "gmp-impl.h"
31: #include "longlong.h"
32:
33: /* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
34: the NSIZE-DSIZE least significant quotient limbs at QP
35: and the DSIZE long remainder at NP. If QEXTRA_LIMBS is
36: non-zero, generate that many fraction bits and append them after the
37: other quotient limbs.
38: Return the most significant limb of the quotient, this is always 0 or 1.
39:
40: Preconditions:
41: 0. NSIZE >= DSIZE.
42: 1. The most significant bit of the divisor must be set.
43: 2. QP must either not overlap with the input operands at all, or
44: QP + DSIZE >= NP must hold true. (This means that it's
45: possible to put the quotient in the high part of NUM, right after the
46: remainder in NUM.
47: 3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero.
48: 4. DSIZE >= 2. */
49:
50:
51: #define PREINVERT_VIABLE \
52: (UDIV_TIME > 2 * UMUL_TIME + 6 /* && ! TARGET_REGISTER_STARVED */)
53:
54: mp_limb_t
55: #if __STDC__
56: mpn_sb_divrem_mn (mp_ptr qp,
57: mp_ptr np, mp_size_t nsize,
58: mp_srcptr dp, mp_size_t dsize)
59: #else
60: mpn_sb_divrem_mn (qp, np, nsize, dp, dsize)
61: mp_ptr qp;
62: mp_ptr np;
63: mp_size_t nsize;
64: mp_srcptr dp;
65: mp_size_t dsize;
66: #endif
67: {
68: mp_limb_t most_significant_q_limb = 0;
69: mp_size_t i;
70: mp_limb_t dx, d1, n0;
71: mp_limb_t dxinv;
72: int have_preinv;
73:
74: ASSERT_ALWAYS (dsize > 2);
75:
76: np += nsize - dsize;
77: dx = dp[dsize - 1];
78: d1 = dp[dsize - 2];
79: n0 = np[dsize - 1];
80:
81: if (n0 >= dx)
82: {
83: if (n0 > dx || mpn_cmp (np, dp, dsize - 1) >= 0)
84: {
85: mpn_sub_n (np, np, dp, dsize);
86: most_significant_q_limb = 1;
87: }
88: }
89:
90: /* If multiplication is much faster than division, preinvert the
91: most significant divisor limb before entering the loop. */
92: if (PREINVERT_VIABLE)
93: {
94: have_preinv = 0;
95: if ((UDIV_TIME - (2 * UMUL_TIME + 6)) * (nsize - dsize) > UDIV_TIME)
96: {
97: invert_limb (dxinv, dx);
98: have_preinv = 1;
99: }
100: }
101:
102: for (i = nsize - dsize - 1; i >= 0; i--)
103: {
104: mp_limb_t q;
105: mp_limb_t nx;
106: mp_limb_t cy_limb;
107:
108: nx = np[dsize - 1];
109: np--;
110:
111: if (nx == dx)
112: {
113: /* This might over-estimate q, but it's probably not worth
114: the extra code here to find out. */
115: q = ~(mp_limb_t) 0;
116:
117: #if 1
118: cy_limb = mpn_submul_1 (np, dp, dsize, q);
119: #else
120: /* This should be faster on many machines */
121: cy_limb = mpn_sub_n (np + 1, np + 1, dp, dsize);
122: cy = mpn_add_n (np, np, dp, dsize);
123: np[dsize] += cy;
124: #endif
125:
126: if (nx != cy_limb)
127: {
128: mpn_add_n (np, np, dp, dsize);
129: q--;
130: }
131:
132: qp[i] = q;
133: }
134: else
135: {
136: mp_limb_t rx, r1, r0, p1, p0;
137:
138: /* "workaround" avoids a problem with gcc 2.7.2.3 i386 register
139: usage when np[dsize-1] is used in an asm statement like
140: umul_ppmm in udiv_qrnnd_preinv. The symptom is seg faults due
141: to registers being clobbered. gcc 2.95 i386 doesn't have the
142: problem. */
143: {
144: mp_limb_t workaround = np[dsize - 1];
145: if (PREINVERT_VIABLE && have_preinv)
146: udiv_qrnnd_preinv (q, r1, nx, workaround, dx, dxinv);
147: else
148: udiv_qrnnd (q, r1, nx, workaround, dx);
149: }
150: umul_ppmm (p1, p0, d1, q);
151:
152: r0 = np[dsize - 2];
153: rx = 0;
154: if (r1 < p1 || (r1 == p1 && r0 < p0))
155: {
156: p1 -= p0 < d1;
157: p0 -= d1;
158: q--;
159: r1 += dx;
160: rx = r1 < dx;
161: }
162:
163: p1 += r0 < p0; /* cannot carry! */
164: rx -= r1 < p1; /* may become 11..1 if q is still too large */
165: r1 -= p1;
166: r0 -= p0;
167:
168: cy_limb = mpn_submul_1 (np, dp, dsize - 2, q);
169:
170: {
171: mp_limb_t cy1, cy2;
172: cy1 = r0 < cy_limb;
173: r0 -= cy_limb;
174: cy2 = r1 < cy1;
175: r1 -= cy1;
176: np[dsize - 1] = r1;
177: np[dsize - 2] = r0;
178: if (cy2 != rx)
179: {
180: mpn_add_n (np, np, dp, dsize);
181: q--;
182: }
183: }
184: qp[i] = q;
185: }
186: }
187:
188: /* ______ ______ ______
189: |__rx__|__r1__|__r0__| partial remainder
190: ______ ______
191: - |__p1__|__p0__| partial product to subtract
192: ______ ______
193: - |______|cylimb|
194:
195: rx is -1, 0 or 1. If rx=1, then q is correct (it should match
196: carry out). If rx=-1 then q is too large. If rx=0, then q might
197: be too large, but it is most likely correct.
198: */
199:
200: return most_significant_q_limb;
201: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sb_divrem_mn.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / generic