Daily bump.
[official-gcc.git] / libquadmath / printf / divrem.c
blob944d1a0c08c081ded89970e7b7ece73974321f39
1 /* mpn_divrem -- Divide natural numbers, producing both remainder and
2 quotient.
4 Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
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 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.
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.
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, write to
20 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21 MA 02111-1307, USA. */
23 #include <config.h>
24 #include "gmp-impl.h"
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.
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. */
42 mp_limb_t
43 #if __STDC__
44 mpn_divrem (mp_ptr qp, mp_size_t qextra_limbs,
45 mp_ptr np, mp_size_t nsize,
46 mp_srcptr dp, mp_size_t dsize)
47 #else
48 mpn_divrem (qp, qextra_limbs, np, nsize, dp, dsize)
49 mp_ptr qp;
50 mp_size_t qextra_limbs;
51 mp_ptr np;
52 mp_size_t nsize;
53 mp_srcptr dp;
54 mp_size_t dsize;
55 #endif
57 mp_limb_t most_significant_q_limb = 0;
59 switch (dsize)
61 case 0:
62 /* We are asked to divide by zero, so go ahead and do it! (To make
63 the compiler not remove this statement, return the value.) */
64 return 1 / dsize;
66 case 1:
68 mp_size_t i;
69 mp_limb_t n1;
70 mp_limb_t d;
72 d = dp[0];
73 n1 = np[nsize - 1];
75 if (n1 >= d)
77 n1 -= d;
78 most_significant_q_limb = 1;
81 qp += qextra_limbs;
82 for (i = nsize - 2; i >= 0; i--)
83 udiv_qrnnd (qp[i], n1, n1, np[i], d);
84 qp -= qextra_limbs;
86 for (i = qextra_limbs - 1; i >= 0; i--)
87 udiv_qrnnd (qp[i], n1, n1, 0, d);
89 np[0] = n1;
91 break;
93 case 2:
95 mp_size_t i;
96 mp_limb_t n1, n0, n2;
97 mp_limb_t d1, d0;
99 np += nsize - 2;
100 d1 = dp[1];
101 d0 = dp[0];
102 n1 = np[1];
103 n0 = np[0];
105 if (n1 >= d1 && (n1 > d1 || n0 >= d0))
107 sub_ddmmss (n1, n0, n1, n0, d1, d0);
108 most_significant_q_limb = 1;
111 for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--)
113 mp_limb_t q;
114 mp_limb_t r;
116 if (i >= qextra_limbs)
117 np--;
118 else
119 np[0] = 0;
121 if (n1 == d1)
123 /* Q should be either 111..111 or 111..110. Need special
124 treatment of this rare case as normal division would
125 give overflow. */
126 q = ~(mp_limb_t) 0;
128 r = n0 + d1;
129 if (r < d1) /* Carry in the addition? */
131 add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
132 qp[i] = q;
133 continue;
135 n1 = d0 - (d0 != 0);
136 n0 = -d0;
138 else
140 udiv_qrnnd (q, r, n1, n0, d1);
141 umul_ppmm (n1, n0, d0, q);
144 n2 = np[0];
145 q_test:
146 if (n1 > r || (n1 == r && n0 > n2))
148 /* The estimated Q was too large. */
149 q--;
151 sub_ddmmss (n1, n0, n1, n0, 0, d0);
152 r += d1;
153 if (r >= d1) /* If not carry, test Q again. */
154 goto q_test;
157 qp[i] = q;
158 sub_ddmmss (n1, n0, r, n2, n1, n0);
160 np[1] = n1;
161 np[0] = n0;
163 break;
165 default:
167 mp_size_t i;
168 mp_limb_t dX, d1, n0;
170 np += nsize - dsize;
171 dX = dp[dsize - 1];
172 d1 = dp[dsize - 2];
173 n0 = np[dsize - 1];
175 if (n0 >= dX)
177 if (n0 > dX || mpn_cmp (np, dp, dsize - 1) >= 0)
179 mpn_sub_n (np, np, dp, dsize);
180 n0 = np[dsize - 1];
181 most_significant_q_limb = 1;
185 for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--)
187 mp_limb_t q;
188 mp_limb_t n1, n2;
189 mp_limb_t cy_limb;
191 if (i >= qextra_limbs)
193 np--;
194 n2 = np[dsize];
196 else
198 n2 = np[dsize - 1];
199 MPN_COPY_DECR (np + 1, np, dsize);
200 np[0] = 0;
203 if (n0 == dX)
204 /* This might over-estimate q, but it's probably not worth
205 the extra code here to find out. */
206 q = ~(mp_limb_t) 0;
207 else
209 mp_limb_t r;
211 udiv_qrnnd (q, r, n0, np[dsize - 1], dX);
212 umul_ppmm (n1, n0, d1, q);
214 while (n1 > r || (n1 == r && n0 > np[dsize - 2]))
216 q--;
217 r += dX;
218 if (r < dX) /* I.e. "carry in previous addition?" */
219 break;
220 n1 -= n0 < d1;
221 n0 -= d1;
225 /* Possible optimization: We already have (q * n0) and (1 * n1)
226 after the calculation of q. Taking advantage of that, we
227 could make this loop make two iterations less. */
229 cy_limb = mpn_submul_1 (np, dp, dsize, q);
231 if (n2 != cy_limb)
233 mpn_add_n (np, np, dp, dsize);
234 q--;
237 qp[i] = q;
238 n0 = np[dsize - 1];
243 return most_significant_q_limb;