xref: /PHP-8.3/ext/bcmath/libbcmath/src/div.c (revision e56ed6e1)
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