1 /* div.c: bcmath library file. */
2 /*
3 Copyright (C) 1991, 1992, 1993, 1994, 1997 Free Software Foundation, Inc.
4 Copyright (C) 2000 Philip A. Nelson
5
6 This library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2 of the License, or (at your option) any later version.
10
11 This library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details. (LICENSE)
15
16 You should have received a copy of the GNU Lesser General Public
17 License along with this library; if not, write to:
18
19 The Free Software Foundation, Inc.
20 59 Temple Place, Suite 330
21 Boston, MA 02111-1307 USA.
22
23 You may contact the author by:
24 e-mail: philnelson@acm.org
25 us-mail: Philip A. Nelson
26 Computer Science Department, 9062
27 Western Washington University
28 Bellingham, WA 98226-9062
29
30 *************************************************************************/
31
32 #include <config.h>
33 #include <stdio.h>
34 #include <stdlib.h>
35 #include <ctype.h>
36 #include <stdarg.h>
37 #include "bcmath.h"
38 #include "private.h"
39
40
41 /* Some utility routines for the divide: First a one digit multiply.
42 NUM (with SIZE digits) is multiplied by DIGIT and the result is
43 placed into RESULT. It is written so that NUM and RESULT can be
44 the same pointers. */
45
46 static void
_one_mult(num,size,digit,result)47 _one_mult (num, size, digit, result)
48 unsigned char *num;
49 int size, digit;
50 unsigned char *result;
51 {
52 int carry, value;
53 unsigned char *nptr, *rptr;
54
55 if (digit == 0)
56 memset (result, 0, size);
57 else
58 {
59 if (digit == 1)
60 memcpy (result, num, size);
61 else
62 {
63 /* Initialize */
64 nptr = (unsigned char *) (num+size-1);
65 rptr = (unsigned char *) (result+size-1);
66 carry = 0;
67
68 while (size-- > 0)
69 {
70 value = *nptr-- * digit + carry;
71 *rptr-- = value % BASE;
72 carry = value / BASE;
73 }
74
75 if (carry != 0) *rptr = carry;
76 }
77 }
78 }
79
80
81 /* The full division routine. This computes N1 / N2. It returns
82 0 if the division is ok and the result is in QUOT. The number of
83 digits after the decimal point is SCALE. It returns -1 if division
84 by zero is tried. The algorithm is found in Knuth Vol 2. p237. */
85
86 int
bc_divide(bc_num n1,bc_num n2,bc_num * quot,int scale)87 bc_divide (bc_num n1, bc_num n2, bc_num *quot, int scale)
88 {
89 bc_num qval;
90 unsigned char *num1, *num2;
91 unsigned char *ptr1, *ptr2, *n2ptr, *qptr;
92 int scale1, val;
93 unsigned int len1, len2, scale2, qdigits, extra, count;
94 unsigned int qdig, qguess, borrow, carry;
95 unsigned char *mval;
96 char zero;
97 unsigned int norm;
98
99 /* Test for divide by zero. */
100 if (bc_is_zero (n2)) return -1;
101
102 /* Test for divide by 1. If it is we must truncate. */
103 if (n2->n_scale == 0)
104 {
105 if (n2->n_len == 1 && *n2->n_value == 1)
106 {
107 qval = bc_new_num (n1->n_len, scale);
108 qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS);
109 memset (&qval->n_value[n1->n_len],0,scale);
110 memcpy (qval->n_value, n1->n_value,
111 n1->n_len + MIN(n1->n_scale,scale));
112 bc_free_num (quot);
113 *quot = qval;
114 }
115 }
116
117 /* Set up the divide. Move the decimal point on n1 by n2's scale.
118 Remember, zeros on the end of num2 are wasted effort for dividing. */
119 scale2 = n2->n_scale;
120 n2ptr = (unsigned char *) n2->n_value+n2->n_len+scale2-1;
121 while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--;
122
123 len1 = n1->n_len + scale2;
124 scale1 = n1->n_scale - scale2;
125 if (scale1 < scale)
126 extra = scale - scale1;
127 else
128 extra = 0;
129 num1 = (unsigned char *) safe_emalloc (1, n1->n_len+n1->n_scale, extra+2);
130 memset (num1, 0, n1->n_len+n1->n_scale+extra+2);
131 memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale);
132
133 len2 = n2->n_len + scale2;
134 num2 = (unsigned char *) safe_emalloc (1, len2, 1);
135 memcpy (num2, n2->n_value, len2);
136 *(num2+len2) = 0;
137 n2ptr = num2;
138 while (*n2ptr == 0)
139 {
140 n2ptr++;
141 len2--;
142 }
143
144 /* Calculate the number of quotient digits. */
145 if (len2 > len1+scale)
146 {
147 qdigits = scale+1;
148 zero = TRUE;
149 }
150 else
151 {
152 zero = FALSE;
153 if (len2>len1)
154 qdigits = scale+1; /* One for the zero integer part. */
155 else
156 qdigits = len1-len2+scale+1;
157 }
158
159 /* Allocate and zero the storage for the quotient. */
160 qval = bc_new_num (qdigits-scale,scale);
161 memset (qval->n_value, 0, qdigits);
162
163 /* Allocate storage for the temporary storage mval. */
164 mval = (unsigned char *) safe_emalloc (1, len2, 1);
165
166 /* Now for the full divide algorithm. */
167 if (!zero)
168 {
169 /* Normalize */
170 norm = 10 / ((int)*n2ptr + 1);
171 if (norm != 1)
172 {
173 _one_mult (num1, len1+scale1+extra+1, norm, num1);
174 _one_mult (n2ptr, len2, norm, n2ptr);
175 }
176
177 /* Initialize divide loop. */
178 qdig = 0;
179 if (len2 > len1)
180 qptr = (unsigned char *) qval->n_value+len2-len1;
181 else
182 qptr = (unsigned char *) qval->n_value;
183
184 /* Loop */
185 while (qdig <= len1+scale-len2)
186 {
187 /* Calculate the quotient digit guess. */
188 if (*n2ptr == num1[qdig])
189 qguess = 9;
190 else
191 qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr;
192
193 /* Test qguess. */
194 if (n2ptr[1]*qguess >
195 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
196 + num1[qdig+2])
197 {
198 qguess--;
199 /* And again. */
200 if (n2ptr[1]*qguess >
201 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
202 + num1[qdig+2])
203 qguess--;
204 }
205
206 /* Multiply and subtract. */
207 borrow = 0;
208 if (qguess != 0)
209 {
210 *mval = 0;
211 _one_mult (n2ptr, len2, qguess, mval+1);
212 ptr1 = (unsigned char *) num1+qdig+len2;
213 ptr2 = (unsigned char *) mval+len2;
214 for (count = 0; count < len2+1; count++)
215 {
216 val = (int) *ptr1 - (int) *ptr2-- - borrow;
217 if (val < 0)
218 {
219 val += 10;
220 borrow = 1;
221 }
222 else
223 borrow = 0;
224 *ptr1-- = val;
225 }
226 }
227
228 /* Test for negative result. */
229 if (borrow == 1)
230 {
231 qguess--;
232 ptr1 = (unsigned char *) num1+qdig+len2;
233 ptr2 = (unsigned char *) n2ptr+len2-1;
234 carry = 0;
235 for (count = 0; count < len2; count++)
236 {
237 val = (int) *ptr1 + (int) *ptr2-- + carry;
238 if (val > 9)
239 {
240 val -= 10;
241 carry = 1;
242 }
243 else
244 carry = 0;
245 *ptr1-- = val;
246 }
247 if (carry == 1) *ptr1 = (*ptr1 + 1) % 10;
248 }
249
250 /* We now know the quotient digit. */
251 *qptr++ = qguess;
252 qdig++;
253 }
254 }
255
256 /* Clean up and return the number. */
257 qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
258 if (bc_is_zero (qval)) qval->n_sign = PLUS;
259 _bc_rm_leading_zeros (qval);
260 bc_free_num (quot);
261 *quot = qval;
262
263 /* Clean up temporary storage. */
264 efree (mval);
265 efree (num1);
266 efree (num2);
267
268 return 0; /* Everything is OK. */
269 }
270
