1 /* mpn_divrem_1 -- mpn by limb division.
3 Copyright 1991, 1993, 1994, 1996, 1998-2000, 2002, 2003 Free Software
6 This file is part of the GNU MP Library.
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of either:
11 * the GNU Lesser General Public License as published by the Free
12 Software Foundation; either version 3 of the License, or (at your
13 option) any later version.
17 * the GNU General Public License as published by the Free Software
18 Foundation; either version 2 of the License, or (at your option) any
21 or both in parallel, as here.
23 The GNU MP Library is distributed in the hope that it will be useful, but
24 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
25 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
28 You should have received copies of the GNU General Public License and the
29 GNU Lesser General Public License along with the GNU MP Library. If not,
30 see https://www.gnu.org/licenses/. */
37 /* The size where udiv_qrnnd_preinv should be used rather than udiv_qrnnd,
38 meaning the quotient size where that should happen, the quotient size
39 being how many udiv divisions will be done.
41 The default is to use preinv always, CPUs where this doesn't suit have
42 tuned thresholds. Note in particular that preinv should certainly be
43 used if that's the only division available (USE_PREINV_ALWAYS). */
45 #ifndef DIVREM_1_NORM_THRESHOLD
46 #define DIVREM_1_NORM_THRESHOLD 0
48 #ifndef DIVREM_1_UNNORM_THRESHOLD
49 #define DIVREM_1_UNNORM_THRESHOLD 0
54 /* If the cpu only has multiply-by-inverse division (eg. alpha), then NORM
55 and UNNORM thresholds are 0 and only the inversion code is included.
57 If multiply-by-inverse is never viable, then NORM and UNNORM thresholds
58 will be MP_SIZE_T_MAX and only the plain division code is included.
60 Otherwise mul-by-inverse is better than plain division above some
61 threshold, and best results are obtained by having code for both present.
63 The main reason for separating the norm and unnorm cases is that not all
64 CPUs give zero for "n0 >> GMP_LIMB_BITS" which would arise in the unnorm
65 code used on an already normalized divisor.
67 If UDIV_NEEDS_NORMALIZATION is false then plain division uses the same
68 non-shifting code for both the norm and unnorm cases, though with
69 different criteria for skipping a division, and with different thresholds
70 of course. And in fact if inversion is never viable, then that simple
71 non-shifting division would be all that's left.
73 The NORM and UNNORM thresholds might not differ much, but if there's
74 going to be separate code for norm and unnorm then it makes sense to have
75 separate thresholds. One thing that's possible is that the
76 mul-by-inverse might be better only for normalized divisors, due to that
77 case not needing variable bit shifts.
79 Notice that the thresholds are tested after the decision to possibly skip
80 one divide step, so they're based on the actual number of divisions done.
82 For the unnorm case, it would be possible to call mpn_lshift to adjust
83 the dividend all in one go (into the quotient space say), rather than
84 limb-by-limb in the loop. This might help if mpn_lshift is a lot faster
85 than what the compiler can generate for EXTRACT. But this is left to CPU
86 specific implementations to consider, especially since EXTRACT isn't on
87 the dependent chain. */
90 mpn_divrem_1 (mp_ptr qp
, mp_size_t qxn
,
91 mp_srcptr up
, mp_size_t un
, mp_limb_t d
)
101 /* FIXME: What's the correct overlap rule when qxn!=0? */
102 ASSERT (MPN_SAME_OR_SEPARATE_P (qp
+qxn
, up
, un
));
110 qp
+= (n
- 1); /* Make qp point at most significant quotient limb */
112 if ((d
& GMP_LIMB_HIGHBIT
) != 0)
116 /* High quotient limb is 0 or 1, skip a divide step. */
118 r
= up
[un
- 1] << GMP_NAIL_BITS
;
127 if (BELOW_THRESHOLD (n
, DIVREM_1_NORM_THRESHOLD
))
130 for (i
= un
- 1; i
>= 0; i
--)
132 n0
= up
[i
] << GMP_NAIL_BITS
;
133 udiv_qrnnd (*qp
, r
, r
, n0
, d
);
137 for (i
= qxn
- 1; i
>= 0; i
--)
139 udiv_qrnnd (*qp
, r
, r
, CNST_LIMB(0), d
);
147 /* Multiply-by-inverse, divisor already normalized. */
149 invert_limb (dinv
, d
);
151 for (i
= un
- 1; i
>= 0; i
--)
153 n0
= up
[i
] << GMP_NAIL_BITS
;
154 udiv_qrnnd_preinv (*qp
, r
, r
, n0
, d
, dinv
);
158 for (i
= qxn
- 1; i
>= 0; i
--)
160 udiv_qrnnd_preinv (*qp
, r
, r
, CNST_LIMB(0), d
, dinv
);
169 /* Most significant bit of divisor == 0. */
172 /* Skip a division if high < divisor (high quotient 0). Testing here
173 before normalizing will still skip as often as possible. */
176 n1
= up
[un
- 1] << GMP_NAIL_BITS
;
179 r
= n1
>> GMP_NAIL_BITS
;
188 if (! UDIV_NEEDS_NORMALIZATION
189 && BELOW_THRESHOLD (n
, DIVREM_1_UNNORM_THRESHOLD
))
192 count_leading_zeros (cnt
, d
);
196 if (UDIV_NEEDS_NORMALIZATION
197 && BELOW_THRESHOLD (n
, DIVREM_1_UNNORM_THRESHOLD
))
202 n1
= up
[un
- 1] << GMP_NAIL_BITS
;
203 r
|= (n1
>> (GMP_LIMB_BITS
- cnt
));
204 for (i
= un
- 2; i
>= 0; i
--)
206 n0
= up
[i
] << GMP_NAIL_BITS
;
207 nshift
= (n1
<< cnt
) | (n0
>> (GMP_NUMB_BITS
- cnt
));
208 udiv_qrnnd (*qp
, r
, r
, nshift
, d
);
213 udiv_qrnnd (*qp
, r
, r
, n1
<< cnt
, d
);
217 for (i
= qxn
- 1; i
>= 0; i
--)
219 udiv_qrnnd (*qp
, r
, r
, CNST_LIMB(0), d
);
227 mp_limb_t dinv
, nshift
;
228 invert_limb (dinv
, d
);
231 n1
= up
[un
- 1] << GMP_NAIL_BITS
;
232 r
|= (n1
>> (GMP_LIMB_BITS
- cnt
));
233 for (i
= un
- 2; i
>= 0; i
--)
235 n0
= up
[i
] << GMP_NAIL_BITS
;
236 nshift
= (n1
<< cnt
) | (n0
>> (GMP_NUMB_BITS
- cnt
));
237 udiv_qrnnd_preinv (*qp
, r
, r
, nshift
, d
, dinv
);
242 udiv_qrnnd_preinv (*qp
, r
, r
, n1
<< cnt
, d
, dinv
);
246 for (i
= qxn
- 1; i
>= 0; i
--)
248 udiv_qrnnd_preinv (*qp
, r
, r
, CNST_LIMB(0), d
, dinv
);