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