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 TSRMLS_DC)88 bc_divide (bc_num n1, bc_num n2, bc_num *quot, int scale TSRMLS_DC)
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 TSRMLS_CC)) 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 if (num1 == NULL) bc_out_of_memory();
132 memset (num1, 0, n1->n_len+n1->n_scale+extra+2);
133 memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale);
134
135 len2 = n2->n_len + scale2;
136 num2 = (unsigned char *) safe_emalloc (1, len2, 1);
137 if (num2 == NULL) bc_out_of_memory();
138 memcpy (num2, n2->n_value, len2);
139 *(num2+len2) = 0;
140 n2ptr = num2;
141 while (*n2ptr == 0)
142 {
143 n2ptr++;
144 len2--;
145 }
146
147 /* Calculate the number of quotient digits. */
148 if (len2 > len1+scale)
149 {
150 qdigits = scale+1;
151 zero = TRUE;
152 }
153 else
154 {
155 zero = FALSE;
156 if (len2>len1)
157 qdigits = scale+1; /* One for the zero integer part. */
158 else
159 qdigits = len1-len2+scale+1;
160 }
161
162 /* Allocate and zero the storage for the quotient. */
163 qval = bc_new_num (qdigits-scale,scale);
164 memset (qval->n_value, 0, qdigits);
165
166 /* Allocate storage for the temporary storage mval. */
167 mval = (unsigned char *) safe_emalloc (1, len2, 1);
168 if (mval == NULL) bc_out_of_memory ();
169
170 /* Now for the full divide algorithm. */
171 if (!zero)
172 {
173 /* Normalize */
174 norm = 10 / ((int)*n2ptr + 1);
175 if (norm != 1)
176 {
177 _one_mult (num1, len1+scale1+extra+1, norm, num1);
178 _one_mult (n2ptr, len2, norm, n2ptr);
179 }
180
181 /* Initialize divide loop. */
182 qdig = 0;
183 if (len2 > len1)
184 qptr = (unsigned char *) qval->n_value+len2-len1;
185 else
186 qptr = (unsigned char *) qval->n_value;
187
188 /* Loop */
189 while (qdig <= len1+scale-len2)
190 {
191 /* Calculate the quotient digit guess. */
192 if (*n2ptr == num1[qdig])
193 qguess = 9;
194 else
195 qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr;
196
197 /* Test qguess. */
198 if (n2ptr[1]*qguess >
199 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
200 + num1[qdig+2])
201 {
202 qguess--;
203 /* And again. */
204 if (n2ptr[1]*qguess >
205 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
206 + num1[qdig+2])
207 qguess--;
208 }
209
210 /* Multiply and subtract. */
211 borrow = 0;
212 if (qguess != 0)
213 {
214 *mval = 0;
215 _one_mult (n2ptr, len2, qguess, mval+1);
216 ptr1 = (unsigned char *) num1+qdig+len2;
217 ptr2 = (unsigned char *) mval+len2;
218 for (count = 0; count < len2+1; count++)
219 {
220 val = (int) *ptr1 - (int) *ptr2-- - borrow;
221 if (val < 0)
222 {
223 val += 10;
224 borrow = 1;
225 }
226 else
227 borrow = 0;
228 *ptr1-- = val;
229 }
230 }
231
232 /* Test for negative result. */
233 if (borrow == 1)
234 {
235 qguess--;
236 ptr1 = (unsigned char *) num1+qdig+len2;
237 ptr2 = (unsigned char *) n2ptr+len2-1;
238 carry = 0;
239 for (count = 0; count < len2; count++)
240 {
241 val = (int) *ptr1 + (int) *ptr2-- + carry;
242 if (val > 9)
243 {
244 val -= 10;
245 carry = 1;
246 }
247 else
248 carry = 0;
249 *ptr1-- = val;
250 }
251 if (carry == 1) *ptr1 = (*ptr1 + 1) % 10;
252 }
253
254 /* We now know the quotient digit. */
255 *qptr++ = qguess;
256 qdig++;
257 }
258 }
259
260 /* Clean up and return the number. */
261 qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
262 if (bc_is_zero (qval TSRMLS_CC)) qval->n_sign = PLUS;
263 _bc_rm_leading_zeros (qval);
264 bc_free_num (quot);
265 *quot = qval;
266
267 /* Clean up temporary storage. */
268 efree (mval);
269 efree (num1);
270 efree (num2);
271
272 return 0; /* Everything is OK. */
273 }
274
275
