1 /* mpn_divrem -- Divide natural numbers, producing both remainder and
2 quotient.
3
4 Copyright (C) 1993-2023 Free Software Foundation, Inc.
5
6 This file is part of the GNU MP Library.
7
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 2.1 of the License, or (at your
11 option) any later version.
12
13 The GNU MP Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MP Library; see the file COPYING.LIB. If not, see
20 <https://www.gnu.org/licenses/>. */
21
22 #include <gmp.h>
23 #include "gmp-impl.h"
24 #include "longlong.h"
25
26 /* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
27 the NSIZE-DSIZE least significant quotient limbs at QP
28 and the DSIZE long remainder at NP. If QEXTRA_LIMBS is
29 non-zero, generate that many fraction bits and append them after the
30 other quotient limbs.
31 Return the most significant limb of the quotient, this is always 0 or 1.
32
33 Preconditions:
34 0. NSIZE >= DSIZE.
35 1. The most significant bit of the divisor must be set.
36 2. QP must either not overlap with the input operands at all, or
37 QP + DSIZE >= NP must hold true. (This means that it's
38 possible to put the quotient in the high part of NUM, right after the
39 remainder in NUM.
40 3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero. */
41
42 mp_limb_t
43 mpn_divrem (mp_ptr qp, mp_size_t qextra_limbs,
44 mp_ptr np, mp_size_t nsize,
45 mp_srcptr dp, mp_size_t dsize)
46 {
47 mp_limb_t most_significant_q_limb = 0;
48
49 switch (dsize)
50 {
51 case 0:
52 /* We are asked to divide by zero, so go ahead and do it! (To make
53 the compiler not remove this statement, return the value.) */
54 return 1 / dsize;
55
56 case 1:
57 {
58 mp_size_t i;
59 mp_limb_t n1;
60 mp_limb_t d;
61
62 d = dp[0];
63 n1 = np[nsize - 1];
64
65 if (n1 >= d)
66 {
67 n1 -= d;
68 most_significant_q_limb = 1;
69 }
70
71 qp += qextra_limbs;
72 for (i = nsize - 2; i >= 0; i--)
73 udiv_qrnnd (qp[i], n1, n1, np[i], d);
74 qp -= qextra_limbs;
75
76 for (i = qextra_limbs - 1; i >= 0; i--)
77 udiv_qrnnd (qp[i], n1, n1, 0, d);
78
79 np[0] = n1;
80 }
81 break;
82
83 case 2:
84 {
85 mp_size_t i;
86 mp_limb_t n1, n0, n2;
87 mp_limb_t d1, d0;
88
89 np += nsize - 2;
90 d1 = dp[1];
91 d0 = dp[0];
92 n1 = np[1];
93 n0 = np[0];
94
95 if (n1 >= d1 && (n1 > d1 || n0 >= d0))
96 {
97 sub_ddmmss (n1, n0, n1, n0, d1, d0);
98 most_significant_q_limb = 1;
99 }
100
101 for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--)
102 {
103 mp_limb_t q;
104 mp_limb_t r;
105
106 if (i >= qextra_limbs)
107 np--;
108 else
109 np[0] = 0;
110
111 if (n1 == d1)
112 {
113 /* Q should be either 111..111 or 111..110. Need special
114 treatment of this rare case as normal division would
115 give overflow. */
116 q = ~(mp_limb_t) 0;
117
118 r = n0 + d1;
119 if (r < d1) /* Carry in the addition? */
120 {
121 add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
122 qp[i] = q;
123 continue;
124 }
125 n1 = d0 - (d0 != 0);
126 n0 = -d0;
127 }
128 else
129 {
130 udiv_qrnnd (q, r, n1, n0, d1);
131 umul_ppmm (n1, n0, d0, q);
132 }
133
134 n2 = np[0];
135 q_test:
136 if (n1 > r || (n1 == r && n0 > n2))
137 {
138 /* The estimated Q was too large. */
139 q--;
140
141 sub_ddmmss (n1, n0, n1, n0, 0, d0);
142 r += d1;
143 if (r >= d1) /* If not carry, test Q again. */
144 goto q_test;
145 }
146
147 qp[i] = q;
148 sub_ddmmss (n1, n0, r, n2, n1, n0);
149 }
150 np[1] = n1;
151 np[0] = n0;
152 }
153 break;
154
155 default:
156 {
157 mp_size_t i;
158 mp_limb_t dX, d1, n0;
159
160 np += nsize - dsize;
161 dX = dp[dsize - 1];
162 d1 = dp[dsize - 2];
163 n0 = np[dsize - 1];
164
165 if (n0 >= dX)
166 {
167 if (n0 > dX || mpn_cmp (np, dp, dsize - 1) >= 0)
168 {
169 mpn_sub_n (np, np, dp, dsize);
170 n0 = np[dsize - 1];
171 most_significant_q_limb = 1;
172 }
173 }
174
175 for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--)
176 {
177 mp_limb_t q;
178 mp_limb_t n1, n2;
179 mp_limb_t cy_limb;
180
181 if (i >= qextra_limbs)
182 {
183 np--;
184 n2 = np[dsize];
185 }
186 else
187 {
188 n2 = np[dsize - 1];
189 MPN_COPY_DECR (np + 1, np, dsize);
190 np[0] = 0;
191 }
192
193 if (n0 == dX)
194 /* This might over-estimate q, but it's probably not worth
195 the extra code here to find out. */
196 q = ~(mp_limb_t) 0;
197 else
198 {
199 mp_limb_t r;
200
201 udiv_qrnnd (q, r, n0, np[dsize - 1], dX);
202 umul_ppmm (n1, n0, d1, q);
203
204 while (n1 > r || (n1 == r && n0 > np[dsize - 2]))
205 {
206 q--;
207 r += dX;
208 if (r < dX) /* I.e. "carry in previous addition?" */
209 break;
210 n1 -= n0 < d1;
211 n0 -= d1;
212 }
213 }
214
215 /* Possible optimization: We already have (q * n0) and (1 * n1)
216 after the calculation of q. Taking advantage of that, we
217 could make this loop make two iterations less. */
218
219 cy_limb = mpn_submul_1 (np, dp, dsize, q);
220
221 if (n2 != cy_limb)
222 {
223 mpn_add_n (np, np, dp, dsize);
224 q--;
225 }
226
227 qp[i] = q;
228 n0 = np[dsize - 1];
229 }
230 }
231 }
232
233 return most_significant_q_limb;
234 }