Annotation of OpenXM_contrib/gmp/mpn/generic/mod_1.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpn_mod_1(dividend_ptr, dividend_size, divisor_limb) --
2: Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
3: Return the single-limb remainder.
4: There are no constraints on the value of the divisor.
5:
1.1.1.2 ! maekawa 6: Copyright (C) 1991, 1993, 1994, 1999 Free Software Foundation, Inc.
1.1 maekawa 7:
8: This file is part of the GNU MP Library.
9:
10: The GNU MP Library is free software; you can redistribute it and/or modify
1.1.1.2 ! maekawa 11: it under the terms of the GNU Lesser General Public License as published by
! 12: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 13: option) any later version.
14:
15: The GNU MP Library is distributed in the hope that it will be useful, but
16: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! maekawa 17: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 18: License for more details.
19:
1.1.1.2 ! maekawa 20: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 21: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
22: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
23: MA 02111-1307, USA. */
24:
25: #include "gmp.h"
26: #include "gmp-impl.h"
27: #include "longlong.h"
28:
29: #ifndef UMUL_TIME
30: #define UMUL_TIME 1
31: #endif
32:
33: #ifndef UDIV_TIME
34: #define UDIV_TIME UMUL_TIME
35: #endif
36:
37: mp_limb_t
38: #if __STDC__
39: mpn_mod_1 (mp_srcptr dividend_ptr, mp_size_t dividend_size,
40: mp_limb_t divisor_limb)
41: #else
42: mpn_mod_1 (dividend_ptr, dividend_size, divisor_limb)
43: mp_srcptr dividend_ptr;
44: mp_size_t dividend_size;
45: mp_limb_t divisor_limb;
46: #endif
47: {
48: mp_size_t i;
49: mp_limb_t n1, n0, r;
50: int dummy;
51:
52: /* Botch: Should this be handled at all? Rely on callers? */
53: if (dividend_size == 0)
54: return 0;
55:
56: /* If multiplication is much faster than division, and the
57: dividend is large, pre-invert the divisor, and use
58: only multiplications in the inner loop. */
59:
60: /* This test should be read:
61: Does it ever help to use udiv_qrnnd_preinv?
62: && Does what we save compensate for the inversion overhead? */
63: if (UDIV_TIME > (2 * UMUL_TIME + 6)
64: && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME)
65: {
66: int normalization_steps;
67:
68: count_leading_zeros (normalization_steps, divisor_limb);
69: if (normalization_steps != 0)
70: {
71: mp_limb_t divisor_limb_inverted;
72:
73: divisor_limb <<= normalization_steps;
1.1.1.2 ! maekawa 74: invert_limb (divisor_limb_inverted, divisor_limb);
1.1 maekawa 75:
76: n1 = dividend_ptr[dividend_size - 1];
77: r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
78:
79: /* Possible optimization:
80: if (r == 0
81: && divisor_limb > ((n1 << normalization_steps)
82: | (dividend_ptr[dividend_size - 2] >> ...)))
83: ...one division less... */
84:
85: for (i = dividend_size - 2; i >= 0; i--)
86: {
87: n0 = dividend_ptr[i];
88: udiv_qrnnd_preinv (dummy, r, r,
89: ((n1 << normalization_steps)
90: | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
91: divisor_limb, divisor_limb_inverted);
92: n1 = n0;
93: }
94: udiv_qrnnd_preinv (dummy, r, r,
95: n1 << normalization_steps,
96: divisor_limb, divisor_limb_inverted);
97: return r >> normalization_steps;
98: }
99: else
100: {
101: mp_limb_t divisor_limb_inverted;
102:
1.1.1.2 ! maekawa 103: invert_limb (divisor_limb_inverted, divisor_limb);
1.1 maekawa 104:
105: i = dividend_size - 1;
106: r = dividend_ptr[i];
107:
108: if (r >= divisor_limb)
109: r = 0;
110: else
111: i--;
112:
113: for (; i >= 0; i--)
114: {
115: n0 = dividend_ptr[i];
116: udiv_qrnnd_preinv (dummy, r, r,
117: n0, divisor_limb, divisor_limb_inverted);
118: }
119: return r;
120: }
121: }
122: else
123: {
124: if (UDIV_NEEDS_NORMALIZATION)
125: {
126: int normalization_steps;
127:
128: count_leading_zeros (normalization_steps, divisor_limb);
129: if (normalization_steps != 0)
130: {
131: divisor_limb <<= normalization_steps;
132:
133: n1 = dividend_ptr[dividend_size - 1];
134: r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
135:
136: /* Possible optimization:
137: if (r == 0
138: && divisor_limb > ((n1 << normalization_steps)
139: | (dividend_ptr[dividend_size - 2] >> ...)))
140: ...one division less... */
141:
142: for (i = dividend_size - 2; i >= 0; i--)
143: {
144: n0 = dividend_ptr[i];
145: udiv_qrnnd (dummy, r, r,
146: ((n1 << normalization_steps)
147: | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
148: divisor_limb);
149: n1 = n0;
150: }
151: udiv_qrnnd (dummy, r, r,
152: n1 << normalization_steps,
153: divisor_limb);
154: return r >> normalization_steps;
155: }
156: }
157: /* No normalization needed, either because udiv_qrnnd doesn't require
158: it, or because DIVISOR_LIMB is already normalized. */
159:
160: i = dividend_size - 1;
161: r = dividend_ptr[i];
162:
163: if (r >= divisor_limb)
164: r = 0;
165: else
166: i--;
167:
168: for (; i >= 0; i--)
169: {
170: n0 = dividend_ptr[i];
171: udiv_qrnnd (dummy, r, r, n0, divisor_limb);
172: }
173: return r;
174: }
175: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>