xref: /PHP-8.3/ext/bcmath/libbcmath/src/recmul.c (revision e56ed6e1)
1 /* recmul.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 <stddef.h>
34 #include <assert.h>
35 #include <stdbool.h>
36 #include "private.h" /* For _bc_rm_leading_zeros() */
37 #include "zend_alloc.h"
38 
39 /* Recursive vs non-recursive multiply crossover ranges. */
40 #if defined(MULDIGITS)
41 #include "muldigits.h"
42 #else
43 #define MUL_BASE_DIGITS 80
44 #endif
45 
46 int mul_base_digits = MUL_BASE_DIGITS;
47 #define MUL_SMALL_DIGITS mul_base_digits/4
48 
49 /* Multiply utility routines */
50 
new_sub_num(size_t length,size_t scale,char * value)51 static bc_num new_sub_num(size_t length, size_t scale, char *value)
52 {
53 	bc_num temp = (bc_num) emalloc(sizeof(bc_struct));
54 
55 	temp->n_sign = PLUS;
56 	temp->n_len = length;
57 	temp->n_scale = scale;
58 	temp->n_refs = 1;
59 	temp->n_ptr = NULL;
60 	temp->n_value = value;
61 	return temp;
62 }
63 
_bc_simp_mul(bc_num n1,size_t n1len,bc_num n2,int n2len,bc_num * prod)64 static void _bc_simp_mul(bc_num n1, size_t n1len, bc_num n2, int n2len, bc_num *prod)
65 {
66 	char *n1ptr, *n2ptr, *pvptr;
67 	char *n1end, *n2end;        /* To the end of n1 and n2. */
68 	int sum = 0;
69 
70 	int prodlen = n1len + n2len + 1;
71 
72 	*prod = bc_new_num (prodlen, 0);
73 
74 	n1end = (char *) (n1->n_value + n1len - 1);
75 	n2end = (char *) (n2->n_value + n2len - 1);
76 	pvptr = (char *) ((*prod)->n_value + prodlen - 1);
77 
78 	/* Here is the loop... */
79 	for (int index = 0; index < prodlen - 1; index++) {
80 		n1ptr = (char *) (n1end - MAX(0, index - n2len + 1));
81 		n2ptr = (char *) (n2end - MIN(index, n2len - 1));
82 		while ((n1ptr >= n1->n_value) && (n2ptr <= n2end)) {
83 			sum += *n1ptr * *n2ptr;
84 			n1ptr--;
85 			n2ptr++;
86 		}
87 		*pvptr-- = sum % BASE;
88 		sum = sum / BASE;
89 	}
90 	*pvptr = sum;
91 }
92 
93 
94 /* A special adder/subtractor for the recursive divide and conquer
95    multiply algorithm.  Note: if sub is called, accum must
96    be larger that what is being subtracted.  Also, accum and val
97    must have n_scale = 0.  (e.g. they must look like integers. *) */
_bc_shift_addsub(bc_num accum,bc_num val,int shift,bool sub)98 static void _bc_shift_addsub(bc_num accum, bc_num val, int shift, bool sub)
99 {
100 	signed char *accp, *valp;
101 	unsigned int carry = 0;
102 	size_t count = val->n_len;
103 
104 	if (val->n_value[0] == 0) {
105 		count--;
106 	}
107 	assert(accum->n_len + accum->n_scale >= shift + count);
108 
109 	/* Set up pointers and others */
110 	accp = (signed char *) (accum->n_value + accum->n_len + accum->n_scale - shift - 1);
111 	valp = (signed char *) (val->n_value + val->n_len - 1);
112 
113 	if (sub) {
114 		/* Subtraction, carry is really borrow. */
115 		while (count--) {
116 			*accp -= *valp-- + carry;
117 			if (*accp < 0) {
118 				carry = 1;
119 				*accp-- += BASE;
120 			} else {
121 				carry = 0;
122 				accp--;
123 			}
124 		}
125 		while (carry) {
126 			*accp -= carry;
127 			if (*accp < 0) {
128 				*accp-- += BASE;
129 			} else {
130 				carry = 0;
131 			}
132 		}
133 	} else {
134 		/* Addition */
135 		while (count--) {
136 			*accp += *valp-- + carry;
137 			if (*accp > (BASE - 1)) {
138 				carry = 1;
139 				*accp-- -= BASE;
140 			} else {
141 				carry = 0;
142 				accp--;
143 			}
144 		}
145 		while (carry) {
146 			*accp += carry;
147 			if (*accp > (BASE - 1)) {
148 				*accp-- -= BASE;
149 			} else {
150 				carry = 0;
151 			}
152 		}
153 	}
154 }
155 
156 /* Recursive divide and conquer multiply algorithm.
157    Based on
158    Let u = u0 + u1*(b^n)
159    Let v = v0 + v1*(b^n)
160    Then uv = (B^2n+B^n)*u1*v1 + B^n*(u1-u0)*(v0-v1) + (B^n+1)*u0*v0
161 
162    B is the base of storage, number of digits in u1,u0 close to equal.
163 */
_bc_rec_mul(bc_num u,size_t ulen,bc_num v,size_t vlen,bc_num * prod)164 static void _bc_rec_mul(bc_num u, size_t ulen, bc_num v, size_t vlen, bc_num *prod)
165 {
166 	bc_num u0, u1, v0, v1;
167 	bc_num m1, m2, m3, d1, d2;
168 	size_t n;
169 	bool m1zero;
170 
171 	/* Base case? */
172 	if ((ulen + vlen) < mul_base_digits
173 		|| ulen < MUL_SMALL_DIGITS
174 		|| vlen < MUL_SMALL_DIGITS
175 	) {
176 		_bc_simp_mul(u, ulen, v, vlen, prod);
177 		return;
178 	}
179 
180 	/* Calculate n -- the u and v split point in digits. */
181 	n = (MAX(ulen, vlen) + 1) / 2;
182 
183 	/* Split u and v. */
184 	if (ulen < n) {
185 		u1 = bc_copy_num(BCG(_zero_));
186 		u0 = new_sub_num(ulen, 0, u->n_value);
187 	} else {
188 		u1 = new_sub_num(ulen - n, 0, u->n_value);
189 		u0 = new_sub_num(n, 0, u->n_value + ulen - n);
190 	}
191 	if (vlen < n) {
192 		v1 = bc_copy_num(BCG(_zero_));
193 		v0 = new_sub_num(vlen, 0, v->n_value);
194 	} else {
195 		v1 = new_sub_num(vlen - n, 0, v->n_value);
196 		v0 = new_sub_num(n, 0, v->n_value + vlen - n);
197 	}
198 	_bc_rm_leading_zeros(u1);
199 	_bc_rm_leading_zeros(u0);
200 	_bc_rm_leading_zeros(v1);
201 	_bc_rm_leading_zeros(v0);
202 
203 	m1zero = bc_is_zero(u1) || bc_is_zero(v1);
204 
205 	/* Calculate sub results ... */
206 
207 	bc_init_num(&d1);
208 	bc_init_num(&d2);
209 	bc_sub(u1, u0, &d1, 0);
210 	bc_sub(v0, v1, &d2, 0);
211 
212 
213 	/* Do recursive multiplies and shifted adds. */
214 	if (m1zero) {
215 		m1 = bc_copy_num(BCG(_zero_));
216 	} else {
217 		_bc_rec_mul(u1, u1->n_len, v1, v1->n_len, &m1);
218 	}
219 
220 	if (bc_is_zero(d1) || bc_is_zero(d2)) {
221 		m2 = bc_copy_num(BCG(_zero_));
222 	} else {
223 		_bc_rec_mul(d1, d1->n_len, d2, d2->n_len, &m2);
224 	}
225 
226 	if (bc_is_zero(u0) || bc_is_zero(v0)) {
227 		m3 = bc_copy_num(BCG(_zero_));
228 	} else {
229 		_bc_rec_mul(u0, u0->n_len, v0, v0->n_len, &m3);
230 	}
231 
232 	/* Initialize product */
233 	*prod = bc_new_num(ulen + vlen + 1, 0);
234 
235 	if (!m1zero) {
236 		_bc_shift_addsub(*prod, m1, 2 * n, false);
237 		_bc_shift_addsub(*prod, m1, n, false);
238 	}
239 	_bc_shift_addsub(*prod, m3, n, false);
240 	_bc_shift_addsub(*prod, m3, 0, false);
241 	_bc_shift_addsub(*prod, m2, n, d1->n_sign != d2->n_sign);
242 
243 	/* Now clean up! */
244 	bc_free_num (&u1);
245 	bc_free_num (&u0);
246 	bc_free_num (&v1);
247 	bc_free_num (&m1);
248 	bc_free_num (&v0);
249 	bc_free_num (&m2);
250 	bc_free_num (&m3);
251 	bc_free_num (&d1);
252 	bc_free_num (&d2);
253 }
254 
255 /* The multiply routine.  N2 times N1 is put int PROD with the scale of
256    the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
257    */
258 
bc_multiply(bc_num n1,bc_num n2,bc_num * prod,size_t scale)259 void bc_multiply(bc_num n1, bc_num n2, bc_num *prod, size_t scale)
260 {
261 	bc_num pval;
262 	size_t len1, len2;
263 	size_t full_scale, prod_scale;
264 
265 	/* Initialize things. */
266 	len1 = n1->n_len + n1->n_scale;
267 	len2 = n2->n_len + n2->n_scale;
268 	full_scale = n1->n_scale + n2->n_scale;
269 	prod_scale = MIN(full_scale, MAX(scale, MAX(n1->n_scale, n2->n_scale)));
270 
271 	/* Do the multiply */
272 	_bc_rec_mul(n1, len1, n2, len2, &pval);
273 
274 	/* Assign to prod and clean up the number. */
275 	pval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS);
276 	pval->n_value = pval->n_ptr;
277 	pval->n_len = len2 + len1 + 1 - full_scale;
278 	pval->n_scale = prod_scale;
279 	_bc_rm_leading_zeros(pval);
280 	if (bc_is_zero(pval)) {
281 		pval->n_sign = PLUS;
282 	}
283 	bc_free_num(prod);
284 	*prod = pval;
285 }
286