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"
37 #include "convert.h"
38 #include "zend_alloc.h"
39
40
41 /* Multiply utility routines */
42
bc_mul_carry_calc(BC_VECTOR * prod_vector,size_t prod_arr_size)43 static inline void bc_mul_carry_calc(BC_VECTOR *prod_vector, size_t prod_arr_size)
44 {
45 for (size_t i = 0; i < prod_arr_size - 1; i++) {
46 prod_vector[i + 1] += prod_vector[i] / BC_VECTOR_BOUNDARY_NUM;
47 prod_vector[i] %= BC_VECTOR_BOUNDARY_NUM;
48 }
49 }
50
51 /*
52 * If the n_values of n1 and n2 are both 4 (32-bit) or 8 (64-bit) digits or less,
53 * the calculation will be performed at high speed without using an array.
54 */
bc_fast_mul(bc_num n1,size_t n1len,bc_num n2,size_t n2len,bc_num * prod)55 static inline void bc_fast_mul(bc_num n1, size_t n1len, bc_num n2, size_t n2len, bc_num *prod)
56 {
57 const char *n1end = n1->n_value + n1len - 1;
58 const char *n2end = n2->n_value + n2len - 1;
59
60 BC_VECTOR n1_vector = bc_partial_convert_to_vector(n1end, n1len);
61 BC_VECTOR n2_vector = bc_partial_convert_to_vector(n2end, n2len);
62 BC_VECTOR prod_vector = n1_vector * n2_vector;
63
64 size_t prodlen = n1len + n2len;
65 *prod = bc_new_num_nonzeroed(prodlen, 0);
66 char *pptr = (*prod)->n_value;
67 char *pend = pptr + prodlen - 1;
68
69 while (pend >= pptr) {
70 *pend-- = prod_vector % BASE;
71 prod_vector /= BASE;
72 }
73 }
74
75 /*
76 * Equivalent of bc_fast_mul for small numbers to perform computations
77 * without using array.
78 */
bc_fast_square(bc_num n1,size_t n1len,bc_num * prod)79 static inline void bc_fast_square(bc_num n1, size_t n1len, bc_num *prod)
80 {
81 const char *n1end = n1->n_value + n1len - 1;
82
83 BC_VECTOR n1_vector = bc_partial_convert_to_vector(n1end, n1len);
84 BC_VECTOR prod_vector = n1_vector * n1_vector;
85
86 size_t prodlen = n1len + n1len;
87 *prod = bc_new_num_nonzeroed(prodlen, 0);
88 char *pptr = (*prod)->n_value;
89 char *pend = pptr + prodlen - 1;
90
91 while (pend >= pptr) {
92 *pend-- = prod_vector % BASE;
93 prod_vector /= BASE;
94 }
95 }
96
97 /* Common part of functions bc_standard_mul and bc_standard_square
98 * that takes a vector and converts it to a bc_num */
bc_mul_finish_from_vector(BC_VECTOR * prod_vector,size_t prod_arr_size,size_t prodlen,bc_num * prod)99 static inline void bc_mul_finish_from_vector(BC_VECTOR *prod_vector, size_t prod_arr_size, size_t prodlen, bc_num *prod) {
100 /*
101 * Move a value exceeding 4/8 digits by carrying to the next digit.
102 * However, the last digit does nothing.
103 */
104 bc_mul_carry_calc(prod_vector, prod_arr_size);
105
106 /* Convert to bc_num */
107 *prod = bc_new_num_nonzeroed(prodlen, 0);
108 char *pptr = (*prod)->n_value;
109 char *pend = pptr + prodlen - 1;
110 size_t i = 0;
111 while (i < prod_arr_size - 1) {
112 #if BC_VECTOR_SIZE == 4
113 bc_write_bcd_representation(prod_vector[i], pend - 3);
114 pend -= 4;
115 #else
116 bc_write_bcd_representation(prod_vector[i] / 10000, pend - 7);
117 bc_write_bcd_representation(prod_vector[i] % 10000, pend - 3);
118 pend -= 8;
119 #endif
120 i++;
121 }
122
123 /*
124 * The last digit may carry over.
125 * Also need to fill it to the end with zeros, so loop until the end of the string.
126 */
127 while (pend >= pptr) {
128 *pend-- = prod_vector[i] % BASE;
129 prod_vector[i] /= BASE;
130 }
131 }
132
133 /*
134 * Converts the BCD of bc_num by 4 (32 bits) or 8 (64 bits) digits to an array of BC_VECTOR.
135 * The array is generated starting with the smaller digits.
136 * e.g. 12345678901234567890 => {34567890, 56789012, 1234}
137 *
138 * Multiply and add these groups of numbers to perform multiplication fast.
139 * How much to shift the digits when adding values can be calculated from the index of the array.
140 */
bc_standard_mul(bc_num n1,size_t n1len,bc_num n2,size_t n2len,bc_num * prod)141 static void bc_standard_mul(bc_num n1, size_t n1len, bc_num n2, size_t n2len, bc_num *prod)
142 {
143 size_t i;
144 const char *n1end = n1->n_value + n1len - 1;
145 const char *n2end = n2->n_value + n2len - 1;
146 size_t prodlen = n1len + n2len;
147
148 size_t n1_arr_size = (n1len + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE;
149 size_t n2_arr_size = (n2len + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE;
150 size_t prod_arr_size = (prodlen + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE;
151
152 /*
153 * let's say that N is the max of n1len and n2len (and a multiple of BC_VECTOR_SIZE for simplicity),
154 * then this sum is <= N/BC_VECTOR_SIZE + N/BC_VECTOR_SIZE + N/BC_VECTOR_SIZE + N/BC_VECTOR_SIZE - 1
155 * which is equal to N - 1 if BC_VECTOR_SIZE is 4, and N/2 - 1 if BC_VECTOR_SIZE is 8.
156 */
157 BC_VECTOR *buf = safe_emalloc(n1_arr_size + n2_arr_size + prod_arr_size, sizeof(BC_VECTOR), 0);
158
159 BC_VECTOR *n1_vector = buf;
160 BC_VECTOR *n2_vector = buf + n1_arr_size;
161 BC_VECTOR *prod_vector = n2_vector + n2_arr_size;
162
163 for (i = 0; i < prod_arr_size; i++) {
164 prod_vector[i] = 0;
165 }
166
167 /* Convert to BC_VECTOR[] */
168 bc_convert_to_vector(n1_vector, n1end, n1len);
169 bc_convert_to_vector(n2_vector, n2end, n2len);
170
171 /* Multiplication and addition */
172 size_t count = 0;
173 for (i = 0; i < n1_arr_size; i++) {
174 /*
175 * This calculation adds the result multiple times to the array entries.
176 * When multiplying large numbers of digits, there is a possibility of
177 * overflow, so digit adjustment is performed beforehand.
178 */
179 if (UNEXPECTED(count >= BC_VECTOR_NO_OVERFLOW_ADD_COUNT)) {
180 bc_mul_carry_calc(prod_vector, prod_arr_size);
181 count = 0;
182 }
183 count++;
184 for (size_t j = 0; j < n2_arr_size; j++) {
185 prod_vector[i + j] += n1_vector[i] * n2_vector[j];
186 }
187 }
188
189 bc_mul_finish_from_vector(prod_vector, prod_arr_size, prodlen, prod);
190
191 efree(buf);
192 }
193
194 /** This is bc_standard_mul implementation for square */
bc_standard_square(bc_num n1,size_t n1len,bc_num * prod)195 static void bc_standard_square(bc_num n1, size_t n1len, bc_num *prod)
196 {
197 size_t i;
198 const char *n1end = n1->n_value + n1len - 1;
199 size_t prodlen = n1len + n1len;
200
201 size_t n1_arr_size = (n1len + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE;
202 size_t prod_arr_size = (prodlen + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE;
203
204 BC_VECTOR *buf = safe_emalloc(n1_arr_size + n1_arr_size + prod_arr_size, sizeof(BC_VECTOR), 0);
205
206 BC_VECTOR *n1_vector = buf;
207 BC_VECTOR *prod_vector = n1_vector + n1_arr_size + n1_arr_size;
208
209 for (i = 0; i < prod_arr_size; i++) {
210 prod_vector[i] = 0;
211 }
212
213 /* Convert to BC_VECTOR[] */
214 bc_convert_to_vector(n1_vector, n1end, n1len);
215
216 /* Multiplication and addition */
217 size_t count = 0;
218 for (i = 0; i < n1_arr_size; i++) {
219 /*
220 * This calculation adds the result multiple times to the array entries.
221 * When multiplying large numbers of digits, there is a possibility of
222 * overflow, so digit adjustment is performed beforehand.
223 */
224 if (UNEXPECTED(count >= BC_VECTOR_NO_OVERFLOW_ADD_COUNT)) {
225 bc_mul_carry_calc(prod_vector, prod_arr_size);
226 count = 0;
227 }
228 count++;
229 for (size_t j = 0; j < n1_arr_size; j++) {
230 prod_vector[i + j] += n1_vector[i] * n1_vector[j];
231 }
232 }
233
234 bc_mul_finish_from_vector(prod_vector, prod_arr_size, prodlen, prod);
235
236 efree(buf);
237 }
238
239 /* The multiply routine. N2 times N1 is put int PROD with the scale of
240 the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
241 */
242
bc_multiply(bc_num n1,bc_num n2,size_t scale)243 bc_num bc_multiply(bc_num n1, bc_num n2, size_t scale)
244 {
245 bc_num prod;
246
247 /* Initialize things. */
248 size_t len1 = n1->n_len + n1->n_scale;
249 size_t len2 = n2->n_len + n2->n_scale;
250 size_t full_scale = n1->n_scale + n2->n_scale;
251 size_t prod_scale = MIN(full_scale, MAX(scale, MAX(n1->n_scale, n2->n_scale)));
252
253 /* Do the multiply */
254 if (len1 <= BC_VECTOR_SIZE && len2 <= BC_VECTOR_SIZE) {
255 bc_fast_mul(n1, len1, n2, len2, &prod);
256 } else {
257 bc_standard_mul(n1, len1, n2, len2, &prod);
258 }
259
260 /* Assign to prod and clean up the number. */
261 prod->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS);
262 prod->n_len -= full_scale;
263 prod->n_scale = prod_scale;
264 _bc_rm_leading_zeros(prod);
265 if (bc_is_zero(prod)) {
266 prod->n_sign = PLUS;
267 }
268 return prod;
269 }
270
bc_square(bc_num n1,size_t scale)271 bc_num bc_square(bc_num n1, size_t scale)
272 {
273 bc_num prod;
274
275 size_t len1 = n1->n_len + n1->n_scale;
276 size_t full_scale = n1->n_scale + n1->n_scale;
277 size_t prod_scale = MIN(full_scale, MAX(scale, n1->n_scale));
278
279 if (len1 <= BC_VECTOR_SIZE) {
280 bc_fast_square(n1, len1, &prod);
281 } else {
282 bc_standard_square(n1, len1, &prod);
283 }
284
285 prod->n_sign = PLUS;
286 prod->n_len -= full_scale;
287 prod->n_scale = prod_scale;
288 _bc_rm_leading_zeros(prod);
289
290 return prod;
291 }
292
