.
[kopensolaris-gnu/glibc.git] / sysdeps / generic / divmod.c
1 /* __mpn_divmod -- Divide natural numbers, producing both remainder and
2    quotient.
3
4 Copyright (C) 1993, 1994 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 Library General Public License as published by
10 the Free Software Foundation; either version 2 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 Library General Public
16 License for more details.
17
18 You should have received a copy of the GNU Library 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., 675 Mass Ave, Cambridge, MA 02139, USA. */
21
22 #include "gmp.h"
23 #include "gmp-impl.h"
24 #include "longlong.h"
25
26 /* Divide num (NUM_PTR/NUM_SIZE) by den (DEN_PTR/DEN_SIZE) and write
27    the NUM_SIZE-DEN_SIZE least significant quotient limbs at QUOT_PTR
28    and the DEN_SIZE long remainder at NUM_PTR.
29    Return the most significant limb of the quotient, this is always 0 or 1.
30
31    Argument constraints:
32    1. The most significant bit of the divisor must be set.
33    2. QUOT_PTR must either not overlap with the input operands at all, or
34       QUOT_PTR + DEN_SIZE >= NUM_PTR must hold true.  (This means that it's
35       possible to put the quotient in the high part of NUM, right after the
36       remainder in NUM.  */
37
38 mp_limb
39 #if __STDC__
40 __mpn_divmod (mp_ptr quot_ptr,
41               mp_ptr num_ptr, mp_size_t num_size,
42               mp_srcptr den_ptr, mp_size_t den_size)
43 #else
44 __mpn_divmod (quot_ptr, num_ptr, num_size, den_ptr, den_size)
45      mp_ptr quot_ptr;
46      mp_ptr num_ptr;
47      mp_size_t num_size;
48      mp_srcptr den_ptr;
49      mp_size_t den_size;
50 #endif
51 {
52   mp_limb most_significant_q_limb = 0;
53
54   switch (den_size)
55     {
56     case 0:
57       /* We are asked to divide by zero, so go ahead and do it!  (To make
58          the compiler not remove this statement, return the value.)  */
59       return 1 / den_size;
60
61     case 1:
62       {
63         mp_size_t i;
64         mp_limb n1, n0;
65         mp_limb d;
66
67         d = den_ptr[0];
68         n1 = num_ptr[num_size - 1];
69
70         if (n1 >= d)
71           {
72             most_significant_q_limb = 1;
73             n1 -= d;
74           }
75
76         for (i = num_size - 2; i >= 0; i--)
77           {
78             n0 = num_ptr[i];
79             udiv_qrnnd (quot_ptr[i], n1, n1, n0, d);
80           }
81
82         num_ptr[0] = n1;
83       }
84       break;
85
86     case 2:
87       {
88         mp_size_t i;
89         mp_limb n1, n0, n2;
90         mp_limb d1, d0;
91
92         num_ptr += num_size - 2;
93         d1 = den_ptr[1];
94         d0 = den_ptr[0];
95         n1 = num_ptr[1];
96         n0 = num_ptr[0];
97
98         if (n1 >= d1 && (n1 > d1 || n0 >= d0))
99           {
100             most_significant_q_limb = 1;
101             sub_ddmmss (n1, n0, n1, n0, d1, d0);
102           }
103
104         for (i = num_size - den_size - 1; i >= 0; i--)
105           {
106             mp_limb q;
107             mp_limb r;
108
109             num_ptr--;
110             if (n1 == d1)
111               {
112                 /* Q should be either 111..111 or 111..110.  Need special
113                    treatment of this rare case as normal division would
114                    give overflow.  */
115                 q = ~(mp_limb) 0;
116
117                 r = n0 + d1;
118                 if (r < d1)     /* Carry in the addition? */
119                   {
120                     add_ssaaaa (n1, n0, r - d0, num_ptr[0], 0, d0);
121                     quot_ptr[i] = q;
122                     continue;
123                   }
124                 n1 = d0 - (d0 != 0);
125                 n0 = -d0;
126               }
127             else
128               {
129                 udiv_qrnnd (q, r, n1, n0, d1);
130                 umul_ppmm (n1, n0, d0, q);
131               }
132
133             n2 = num_ptr[0];
134           q_test:
135             if (n1 > r || (n1 == r && n0 > n2))
136               {
137                 /* The estimated Q was too large.  */
138                 q--;
139
140                 sub_ddmmss (n1, n0, n1, n0, 0, d0);
141                 r += d1;
142                 if (r >= d1)    /* If not carry, test Q again.  */
143                   goto q_test;
144               }
145
146             quot_ptr[i] = q;
147             sub_ddmmss (n1, n0, r, n2, n1, n0);
148           }
149         num_ptr[1] = n1;
150         num_ptr[0] = n0;
151       }
152       break;
153
154     default:
155       {
156         mp_size_t i;
157         mp_limb dX, d1, n0;
158
159         num_ptr += num_size;
160         den_ptr += den_size;
161         dX = den_ptr[-1];
162         d1 = den_ptr[-2];
163         n0 = num_ptr[-1];
164
165         if (n0 >= dX)
166           {
167             if (n0 > dX
168                 || __mpn_cmp (num_ptr - den_size, den_ptr - den_size,
169                               den_size - 1) >= 0)
170               {
171                 __mpn_sub_n (num_ptr - den_size,
172                              num_ptr - den_size, den_ptr - den_size,
173                              den_size);
174                 most_significant_q_limb = 1;
175               }
176
177             n0 = num_ptr[-1];
178           }
179
180         for (i = num_size - den_size - 1; i >= 0; i--)
181           {
182             mp_limb q;
183             mp_limb n1;
184             mp_limb cy_limb;
185
186             num_ptr--;
187             if (n0 == dX)
188               /* This might over-estimate q, but it's probably not worth
189                  the extra code here to find out.  */
190               q = ~(mp_limb) 0;
191             else
192               {
193                 mp_limb r;
194
195                 udiv_qrnnd (q, r, n0, num_ptr[-1], dX);
196                 umul_ppmm (n1, n0, d1, q);
197
198                 while (n1 > r || (n1 == r && n0 > num_ptr[-2]))
199                   {
200                     q--;
201                     r += dX;
202                     if (r < dX) /* I.e. "carry in previous addition?"  */
203                       break;
204                     n1 -= n0 < d1;
205                     n0 -= d1;
206                   }
207               }
208
209             /* Possible optimization: We already have (q * n0) and (1 * n1)
210                after the calculation of q.  Taking advantage of that, we
211                could make this loop make two iterations less.  */
212
213             cy_limb = __mpn_submul_1 (num_ptr - den_size,
214                                       den_ptr - den_size, den_size, q);
215
216             if (num_ptr[0] != cy_limb)
217               {
218                 mp_limb cy;
219                 cy = __mpn_add_n (num_ptr - den_size,
220                                   num_ptr - den_size,
221                                   den_ptr - den_size, den_size);
222                 if (cy == 0)
223                   abort ();
224                 q--;
225               }
226
227             quot_ptr[i] = q;
228             n0 = num_ptr[-1];
229           }
230       }
231     }
232
233   return most_significant_q_limb;
234 }