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