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.
6
7 QUOT_PTR and DIVIDEND_PTR might point to the same limb.
8
9 Copyright (C) 1991-2023 Free Software Foundation, Inc.
10
11 This file is part of the GNU MP Library.
12
13 The GNU MP Library is free software; you can redistribute it and/or modify
14 it under the terms of the GNU Lesser General Public License as published by
15 the Free Software Foundation; either version 2.1 of the License, or (at your
16 option) any later version.
17
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 Lesser General Public
21 License for more details.
22
23 You should have received a copy of the GNU Lesser General Public License
24 along with the GNU MP Library; see the file COPYING.LIB. If not, see
25 <https://www.gnu.org/licenses/>. */
26
27 #include <gmp.h>
28 #include "gmp-impl.h"
29 #include "longlong.h"
30
31 #ifndef UMUL_TIME
32 #define UMUL_TIME 1
33 #endif
34
35 #ifndef UDIV_TIME
36 #define UDIV_TIME UMUL_TIME
37 #endif
38
39 /* FIXME: We should be using invert_limb (or invert_normalized_limb)
40 here (not udiv_qrnnd). */
41
42 mp_limb_t
43 mpn_divmod_1 (mp_ptr quot_ptr,
44 mp_srcptr dividend_ptr, mp_size_t dividend_size,
45 mp_limb_t divisor_limb)
46 {
47 mp_size_t i;
48 mp_limb_t n1, n0, r;
49 mp_limb_t dummy __attribute__ ((unused));
50
51 /* ??? Should this be handled at all? Rely on callers? */
52 if (dividend_size == 0)
53 return 0;
54
55 /* If multiplication is much faster than division, and the
56 dividend is large, pre-invert the divisor, and use
57 only multiplications in the inner loop. */
58
59 /* This test should be read:
60 Does it ever help to use udiv_qrnnd_preinv?
61 && Does what we save compensate for the inversion overhead? */
62 if (UDIV_TIME > (2 * UMUL_TIME + 6)
63 && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME)
64 {
65 int normalization_steps;
66
67 count_leading_zeros (normalization_steps, divisor_limb);
68 if (normalization_steps != 0)
69 {
70 mp_limb_t divisor_limb_inverted;
71
72 divisor_limb <<= normalization_steps;
73
74 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
75 result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
76 most significant bit (with weight 2**N) implicit. */
77
78 /* Special case for DIVISOR_LIMB == 100...000. */
79 if (divisor_limb << 1 == 0)
80 divisor_limb_inverted = ~(mp_limb_t) 0;
81 else
82 udiv_qrnnd (divisor_limb_inverted, dummy,
83 -divisor_limb, 0, divisor_limb);
84
85 n1 = dividend_ptr[dividend_size - 1];
86 r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
87
88 /* Possible optimization:
89 if (r == 0
90 && divisor_limb > ((n1 << normalization_steps)
91 | (dividend_ptr[dividend_size - 2] >> ...)))
92 ...one division less... */
93
94 for (i = dividend_size - 2; i >= 0; i--)
95 {
96 n0 = dividend_ptr[i];
97 udiv_qrnnd_preinv (quot_ptr[i + 1], r, r,
98 ((n1 << normalization_steps)
99 | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
100 divisor_limb, divisor_limb_inverted);
101 n1 = n0;
102 }
103 udiv_qrnnd_preinv (quot_ptr[0], r, r,
104 n1 << normalization_steps,
105 divisor_limb, divisor_limb_inverted);
106 return r >> normalization_steps;
107 }
108 else
109 {
110 mp_limb_t divisor_limb_inverted;
111
112 /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
113 result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
114 most significant bit (with weight 2**N) implicit. */
115
116 /* Special case for DIVISOR_LIMB == 100...000. */
117 if (divisor_limb << 1 == 0)
118 divisor_limb_inverted = ~(mp_limb_t) 0;
119 else
120 udiv_qrnnd (divisor_limb_inverted, dummy,
121 -divisor_limb, 0, divisor_limb);
122
123 i = dividend_size - 1;
124 r = dividend_ptr[i];
125
126 if (r >= divisor_limb)
127 r = 0;
128 else
129 {
130 quot_ptr[i] = 0;
131 i--;
132 }
133
134 for (; i >= 0; i--)
135 {
136 n0 = dividend_ptr[i];
137 udiv_qrnnd_preinv (quot_ptr[i], r, r,
138 n0, divisor_limb, divisor_limb_inverted);
139 }
140 return r;
141 }
142 }
143 else
144 {
145 if (UDIV_NEEDS_NORMALIZATION)
146 {
147 int normalization_steps;
148
149 count_leading_zeros (normalization_steps, divisor_limb);
150 if (normalization_steps != 0)
151 {
152 divisor_limb <<= normalization_steps;
153
154 n1 = dividend_ptr[dividend_size - 1];
155 r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
156
157 /* Possible optimization:
158 if (r == 0
159 && divisor_limb > ((n1 << normalization_steps)
160 | (dividend_ptr[dividend_size - 2] >> ...)))
161 ...one division less... */
162
163 for (i = dividend_size - 2; i >= 0; i--)
164 {
165 n0 = dividend_ptr[i];
166 udiv_qrnnd (quot_ptr[i + 1], r, r,
167 ((n1 << normalization_steps)
168 | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
169 divisor_limb);
170 n1 = n0;
171 }
172 udiv_qrnnd (quot_ptr[0], r, r,
173 n1 << normalization_steps,
174 divisor_limb);
175 return r >> normalization_steps;
176 }
177 }
178 /* No normalization needed, either because udiv_qrnnd doesn't require
179 it, or because DIVISOR_LIMB is already normalized. */
180
181 i = dividend_size - 1;
182 r = dividend_ptr[i];
183
184 if (r >= divisor_limb)
185 r = 0;
186 else
187 {
188 quot_ptr[i] = 0;
189 i--;
190 }
191
192 for (; i >= 0; i--)
193 {
194 n0 = dividend_ptr[i];
195 udiv_qrnnd (quot_ptr[i], r, r, n0, divisor_limb);
196 }
197 return r;
198 }
199 }