1/* mpn_divrem -- Divide natural numbers, producing both remainder and
2   quotient.
3
4Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
5
6This file is part of the GNU MP Library.
7
8The GNU MP Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 2.1 of the License, or (at your
11option) any later version.
12
13The GNU MP Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
20the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21MA 02111-1307, USA. */
22
23#include <config.h>
24#include "gmp-impl.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
42mp_limb_t
43#if __STDC__
44mpn_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
48mpn_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
56{
57  mp_limb_t most_significant_q_limb = 0;
58
59  switch (dsize)
60    {
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;
65
66    case 1:
67      {
68	mp_size_t i;
69	mp_limb_t n1;
70	mp_limb_t d;
71
72	d = dp[0];
73	n1 = np[nsize - 1];
74
75	if (n1 >= d)
76	  {
77	    n1 -= d;
78	    most_significant_q_limb = 1;
79	  }
80
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;
85
86	for (i = qextra_limbs - 1; i >= 0; i--)
87	  udiv_qrnnd (qp[i], n1, n1, 0, d);
88
89	np[0] = n1;
90      }
91      break;
92
93    case 2:
94      {
95	mp_size_t i;
96	mp_limb_t n1, n0, n2;
97	mp_limb_t d1, d0;
98
99	np += nsize - 2;
100	d1 = dp[1];
101	d0 = dp[0];
102	n1 = np[1];
103	n0 = np[0];
104
105	if (n1 >= d1 && (n1 > d1 || n0 >= d0))
106	  {
107	    sub_ddmmss (n1, n0, n1, n0, d1, d0);
108	    most_significant_q_limb = 1;
109	  }
110
111	for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--)
112	  {
113	    mp_limb_t q;
114	    mp_limb_t r;
115
116	    if (i >= qextra_limbs)
117	      np--;
118	    else
119	      np[0] = 0;
120
121	    if (n1 == d1)
122	      {
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;
127
128		r = n0 + d1;
129		if (r < d1)	/* Carry in the addition? */
130		  {
131		    add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
132		    qp[i] = q;
133		    continue;
134		  }
135		n1 = d0 - (d0 != 0);
136		n0 = -d0;
137	      }
138	    else
139	      {
140		udiv_qrnnd (q, r, n1, n0, d1);
141		umul_ppmm (n1, n0, d0, q);
142	      }
143
144	    n2 = np[0];
145	  q_test:
146	    if (n1 > r || (n1 == r && n0 > n2))
147	      {
148		/* The estimated Q was too large.  */
149		q--;
150
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;
155	      }
156
157	    qp[i] = q;
158	    sub_ddmmss (n1, n0, r, n2, n1, n0);
159	  }
160	np[1] = n1;
161	np[0] = n0;
162      }
163      break;
164
165    default:
166      {
167	mp_size_t i;
168	mp_limb_t dX, d1, n0;
169
170	np += nsize - dsize;
171	dX = dp[dsize - 1];
172	d1 = dp[dsize - 2];
173	n0 = np[dsize - 1];
174
175	if (n0 >= dX)
176	  {
177	    if (n0 > dX || mpn_cmp (np, dp, dsize - 1) >= 0)
178	      {
179		mpn_sub_n (np, np, dp, dsize);
180		n0 = np[dsize - 1];
181		most_significant_q_limb = 1;
182	      }
183	  }
184
185	for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--)
186	  {
187	    mp_limb_t q;
188	    mp_limb_t n1, n2;
189	    mp_limb_t cy_limb;
190
191	    if (i >= qextra_limbs)
192	      {
193		np--;
194		n2 = np[dsize];
195	      }
196	    else
197	      {
198		n2 = np[dsize - 1];
199		MPN_COPY_DECR (np + 1, np, dsize);
200		np[0] = 0;
201	      }
202
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
208	      {
209		mp_limb_t r;
210
211		udiv_qrnnd (q, r, n0, np[dsize - 1], dX);
212		umul_ppmm (n1, n0, d1, q);
213
214		while (n1 > r || (n1 == r && n0 > np[dsize - 2]))
215		  {
216		    q--;
217		    r += dX;
218		    if (r < dX)	/* I.e. "carry in previous addition?"  */
219		      break;
220		    n1 -= n0 < d1;
221		    n0 -= d1;
222		  }
223	      }
224
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.  */
228
229	    cy_limb = mpn_submul_1 (np, dp, dsize, q);
230
231	    if (n2 != cy_limb)
232	      {
233		mpn_add_n (np, np, dp, dsize);
234		q--;
235	      }
236
237	    qp[i] = q;
238	    n0 = np[dsize - 1];
239	  }
240      }
241    }
242
243  return most_significant_q_limb;
244}
245