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
48 : _one_mult (num, size, digit, result)
49 : unsigned char *num;
50 : int size, digit;
51 : unsigned char *result;
52 217 : {
53 : int carry, value;
54 : unsigned char *nptr, *rptr;
55 :
56 217 : if (digit == 0)
57 0 : memset (result, 0, size);
58 : else
59 : {
60 217 : if (digit == 1)
61 21 : memcpy (result, num, size);
62 : else
63 : {
64 : /* Initialize */
65 196 : nptr = (unsigned char *) (num+size-1);
66 196 : rptr = (unsigned char *) (result+size-1);
67 196 : carry = 0;
68 :
69 2926 : while (size-- > 0)
70 : {
71 2534 : value = *nptr-- * digit + carry;
72 2534 : *rptr-- = value % BASE;
73 2534 : carry = value / BASE;
74 : }
75 :
76 196 : if (carry != 0) *rptr = carry;
77 : }
78 : }
79 217 : }
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
88 : bc_divide (bc_num n1, bc_num n2, bc_num *quot, int scale TSRMLS_DC)
89 24 : {
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 24 : if (bc_is_zero (n2 TSRMLS_CC)) return -1;
102 :
103 : /* Test for divide by 1. If it is we must truncate. */
104 23 : if (n2->n_scale == 0)
105 : {
106 11 : if (n2->n_len == 1 && *n2->n_value == 1)
107 : {
108 3 : qval = bc_new_num (n1->n_len, scale);
109 3 : qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS);
110 3 : memset (&qval->n_value[n1->n_len],0,scale);
111 3 : memcpy (qval->n_value, n1->n_value,
112 : n1->n_len + MIN(n1->n_scale,scale));
113 3 : bc_free_num (quot);
114 3 : *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 23 : scale2 = n2->n_scale;
121 23 : n2ptr = (unsigned char *) n2->n_value+n2->n_len+scale2-1;
122 23 : while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--;
123 :
124 23 : len1 = n1->n_len + scale2;
125 23 : scale1 = n1->n_scale - scale2;
126 23 : if (scale1 < scale)
127 17 : extra = scale - scale1;
128 : else
129 6 : extra = 0;
130 23 : num1 = (unsigned char *) safe_emalloc (1, n1->n_len+n1->n_scale, extra+2);
131 23 : if (num1 == NULL) bc_out_of_memory();
132 23 : memset (num1, 0, n1->n_len+n1->n_scale+extra+2);
133 23 : memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale);
134 :
135 23 : len2 = n2->n_len + scale2;
136 23 : num2 = (unsigned char *) safe_emalloc (1, len2, 1);
137 23 : if (num2 == NULL) bc_out_of_memory();
138 23 : memcpy (num2, n2->n_value, len2);
139 23 : *(num2+len2) = 0;
140 23 : n2ptr = num2;
141 46 : while (*n2ptr == 0)
142 : {
143 0 : n2ptr++;
144 0 : len2--;
145 : }
146 :
147 : /* Calculate the number of quotient digits. */
148 23 : if (len2 > len1+scale)
149 : {
150 0 : qdigits = scale+1;
151 0 : zero = TRUE;
152 : }
153 : else
154 : {
155 23 : zero = FALSE;
156 23 : if (len2>len1)
157 0 : qdigits = scale+1; /* One for the zero integer part. */
158 : else
159 23 : qdigits = len1-len2+scale+1;
160 : }
161 :
162 : /* Allocate and zero the storage for the quotient. */
163 23 : qval = bc_new_num (qdigits-scale,scale);
164 23 : memset (qval->n_value, 0, qdigits);
165 :
166 : /* Allocate storage for the temporary storage mval. */
167 23 : mval = (unsigned char *) safe_emalloc (1, len2, 1);
168 23 : if (mval == NULL) bc_out_of_memory ();
169 :
170 : /* Now for the full divide algorithm. */
171 23 : if (!zero)
172 : {
173 : /* Normalize */
174 23 : norm = 10 / ((int)*n2ptr + 1);
175 23 : if (norm != 1)
176 : {
177 19 : _one_mult (num1, len1+scale1+extra+1, norm, num1);
178 19 : _one_mult (n2ptr, len2, norm, n2ptr);
179 : }
180 :
181 : /* Initialize divide loop. */
182 23 : qdig = 0;
183 23 : if (len2 > len1)
184 0 : qptr = (unsigned char *) qval->n_value+len2-len1;
185 : else
186 23 : qptr = (unsigned char *) qval->n_value;
187 :
188 : /* Loop */
189 257 : while (qdig <= len1+scale-len2)
190 : {
191 : /* Calculate the quotient digit guess. */
192 211 : if (*n2ptr == num1[qdig])
193 10 : qguess = 9;
194 : else
195 201 : qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr;
196 :
197 : /* Test qguess. */
198 211 : if (n2ptr[1]*qguess >
199 : (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
200 : + num1[qdig+2])
201 : {
202 46 : qguess--;
203 : /* And again. */
204 46 : if (n2ptr[1]*qguess >
205 : (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
206 : + num1[qdig+2])
207 0 : qguess--;
208 : }
209 :
210 : /* Multiply and subtract. */
211 211 : borrow = 0;
212 211 : if (qguess != 0)
213 : {
214 179 : *mval = 0;
215 179 : _one_mult (n2ptr, len2, qguess, mval+1);
216 179 : ptr1 = (unsigned char *) num1+qdig+len2;
217 179 : ptr2 = (unsigned char *) mval+len2;
218 2634 : for (count = 0; count < len2+1; count++)
219 : {
220 2455 : val = (int) *ptr1 - (int) *ptr2-- - borrow;
221 2455 : if (val < 0)
222 : {
223 1010 : val += 10;
224 1010 : borrow = 1;
225 : }
226 : else
227 1445 : borrow = 0;
228 2455 : *ptr1-- = val;
229 : }
230 : }
231 :
232 : /* Test for negative result. */
233 211 : if (borrow == 1)
234 : {
235 2 : qguess--;
236 2 : ptr1 = (unsigned char *) num1+qdig+len2;
237 2 : ptr2 = (unsigned char *) n2ptr+len2-1;
238 2 : carry = 0;
239 20 : for (count = 0; count < len2; count++)
240 : {
241 18 : val = (int) *ptr1 + (int) *ptr2-- + carry;
242 18 : if (val > 9)
243 : {
244 9 : val -= 10;
245 9 : carry = 1;
246 : }
247 : else
248 9 : carry = 0;
249 18 : *ptr1-- = val;
250 : }
251 2 : if (carry == 1) *ptr1 = (*ptr1 + 1) % 10;
252 : }
253 :
254 : /* We now know the quotient digit. */
255 211 : *qptr++ = qguess;
256 211 : qdig++;
257 : }
258 : }
259 :
260 : /* Clean up and return the number. */
261 23 : qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
262 23 : if (bc_is_zero (qval TSRMLS_CC)) qval->n_sign = PLUS;
263 23 : _bc_rm_leading_zeros (qval);
264 23 : bc_free_num (quot);
265 23 : *quot = qval;
266 :
267 : /* Clean up temporary storage. */
268 23 : efree (mval);
269 23 : efree (num1);
270 23 : efree (num2);
271 :
272 23 : return 0; /* Everything is OK. */
273 : }
274 :
|
