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