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