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. (COPYING.LIB)
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 <assert.h>
35 #include <stdlib.h>
36 #include <ctype.h>
37 #include <stdarg.h>
38 #include "bcmath.h"
39 #include "private.h"
40
41
42 /* Some utility routines for the divide: First a one digit multiply.
43 NUM (with SIZE digits) is multiplied by DIGIT and the result is
44 placed into RESULT. It is written so that NUM and RESULT can be
45 the same pointers. */
46
47 static void
_one_mult(num,size,digit,result)48 _one_mult (num, size, digit, result)
49 unsigned char *num;
50 int size, digit;
51 unsigned char *result;
52 {
53 int carry, value;
54 unsigned char *nptr, *rptr;
55
56 if (digit == 0)
57 memset (result, 0, size);
58 else
59 {
60 if (digit == 1)
61 memcpy (result, num, size);
62 else
63 {
64 /* Initialize */
65 nptr = (unsigned char *) (num+size-1);
66 rptr = (unsigned char *) (result+size-1);
67 carry = 0;
68
69 while (size-- > 0)
70 {
71 value = *nptr-- * digit + carry;
72 *rptr-- = value % BASE;
73 carry = value / BASE;
74 }
75
76 if (carry != 0) *rptr = carry;
77 }
78 }
79 }
80
81
82 /* The full division routine. This computes N1 / N2. It returns
83 0 if the division is ok and the result is in QUOT. The number of
84 digits after the decimal point is SCALE. It returns -1 if division
85 by zero is tried. The algorithm is found in Knuth Vol 2. p237. */
86
87 int
bc_divide(bc_num n1,bc_num n2,bc_num * quot,int scale)88 bc_divide (bc_num n1, bc_num n2, bc_num *quot, int scale)
89 {
90 bc_num qval;
91 unsigned char *num1, *num2;
92 unsigned char *ptr1, *ptr2, *n2ptr, *qptr;
93 int scale1, val;
94 unsigned int len1, len2, scale2, qdigits, extra, count;
95 unsigned int qdig, qguess, borrow, carry;
96 unsigned char *mval;
97 char zero;
98 unsigned int norm;
99
100 /* Test for divide by zero. */
101 if (bc_is_zero (n2)) return -1;
102
103 /* Test for divide by 1. If it is we must truncate. */
104 if (n2->n_scale == 0)
105 {
106 if (n2->n_len == 1 && *n2->n_value == 1)
107 {
108 qval = bc_new_num (n1->n_len, scale);
109 qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS);
110 memset (&qval->n_value[n1->n_len],0,scale);
111 memcpy (qval->n_value, n1->n_value,
112 n1->n_len + MIN(n1->n_scale,scale));
113 bc_free_num (quot);
114 *quot = qval;
115 }
116 }
117
118 /* Set up the divide. Move the decimal point on n1 by n2's scale.
119 Remember, zeros on the end of num2 are wasted effort for dividing. */
120 scale2 = n2->n_scale;
121 n2ptr = (unsigned char *) n2->n_value+n2->n_len+scale2-1;
122 while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--;
123
124 len1 = n1->n_len + scale2;
125 scale1 = n1->n_scale - scale2;
126 if (scale1 < scale)
127 extra = scale - scale1;
128 else
129 extra = 0;
130 num1 = (unsigned char *) safe_emalloc (1, n1->n_len+n1->n_scale, extra+2);
131 memset (num1, 0, n1->n_len+n1->n_scale+extra+2);
132 memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale);
133
134 len2 = n2->n_len + scale2;
135 num2 = (unsigned char *) safe_emalloc (1, len2, 1);
136 memcpy (num2, n2->n_value, len2);
137 *(num2+len2) = 0;
138 n2ptr = num2;
139 while (*n2ptr == 0)
140 {
141 n2ptr++;
142 len2--;
143 }
144
145 /* Calculate the number of quotient digits. */
146 if (len2 > len1+scale)
147 {
148 qdigits = scale+1;
149 zero = TRUE;
150 }
151 else
152 {
153 zero = FALSE;
154 if (len2>len1)
155 qdigits = scale+1; /* One for the zero integer part. */
156 else
157 qdigits = len1-len2+scale+1;
158 }
159
160 /* Allocate and zero the storage for the quotient. */
161 qval = bc_new_num (qdigits-scale,scale);
162 memset (qval->n_value, 0, qdigits);
163
164 /* Allocate storage for the temporary storage mval. */
165 mval = (unsigned char *) safe_emalloc (1, len2, 1);
166
167 /* Now for the full divide algorithm. */
168 if (!zero)
169 {
170 /* Normalize */
171 norm = 10 / ((int)*n2ptr + 1);
172 if (norm != 1)
173 {
174 _one_mult (num1, len1+scale1+extra+1, norm, num1);
175 _one_mult (n2ptr, len2, norm, n2ptr);
176 }
177
178 /* Initialize divide loop. */
179 qdig = 0;
180 if (len2 > len1)
181 qptr = (unsigned char *) qval->n_value+len2-len1;
182 else
183 qptr = (unsigned char *) qval->n_value;
184
185 /* Loop */
186 while (qdig <= len1+scale-len2)
187 {
188 /* Calculate the quotient digit guess. */
189 if (*n2ptr == num1[qdig])
190 qguess = 9;
191 else
192 qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr;
193
194 /* Test qguess. */
195 if (n2ptr[1]*qguess >
196 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
197 + num1[qdig+2])
198 {
199 qguess--;
200 /* And again. */
201 if (n2ptr[1]*qguess >
202 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
203 + num1[qdig+2])
204 qguess--;
205 }
206
207 /* Multiply and subtract. */
208 borrow = 0;
209 if (qguess != 0)
210 {
211 *mval = 0;
212 _one_mult (n2ptr, len2, qguess, mval+1);
213 ptr1 = (unsigned char *) num1+qdig+len2;
214 ptr2 = (unsigned char *) mval+len2;
215 for (count = 0; count < len2+1; count++)
216 {
217 val = (int) *ptr1 - (int) *ptr2-- - borrow;
218 if (val < 0)
219 {
220 val += 10;
221 borrow = 1;
222 }
223 else
224 borrow = 0;
225 *ptr1-- = val;
226 }
227 }
228
229 /* Test for negative result. */
230 if (borrow == 1)
231 {
232 qguess--;
233 ptr1 = (unsigned char *) num1+qdig+len2;
234 ptr2 = (unsigned char *) n2ptr+len2-1;
235 carry = 0;
236 for (count = 0; count < len2; count++)
237 {
238 val = (int) *ptr1 + (int) *ptr2-- + carry;
239 if (val > 9)
240 {
241 val -= 10;
242 carry = 1;
243 }
244 else
245 carry = 0;
246 *ptr1-- = val;
247 }
248 if (carry == 1) *ptr1 = (*ptr1 + 1) % 10;
249 }
250
251 /* We now know the quotient digit. */
252 *qptr++ = qguess;
253 qdig++;
254 }
255 }
256
257 /* Clean up and return the number. */
258 qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
259 if (bc_is_zero (qval)) qval->n_sign = PLUS;
260 _bc_rm_leading_zeros (qval);
261 bc_free_num (quot);
262 *quot = qval;
263
264 /* Clean up temporary storage. */
265 efree (mval);
266 efree (num1);
267 efree (num2);
268
269 return 0; /* Everything is OK. */
270 }
271
