1 /* __mpn_divmod_1(quot_ptr, dividend_ptr, dividend_size, divisor_limb) --
2 Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
3 Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
4 Return the single-limb remainder.
5 There are no constraints on the value of the divisor.
7 QUOT_PTR and DIVIDEND_PTR might point to the same limb.
9 Copyright (C) 1991, 1993, 1994 Free Software Foundation, Inc.
11 This file is part of the GNU MP Library.
13 The GNU MP Library is free software; you can redistribute it and/or modify
14 it under the terms of the GNU Library General Public License as published by
15 the Free Software Foundation; either version 2 of the License, or (at your
16 option) any later version.
18 The GNU MP Library is distributed in the hope that it will be useful, but
19 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
20 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
21 License for more details.
23 You should have received a copy of the GNU Library General Public License
24 along with the GNU MP Library; see the file COPYING.LIB. If not, write to
25 the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
36 #define UDIV_TIME UMUL_TIME
39 /* FIXME: We should be using invert_limb (or invert_normalized_limb)
40 here (not udiv_qrnnd). */
44 __mpn_divmod_1 (mp_ptr quot_ptr
,
45 mp_srcptr dividend_ptr
, mp_size_t dividend_size
,
48 __mpn_divmod_1 (quot_ptr
, dividend_ptr
, dividend_size
, divisor_limb
)
50 mp_srcptr dividend_ptr
;
51 mp_size_t dividend_size
;
59 /* ??? Should this be handled at all? Rely on callers? */
60 if (dividend_size
== 0)
63 /* If multiplication is much faster than division, and the
64 dividend is large, pre-invert the divisor, and use
65 only multiplications in the inner loop. */
67 /* This test should be read:
68 Does it ever help to use udiv_qrnnd_preinv?
69 && Does what we save compensate for the inversion overhead? */
70 if (UDIV_TIME
> (2 * UMUL_TIME
+ 6)
71 && (UDIV_TIME
- (2 * UMUL_TIME
+ 6)) * dividend_size
> UDIV_TIME
)
73 int normalization_steps
;
75 count_leading_zeros (normalization_steps
, divisor_limb
);
76 if (normalization_steps
!= 0)
78 mp_limb divisor_limb_inverted
;
80 divisor_limb
<<= normalization_steps
;
82 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
83 result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
84 most significant bit (with weight 2**N) implicit. */
86 /* Special case for DIVISOR_LIMB == 100...000. */
87 if (divisor_limb
<< 1 == 0)
88 divisor_limb_inverted
= ~(mp_limb
) 0;
90 udiv_qrnnd (divisor_limb_inverted
, dummy
,
91 -divisor_limb
, 0, divisor_limb
);
93 n1
= dividend_ptr
[dividend_size
- 1];
94 r
= n1
>> (BITS_PER_MP_LIMB
- normalization_steps
);
96 /* Possible optimization:
98 && divisor_limb > ((n1 << normalization_steps)
99 | (dividend_ptr[dividend_size - 2] >> ...)))
100 ...one division less... */
102 for (i
= dividend_size
- 2; i
>= 0; i
--)
104 n0
= dividend_ptr
[i
];
105 udiv_qrnnd_preinv (quot_ptr
[i
+ 1], r
, r
,
106 ((n1
<< normalization_steps
)
107 | (n0
>> (BITS_PER_MP_LIMB
- normalization_steps
))),
108 divisor_limb
, divisor_limb_inverted
);
111 udiv_qrnnd_preinv (quot_ptr
[0], r
, r
,
112 n1
<< normalization_steps
,
113 divisor_limb
, divisor_limb_inverted
);
114 return r
>> normalization_steps
;
118 mp_limb divisor_limb_inverted
;
120 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
121 result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
122 most significant bit (with weight 2**N) implicit. */
124 /* Special case for DIVISOR_LIMB == 100...000. */
125 if (divisor_limb
<< 1 == 0)
126 divisor_limb_inverted
= ~(mp_limb
) 0;
128 udiv_qrnnd (divisor_limb_inverted
, dummy
,
129 -divisor_limb
, 0, divisor_limb
);
131 i
= dividend_size
- 1;
134 if (r
>= divisor_limb
)
144 n0
= dividend_ptr
[i
];
145 udiv_qrnnd_preinv (quot_ptr
[i
], r
, r
,
146 n0
, divisor_limb
, divisor_limb_inverted
);
153 if (UDIV_NEEDS_NORMALIZATION
)
155 int normalization_steps
;
157 count_leading_zeros (normalization_steps
, divisor_limb
);
158 if (normalization_steps
!= 0)
160 divisor_limb
<<= normalization_steps
;
162 n1
= dividend_ptr
[dividend_size
- 1];
163 r
= n1
>> (BITS_PER_MP_LIMB
- normalization_steps
);
165 /* Possible optimization:
167 && divisor_limb > ((n1 << normalization_steps)
168 | (dividend_ptr[dividend_size - 2] >> ...)))
169 ...one division less... */
171 for (i
= dividend_size
- 2; i
>= 0; i
--)
173 n0
= dividend_ptr
[i
];
174 udiv_qrnnd (quot_ptr
[i
+ 1], r
, r
,
175 ((n1
<< normalization_steps
)
176 | (n0
>> (BITS_PER_MP_LIMB
- normalization_steps
))),
180 udiv_qrnnd (quot_ptr
[0], r
, r
,
181 n1
<< normalization_steps
,
183 return r
>> normalization_steps
;
186 /* No normalization needed, either because udiv_qrnnd doesn't require
187 it, or because DIVISOR_LIMB is already normalized. */
189 i
= dividend_size
- 1;
192 if (r
>= divisor_limb
)
202 n0
= dividend_ptr
[i
];
203 udiv_qrnnd (quot_ptr
[i
], r
, r
, n0
, divisor_limb
);