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.
6 Copyright (C) 1991, 1993, 1994, Free Software Foundation, Inc.
8 This file is part of the GNU MP Library.
10 The GNU MP Library is free software; you can redistribute it and/or modify
11 it under the terms of the GNU Library General Public License as published by
12 the Free Software Foundation; either version 2 of the License, or (at your
13 option) any later version.
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
17 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
18 License for more details.
20 You should have received a copy of the GNU Library General Public License
21 along with the GNU MP Library; see the file COPYING.LIB. If not, write to
22 the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
33 #define UDIV_TIME UMUL_TIME
36 /* FIXME: We should be using invert_limb (or invert_normalized_limb)
37 here (not udiv_qrnnd). */
41 __mpn_mod_1 (mp_srcptr dividend_ptr
, mp_size_t dividend_size
,
44 __mpn_mod_1 (dividend_ptr
, dividend_size
, divisor_limb
)
45 mp_srcptr dividend_ptr
;
46 mp_size_t dividend_size
;
54 /* Botch: Should this be handled at all? Rely on callers? */
55 if (dividend_size
== 0)
58 /* If multiplication is much faster than division, and the
59 dividend is large, pre-invert the divisor, and use
60 only multiplications in the inner loop. */
62 /* This test should be read:
63 Does it ever help to use udiv_qrnnd_preinv?
64 && Does what we save compensate for the inversion overhead? */
65 if (UDIV_TIME
> (2 * UMUL_TIME
+ 6)
66 && (UDIV_TIME
- (2 * UMUL_TIME
+ 6)) * dividend_size
> UDIV_TIME
)
68 int normalization_steps
;
70 count_leading_zeros (normalization_steps
, divisor_limb
);
71 if (normalization_steps
!= 0)
73 mp_limb divisor_limb_inverted
;
75 divisor_limb
<<= normalization_steps
;
77 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
78 result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
79 most significant bit (with weight 2**N) implicit. */
81 /* Special case for DIVISOR_LIMB == 100...000. */
82 if (divisor_limb
<< 1 == 0)
83 divisor_limb_inverted
= ~(mp_limb
) 0;
85 udiv_qrnnd (divisor_limb_inverted
, dummy
,
86 -divisor_limb
, 0, divisor_limb
);
88 n1
= dividend_ptr
[dividend_size
- 1];
89 r
= n1
>> (BITS_PER_MP_LIMB
- normalization_steps
);
91 /* Possible optimization:
93 && divisor_limb > ((n1 << normalization_steps)
94 | (dividend_ptr[dividend_size - 2] >> ...)))
95 ...one division less... */
97 for (i
= dividend_size
- 2; i
>= 0; i
--)
100 udiv_qrnnd_preinv (dummy
, r
, r
,
101 ((n1
<< normalization_steps
)
102 | (n0
>> (BITS_PER_MP_LIMB
- normalization_steps
))),
103 divisor_limb
, divisor_limb_inverted
);
106 udiv_qrnnd_preinv (dummy
, r
, r
,
107 n1
<< normalization_steps
,
108 divisor_limb
, divisor_limb_inverted
);
109 return r
>> normalization_steps
;
113 mp_limb divisor_limb_inverted
;
115 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
116 result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
117 most significant bit (with weight 2**N) implicit. */
119 /* Special case for DIVISOR_LIMB == 100...000. */
120 if (divisor_limb
<< 1 == 0)
121 divisor_limb_inverted
= ~(mp_limb
) 0;
123 udiv_qrnnd (divisor_limb_inverted
, dummy
,
124 -divisor_limb
, 0, divisor_limb
);
126 i
= dividend_size
- 1;
129 if (r
>= divisor_limb
)
136 n0
= dividend_ptr
[i
];
137 udiv_qrnnd_preinv (dummy
, r
, r
,
138 n0
, divisor_limb
, divisor_limb_inverted
);
145 if (UDIV_NEEDS_NORMALIZATION
)
147 int normalization_steps
;
149 count_leading_zeros (normalization_steps
, divisor_limb
);
150 if (normalization_steps
!= 0)
152 divisor_limb
<<= normalization_steps
;
154 n1
= dividend_ptr
[dividend_size
- 1];
155 r
= n1
>> (BITS_PER_MP_LIMB
- normalization_steps
);
157 /* Possible optimization:
159 && divisor_limb > ((n1 << normalization_steps)
160 | (dividend_ptr[dividend_size - 2] >> ...)))
161 ...one division less... */
163 for (i
= dividend_size
- 2; i
>= 0; i
--)
165 n0
= dividend_ptr
[i
];
166 udiv_qrnnd (dummy
, r
, r
,
167 ((n1
<< normalization_steps
)
168 | (n0
>> (BITS_PER_MP_LIMB
- normalization_steps
))),
172 udiv_qrnnd (dummy
, r
, r
,
173 n1
<< normalization_steps
,
175 return r
>> normalization_steps
;
178 /* No normalization needed, either because udiv_qrnnd doesn't require
179 it, or because DIVISOR_LIMB is already normalized. */
181 i
= dividend_size
- 1;
184 if (r
>= divisor_limb
)
191 n0
= dividend_ptr
[i
];
192 udiv_qrnnd (dummy
, r
, r
, n0
, divisor_limb
);