xref: /freebsd/contrib/bc/src/num.c (revision cfd6422a5217410fbd66f7a7a8a64d9d85e61229)
1 /*
2  * *****************************************************************************
3  *
4  * SPDX-License-Identifier: BSD-2-Clause
5  *
6  * Copyright (c) 2018-2020 Gavin D. Howard and contributors.
7  *
8  * Redistribution and use in source and binary forms, with or without
9  * modification, are permitted provided that the following conditions are met:
10  *
11  * * Redistributions of source code must retain the above copyright notice, this
12  *   list of conditions and the following disclaimer.
13  *
14  * * Redistributions in binary form must reproduce the above copyright notice,
15  *   this list of conditions and the following disclaimer in the documentation
16  *   and/or other materials provided with the distribution.
17  *
18  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
22  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
28  * POSSIBILITY OF SUCH DAMAGE.
29  *
30  * *****************************************************************************
31  *
32  * Code for the number type.
33  *
34  */
35 
36 #include <assert.h>
37 #include <ctype.h>
38 #include <stdbool.h>
39 #include <stdlib.h>
40 #include <string.h>
41 #include <setjmp.h>
42 #include <limits.h>
43 
44 #include <num.h>
45 #include <rand.h>
46 #include <vm.h>
47 
48 static void bc_num_m(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale);
49 
50 static inline ssize_t bc_num_neg(size_t n, bool neg) {
51 	return (((ssize_t) n) ^ -((ssize_t) neg)) + neg;
52 }
53 
54 ssize_t bc_num_cmpZero(const BcNum *n) {
55 	return bc_num_neg((n)->len != 0, BC_NUM_NEG(n));
56 }
57 
58 static inline size_t bc_num_int(const BcNum *n) {
59 	return n->len ? n->len - BC_NUM_RDX_VAL(n) : 0;
60 }
61 
62 static void bc_num_expand(BcNum *restrict n, size_t req) {
63 
64 	assert(n != NULL);
65 
66 	req = req >= BC_NUM_DEF_SIZE ? req : BC_NUM_DEF_SIZE;
67 
68 	if (req > n->cap) {
69 
70 		BC_SIG_LOCK;
71 
72 		n->num = bc_vm_realloc(n->num, BC_NUM_SIZE(req));
73 		n->cap = req;
74 
75 		BC_SIG_UNLOCK;
76 	}
77 }
78 
79 static void bc_num_setToZero(BcNum *restrict n, size_t scale) {
80 	assert(n != NULL);
81 	n->scale = scale;
82 	n->len = n->rdx = 0;
83 }
84 
85 void bc_num_zero(BcNum *restrict n) {
86 	bc_num_setToZero(n, 0);
87 }
88 
89 void bc_num_one(BcNum *restrict n) {
90 	bc_num_zero(n);
91 	n->len = 1;
92 	n->num[0] = 1;
93 }
94 
95 static void bc_num_clean(BcNum *restrict n) {
96 
97 	while (BC_NUM_NONZERO(n) && !n->num[n->len - 1]) n->len -= 1;
98 
99 	if (BC_NUM_ZERO(n)) n->rdx = 0;
100 	else {
101 		size_t rdx = BC_NUM_RDX_VAL(n);
102 		if (n->len < rdx) n->len = rdx;
103 	}
104 }
105 
106 static size_t bc_num_log10(size_t i) {
107 	size_t len;
108 	for (len = 1; i; i /= BC_BASE, ++len);
109 	assert(len - 1 <= BC_BASE_DIGS + 1);
110 	return len - 1;
111 }
112 
113 static inline size_t bc_num_zeroDigits(const BcDig *n) {
114 	assert(*n >= 0);
115 	assert(((size_t) *n) < BC_BASE_POW);
116 	return BC_BASE_DIGS - bc_num_log10((size_t) *n);
117 }
118 
119 static size_t bc_num_intDigits(const BcNum *n) {
120 	size_t digits = bc_num_int(n) * BC_BASE_DIGS;
121 	if (digits > 0) digits -= bc_num_zeroDigits(n->num + n->len - 1);
122 	return digits;
123 }
124 
125 static size_t bc_num_nonzeroLen(const BcNum *restrict n) {
126 	size_t i, len = n->len;
127 	assert(len == BC_NUM_RDX_VAL(n));
128 	for (i = len - 1; i < len && !n->num[i]; --i);
129 	assert(i + 1 > 0);
130 	return i + 1;
131 }
132 
133 static BcDig bc_num_addDigits(BcDig a, BcDig b, bool *carry) {
134 
135 	assert(((BcBigDig) BC_BASE_POW) * 2 == ((BcDig) BC_BASE_POW) * 2);
136 	assert(a < BC_BASE_POW);
137 	assert(b < BC_BASE_POW);
138 
139 	a += b + *carry;
140 	*carry = (a >= BC_BASE_POW);
141 	if (*carry) a -= BC_BASE_POW;
142 
143 	assert(a >= 0);
144 	assert(a < BC_BASE_POW);
145 
146 	return a;
147 }
148 
149 static BcDig bc_num_subDigits(BcDig a, BcDig b, bool *carry) {
150 
151 	assert(a < BC_BASE_POW);
152 	assert(b < BC_BASE_POW);
153 
154 	b += *carry;
155 	*carry = (a < b);
156 	if (*carry) a += BC_BASE_POW;
157 
158 	assert(a - b >= 0);
159 	assert(a - b < BC_BASE_POW);
160 
161 	return a - b;
162 }
163 
164 static void bc_num_addArrays(BcDig *restrict a, const BcDig *restrict b,
165                              size_t len)
166 {
167 	size_t i;
168 	bool carry = false;
169 
170 	for (i = 0; i < len; ++i) a[i] = bc_num_addDigits(a[i], b[i], &carry);
171 
172 	for (; carry; ++i) a[i] = bc_num_addDigits(a[i], 0, &carry);
173 }
174 
175 static void bc_num_subArrays(BcDig *restrict a, const BcDig *restrict b,
176                              size_t len)
177 {
178 	size_t i;
179 	bool carry = false;
180 
181 	for (i = 0; i < len; ++i) a[i] = bc_num_subDigits(a[i], b[i], &carry);
182 
183 	for (; carry; ++i) a[i] = bc_num_subDigits(a[i], 0, &carry);
184 }
185 
186 static void bc_num_mulArray(const BcNum *restrict a, BcBigDig b,
187                             BcNum *restrict c)
188 {
189 	size_t i;
190 	BcBigDig carry = 0;
191 
192 	assert(b <= BC_BASE_POW);
193 
194 	if (a->len + 1 > c->cap) bc_num_expand(c, a->len + 1);
195 
196 	memset(c->num, 0, BC_NUM_SIZE(c->cap));
197 
198 	for (i = 0; i < a->len; ++i) {
199 		BcBigDig in = ((BcBigDig) a->num[i]) * b + carry;
200 		c->num[i] = in % BC_BASE_POW;
201 		carry = in / BC_BASE_POW;
202 	}
203 
204 	assert(carry < BC_BASE_POW);
205 	c->num[i] = (BcDig) carry;
206 	c->len = a->len;
207 	c->len += (carry != 0);
208 
209 	bc_num_clean(c);
210 
211 	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
212 	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
213 	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
214 }
215 
216 static void bc_num_divArray(const BcNum *restrict a, BcBigDig b,
217                             BcNum *restrict c, BcBigDig *rem)
218 {
219 	size_t i;
220 	BcBigDig carry = 0;
221 
222 	assert(c->cap >= a->len);
223 
224 	for (i = a->len - 1; i < a->len; --i) {
225 		BcBigDig in = ((BcBigDig) a->num[i]) + carry * BC_BASE_POW;
226 		assert(in / b < BC_BASE_POW);
227 		c->num[i] = (BcDig) (in / b);
228 		carry = in % b;
229 	}
230 
231 	c->len = a->len;
232 	bc_num_clean(c);
233 	*rem = carry;
234 
235 	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
236 	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
237 	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
238 }
239 
240 static ssize_t bc_num_compare(const BcDig *restrict a, const BcDig *restrict b,
241                               size_t len)
242 {
243 	size_t i;
244 	BcDig c = 0;
245 	for (i = len - 1; i < len && !(c = a[i] - b[i]); --i);
246 	return bc_num_neg(i + 1, c < 0);
247 }
248 
249 ssize_t bc_num_cmp(const BcNum *a, const BcNum *b) {
250 
251 	size_t i, min, a_int, b_int, diff, ardx, brdx;
252 	BcDig *max_num, *min_num;
253 	bool a_max, neg = false;
254 	ssize_t cmp;
255 
256 	assert(a != NULL && b != NULL);
257 
258 	if (a == b) return 0;
259 	if (BC_NUM_ZERO(a)) return bc_num_neg(b->len != 0, !BC_NUM_NEG(b));
260 	if (BC_NUM_ZERO(b)) return bc_num_cmpZero(a);
261 	if (BC_NUM_NEG(a)) {
262 		if (BC_NUM_NEG(b)) neg = true;
263 		else return -1;
264 	}
265 	else if (BC_NUM_NEG(b)) return 1;
266 
267 	a_int = bc_num_int(a);
268 	b_int = bc_num_int(b);
269 	a_int -= b_int;
270 
271 	if (a_int) return neg ? -((ssize_t) a_int) : (ssize_t) a_int;
272 
273 	ardx = BC_NUM_RDX_VAL(a);
274 	brdx = BC_NUM_RDX_VAL(b);
275 	a_max = (ardx > brdx);
276 
277 	if (a_max) {
278 		min = brdx;
279 		diff = ardx - brdx;
280 		max_num = a->num + diff;
281 		min_num = b->num;
282 	}
283 	else {
284 		min = ardx;
285 		diff = brdx - ardx;
286 		max_num = b->num + diff;
287 		min_num = a->num;
288 	}
289 
290 	cmp = bc_num_compare(max_num, min_num, b_int + min);
291 
292 	if (cmp) return bc_num_neg((size_t) cmp, !a_max == !neg);
293 
294 	for (max_num -= diff, i = diff - 1; i < diff; --i) {
295 		if (max_num[i]) return bc_num_neg(1, !a_max == !neg);
296 	}
297 
298 	return 0;
299 }
300 
301 void bc_num_truncate(BcNum *restrict n, size_t places) {
302 
303 	size_t nrdx, places_rdx;
304 
305 	if (!places) return;
306 
307 	nrdx = BC_NUM_RDX_VAL(n);
308 	places_rdx = nrdx ? nrdx - BC_NUM_RDX(n->scale - places) : 0;
309 	assert(places <= n->scale && (BC_NUM_ZERO(n) || places_rdx <= n->len));
310 
311 	n->scale -= places;
312 	BC_NUM_RDX_SET(n, nrdx - places_rdx);
313 
314 	if (BC_NUM_NONZERO(n)) {
315 
316 		size_t pow;
317 
318 		pow = n->scale % BC_BASE_DIGS;
319 		pow = pow ? BC_BASE_DIGS - pow : 0;
320 		pow = bc_num_pow10[pow];
321 
322 		n->len -= places_rdx;
323 		memmove(n->num, n->num + places_rdx, BC_NUM_SIZE(n->len));
324 
325 		// Clear the lower part of the last digit.
326 		if (BC_NUM_NONZERO(n)) n->num[0] -= n->num[0] % (BcDig) pow;
327 
328 		bc_num_clean(n);
329 	}
330 }
331 
332 void bc_num_extend(BcNum *restrict n, size_t places) {
333 
334 	size_t nrdx, places_rdx;
335 
336 	if (!places) return;
337 	if (BC_NUM_ZERO(n)) {
338 		n->scale += places;
339 		return;
340 	}
341 
342 	nrdx = BC_NUM_RDX_VAL(n);
343 	places_rdx = BC_NUM_RDX(places + n->scale) - nrdx;
344 
345 	if (places_rdx) {
346 		bc_num_expand(n, bc_vm_growSize(n->len, places_rdx));
347 		memmove(n->num + places_rdx, n->num, BC_NUM_SIZE(n->len));
348 		memset(n->num, 0, BC_NUM_SIZE(places_rdx));
349 	}
350 
351 	BC_NUM_RDX_SET(n, nrdx + places_rdx);
352 	n->scale += places;
353 	n->len += places_rdx;
354 
355 	assert(BC_NUM_RDX_VAL(n) == BC_NUM_RDX(n->scale));
356 }
357 
358 static void bc_num_retireMul(BcNum *restrict n, size_t scale,
359                              bool neg1, bool neg2)
360 {
361 	if (n->scale < scale) bc_num_extend(n, scale - n->scale);
362 	else bc_num_truncate(n, n->scale - scale);
363 
364 	bc_num_clean(n);
365 	if (BC_NUM_NONZERO(n)) n->rdx = BC_NUM_NEG_VAL(n, !neg1 != !neg2);
366 }
367 
368 static void bc_num_split(const BcNum *restrict n, size_t idx,
369                          BcNum *restrict a, BcNum *restrict b)
370 {
371 	assert(BC_NUM_ZERO(a));
372 	assert(BC_NUM_ZERO(b));
373 
374 	if (idx < n->len) {
375 
376 		b->len = n->len - idx;
377 		a->len = idx;
378 		a->scale = b->scale = 0;
379 		BC_NUM_RDX_SET(a, 0);
380 		BC_NUM_RDX_SET(b, 0);
381 
382 		assert(a->cap >= a->len);
383 		assert(b->cap >= b->len);
384 
385 		memcpy(b->num, n->num + idx, BC_NUM_SIZE(b->len));
386 		memcpy(a->num, n->num, BC_NUM_SIZE(idx));
387 
388 		bc_num_clean(b);
389 	}
390 	else bc_num_copy(a, n);
391 
392 	bc_num_clean(a);
393 }
394 
395 static size_t bc_num_shiftZero(BcNum *restrict n) {
396 
397 	size_t i;
398 
399 	assert(!BC_NUM_RDX_VAL(n) || BC_NUM_ZERO(n));
400 
401 	for (i = 0; i < n->len && !n->num[i]; ++i);
402 
403 	n->len -= i;
404 	n->num += i;
405 
406 	return i;
407 }
408 
409 static void bc_num_unshiftZero(BcNum *restrict n, size_t places_rdx) {
410 	n->len += places_rdx;
411 	n->num -= places_rdx;
412 }
413 
414 static void bc_num_shift(BcNum *restrict n, BcBigDig dig) {
415 
416 	size_t i, len = n->len;
417 	BcBigDig carry = 0, pow;
418 	BcDig *ptr = n->num;
419 
420 	assert(dig < BC_BASE_DIGS);
421 
422 	pow = bc_num_pow10[dig];
423 	dig = bc_num_pow10[BC_BASE_DIGS - dig];
424 
425 	for (i = len - 1; i < len; --i) {
426 		BcBigDig in, temp;
427 		in = ((BcBigDig) ptr[i]);
428 		temp = carry * dig;
429 		carry = in % pow;
430 		ptr[i] = ((BcDig) (in / pow)) + (BcDig) temp;
431 	}
432 
433 	assert(!carry);
434 }
435 
436 static void bc_num_shiftLeft(BcNum *restrict n, size_t places) {
437 
438 	BcBigDig dig;
439 	size_t places_rdx;
440 	bool shift;
441 
442 	if (!places) return;
443 	if (places > n->scale) {
444 		size_t size = bc_vm_growSize(BC_NUM_RDX(places - n->scale), n->len);
445 		if (size > SIZE_MAX - 1) bc_vm_err(BC_ERR_MATH_OVERFLOW);
446 	}
447 	if (BC_NUM_ZERO(n)) {
448 		if (n->scale >= places) n->scale -= places;
449 		else n->scale = 0;
450 		return;
451 	}
452 
453 	dig = (BcBigDig) (places % BC_BASE_DIGS);
454 	shift = (dig != 0);
455 	places_rdx = BC_NUM_RDX(places);
456 
457 	if (n->scale) {
458 
459 		size_t nrdx = BC_NUM_RDX_VAL(n);
460 
461 		if (nrdx >= places_rdx) {
462 
463 			size_t mod = n->scale % BC_BASE_DIGS, revdig;
464 
465 			mod = mod ? mod : BC_BASE_DIGS;
466 			revdig = dig ? BC_BASE_DIGS - dig : 0;
467 
468 			if (mod + revdig > BC_BASE_DIGS) places_rdx = 1;
469 			else places_rdx = 0;
470 		}
471 		else places_rdx -= nrdx;
472 	}
473 
474 	if (places_rdx) {
475 		bc_num_expand(n, bc_vm_growSize(n->len, places_rdx));
476 		memmove(n->num + places_rdx, n->num, BC_NUM_SIZE(n->len));
477 		memset(n->num, 0, BC_NUM_SIZE(places_rdx));
478 		n->len += places_rdx;
479 	}
480 
481 	if (places > n->scale) {
482 		n->scale = 0;
483 		BC_NUM_RDX_SET(n, 0);
484 	}
485 	else {
486 		n->scale -= places;
487 		BC_NUM_RDX_SET(n, BC_NUM_RDX(n->scale));
488 	}
489 
490 	if (shift) bc_num_shift(n, BC_BASE_DIGS - dig);
491 
492 	bc_num_clean(n);
493 }
494 
495 void bc_num_shiftRight(BcNum *restrict n, size_t places) {
496 
497 	BcBigDig dig;
498 	size_t places_rdx, scale, scale_mod, int_len, expand;
499 	bool shift;
500 
501 	if (!places) return;
502 	if (BC_NUM_ZERO(n)) {
503 		n->scale += places;
504 		bc_num_expand(n, BC_NUM_RDX(n->scale));
505 		return;
506 	}
507 
508 	dig = (BcBigDig) (places % BC_BASE_DIGS);
509 	shift = (dig != 0);
510 	scale = n->scale;
511 	scale_mod = scale % BC_BASE_DIGS;
512 	scale_mod = scale_mod ? scale_mod : BC_BASE_DIGS;
513 	int_len = bc_num_int(n);
514 	places_rdx = BC_NUM_RDX(places);
515 
516 	if (scale_mod + dig > BC_BASE_DIGS) {
517 		expand = places_rdx - 1;
518 		places_rdx = 1;
519 	}
520 	else {
521 		expand = places_rdx;
522 		places_rdx = 0;
523 	}
524 
525 	if (expand > int_len) expand -= int_len;
526 	else expand = 0;
527 
528 	bc_num_extend(n, places_rdx * BC_BASE_DIGS);
529 	bc_num_expand(n, bc_vm_growSize(expand, n->len));
530 	memset(n->num + n->len, 0, BC_NUM_SIZE(expand));
531 	n->len += expand;
532 	n->scale = 0;
533 	BC_NUM_RDX_SET(n, 0);
534 
535 	if (shift) bc_num_shift(n, dig);
536 
537 	n->scale = scale + places;
538 	BC_NUM_RDX_SET(n, BC_NUM_RDX(n->scale));
539 
540 	bc_num_clean(n);
541 
542 	assert(BC_NUM_RDX_VAL(n) <= n->len && n->len <= n->cap);
543 	assert(BC_NUM_RDX_VAL(n) == BC_NUM_RDX(n->scale));
544 }
545 
546 static void bc_num_inv(BcNum *a, BcNum *b, size_t scale) {
547 
548 	BcNum one;
549 	BcDig num[2];
550 
551 	assert(BC_NUM_NONZERO(a));
552 
553 	bc_num_setup(&one, num, sizeof(num) / sizeof(BcDig));
554 	bc_num_one(&one);
555 
556 	bc_num_div(&one, a, b, scale);
557 }
558 
559 #if BC_ENABLE_EXTRA_MATH
560 static void bc_num_intop(const BcNum *a, const BcNum *b, BcNum *restrict c,
561                          BcBigDig *v)
562 {
563 	if (BC_ERR(BC_NUM_RDX_VAL(b))) bc_vm_err(BC_ERR_MATH_NON_INTEGER);
564 	bc_num_copy(c, a);
565 	bc_num_bigdig(b, v);
566 }
567 #endif // BC_ENABLE_EXTRA_MATH
568 
569 static void bc_num_as(BcNum *a, BcNum *b, BcNum *restrict c, size_t sub) {
570 
571 	BcDig *ptr_c, *ptr_l, *ptr_r;
572 	size_t i, min_rdx, max_rdx, diff, a_int, b_int, min_len, max_len, max_int;
573 	size_t len_l, len_r, ardx, brdx;
574 	bool b_neg, do_sub, do_rev_sub, carry, c_neg;
575 
576 	// Because this function doesn't need to use scale (per the bc spec),
577 	// I am hijacking it to say whether it's doing an add or a subtract.
578 	// Convert substraction to addition of negative second operand.
579 
580 	if (BC_NUM_ZERO(b)) {
581 		bc_num_copy(c, a);
582 		return;
583 	}
584 	if (BC_NUM_ZERO(a)) {
585 		bc_num_copy(c, b);
586 		c->rdx = BC_NUM_NEG_VAL(c, BC_NUM_NEG(b) != sub);
587 		return;
588 	}
589 
590 	// Invert sign of b if it is to be subtracted. This operation must
591 	// preced the tests for any of the operands being zero.
592 	b_neg = (BC_NUM_NEG(b) != sub);
593 
594 	// Actually add the numbers if their signs are equal, else subtract.
595 	do_sub = (BC_NUM_NEG(a) != b_neg);
596 
597 	a_int = bc_num_int(a);
598 	b_int = bc_num_int(b);
599 	max_int = BC_MAX(a_int, b_int);
600 
601 	ardx = BC_NUM_RDX_VAL(a);
602 	brdx = BC_NUM_RDX_VAL(b);
603 	min_rdx = BC_MIN(ardx, brdx);
604 	max_rdx = BC_MAX(ardx, brdx);
605 	diff = max_rdx - min_rdx;
606 
607 	max_len = max_int + max_rdx;
608 
609 	if (do_sub) {
610 
611 		// Check whether b has to be subtracted from a or a from b.
612 		if (a_int != b_int) do_rev_sub = (a_int < b_int);
613 		else if (ardx > brdx)
614 			do_rev_sub = (bc_num_compare(a->num + diff, b->num, b->len) < 0);
615 		else
616 			do_rev_sub = (bc_num_compare(a->num, b->num + diff, a->len) <= 0);
617 	}
618 	else {
619 
620 		// The result array of the addition might come out one element
621 		// longer than the bigger of the operand arrays.
622 		max_len += 1;
623 		do_rev_sub = (a_int < b_int);
624 	}
625 
626 	assert(max_len <= c->cap);
627 
628 	if (do_rev_sub) {
629 		ptr_l = b->num;
630 		ptr_r = a->num;
631 		len_l = b->len;
632 		len_r = a->len;
633 	}
634 	else {
635 		ptr_l = a->num;
636 		ptr_r = b->num;
637 		len_l = a->len;
638 		len_r = b->len;
639 	}
640 
641 	ptr_c = c->num;
642 	carry = false;
643 
644 	if (diff) {
645 
646 		// If the rdx values of the operands do not match, the result will
647 		// have low end elements that are the positive or negative trailing
648 		// elements of the operand with higher rdx value.
649 		if ((ardx > brdx) != do_rev_sub) {
650 
651 			// !do_rev_sub && ardx > brdx || do_rev_sub && brdx > ardx
652 			// The left operand has BcDig values that need to be copied,
653 			// either from a or from b (in case of a reversed subtraction).
654 			memcpy(ptr_c, ptr_l, BC_NUM_SIZE(diff));
655 			ptr_l += diff;
656 			len_l -= diff;
657 		}
658 		else {
659 
660 			// The right operand has BcDig values that need to be copied
661 			// or subtracted from zero (in case of a subtraction).
662 			if (do_sub) {
663 
664 				// do_sub (do_rev_sub && ardx > brdx ||
665 				// !do_rev_sub && brdx > ardx)
666 				for (i = 0; i < diff; i++)
667 					ptr_c[i] = bc_num_subDigits(0, ptr_r[i], &carry);
668 			}
669 			else {
670 
671 				// !do_sub && brdx > ardx
672 				memcpy(ptr_c, ptr_r, BC_NUM_SIZE(diff));
673 			}
674 
675 			ptr_r += diff;
676 			len_r -= diff;
677 		}
678 
679 		ptr_c += diff;
680 	}
681 
682 	min_len = BC_MIN(len_l, len_r);
683 
684 	// After dealing with possible low array elements that depend on only one
685 	// operand, the actual add or subtract can be performed as if the rdx of
686 	// both operands was the same.
687 	// Inlining takes care of eliminating constant zero arguments to
688 	// addDigit/subDigit (checked in disassembly of resulting bc binary
689 	// compiled with gcc and clang).
690 	if (do_sub) {
691 		for (i = 0; i < min_len; ++i)
692 			ptr_c[i] = bc_num_subDigits(ptr_l[i], ptr_r[i], &carry);
693 		for (; i < len_l; ++i) ptr_c[i] = bc_num_subDigits(ptr_l[i], 0, &carry);
694 	}
695 	else {
696 		for (i = 0; i < min_len; ++i)
697 			ptr_c[i] = bc_num_addDigits(ptr_l[i], ptr_r[i], &carry);
698 		for (; i < len_l; ++i) ptr_c[i] = bc_num_addDigits(ptr_l[i], 0, &carry);
699 		ptr_c[i] = bc_num_addDigits(0, 0, &carry);
700 	}
701 
702 	assert(carry == false);
703 
704 	// The result has the same sign as a, unless the operation was a
705 	// reverse subtraction (b - a).
706 	c_neg = BC_NUM_NEG(a) != (do_sub && do_rev_sub);
707 	BC_NUM_RDX_SET_NEG(c, max_rdx, c_neg);
708 	c->len = max_len;
709 	c->scale = BC_MAX(a->scale, b->scale);
710 
711 	bc_num_clean(c);
712 }
713 
714 static void bc_num_m_simp(const BcNum *a, const BcNum *b, BcNum *restrict c)
715 {
716 	size_t i, alen = a->len, blen = b->len, clen;
717 	BcDig *ptr_a = a->num, *ptr_b = b->num, *ptr_c;
718 	BcBigDig sum = 0, carry = 0;
719 
720 	assert(sizeof(sum) >= sizeof(BcDig) * 2);
721 	assert(!BC_NUM_RDX_VAL(a) && !BC_NUM_RDX_VAL(b));
722 
723 	clen = bc_vm_growSize(alen, blen);
724 	bc_num_expand(c, bc_vm_growSize(clen, 1));
725 
726 	ptr_c = c->num;
727 	memset(ptr_c, 0, BC_NUM_SIZE(c->cap));
728 
729 	for (i = 0; i < clen; ++i) {
730 
731 		ssize_t sidx = (ssize_t) (i - blen + 1);
732 		size_t j = (size_t) BC_MAX(0, sidx), k = BC_MIN(i, blen - 1);
733 
734 		for (; j < alen && k < blen; ++j, --k) {
735 
736 			sum += ((BcBigDig) ptr_a[j]) * ((BcBigDig) ptr_b[k]);
737 
738 			if (sum >= ((BcBigDig) BC_BASE_POW) * BC_BASE_POW) {
739 				carry += sum / BC_BASE_POW;
740 				sum %= BC_BASE_POW;
741 			}
742 		}
743 
744 		if (sum >= BC_BASE_POW) {
745 			carry += sum / BC_BASE_POW;
746 			sum %= BC_BASE_POW;
747 		}
748 
749 		ptr_c[i] = (BcDig) sum;
750 		assert(ptr_c[i] < BC_BASE_POW);
751 		sum = carry;
752 		carry = 0;
753 	}
754 
755 	// This should always be true because there should be no carry on the last
756 	// digit; multiplication never goes above the sum of both lengths.
757 	assert(!sum);
758 
759 	c->len = clen;
760 }
761 
762 static void bc_num_shiftAddSub(BcNum *restrict n, const BcNum *restrict a,
763                                size_t shift, BcNumShiftAddOp op)
764 {
765 	assert(n->len >= shift + a->len);
766 	assert(!BC_NUM_RDX_VAL(n) && !BC_NUM_RDX_VAL(a));
767 	op(n->num + shift, a->num, a->len);
768 }
769 
770 static void bc_num_k(BcNum *a, BcNum *b, BcNum *restrict c) {
771 
772 	size_t max, max2, total;
773 	BcNum l1, h1, l2, h2, m2, m1, z0, z1, z2, temp;
774 	BcDig *digs, *dig_ptr;
775 	BcNumShiftAddOp op;
776 	bool aone = BC_NUM_ONE(a);
777 
778 	assert(BC_NUM_ZERO(c));
779 
780 	if (BC_NUM_ZERO(a) || BC_NUM_ZERO(b)) return;
781 	if (aone || BC_NUM_ONE(b)) {
782 		bc_num_copy(c, aone ? b : a);
783 		if ((aone && BC_NUM_NEG(a)) || BC_NUM_NEG(b)) BC_NUM_NEG_TGL(c);
784 		return;
785 	}
786 	if (a->len < BC_NUM_KARATSUBA_LEN || b->len < BC_NUM_KARATSUBA_LEN) {
787 		bc_num_m_simp(a, b, c);
788 		return;
789 	}
790 
791 	max = BC_MAX(a->len, b->len);
792 	max = BC_MAX(max, BC_NUM_DEF_SIZE);
793 	max2 = (max + 1) / 2;
794 
795 	total = bc_vm_arraySize(BC_NUM_KARATSUBA_ALLOCS, max);
796 
797 	BC_SIG_LOCK;
798 
799 	digs = dig_ptr = bc_vm_malloc(BC_NUM_SIZE(total));
800 
801 	bc_num_setup(&l1, dig_ptr, max);
802 	dig_ptr += max;
803 	bc_num_setup(&h1, dig_ptr, max);
804 	dig_ptr += max;
805 	bc_num_setup(&l2, dig_ptr, max);
806 	dig_ptr += max;
807 	bc_num_setup(&h2, dig_ptr, max);
808 	dig_ptr += max;
809 	bc_num_setup(&m1, dig_ptr, max);
810 	dig_ptr += max;
811 	bc_num_setup(&m2, dig_ptr, max);
812 	max = bc_vm_growSize(max, 1);
813 	bc_num_init(&z0, max);
814 	bc_num_init(&z1, max);
815 	bc_num_init(&z2, max);
816 	max = bc_vm_growSize(max, max) + 1;
817 	bc_num_init(&temp, max);
818 
819 	BC_SETJMP_LOCKED(err);
820 
821 	BC_SIG_UNLOCK;
822 
823 	bc_num_split(a, max2, &l1, &h1);
824 	bc_num_split(b, max2, &l2, &h2);
825 
826 	bc_num_expand(c, max);
827 	c->len = max;
828 	memset(c->num, 0, BC_NUM_SIZE(c->len));
829 
830 	bc_num_sub(&h1, &l1, &m1, 0);
831 	bc_num_sub(&l2, &h2, &m2, 0);
832 
833 	if (BC_NUM_NONZERO(&h1) && BC_NUM_NONZERO(&h2)) {
834 
835 		assert(BC_NUM_RDX_VALID_NP(h1));
836 		assert(BC_NUM_RDX_VALID_NP(h2));
837 
838 		bc_num_m(&h1, &h2, &z2, 0);
839 		bc_num_clean(&z2);
840 
841 		bc_num_shiftAddSub(c, &z2, max2 * 2, bc_num_addArrays);
842 		bc_num_shiftAddSub(c, &z2, max2, bc_num_addArrays);
843 	}
844 
845 	if (BC_NUM_NONZERO(&l1) && BC_NUM_NONZERO(&l2)) {
846 
847 		assert(BC_NUM_RDX_VALID_NP(l1));
848 		assert(BC_NUM_RDX_VALID_NP(l2));
849 
850 		bc_num_m(&l1, &l2, &z0, 0);
851 		bc_num_clean(&z0);
852 
853 		bc_num_shiftAddSub(c, &z0, max2, bc_num_addArrays);
854 		bc_num_shiftAddSub(c, &z0, 0, bc_num_addArrays);
855 	}
856 
857 	if (BC_NUM_NONZERO(&m1) && BC_NUM_NONZERO(&m2)) {
858 
859 		assert(BC_NUM_RDX_VALID_NP(m1));
860 		assert(BC_NUM_RDX_VALID_NP(m1));
861 
862 		bc_num_m(&m1, &m2, &z1, 0);
863 		bc_num_clean(&z1);
864 
865 		op = (BC_NUM_NEG_NP(m1) != BC_NUM_NEG_NP(m2)) ?
866 		     bc_num_subArrays : bc_num_addArrays;
867 		bc_num_shiftAddSub(c, &z1, max2, op);
868 	}
869 
870 err:
871 	BC_SIG_MAYLOCK;
872 	free(digs);
873 	bc_num_free(&temp);
874 	bc_num_free(&z2);
875 	bc_num_free(&z1);
876 	bc_num_free(&z0);
877 	BC_LONGJMP_CONT;
878 }
879 
880 static void bc_num_m(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
881 
882 	BcNum cpa, cpb;
883 	size_t ascale, bscale, ardx, brdx, azero = 0, bzero = 0, zero, len, rscale;
884 
885 	assert(BC_NUM_RDX_VALID(a));
886 	assert(BC_NUM_RDX_VALID(b));
887 
888 	bc_num_zero(c);
889 	ascale = a->scale;
890 	bscale = b->scale;
891 	scale = BC_MAX(scale, ascale);
892 	scale = BC_MAX(scale, bscale);
893 
894 	rscale = ascale + bscale;
895 	scale = BC_MIN(rscale, scale);
896 
897 	if ((a->len == 1 || b->len == 1) && !a->rdx && !b->rdx) {
898 
899 		BcNum *operand;
900 		BcBigDig dig;
901 
902 		if (a->len == 1) {
903 			dig = (BcBigDig) a->num[0];
904 			operand = b;
905 		}
906 		else {
907 			dig = (BcBigDig) b->num[0];
908 			operand = a;
909 		}
910 
911 		bc_num_mulArray(operand, dig, c);
912 
913 		if (BC_NUM_NONZERO(c))
914 			c->rdx = BC_NUM_NEG_VAL(c, BC_NUM_NEG(a) != BC_NUM_NEG(b));
915 
916 		return;
917 	}
918 
919 	assert(BC_NUM_RDX_VALID(a));
920 	assert(BC_NUM_RDX_VALID(b));
921 
922 	BC_SIG_LOCK;
923 
924 	bc_num_init(&cpa, a->len + BC_NUM_RDX_VAL(a));
925 	bc_num_init(&cpb, b->len + BC_NUM_RDX_VAL(b));
926 
927 	BC_SETJMP_LOCKED(err);
928 
929 	BC_SIG_UNLOCK;
930 
931 	bc_num_copy(&cpa, a);
932 	bc_num_copy(&cpb, b);
933 
934 	assert(BC_NUM_RDX_VALID_NP(cpa));
935 	assert(BC_NUM_RDX_VALID_NP(cpb));
936 
937 	BC_NUM_NEG_CLR_NP(cpa);
938 	BC_NUM_NEG_CLR_NP(cpb);
939 
940 	assert(BC_NUM_RDX_VALID_NP(cpa));
941 	assert(BC_NUM_RDX_VALID_NP(cpb));
942 
943 	ardx = BC_NUM_RDX_VAL_NP(cpa) * BC_BASE_DIGS;
944 	bc_num_shiftLeft(&cpa, ardx);
945 
946 	brdx = BC_NUM_RDX_VAL_NP(cpb) * BC_BASE_DIGS;
947 	bc_num_shiftLeft(&cpb, brdx);
948 
949 	// We need to reset the jump here because azero and bzero are used in the
950 	// cleanup, and local variables are not guaranteed to be the same after a
951 	// jump.
952 	BC_SIG_LOCK;
953 
954 	BC_UNSETJMP;
955 
956 	azero = bc_num_shiftZero(&cpa);
957 	bzero = bc_num_shiftZero(&cpb);
958 
959 	BC_SETJMP_LOCKED(err);
960 
961 	BC_SIG_UNLOCK;
962 
963 	bc_num_clean(&cpa);
964 	bc_num_clean(&cpb);
965 
966 	bc_num_k(&cpa, &cpb, c);
967 
968 	zero = bc_vm_growSize(azero, bzero);
969 	len = bc_vm_growSize(c->len, zero);
970 
971 	bc_num_expand(c, len);
972 	bc_num_shiftLeft(c, (len - c->len) * BC_BASE_DIGS);
973 	bc_num_shiftRight(c, ardx + brdx);
974 
975 	bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));
976 
977 err:
978 	BC_SIG_MAYLOCK;
979 	bc_num_unshiftZero(&cpb, bzero);
980 	bc_num_free(&cpb);
981 	bc_num_unshiftZero(&cpa, azero);
982 	bc_num_free(&cpa);
983 	BC_LONGJMP_CONT;
984 }
985 
986 static bool bc_num_nonZeroDig(BcDig *restrict a, size_t len) {
987 	size_t i;
988 	bool nonzero = false;
989 	for (i = len - 1; !nonzero && i < len; --i) nonzero = (a[i] != 0);
990 	return nonzero;
991 }
992 
993 static ssize_t bc_num_divCmp(const BcDig *a, const BcNum *b, size_t len) {
994 
995 	ssize_t cmp;
996 
997 	if (b->len > len && a[len]) cmp = bc_num_compare(a, b->num, len + 1);
998 	else if (b->len <= len) {
999 		if (a[len]) cmp = 1;
1000 		else cmp = bc_num_compare(a, b->num, len);
1001 	}
1002 	else cmp = -1;
1003 
1004 	return cmp;
1005 }
1006 
1007 static void bc_num_divExtend(BcNum *restrict a, BcNum *restrict b,
1008                              BcBigDig divisor)
1009 {
1010 	size_t pow;
1011 
1012 	assert(divisor < BC_BASE_POW);
1013 
1014 	pow = BC_BASE_DIGS - bc_num_log10((size_t) divisor);
1015 
1016 	bc_num_shiftLeft(a, pow);
1017 	bc_num_shiftLeft(b, pow);
1018 }
1019 
1020 static void bc_num_d_long(BcNum *restrict a, BcNum *restrict b,
1021                           BcNum *restrict c, size_t scale)
1022 {
1023 	BcBigDig divisor;
1024 	size_t len, end, i, rdx;
1025 	BcNum cpb;
1026 	bool nonzero = false;
1027 
1028 	assert(b->len < a->len);
1029 	len = b->len;
1030 	end = a->len - len;
1031 	assert(len >= 1);
1032 
1033 	bc_num_expand(c, a->len);
1034 	memset(c->num, 0, c->cap * sizeof(BcDig));
1035 
1036 	BC_NUM_RDX_SET(c, BC_NUM_RDX_VAL(a));
1037 	c->scale = a->scale;
1038 	c->len = a->len;
1039 
1040 	divisor = (BcBigDig) b->num[len - 1];
1041 
1042 	if (len > 1 && bc_num_nonZeroDig(b->num, len - 1)) {
1043 
1044 		nonzero = (divisor > 1 << ((10 * BC_BASE_DIGS) / 6 + 1));
1045 
1046 		if (!nonzero) {
1047 
1048 			bc_num_divExtend(a, b, divisor);
1049 
1050 			len = BC_MAX(a->len, b->len);
1051 			bc_num_expand(a, len + 1);
1052 
1053 			if (len + 1 > a->len) a->len = len + 1;
1054 
1055 			len = b->len;
1056 			end = a->len - len;
1057 			divisor = (BcBigDig) b->num[len - 1];
1058 
1059 			nonzero = bc_num_nonZeroDig(b->num, len - 1);
1060 		}
1061 	}
1062 
1063 	divisor += nonzero;
1064 
1065 	bc_num_expand(c, a->len);
1066 	memset(c->num, 0, BC_NUM_SIZE(c->cap));
1067 
1068 	assert(c->scale >= scale);
1069 	rdx = BC_NUM_RDX_VAL(c) - BC_NUM_RDX(scale);
1070 
1071 	BC_SIG_LOCK;
1072 
1073 	bc_num_init(&cpb, len + 1);
1074 
1075 	BC_SETJMP_LOCKED(err);
1076 
1077 	BC_SIG_UNLOCK;
1078 
1079 	i = end - 1;
1080 
1081 	for (; i < end && i >= rdx && BC_NUM_NONZERO(a); --i) {
1082 
1083 		ssize_t cmp;
1084 		BcDig *n;
1085 		BcBigDig result;
1086 
1087 		n = a->num + i;
1088 		assert(n >= a->num);
1089 		result = 0;
1090 
1091 		cmp = bc_num_divCmp(n, b, len);
1092 
1093 		while (cmp >= 0) {
1094 
1095 			BcBigDig n1, dividend, q;
1096 
1097 			n1 = (BcBigDig) n[len];
1098 			dividend = n1 * BC_BASE_POW + (BcBigDig) n[len - 1];
1099 			q = (dividend / divisor);
1100 
1101 			if (q <= 1) {
1102 				q = 1;
1103 				bc_num_subArrays(n, b->num, len);
1104 			}
1105 			else {
1106 
1107 				assert(q <= BC_BASE_POW);
1108 
1109 				bc_num_mulArray(b, (BcBigDig) q, &cpb);
1110 				bc_num_subArrays(n, cpb.num, cpb.len);
1111 			}
1112 
1113 			result += q;
1114 			assert(result <= BC_BASE_POW);
1115 
1116 			if (nonzero) cmp = bc_num_divCmp(n, b, len);
1117 			else cmp = -1;
1118 		}
1119 
1120 		assert(result < BC_BASE_POW);
1121 
1122 		c->num[i] = (BcDig) result;
1123 	}
1124 
1125 err:
1126 	BC_SIG_MAYLOCK;
1127 	bc_num_free(&cpb);
1128 	BC_LONGJMP_CONT;
1129 }
1130 
1131 static void bc_num_d(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1132 
1133 	size_t len, cpardx;
1134 	BcNum cpa, cpb;
1135 
1136 	if (BC_NUM_ZERO(b)) bc_vm_err(BC_ERR_MATH_DIVIDE_BY_ZERO);
1137 	if (BC_NUM_ZERO(a)) {
1138 		bc_num_setToZero(c, scale);
1139 		return;
1140 	}
1141 	if (BC_NUM_ONE(b)) {
1142 		bc_num_copy(c, a);
1143 		bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));
1144 		return;
1145 	}
1146 	if (!BC_NUM_RDX_VAL(a) && !BC_NUM_RDX_VAL(b) && b->len == 1 && !scale) {
1147 		BcBigDig rem;
1148 		bc_num_divArray(a, (BcBigDig) b->num[0], c, &rem);
1149 		bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));
1150 		return;
1151 	}
1152 
1153 	len = bc_num_divReq(a, b, scale);
1154 
1155 	BC_SIG_LOCK;
1156 
1157 	bc_num_init(&cpa, len);
1158 	bc_num_copy(&cpa, a);
1159 	bc_num_createCopy(&cpb, b);
1160 
1161 	BC_SETJMP_LOCKED(err);
1162 
1163 	BC_SIG_UNLOCK;
1164 
1165 	len = b->len;
1166 
1167 	if (len > cpa.len) {
1168 		bc_num_expand(&cpa, bc_vm_growSize(len, 2));
1169 		bc_num_extend(&cpa, (len - cpa.len) * BC_BASE_DIGS);
1170 	}
1171 
1172 	cpardx = BC_NUM_RDX_VAL_NP(cpa);
1173 	cpa.scale = cpardx * BC_BASE_DIGS;
1174 
1175 	bc_num_extend(&cpa, b->scale);
1176 	cpardx = BC_NUM_RDX_VAL_NP(cpa) - BC_NUM_RDX(b->scale);
1177 	BC_NUM_RDX_SET_NP(cpa, cpardx);
1178 	cpa.scale = cpardx * BC_BASE_DIGS;
1179 
1180 	if (scale > cpa.scale) {
1181 		bc_num_extend(&cpa, scale);
1182 		cpardx = BC_NUM_RDX_VAL_NP(cpa);
1183 		cpa.scale = cpardx * BC_BASE_DIGS;
1184 	}
1185 
1186 	if (cpa.cap == cpa.len) bc_num_expand(&cpa, bc_vm_growSize(cpa.len, 1));
1187 
1188 	// We want an extra zero in front to make things simpler.
1189 	cpa.num[cpa.len++] = 0;
1190 
1191 	if (cpardx == cpa.len) cpa.len = bc_num_nonzeroLen(&cpa);
1192 	if (BC_NUM_RDX_VAL_NP(cpb) == cpb.len) cpb.len = bc_num_nonzeroLen(&cpb);
1193 	cpb.scale = 0;
1194 	BC_NUM_RDX_SET_NP(cpb, 0);
1195 
1196 	bc_num_d_long(&cpa, &cpb, c, scale);
1197 
1198 	bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));
1199 
1200 err:
1201 	BC_SIG_MAYLOCK;
1202 	bc_num_free(&cpb);
1203 	bc_num_free(&cpa);
1204 	BC_LONGJMP_CONT;
1205 }
1206 
1207 static void bc_num_r(BcNum *a, BcNum *b, BcNum *restrict c,
1208                      BcNum *restrict d, size_t scale, size_t ts)
1209 {
1210 	BcNum temp;
1211 	bool neg;
1212 
1213 	if (BC_NUM_ZERO(b)) bc_vm_err(BC_ERR_MATH_DIVIDE_BY_ZERO);
1214 	if (BC_NUM_ZERO(a)) {
1215 		bc_num_setToZero(c, ts);
1216 		bc_num_setToZero(d, ts);
1217 		return;
1218 	}
1219 
1220 	BC_SIG_LOCK;
1221 
1222 	bc_num_init(&temp, d->cap);
1223 
1224 	BC_SETJMP_LOCKED(err);
1225 
1226 	BC_SIG_UNLOCK;
1227 
1228 	bc_num_d(a, b, c, scale);
1229 
1230 	if (scale) scale = ts + 1;
1231 
1232 	assert(BC_NUM_RDX_VALID(c));
1233 	assert(BC_NUM_RDX_VALID(b));
1234 
1235 	bc_num_m(c, b, &temp, scale);
1236 	bc_num_sub(a, &temp, d, scale);
1237 
1238 	if (ts > d->scale && BC_NUM_NONZERO(d)) bc_num_extend(d, ts - d->scale);
1239 
1240 	neg = BC_NUM_NEG(d);
1241 	bc_num_retireMul(d, ts, BC_NUM_NEG(a), BC_NUM_NEG(b));
1242 	d->rdx = BC_NUM_NEG_VAL(d, BC_NUM_NONZERO(d) ? neg : false);
1243 
1244 err:
1245 	BC_SIG_MAYLOCK;
1246 	bc_num_free(&temp);
1247 	BC_LONGJMP_CONT;
1248 }
1249 
1250 static void bc_num_rem(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1251 
1252 	BcNum c1;
1253 	size_t ts;
1254 
1255 	ts = bc_vm_growSize(scale, b->scale);
1256 	ts = BC_MAX(ts, a->scale);
1257 
1258 	BC_SIG_LOCK;
1259 
1260 	bc_num_init(&c1, bc_num_mulReq(a, b, ts));
1261 
1262 	BC_SETJMP_LOCKED(err);
1263 
1264 	BC_SIG_UNLOCK;
1265 
1266 	bc_num_r(a, b, &c1, c, scale, ts);
1267 
1268 err:
1269 	BC_SIG_MAYLOCK;
1270 	bc_num_free(&c1);
1271 	BC_LONGJMP_CONT;
1272 }
1273 
1274 static void bc_num_p(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1275 
1276 	BcNum copy;
1277 	BcBigDig pow = 0;
1278 	size_t i, powrdx, resrdx;
1279 	bool neg, zero;
1280 
1281 	if (BC_ERR(BC_NUM_RDX_VAL(b))) bc_vm_err(BC_ERR_MATH_NON_INTEGER);
1282 
1283 	if (BC_NUM_ZERO(b)) {
1284 		bc_num_one(c);
1285 		return;
1286 	}
1287 	if (BC_NUM_ZERO(a)) {
1288 		if (BC_NUM_NEG(b)) bc_vm_err(BC_ERR_MATH_DIVIDE_BY_ZERO);
1289 		bc_num_setToZero(c, scale);
1290 		return;
1291 	}
1292 	if (BC_NUM_ONE(b)) {
1293 		if (!BC_NUM_NEG(b)) bc_num_copy(c, a);
1294 		else bc_num_inv(a, c, scale);
1295 		return;
1296 	}
1297 
1298 	BC_SIG_LOCK;
1299 
1300 	neg = BC_NUM_NEG(b);
1301 	BC_NUM_NEG_CLR(b);
1302 	bc_num_bigdig(b, &pow);
1303 	b->rdx = BC_NUM_NEG_VAL(b, neg);
1304 
1305 	bc_num_createCopy(&copy, a);
1306 
1307 	BC_SETJMP_LOCKED(err);
1308 
1309 	BC_SIG_UNLOCK;
1310 
1311 	if (!neg) {
1312 		size_t max = BC_MAX(scale, a->scale), scalepow = a->scale * pow;
1313 		scale = BC_MIN(scalepow, max);
1314 	}
1315 
1316 	for (powrdx = a->scale; !(pow & 1); pow >>= 1) {
1317 		powrdx <<= 1;
1318 		assert(BC_NUM_RDX_VALID_NP(copy));
1319 		bc_num_mul(&copy, &copy, &copy, powrdx);
1320 	}
1321 
1322 	bc_num_copy(c, &copy);
1323 	resrdx = powrdx;
1324 
1325 	while (pow >>= 1) {
1326 
1327 		powrdx <<= 1;
1328 		assert(BC_NUM_RDX_VALID_NP(copy));
1329 		bc_num_mul(&copy, &copy, &copy, powrdx);
1330 
1331 		if (pow & 1) {
1332 			resrdx += powrdx;
1333 			assert(BC_NUM_RDX_VALID(c));
1334 			assert(BC_NUM_RDX_VALID_NP(copy));
1335 			bc_num_mul(c, &copy, c, resrdx);
1336 		}
1337 	}
1338 
1339 	if (neg) bc_num_inv(c, c, scale);
1340 
1341 	if (c->scale > scale) bc_num_truncate(c, c->scale - scale);
1342 
1343 	// We can't use bc_num_clean() here.
1344 	for (zero = true, i = 0; zero && i < c->len; ++i) zero = !c->num[i];
1345 	if (zero) bc_num_setToZero(c, scale);
1346 
1347 err:
1348 	BC_SIG_MAYLOCK;
1349 	bc_num_free(&copy);
1350 	BC_LONGJMP_CONT;
1351 }
1352 
1353 #if BC_ENABLE_EXTRA_MATH
1354 static void bc_num_place(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1355 
1356 	BcBigDig val = 0;
1357 
1358 	BC_UNUSED(scale);
1359 
1360 	bc_num_intop(a, b, c, &val);
1361 
1362 	if (val < c->scale) bc_num_truncate(c, c->scale - val);
1363 	else if (val > c->scale) bc_num_extend(c, val - c->scale);
1364 }
1365 
1366 static void bc_num_left(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1367 
1368 	BcBigDig val = 0;
1369 
1370 	BC_UNUSED(scale);
1371 
1372 	bc_num_intop(a, b, c, &val);
1373 
1374 	bc_num_shiftLeft(c, (size_t) val);
1375 }
1376 
1377 static void bc_num_right(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1378 
1379 	BcBigDig val = 0;
1380 
1381 	BC_UNUSED(scale);
1382 
1383 	bc_num_intop(a, b, c, &val);
1384 
1385 	if (BC_NUM_ZERO(c)) return;
1386 
1387 	bc_num_shiftRight(c, (size_t) val);
1388 }
1389 #endif // BC_ENABLE_EXTRA_MATH
1390 
1391 static void bc_num_binary(BcNum *a, BcNum *b, BcNum *c, size_t scale,
1392                           BcNumBinaryOp op, size_t req)
1393 {
1394 	BcNum *ptr_a, *ptr_b, num2;
1395 	bool init = false;
1396 
1397 	assert(a != NULL && b != NULL && c != NULL && op != NULL);
1398 
1399 	assert(BC_NUM_RDX_VALID(a));
1400 	assert(BC_NUM_RDX_VALID(b));
1401 
1402 	BC_SIG_LOCK;
1403 
1404 	if (c == a) {
1405 
1406 		ptr_a = &num2;
1407 
1408 		memcpy(ptr_a, c, sizeof(BcNum));
1409 		init = true;
1410 	}
1411 	else {
1412 		ptr_a = a;
1413 	}
1414 
1415 	if (c == b) {
1416 
1417 		ptr_b = &num2;
1418 
1419 		if (c != a) {
1420 			memcpy(ptr_b, c, sizeof(BcNum));
1421 			init = true;
1422 		}
1423 	}
1424 	else {
1425 		ptr_b = b;
1426 	}
1427 
1428 	if (init) {
1429 
1430 		bc_num_init(c, req);
1431 
1432 		BC_SETJMP_LOCKED(err);
1433 		BC_SIG_UNLOCK;
1434 	}
1435 	else {
1436 		BC_SIG_UNLOCK;
1437 		bc_num_expand(c, req);
1438 	}
1439 
1440 	op(ptr_a, ptr_b, c, scale);
1441 
1442 	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
1443 	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
1444 	assert(BC_NUM_RDX_VALID(c));
1445 	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
1446 
1447 err:
1448 	if (init) {
1449 		BC_SIG_MAYLOCK;
1450 		bc_num_free(&num2);
1451 		BC_LONGJMP_CONT;
1452 	}
1453 }
1454 
1455 #if !defined(NDEBUG) || BC_ENABLE_LIBRARY
1456 bool bc_num_strValid(const char *restrict val) {
1457 
1458 	bool radix = false;
1459 	size_t i, len = strlen(val);
1460 
1461 	if (!len) return true;
1462 
1463 	for (i = 0; i < len; ++i) {
1464 
1465 		BcDig c = val[i];
1466 
1467 		if (c == '.') {
1468 
1469 			if (radix) return false;
1470 
1471 			radix = true;
1472 			continue;
1473 		}
1474 
1475 		if (!(isdigit(c) || isupper(c))) return false;
1476 	}
1477 
1478 	return true;
1479 }
1480 #endif // !defined(NDEBUG) || BC_ENABLE_LIBRARY
1481 
1482 static BcBigDig bc_num_parseChar(char c, size_t base_t) {
1483 
1484 	if (isupper(c)) {
1485 		c = BC_NUM_NUM_LETTER(c);
1486 		c = ((size_t) c) >= base_t ? (char) base_t - 1 : c;
1487 	}
1488 	else c -= '0';
1489 
1490 	return (BcBigDig) (uchar) c;
1491 }
1492 
1493 static void bc_num_parseDecimal(BcNum *restrict n, const char *restrict val) {
1494 
1495 	size_t len, i, temp, mod;
1496 	const char *ptr;
1497 	bool zero = true, rdx;
1498 
1499 	for (i = 0; val[i] == '0'; ++i);
1500 
1501 	val += i;
1502 	assert(!val[0] || isalnum(val[0]) || val[0] == '.');
1503 
1504 	// All 0's. We can just return, since this
1505 	// procedure expects a virgin (already 0) BcNum.
1506 	if (!val[0]) return;
1507 
1508 	len = strlen(val);
1509 
1510 	ptr = strchr(val, '.');
1511 	rdx = (ptr != NULL);
1512 
1513 	for (i = 0; i < len && (zero = (val[i] == '0' || val[i] == '.')); ++i);
1514 
1515 	n->scale = (size_t) (rdx * (((uintptr_t) (val + len)) -
1516 	                            (((uintptr_t) ptr) + 1)));
1517 
1518 	BC_NUM_RDX_SET(n, BC_NUM_RDX(n->scale));
1519 	i = len - (ptr == val ? 0 : i) - rdx;
1520 	temp = BC_NUM_ROUND_POW(i);
1521 	mod = n->scale % BC_BASE_DIGS;
1522 	i = mod ? BC_BASE_DIGS - mod : 0;
1523 	n->len = ((temp + i) / BC_BASE_DIGS);
1524 
1525 	bc_num_expand(n, n->len);
1526 	memset(n->num, 0, BC_NUM_SIZE(n->len));
1527 
1528 	if (zero) {
1529 		// I think I can set rdx directly to zero here because n should be a
1530 		// new number with sign set to false.
1531 		n->len = n->rdx = 0;
1532 	}
1533 	else {
1534 
1535 		BcBigDig exp, pow;
1536 
1537 		assert(i <= BC_NUM_BIGDIG_MAX);
1538 
1539 		exp = (BcBigDig) i;
1540 		pow = bc_num_pow10[exp];
1541 
1542 		for (i = len - 1; i < len; --i, ++exp) {
1543 
1544 			char c = val[i];
1545 
1546 			if (c == '.') exp -= 1;
1547 			else {
1548 
1549 				size_t idx = exp / BC_BASE_DIGS;
1550 
1551 				if (isupper(c)) c = '9';
1552 				n->num[idx] += (((BcBigDig) c) - '0') * pow;
1553 
1554 				if ((exp + 1) % BC_BASE_DIGS == 0) pow = 1;
1555 				else pow *= BC_BASE;
1556 			}
1557 		}
1558 	}
1559 }
1560 
1561 static void bc_num_parseBase(BcNum *restrict n, const char *restrict val,
1562                              BcBigDig base)
1563 {
1564 	BcNum temp, mult1, mult2, result1, result2, *m1, *m2, *ptr;
1565 	char c = 0;
1566 	bool zero = true;
1567 	BcBigDig v;
1568 	size_t i, digs, len = strlen(val);
1569 
1570 	for (i = 0; zero && i < len; ++i) zero = (val[i] == '.' || val[i] == '0');
1571 	if (zero) return;
1572 
1573 	BC_SIG_LOCK;
1574 
1575 	bc_num_init(&temp, BC_NUM_BIGDIG_LOG10);
1576 	bc_num_init(&mult1, BC_NUM_BIGDIG_LOG10);
1577 
1578 	BC_SETJMP_LOCKED(int_err);
1579 
1580 	BC_SIG_UNLOCK;
1581 
1582 	for (i = 0; i < len && (c = val[i]) && c != '.'; ++i) {
1583 
1584 		v = bc_num_parseChar(c, base);
1585 
1586 		bc_num_mulArray(n, base, &mult1);
1587 		bc_num_bigdig2num(&temp, v);
1588 		bc_num_add(&mult1, &temp, n, 0);
1589 	}
1590 
1591 	if (i == len && !(c = val[i])) goto int_err;
1592 
1593 	assert(c == '.');
1594 
1595 	BC_SIG_LOCK;
1596 
1597 	BC_UNSETJMP;
1598 
1599 	bc_num_init(&mult2, BC_NUM_BIGDIG_LOG10);
1600 	bc_num_init(&result1, BC_NUM_DEF_SIZE);
1601 	bc_num_init(&result2, BC_NUM_DEF_SIZE);
1602 	bc_num_one(&mult1);
1603 
1604 	BC_SETJMP_LOCKED(err);
1605 
1606 	BC_SIG_UNLOCK;
1607 
1608 	m1 = &mult1;
1609 	m2 = &mult2;
1610 
1611 	for (i += 1, digs = 0; i < len && (c = val[i]); ++i, ++digs) {
1612 
1613 		size_t rdx;
1614 
1615 		v = bc_num_parseChar(c, base);
1616 
1617 		bc_num_mulArray(&result1, base, &result2);
1618 
1619 		bc_num_bigdig2num(&temp, v);
1620 		bc_num_add(&result2, &temp, &result1, 0);
1621 		bc_num_mulArray(m1, base, m2);
1622 
1623 		rdx = BC_NUM_RDX_VAL(m2);
1624 
1625 		if (m2->len < rdx) m2->len = rdx;
1626 
1627 		ptr = m1;
1628 		m1 = m2;
1629 		m2 = ptr;
1630 	}
1631 
1632 	// This one cannot be a divide by 0 because mult starts out at 1, then is
1633 	// multiplied by base, and base cannot be 0, so mult cannot be 0.
1634 	bc_num_div(&result1, m1, &result2, digs * 2);
1635 	bc_num_truncate(&result2, digs);
1636 	bc_num_add(n, &result2, n, digs);
1637 
1638 	if (BC_NUM_NONZERO(n)) {
1639 		if (n->scale < digs) bc_num_extend(n, digs - n->scale);
1640 	}
1641 	else bc_num_zero(n);
1642 
1643 err:
1644 	BC_SIG_MAYLOCK;
1645 	bc_num_free(&result2);
1646 	bc_num_free(&result1);
1647 	bc_num_free(&mult2);
1648 int_err:
1649 	BC_SIG_MAYLOCK;
1650 	bc_num_free(&mult1);
1651 	bc_num_free(&temp);
1652 	BC_LONGJMP_CONT;
1653 }
1654 
1655 static inline void bc_num_printNewline(void) {
1656 #if !BC_ENABLE_LIBRARY
1657 	if (vm.nchars >= vm.line_len - 1) {
1658 		bc_vm_putchar('\\');
1659 		bc_vm_putchar('\n');
1660 	}
1661 #endif // !BC_ENABLE_LIBRARY
1662 }
1663 
1664 static void bc_num_putchar(int c) {
1665 	if (c != '\n') bc_num_printNewline();
1666 	bc_vm_putchar(c);
1667 }
1668 
1669 #if DC_ENABLED && !BC_ENABLE_LIBRARY
1670 static void bc_num_printChar(size_t n, size_t len, bool rdx) {
1671 	BC_UNUSED(rdx);
1672 	BC_UNUSED(len);
1673 	assert(len == 1);
1674 	bc_vm_putchar((uchar) n);
1675 }
1676 #endif // DC_ENABLED && !BC_ENABLE_LIBRARY
1677 
1678 static void bc_num_printDigits(size_t n, size_t len, bool rdx) {
1679 
1680 	size_t exp, pow;
1681 
1682 	bc_num_putchar(rdx ? '.' : ' ');
1683 
1684 	for (exp = 0, pow = 1; exp < len - 1; ++exp, pow *= BC_BASE);
1685 
1686 	for (exp = 0; exp < len; pow /= BC_BASE, ++exp) {
1687 		size_t dig = n / pow;
1688 		n -= dig * pow;
1689 		bc_num_putchar(((uchar) dig) + '0');
1690 	}
1691 }
1692 
1693 static void bc_num_printHex(size_t n, size_t len, bool rdx) {
1694 
1695 	BC_UNUSED(len);
1696 
1697 	assert(len == 1);
1698 
1699 	if (rdx) bc_num_putchar('.');
1700 
1701 	bc_num_putchar(bc_num_hex_digits[n]);
1702 }
1703 
1704 static void bc_num_printDecimal(const BcNum *restrict n) {
1705 
1706 	size_t i, j, rdx = BC_NUM_RDX_VAL(n);
1707 	bool zero = true;
1708 	size_t buffer[BC_BASE_DIGS];
1709 
1710 	if (BC_NUM_NEG(n)) bc_num_putchar('-');
1711 
1712 	for (i = n->len - 1; i < n->len; --i) {
1713 
1714 		BcDig n9 = n->num[i];
1715 		size_t temp;
1716 		bool irdx = (i == rdx - 1);
1717 
1718 		zero = (zero & !irdx);
1719 		temp = n->scale % BC_BASE_DIGS;
1720 		temp = i || !temp ? 0 : BC_BASE_DIGS - temp;
1721 
1722 		memset(buffer, 0, BC_BASE_DIGS * sizeof(size_t));
1723 
1724 		for (j = 0; n9 && j < BC_BASE_DIGS; ++j) {
1725 			buffer[j] = ((size_t) n9) % BC_BASE;
1726 			n9 /= BC_BASE;
1727 		}
1728 
1729 		for (j = BC_BASE_DIGS - 1; j < BC_BASE_DIGS && j >= temp; --j) {
1730 			bool print_rdx = (irdx & (j == BC_BASE_DIGS - 1));
1731 			zero = (zero && buffer[j] == 0);
1732 			if (!zero) bc_num_printHex(buffer[j], 1, print_rdx);
1733 		}
1734 	}
1735 }
1736 
1737 #if BC_ENABLE_EXTRA_MATH
1738 static void bc_num_printExponent(const BcNum *restrict n, bool eng) {
1739 
1740 	size_t places, mod, nrdx = BC_NUM_RDX_VAL(n);
1741 	bool neg = (n->len <= nrdx);
1742 	BcNum temp, exp;
1743 	BcDig digs[BC_NUM_BIGDIG_LOG10];
1744 
1745 	BC_SIG_LOCK;
1746 
1747 	bc_num_createCopy(&temp, n);
1748 
1749 	BC_SETJMP_LOCKED(exit);
1750 
1751 	BC_SIG_UNLOCK;
1752 
1753 	if (neg) {
1754 
1755 		size_t i, idx = bc_num_nonzeroLen(n) - 1;
1756 
1757 		places = 1;
1758 
1759 		for (i = BC_BASE_DIGS - 1; i < BC_BASE_DIGS; --i) {
1760 			if (bc_num_pow10[i] > (BcBigDig) n->num[idx]) places += 1;
1761 			else break;
1762 		}
1763 
1764 		places += (nrdx - (idx + 1)) * BC_BASE_DIGS;
1765 		mod = places % 3;
1766 
1767 		if (eng && mod != 0) places += 3 - mod;
1768 		bc_num_shiftLeft(&temp, places);
1769 	}
1770 	else {
1771 		places = bc_num_intDigits(n) - 1;
1772 		mod = places % 3;
1773 		if (eng && mod != 0) places -= 3 - (3 - mod);
1774 		bc_num_shiftRight(&temp, places);
1775 	}
1776 
1777 	bc_num_printDecimal(&temp);
1778 	bc_num_putchar('e');
1779 
1780 	if (!places) {
1781 		bc_num_printHex(0, 1, false);
1782 		goto exit;
1783 	}
1784 
1785 	if (neg) bc_num_putchar('-');
1786 
1787 	bc_num_setup(&exp, digs, BC_NUM_BIGDIG_LOG10);
1788 	bc_num_bigdig2num(&exp, (BcBigDig) places);
1789 
1790 	bc_num_printDecimal(&exp);
1791 
1792 exit:
1793 	BC_SIG_MAYLOCK;
1794 	bc_num_free(&temp);
1795 	BC_LONGJMP_CONT;
1796 }
1797 #endif // BC_ENABLE_EXTRA_MATH
1798 
1799 static void bc_num_printFixup(BcNum *restrict n, BcBigDig rem,
1800                               BcBigDig pow, size_t idx)
1801 {
1802 	size_t i, len = n->len - idx;
1803 	BcBigDig acc;
1804 	BcDig *a = n->num + idx;
1805 
1806 	if (len < 2) return;
1807 
1808 	for (i = len - 1; i > 0; --i) {
1809 
1810 		acc = ((BcBigDig) a[i]) * rem + ((BcBigDig) a[i - 1]);
1811 		a[i - 1] = (BcDig) (acc % pow);
1812 		acc /= pow;
1813 		acc += (BcBigDig) a[i];
1814 
1815 		if (acc >= BC_BASE_POW) {
1816 
1817 			if (i == len - 1) {
1818 				len = bc_vm_growSize(len, 1);
1819 				bc_num_expand(n, bc_vm_growSize(len, idx));
1820 				a = n->num + idx;
1821 				a[len - 1] = 0;
1822 			}
1823 
1824 			a[i + 1] += acc / BC_BASE_POW;
1825 			acc %= BC_BASE_POW;
1826 		}
1827 
1828 		assert(acc < BC_BASE_POW);
1829 		a[i] = (BcDig) acc;
1830 	}
1831 
1832 	n->len = len + idx;
1833 }
1834 
1835 static void bc_num_printPrepare(BcNum *restrict n, BcBigDig rem,
1836                                 BcBigDig pow)
1837 {
1838 	size_t i;
1839 
1840 	for (i = 0; i < n->len; ++i) bc_num_printFixup(n, rem, pow, i);
1841 
1842 	for (i = 0; i < n->len; ++i) {
1843 
1844 		assert(pow == ((BcBigDig) ((BcDig) pow)));
1845 
1846 		if (n->num[i] >= (BcDig) pow) {
1847 
1848 			if (i + 1 == n->len) {
1849 				n->len = bc_vm_growSize(n->len, 1);
1850 				bc_num_expand(n, n->len);
1851 				n->num[i + 1] = 0;
1852 			}
1853 
1854 			assert(pow < BC_BASE_POW);
1855 			n->num[i + 1] += n->num[i] / ((BcDig) pow);
1856 			n->num[i] %= (BcDig) pow;
1857 		}
1858 	}
1859 }
1860 
1861 static void bc_num_printNum(BcNum *restrict n, BcBigDig base,
1862                             size_t len, BcNumDigitOp print)
1863 {
1864 	BcVec stack;
1865 	BcNum intp, fracp1, fracp2, digit, flen1, flen2, *n1, *n2, *temp;
1866 	BcBigDig dig = 0, *ptr, acc, exp;
1867 	size_t i, j, nrdx;
1868 	bool radix;
1869 	BcDig digit_digs[BC_NUM_BIGDIG_LOG10 + 1];
1870 
1871 	assert(base > 1);
1872 
1873 	if (BC_NUM_ZERO(n)) {
1874 		print(0, len, false);
1875 		return;
1876 	}
1877 
1878 	// This function uses an algorithm that Stefan Esser <se@freebsd.org> came
1879 	// up with to print the integer part of a number. What it does is convert
1880 	// intp into a number of the specified base, but it does it directly,
1881 	// instead of just doing a series of divisions and printing the remainders
1882 	// in reverse order.
1883 	//
1884 	// Let me explain in a bit more detail:
1885 	//
1886 	// The algorithm takes the current least significant digit (after intp has
1887 	// been converted to an integer) and the next to least significant digit,
1888 	// and it converts the least significant digit into one of the specified
1889 	// base, putting any overflow into the next to least significant digit. It
1890 	// iterates through the whole number, from least significant to most
1891 	// significant, doing this conversion. At the end of that iteration, the
1892 	// least significant digit is converted, but the others are not, so it
1893 	// iterates again, starting at the next to least significant digit. It keeps
1894 	// doing that conversion, skipping one more digit than the last time, until
1895 	// all digits have been converted. Then it prints them in reverse order.
1896 	//
1897 	// That is the gist of the algorithm. It leaves out several things, such as
1898 	// the fact that digits are not always converted into the specified base,
1899 	// but into something close, basically a power of the specified base. In
1900 	// Stefan's words, "You could consider BcDigs to be of base 10^BC_BASE_DIGS
1901 	// in the normal case and obase^N for the largest value of N that satisfies
1902 	// obase^N <= 10^BC_BASE_DIGS. [This means that] the result is not in base
1903 	// "obase", but in base "obase^N", which happens to be printable as a number
1904 	// of base "obase" without consideration for neighbouring BcDigs." This fact
1905 	// is what necessitates the existence of the loop later in this function.
1906 	//
1907 	// The conversion happens in bc_num_printPrepare() where the outer loop
1908 	// happens and bc_num_printFixup() where the inner loop, or actual
1909 	// conversion, happens.
1910 
1911 	nrdx = BC_NUM_RDX_VAL(n);
1912 
1913 	BC_SIG_LOCK;
1914 
1915 	bc_vec_init(&stack, sizeof(BcBigDig), NULL);
1916 	bc_num_init(&fracp1, nrdx);
1917 
1918 	bc_num_createCopy(&intp, n);
1919 
1920 	BC_SETJMP_LOCKED(err);
1921 
1922 	BC_SIG_UNLOCK;
1923 
1924 	bc_num_truncate(&intp, intp.scale);
1925 
1926 	bc_num_sub(n, &intp, &fracp1, 0);
1927 
1928 	if (base != vm.last_base) {
1929 
1930 		vm.last_pow = 1;
1931 		vm.last_exp = 0;
1932 
1933 		while (vm.last_pow * base <= BC_BASE_POW) {
1934 			vm.last_pow *= base;
1935 			vm.last_exp += 1;
1936 		}
1937 
1938 		vm.last_rem = BC_BASE_POW - vm.last_pow;
1939 		vm.last_base = base;
1940 	}
1941 
1942 	exp = vm.last_exp;
1943 
1944 	if (vm.last_rem != 0) bc_num_printPrepare(&intp, vm.last_rem, vm.last_pow);
1945 
1946 	for (i = 0; i < intp.len; ++i) {
1947 
1948 		acc = (BcBigDig) intp.num[i];
1949 
1950 		for (j = 0; j < exp && (i < intp.len - 1 || acc != 0); ++j)
1951 		{
1952 			if (j != exp - 1) {
1953 				dig = acc % base;
1954 				acc /= base;
1955 			}
1956 			else {
1957 				dig = acc;
1958 				acc = 0;
1959 			}
1960 
1961 			assert(dig < base);
1962 
1963 			bc_vec_push(&stack, &dig);
1964 		}
1965 
1966 		assert(acc == 0);
1967 	}
1968 
1969 	for (i = 0; i < stack.len; ++i) {
1970 		ptr = bc_vec_item_rev(&stack, i);
1971 		assert(ptr != NULL);
1972 		print(*ptr, len, false);
1973 	}
1974 
1975 	if (!n->scale) goto err;
1976 
1977 	BC_SIG_LOCK;
1978 
1979 	BC_UNSETJMP;
1980 
1981 	bc_num_init(&fracp2, nrdx);
1982 	bc_num_setup(&digit, digit_digs, sizeof(digit_digs) / sizeof(BcDig));
1983 	bc_num_init(&flen1, BC_NUM_BIGDIG_LOG10);
1984 	bc_num_init(&flen2, BC_NUM_BIGDIG_LOG10);
1985 
1986 	BC_SETJMP_LOCKED(frac_err);
1987 
1988 	BC_SIG_UNLOCK;
1989 
1990 	bc_num_one(&flen1);
1991 
1992 	radix = true;
1993 	n1 = &flen1;
1994 	n2 = &flen2;
1995 
1996 	fracp2.scale = n->scale;
1997 	BC_NUM_RDX_SET_NP(fracp2, BC_NUM_RDX(fracp2.scale));
1998 
1999 	while (bc_num_intDigits(n1) < n->scale + 1) {
2000 
2001 		bc_num_expand(&fracp2, fracp1.len + 1);
2002 		bc_num_mulArray(&fracp1, base, &fracp2);
2003 
2004 		nrdx = BC_NUM_RDX_VAL_NP(fracp2);
2005 
2006 		if (fracp2.len < nrdx) fracp2.len = nrdx;
2007 
2008 		// fracp is guaranteed to be non-negative and small enough.
2009 		bc_num_bigdig2(&fracp2, &dig);
2010 
2011 		bc_num_bigdig2num(&digit, dig);
2012 		bc_num_sub(&fracp2, &digit, &fracp1, 0);
2013 
2014 		print(dig, len, radix);
2015 		bc_num_mulArray(n1, base, n2);
2016 
2017 		radix = false;
2018 		temp = n1;
2019 		n1 = n2;
2020 		n2 = temp;
2021 	}
2022 
2023 frac_err:
2024 	BC_SIG_MAYLOCK;
2025 	bc_num_free(&flen2);
2026 	bc_num_free(&flen1);
2027 	bc_num_free(&fracp2);
2028 err:
2029 	BC_SIG_MAYLOCK;
2030 	bc_num_free(&fracp1);
2031 	bc_num_free(&intp);
2032 	bc_vec_free(&stack);
2033 	BC_LONGJMP_CONT;
2034 }
2035 
2036 static void bc_num_printBase(BcNum *restrict n, BcBigDig base) {
2037 
2038 	size_t width;
2039 	BcNumDigitOp print;
2040 	bool neg = BC_NUM_NEG(n);
2041 
2042 	if (neg) bc_num_putchar('-');
2043 
2044 	BC_NUM_NEG_CLR(n);
2045 
2046 	if (base <= BC_NUM_MAX_POSIX_IBASE) {
2047 		width = 1;
2048 		print = bc_num_printHex;
2049 	}
2050 	else {
2051 		assert(base <= BC_BASE_POW);
2052 		width = bc_num_log10(base - 1);
2053 		print = bc_num_printDigits;
2054 	}
2055 
2056 	bc_num_printNum(n, base, width, print);
2057 	n->rdx = BC_NUM_NEG_VAL(n, neg);
2058 }
2059 
2060 #if DC_ENABLED && !BC_ENABLE_LIBRARY
2061 void bc_num_stream(BcNum *restrict n, BcBigDig base) {
2062 	bc_num_printNum(n, base, 1, bc_num_printChar);
2063 }
2064 #endif // DC_ENABLED && !BC_ENABLE_LIBRARY
2065 
2066 void bc_num_setup(BcNum *restrict n, BcDig *restrict num, size_t cap) {
2067 	assert(n != NULL);
2068 	n->num = num;
2069 	n->cap = cap;
2070 	bc_num_zero(n);
2071 }
2072 
2073 void bc_num_init(BcNum *restrict n, size_t req) {
2074 
2075 	BcDig *num;
2076 
2077 	BC_SIG_ASSERT_LOCKED;
2078 
2079 	assert(n != NULL);
2080 
2081 	req = req >= BC_NUM_DEF_SIZE ? req : BC_NUM_DEF_SIZE;
2082 
2083 	if (req == BC_NUM_DEF_SIZE && vm.temps.len) {
2084 		BcNum *nptr = bc_vec_top(&vm.temps);
2085 		num = nptr->num;
2086 		bc_vec_pop(&vm.temps);
2087 	}
2088 	else num = bc_vm_malloc(BC_NUM_SIZE(req));
2089 
2090 	bc_num_setup(n, num, req);
2091 }
2092 
2093 void bc_num_clear(BcNum *restrict n) {
2094 	n->num = NULL;
2095 	n->cap = 0;
2096 }
2097 
2098 void bc_num_free(void *num) {
2099 
2100 	BcNum *n = (BcNum*) num;
2101 
2102 	BC_SIG_ASSERT_LOCKED;
2103 
2104 	assert(n != NULL);
2105 
2106 	if (n->cap == BC_NUM_DEF_SIZE) bc_vec_push(&vm.temps, n);
2107 	else free(n->num);
2108 }
2109 
2110 void bc_num_copy(BcNum *d, const BcNum *s) {
2111 	assert(d != NULL && s != NULL);
2112 	if (d == s) return;
2113 	bc_num_expand(d, s->len);
2114 	d->len = s->len;
2115 	// I can just copy directly here.
2116 	d->rdx = s->rdx;
2117 	d->scale = s->scale;
2118 	memcpy(d->num, s->num, BC_NUM_SIZE(d->len));
2119 }
2120 
2121 void bc_num_createCopy(BcNum *d, const BcNum *s) {
2122 	BC_SIG_ASSERT_LOCKED;
2123 	bc_num_init(d, s->len);
2124 	bc_num_copy(d, s);
2125 }
2126 
2127 void bc_num_createFromBigdig(BcNum *n, BcBigDig val) {
2128 	BC_SIG_ASSERT_LOCKED;
2129 	bc_num_init(n, BC_NUM_BIGDIG_LOG10);
2130 	bc_num_bigdig2num(n, val);
2131 }
2132 
2133 size_t bc_num_scale(const BcNum *restrict n) {
2134 	return n->scale;
2135 }
2136 
2137 size_t bc_num_len(const BcNum *restrict n) {
2138 
2139 	size_t len = n->len;
2140 
2141 	if (BC_NUM_ZERO(n)) return 0;
2142 
2143 	if (BC_NUM_RDX_VAL(n) == len) {
2144 
2145 		size_t zero, scale;
2146 
2147 		len = bc_num_nonzeroLen(n);
2148 
2149 		scale = n->scale % BC_BASE_DIGS;
2150 		scale = scale ? scale : BC_BASE_DIGS;
2151 
2152 		zero = bc_num_zeroDigits(n->num + len - 1);
2153 
2154 		len = len * BC_BASE_DIGS - zero - (BC_BASE_DIGS - scale);
2155 	}
2156 	else len = bc_num_intDigits(n) + n->scale;
2157 
2158 	return len;
2159 }
2160 
2161 void bc_num_parse(BcNum *restrict n, const char *restrict val, BcBigDig base) {
2162 
2163 	assert(n != NULL && val != NULL && base);
2164 	assert(base >= BC_NUM_MIN_BASE && base <= vm.maxes[BC_PROG_GLOBALS_IBASE]);
2165 	assert(bc_num_strValid(val));
2166 
2167 	if (!val[1]) {
2168 		BcBigDig dig = bc_num_parseChar(val[0], BC_NUM_MAX_LBASE);
2169 		bc_num_bigdig2num(n, dig);
2170 	}
2171 	else if (base == BC_BASE) bc_num_parseDecimal(n, val);
2172 	else bc_num_parseBase(n, val, base);
2173 
2174 	assert(BC_NUM_RDX_VALID(n));
2175 }
2176 
2177 void bc_num_print(BcNum *restrict n, BcBigDig base, bool newline) {
2178 
2179 	assert(n != NULL);
2180 	assert(BC_ENABLE_EXTRA_MATH || base >= BC_NUM_MIN_BASE);
2181 
2182 	bc_num_printNewline();
2183 
2184 	if (BC_NUM_ZERO(n)) bc_num_printHex(0, 1, false);
2185 	else if (base == BC_BASE) bc_num_printDecimal(n);
2186 #if BC_ENABLE_EXTRA_MATH
2187 	else if (base == 0 || base == 1) bc_num_printExponent(n, base != 0);
2188 #endif // BC_ENABLE_EXTRA_MATH
2189 	else bc_num_printBase(n, base);
2190 
2191 	if (newline) bc_num_putchar('\n');
2192 }
2193 
2194 void bc_num_bigdig2(const BcNum *restrict n, BcBigDig *result) {
2195 
2196 	// This function returns no errors because it's guaranteed to succeed if
2197 	// its preconditions are met. Those preconditions include both parameters
2198 	// being non-NULL, n being non-negative, and n being less than vm.max. If
2199 	// all of that is true, then we can just convert without worrying about
2200 	// negative errors or overflow.
2201 
2202 	BcBigDig r = 0;
2203 	size_t nrdx = BC_NUM_RDX_VAL(n);
2204 
2205 	assert(n != NULL && result != NULL);
2206 	assert(!BC_NUM_NEG(n));
2207 	assert(bc_num_cmp(n, &vm.max) < 0);
2208 	assert(n->len - nrdx <= 3);
2209 
2210 	// There is a small speed win from unrolling the loop here, and since it
2211 	// only adds 53 bytes, I decided that it was worth it.
2212 	switch (n->len - nrdx) {
2213 
2214 		case 3:
2215 		{
2216 			r = (BcBigDig) n->num[nrdx + 2];
2217 		}
2218 		// Fallthrough.
2219 		BC_FALLTHROUGH
2220 
2221 		case 2:
2222 		{
2223 			r = r * BC_BASE_POW + (BcBigDig) n->num[nrdx + 1];
2224 		}
2225 		// Fallthrough.
2226 		BC_FALLTHROUGH
2227 
2228 		case 1:
2229 		{
2230 			r = r * BC_BASE_POW + (BcBigDig) n->num[nrdx];
2231 		}
2232 	}
2233 
2234 	*result = r;
2235 }
2236 
2237 void bc_num_bigdig(const BcNum *restrict n, BcBigDig *result) {
2238 
2239 	assert(n != NULL && result != NULL);
2240 
2241 	if (BC_ERR(BC_NUM_NEG(n))) bc_vm_err(BC_ERR_MATH_NEGATIVE);
2242 	if (BC_ERR(bc_num_cmp(n, &vm.max) >= 0))
2243 		bc_vm_err(BC_ERR_MATH_OVERFLOW);
2244 
2245 	bc_num_bigdig2(n, result);
2246 }
2247 
2248 void bc_num_bigdig2num(BcNum *restrict n, BcBigDig val) {
2249 
2250 	BcDig *ptr;
2251 	size_t i;
2252 
2253 	assert(n != NULL);
2254 
2255 	bc_num_zero(n);
2256 
2257 	if (!val) return;
2258 
2259 	bc_num_expand(n, BC_NUM_BIGDIG_LOG10);
2260 
2261 	for (ptr = n->num, i = 0; val; ++i, val /= BC_BASE_POW)
2262 		ptr[i] = val % BC_BASE_POW;
2263 
2264 	n->len = i;
2265 }
2266 
2267 #if BC_ENABLE_EXTRA_MATH && BC_ENABLE_RAND
2268 void bc_num_rng(const BcNum *restrict n, BcRNG *rng) {
2269 
2270 	BcNum temp, temp2, intn, frac;
2271 	BcRand state1, state2, inc1, inc2;
2272 	size_t nrdx = BC_NUM_RDX_VAL(n);
2273 
2274 	BC_SIG_LOCK;
2275 
2276 	bc_num_init(&temp, n->len);
2277 	bc_num_init(&temp2, n->len);
2278 	bc_num_init(&frac, nrdx);
2279 	bc_num_init(&intn, bc_num_int(n));
2280 
2281 	BC_SETJMP_LOCKED(err);
2282 
2283 	BC_SIG_UNLOCK;
2284 
2285 	assert(BC_NUM_RDX_VALID_NP(vm.max));
2286 
2287 	memcpy(frac.num, n->num, BC_NUM_SIZE(nrdx));
2288 	frac.len = nrdx;
2289 	BC_NUM_RDX_SET_NP(frac, nrdx);
2290 	frac.scale = n->scale;
2291 
2292 	assert(BC_NUM_RDX_VALID_NP(frac));
2293 	assert(BC_NUM_RDX_VALID_NP(vm.max2));
2294 
2295 	bc_num_mul(&frac, &vm.max2, &temp, 0);
2296 
2297 	bc_num_truncate(&temp, temp.scale);
2298 	bc_num_copy(&frac, &temp);
2299 
2300 	memcpy(intn.num, n->num + nrdx, BC_NUM_SIZE(bc_num_int(n)));
2301 	intn.len = bc_num_int(n);
2302 
2303 	// This assert is here because it has to be true. It is also here to justify
2304 	// the use of BC_ERR_SIGNAL_ONLY() on each of the divmod's and mod's below.
2305 	assert(BC_NUM_NONZERO(&vm.max));
2306 
2307 	if (BC_NUM_NONZERO(&frac)) {
2308 
2309 		bc_num_divmod(&frac, &vm.max, &temp, &temp2, 0);
2310 
2311 		// frac is guaranteed to be smaller than vm.max * vm.max (pow).
2312 		// This means that when dividing frac by vm.max, as above, the
2313 		// quotient and remainder are both guaranteed to be less than vm.max,
2314 		// which means we can use bc_num_bigdig2() here and not worry about
2315 		// overflow.
2316 		bc_num_bigdig2(&temp2, (BcBigDig*) &state1);
2317 		bc_num_bigdig2(&temp, (BcBigDig*) &state2);
2318 	}
2319 	else state1 = state2 = 0;
2320 
2321 	if (BC_NUM_NONZERO(&intn)) {
2322 
2323 		bc_num_divmod(&intn, &vm.max, &temp, &temp2, 0);
2324 
2325 		// Because temp2 is the mod of vm.max, from above, it is guaranteed
2326 		// to be small enough to use bc_num_bigdig2().
2327 		bc_num_bigdig2(&temp2, (BcBigDig*) &inc1);
2328 
2329 		if (bc_num_cmp(&temp, &vm.max) >= 0) {
2330 			bc_num_copy(&temp2, &temp);
2331 			bc_num_mod(&temp2, &vm.max, &temp, 0);
2332 		}
2333 
2334 		// The if statement above ensures that temp is less than vm.max, which
2335 		// means that we can use bc_num_bigdig2() here.
2336 		bc_num_bigdig2(&temp, (BcBigDig*) &inc2);
2337 	}
2338 	else inc1 = inc2 = 0;
2339 
2340 	bc_rand_seed(rng, state1, state2, inc1, inc2);
2341 
2342 err:
2343 	BC_SIG_MAYLOCK;
2344 	bc_num_free(&intn);
2345 	bc_num_free(&frac);
2346 	bc_num_free(&temp2);
2347 	bc_num_free(&temp);
2348 	BC_LONGJMP_CONT;
2349 }
2350 
2351 void bc_num_createFromRNG(BcNum *restrict n, BcRNG *rng) {
2352 
2353 	BcRand s1, s2, i1, i2;
2354 	BcNum conv, temp1, temp2, temp3;
2355 	BcDig temp1_num[BC_RAND_NUM_SIZE], temp2_num[BC_RAND_NUM_SIZE];
2356 	BcDig conv_num[BC_NUM_BIGDIG_LOG10];
2357 
2358 	BC_SIG_LOCK;
2359 
2360 	bc_num_init(&temp3, 2 * BC_RAND_NUM_SIZE);
2361 
2362 	BC_SETJMP_LOCKED(err);
2363 
2364 	BC_SIG_UNLOCK;
2365 
2366 	bc_num_setup(&temp1, temp1_num, sizeof(temp1_num) / sizeof(BcDig));
2367 	bc_num_setup(&temp2, temp2_num, sizeof(temp2_num) / sizeof(BcDig));
2368 	bc_num_setup(&conv, conv_num, sizeof(conv_num) / sizeof(BcDig));
2369 
2370 	// This assert is here because it has to be true. It is also here to justify
2371 	// the assumption that vm.max2 is not zero.
2372 	assert(BC_NUM_NONZERO(&vm.max));
2373 
2374 	// Because this is true, we can just use BC_ERR_SIGNAL_ONLY() below when
2375 	// dividing by vm.max2.
2376 	assert(BC_NUM_NONZERO(&vm.max2));
2377 
2378 	bc_rand_getRands(rng, &s1, &s2, &i1, &i2);
2379 
2380 	bc_num_bigdig2num(&conv, (BcBigDig) s2);
2381 
2382 	assert(BC_NUM_RDX_VALID_NP(conv));
2383 
2384 	bc_num_mul(&conv, &vm.max, &temp1, 0);
2385 
2386 	bc_num_bigdig2num(&conv, (BcBigDig) s1);
2387 
2388 	bc_num_add(&conv, &temp1, &temp2, 0);
2389 
2390 	bc_num_div(&temp2, &vm.max2, &temp3, BC_RAND_STATE_BITS);
2391 
2392 	bc_num_bigdig2num(&conv, (BcBigDig) i2);
2393 
2394 	assert(BC_NUM_RDX_VALID_NP(conv));
2395 
2396 	bc_num_mul(&conv, &vm.max, &temp1, 0);
2397 
2398 	bc_num_bigdig2num(&conv, (BcBigDig) i1);
2399 
2400 	bc_num_add(&conv, &temp1, &temp2, 0);
2401 
2402 	bc_num_add(&temp2, &temp3, n, 0);
2403 
2404 	assert(BC_NUM_RDX_VALID(n));
2405 
2406 err:
2407 	BC_SIG_MAYLOCK;
2408 	bc_num_free(&temp3);
2409 	BC_LONGJMP_CONT;
2410 }
2411 
2412 void bc_num_irand(const BcNum *restrict a, BcNum *restrict b,
2413                   BcRNG *restrict rng)
2414 {
2415 	BcRand r;
2416 	BcBigDig modl;
2417 	BcNum pow, pow2, cp, cp2, mod, temp1, temp2, rand;
2418 	BcNum *p1, *p2, *t1, *t2, *c1, *c2, *tmp;
2419 	BcDig rand_num[BC_NUM_BIGDIG_LOG10];
2420 	bool carry;
2421 	ssize_t cmp;
2422 
2423 	assert(a != b);
2424 
2425 	if (BC_ERR(BC_NUM_NEG(a))) bc_vm_err(BC_ERR_MATH_NEGATIVE);
2426 	if (BC_ERR(BC_NUM_RDX_VAL(a))) bc_vm_err(BC_ERR_MATH_NON_INTEGER);
2427 	if (BC_NUM_ZERO(a) || BC_NUM_ONE(a)) return;
2428 
2429 	cmp = bc_num_cmp(a, &vm.max);
2430 
2431 	if (cmp <= 0) {
2432 
2433 		BcRand bits = 0;
2434 
2435 		if (cmp < 0) bc_num_bigdig2(a, (BcBigDig*) &bits);
2436 
2437 		// This condition means that bits is a power of 2. In that case, we
2438 		// can just grab a full-size int and mask out the unneeded bits.
2439 		// Also, this condition says that 0 is a power of 2, which works for
2440 		// us, since a value of 0 means a == rng->max. The bitmask will mask
2441 		// nothing in that case as well.
2442 		if (!(bits & (bits - 1))) r = bc_rand_int(rng) & (bits - 1);
2443 		else r = bc_rand_bounded(rng, bits);
2444 
2445 		// We made sure that r is less than vm.max,
2446 		// so we can use bc_num_bigdig2() here.
2447 		bc_num_bigdig2num(b, r);
2448 
2449 		return;
2450 	}
2451 
2452 	// In the case where a is less than rng->max, we have to make sure we have
2453 	// an exclusive bound. This ensures that it happens. (See below.)
2454 	carry = (cmp < 0);
2455 
2456 	BC_SIG_LOCK;
2457 
2458 	bc_num_createCopy(&cp, a);
2459 
2460 	bc_num_init(&cp2, cp.len);
2461 	bc_num_init(&mod, BC_NUM_BIGDIG_LOG10);
2462 	bc_num_init(&temp1, BC_NUM_DEF_SIZE);
2463 	bc_num_init(&temp2, BC_NUM_DEF_SIZE);
2464 	bc_num_init(&pow2, BC_NUM_DEF_SIZE);
2465 	bc_num_init(&pow, BC_NUM_DEF_SIZE);
2466 	bc_num_one(&pow);
2467 	bc_num_setup(&rand, rand_num, sizeof(rand_num) / sizeof(BcDig));
2468 
2469 	BC_SETJMP_LOCKED(err);
2470 
2471 	BC_SIG_UNLOCK;
2472 
2473 	p1 = &pow;
2474 	p2 = &pow2;
2475 	t1 = &temp1;
2476 	t2 = &temp2;
2477 	c1 = &cp;
2478 	c2 = &cp2;
2479 
2480 	// This assert is here because it has to be true. It is also here to justify
2481 	// the use of BC_ERR_SIGNAL_ONLY() on each of the divmod's and mod's below.
2482 	assert(BC_NUM_NONZERO(&vm.max));
2483 
2484 	while (BC_NUM_NONZERO(c1)) {
2485 
2486 		bc_num_divmod(c1, &vm.max, c2, &mod, 0);
2487 
2488 		// Because mod is the mod of vm.max, it is guaranteed to be smaller,
2489 		// which means we can use bc_num_bigdig2() here.
2490 		bc_num_bigdig(&mod, &modl);
2491 
2492 		if (bc_num_cmp(c1, &vm.max) < 0) {
2493 
2494 			// In this case, if there is no carry, then we know we can generate
2495 			// an integer *equal* to modl. Thus, we add one if there is no
2496 			// carry. Otherwise, we add zero, and we are still bounded properly.
2497 			// Since the last portion is guaranteed to be greater than 1, we
2498 			// know modl isn't 0 unless there is no carry.
2499 			modl += !carry;
2500 
2501 			if (modl == 1) r = 0;
2502 			else if (!modl) r = bc_rand_int(rng);
2503 			else r = bc_rand_bounded(rng, (BcRand) modl);
2504 		}
2505 		else {
2506 			if (modl) modl -= carry;
2507 			r = bc_rand_int(rng);
2508 			carry = (r >= (BcRand) modl);
2509 		}
2510 
2511 		bc_num_bigdig2num(&rand, r);
2512 
2513 		assert(BC_NUM_RDX_VALID_NP(rand));
2514 		assert(BC_NUM_RDX_VALID(p1));
2515 
2516 		bc_num_mul(&rand, p1, p2, 0);
2517 		bc_num_add(p2, t1, t2, 0);
2518 
2519 		if (BC_NUM_NONZERO(c2)) {
2520 
2521 			assert(BC_NUM_RDX_VALID_NP(vm.max));
2522 			assert(BC_NUM_RDX_VALID(p1));
2523 
2524 			bc_num_mul(&vm.max, p1, p2, 0);
2525 
2526 			tmp = p1;
2527 			p1 = p2;
2528 			p2 = tmp;
2529 
2530 			tmp = c1;
2531 			c1 = c2;
2532 			c2 = tmp;
2533 		}
2534 		else c1 = c2;
2535 
2536 		tmp = t1;
2537 		t1 = t2;
2538 		t2 = tmp;
2539 	}
2540 
2541 	bc_num_copy(b, t1);
2542 	bc_num_clean(b);
2543 
2544 	assert(BC_NUM_RDX_VALID(b));
2545 
2546 err:
2547 	BC_SIG_MAYLOCK;
2548 	bc_num_free(&pow);
2549 	bc_num_free(&pow2);
2550 	bc_num_free(&temp2);
2551 	bc_num_free(&temp1);
2552 	bc_num_free(&mod);
2553 	bc_num_free(&cp2);
2554 	bc_num_free(&cp);
2555 	BC_LONGJMP_CONT;
2556 }
2557 #endif // BC_ENABLE_EXTRA_MATH && BC_ENABLE_RAND
2558 
2559 size_t bc_num_addReq(const BcNum *a, const BcNum *b, size_t scale) {
2560 
2561 	size_t aint, bint, ardx, brdx;
2562 
2563 	BC_UNUSED(scale);
2564 
2565 	ardx = BC_NUM_RDX_VAL(a);
2566 	aint = bc_num_int(a);
2567 	assert(aint <= a->len && ardx <= a->len);
2568 
2569 	brdx = BC_NUM_RDX_VAL(b);
2570 	bint = bc_num_int(b);
2571 	assert(bint <= b->len && brdx <= b->len);
2572 
2573 	ardx = BC_MAX(ardx, brdx);
2574 	aint = BC_MAX(aint, bint);
2575 
2576 	return bc_vm_growSize(bc_vm_growSize(ardx, aint), 1);
2577 }
2578 
2579 size_t bc_num_mulReq(const BcNum *a, const BcNum *b, size_t scale) {
2580 	size_t max, rdx;
2581 	rdx = bc_vm_growSize(BC_NUM_RDX_VAL(a), BC_NUM_RDX_VAL(b));
2582 	max = BC_NUM_RDX(scale);
2583 	max = bc_vm_growSize(BC_MAX(max, rdx), 1);
2584 	rdx = bc_vm_growSize(bc_vm_growSize(bc_num_int(a), bc_num_int(b)), max);
2585 	return rdx;
2586 }
2587 
2588 size_t bc_num_divReq(const BcNum *a, const BcNum *b, size_t scale) {
2589 	size_t max, rdx;
2590 	rdx = bc_vm_growSize(BC_NUM_RDX_VAL(a), BC_NUM_RDX_VAL(b));
2591 	max = BC_NUM_RDX(scale);
2592 	max = bc_vm_growSize(BC_MAX(max, rdx), 1);
2593 	rdx = bc_vm_growSize(bc_num_int(a), max);
2594 	return rdx;
2595 }
2596 
2597 size_t bc_num_powReq(const BcNum *a, const BcNum *b, size_t scale) {
2598 	BC_UNUSED(scale);
2599 	return bc_vm_growSize(bc_vm_growSize(a->len, b->len), 1);
2600 }
2601 
2602 #if BC_ENABLE_EXTRA_MATH
2603 size_t bc_num_placesReq(const BcNum *a, const BcNum *b, size_t scale) {
2604 	BC_UNUSED(scale);
2605 	return a->len + b->len - BC_NUM_RDX_VAL(a) - BC_NUM_RDX_VAL(b);
2606 }
2607 #endif // BC_ENABLE_EXTRA_MATH
2608 
2609 void bc_num_add(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2610 	assert(BC_NUM_RDX_VALID(a));
2611 	assert(BC_NUM_RDX_VALID(b));
2612 	bc_num_binary(a, b, c, false, bc_num_as, bc_num_addReq(a, b, scale));
2613 }
2614 
2615 void bc_num_sub(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2616 	assert(BC_NUM_RDX_VALID(a));
2617 	assert(BC_NUM_RDX_VALID(b));
2618 	bc_num_binary(a, b, c, true, bc_num_as, bc_num_addReq(a, b, scale));
2619 }
2620 
2621 void bc_num_mul(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2622 	assert(BC_NUM_RDX_VALID(a));
2623 	assert(BC_NUM_RDX_VALID(b));
2624 	bc_num_binary(a, b, c, scale, bc_num_m, bc_num_mulReq(a, b, scale));
2625 }
2626 
2627 void bc_num_div(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2628 	assert(BC_NUM_RDX_VALID(a));
2629 	assert(BC_NUM_RDX_VALID(b));
2630 	bc_num_binary(a, b, c, scale, bc_num_d, bc_num_divReq(a, b, scale));
2631 }
2632 
2633 void bc_num_mod(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2634 	assert(BC_NUM_RDX_VALID(a));
2635 	assert(BC_NUM_RDX_VALID(b));
2636 	bc_num_binary(a, b, c, scale, bc_num_rem, bc_num_divReq(a, b, scale));
2637 }
2638 
2639 void bc_num_pow(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2640 	assert(BC_NUM_RDX_VALID(a));
2641 	assert(BC_NUM_RDX_VALID(b));
2642 	bc_num_binary(a, b, c, scale, bc_num_p, bc_num_powReq(a, b, scale));
2643 }
2644 
2645 #if BC_ENABLE_EXTRA_MATH
2646 void bc_num_places(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2647 	assert(BC_NUM_RDX_VALID(a));
2648 	assert(BC_NUM_RDX_VALID(b));
2649 	bc_num_binary(a, b, c, scale, bc_num_place, bc_num_placesReq(a, b, scale));
2650 }
2651 
2652 void bc_num_lshift(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2653 	assert(BC_NUM_RDX_VALID(a));
2654 	assert(BC_NUM_RDX_VALID(b));
2655 	bc_num_binary(a, b, c, scale, bc_num_left, bc_num_placesReq(a, b, scale));
2656 }
2657 
2658 void bc_num_rshift(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2659 	assert(BC_NUM_RDX_VALID(a));
2660 	assert(BC_NUM_RDX_VALID(b));
2661 	bc_num_binary(a, b, c, scale, bc_num_right, bc_num_placesReq(a, b, scale));
2662 }
2663 #endif // BC_ENABLE_EXTRA_MATH
2664 
2665 void bc_num_sqrt(BcNum *restrict a, BcNum *restrict b, size_t scale) {
2666 
2667 	BcNum num1, num2, half, f, fprime, *x0, *x1, *temp;
2668 	size_t pow, len, rdx, req, digs, digs1, digs2, resscale;
2669 	BcDig half_digs[1];
2670 
2671 	assert(a != NULL && b != NULL && a != b);
2672 
2673 	if (BC_ERR(BC_NUM_NEG(a))) bc_vm_err(BC_ERR_MATH_NEGATIVE);
2674 
2675 	if (a->scale > scale) scale = a->scale;
2676 
2677 	len = bc_vm_growSize(bc_num_intDigits(a), 1);
2678 	rdx = BC_NUM_RDX(scale);
2679 	req = bc_vm_growSize(BC_MAX(rdx, BC_NUM_RDX_VAL(a)), len >> 1);
2680 
2681 	BC_SIG_LOCK;
2682 
2683 	bc_num_init(b, bc_vm_growSize(req, 1));
2684 
2685 	BC_SIG_UNLOCK;
2686 
2687 	assert(a != NULL && b != NULL && a != b);
2688 	assert(a->num != NULL && b->num != NULL);
2689 
2690 	if (BC_NUM_ZERO(a)) {
2691 		bc_num_setToZero(b, scale);
2692 		return;
2693 	}
2694 	if (BC_NUM_ONE(a)) {
2695 		bc_num_one(b);
2696 		bc_num_extend(b, scale);
2697 		return;
2698 	}
2699 
2700 	rdx = BC_NUM_RDX(scale);
2701 	rdx = BC_MAX(rdx, BC_NUM_RDX_VAL(a));
2702 	len = bc_vm_growSize(a->len, rdx);
2703 
2704 	BC_SIG_LOCK;
2705 
2706 	bc_num_init(&num1, len);
2707 	bc_num_init(&num2, len);
2708 	bc_num_setup(&half, half_digs, sizeof(half_digs) / sizeof(BcDig));
2709 
2710 	bc_num_one(&half);
2711 	half.num[0] = BC_BASE_POW / 2;
2712 	half.len = 1;
2713 	BC_NUM_RDX_SET_NP(half, 1);
2714 	half.scale = 1;
2715 
2716 	bc_num_init(&f, len);
2717 	bc_num_init(&fprime, len);
2718 
2719 	BC_SETJMP_LOCKED(err);
2720 
2721 	BC_SIG_UNLOCK;
2722 
2723 	x0 = &num1;
2724 	x1 = &num2;
2725 
2726 	bc_num_one(x0);
2727 	pow = bc_num_intDigits(a);
2728 
2729 	if (pow) {
2730 
2731 		if (pow & 1) x0->num[0] = 2;
2732 		else x0->num[0] = 6;
2733 
2734 		pow -= 2 - (pow & 1);
2735 		bc_num_shiftLeft(x0, pow / 2);
2736 	}
2737 
2738 	// I can set the rdx here directly because neg should be false.
2739 	x0->scale = x0->rdx = digs = digs1 = digs2 = 0;
2740 	resscale = (scale + BC_BASE_DIGS) + 2;
2741 
2742 	while (bc_num_cmp(x1, x0)) {
2743 
2744 		assert(BC_NUM_NONZERO(x0));
2745 
2746 		bc_num_div(a, x0, &f, resscale);
2747 		bc_num_add(x0, &f, &fprime, resscale);
2748 
2749 		assert(BC_NUM_RDX_VALID_NP(fprime));
2750 		assert(BC_NUM_RDX_VALID_NP(half));
2751 
2752 		bc_num_mul(&fprime, &half, x1, resscale);
2753 
2754 		temp = x0;
2755 		x0 = x1;
2756 		x1 = temp;
2757 	}
2758 
2759 	bc_num_copy(b, x0);
2760 	if (b->scale > scale) bc_num_truncate(b, b->scale - scale);
2761 
2762 	assert(!BC_NUM_NEG(b) || BC_NUM_NONZERO(b));
2763 	assert(BC_NUM_RDX_VALID(b));
2764 	assert(BC_NUM_RDX_VAL(b) <= b->len || !b->len);
2765 	assert(!b->len || b->num[b->len - 1] || BC_NUM_RDX_VAL(b) == b->len);
2766 
2767 err:
2768 	BC_SIG_MAYLOCK;
2769 	bc_num_free(&fprime);
2770 	bc_num_free(&f);
2771 	bc_num_free(&num2);
2772 	bc_num_free(&num1);
2773 	BC_LONGJMP_CONT;
2774 }
2775 
2776 void bc_num_divmod(BcNum *a, BcNum *b, BcNum *c, BcNum *d, size_t scale) {
2777 
2778 	size_t ts, len;
2779 	BcNum *ptr_a, num2;
2780 	bool init = false;
2781 
2782 	ts = BC_MAX(scale + b->scale, a->scale);
2783 	len = bc_num_mulReq(a, b, ts);
2784 
2785 	assert(a != NULL && b != NULL && c != NULL && d != NULL);
2786 	assert(c != d && a != d && b != d && b != c);
2787 
2788 	if (c == a) {
2789 
2790 		memcpy(&num2, c, sizeof(BcNum));
2791 		ptr_a = &num2;
2792 
2793 		BC_SIG_LOCK;
2794 
2795 		bc_num_init(c, len);
2796 
2797 		init = true;
2798 
2799 		BC_SETJMP_LOCKED(err);
2800 
2801 		BC_SIG_UNLOCK;
2802 	}
2803 	else {
2804 		ptr_a = a;
2805 		bc_num_expand(c, len);
2806 	}
2807 
2808 	if (BC_NUM_NONZERO(a) && !BC_NUM_RDX_VAL(a) &&
2809 	    !BC_NUM_RDX_VAL(b) && b->len == 1 && !scale)
2810 	{
2811 		BcBigDig rem;
2812 
2813 		bc_num_divArray(ptr_a, (BcBigDig) b->num[0], c, &rem);
2814 
2815 		assert(rem < BC_BASE_POW);
2816 
2817 		d->num[0] = (BcDig) rem;
2818 		d->len = (rem != 0);
2819 	}
2820 	else bc_num_r(ptr_a, b, c, d, scale, ts);
2821 
2822 	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
2823 	assert(BC_NUM_RDX_VALID(c));
2824 	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
2825 	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
2826 	assert(!BC_NUM_NEG(d) || BC_NUM_NONZERO(d));
2827 	assert(BC_NUM_RDX_VALID(d));
2828 	assert(BC_NUM_RDX_VAL(d) <= d->len || !d->len);
2829 	assert(!d->len || d->num[d->len - 1] || BC_NUM_RDX_VAL(d) == d->len);
2830 
2831 err:
2832 	if (init) {
2833 		BC_SIG_MAYLOCK;
2834 		bc_num_free(&num2);
2835 		BC_LONGJMP_CONT;
2836 	}
2837 }
2838 
2839 #if DC_ENABLED
2840 void bc_num_modexp(BcNum *a, BcNum *b, BcNum *c, BcNum *restrict d) {
2841 
2842 	BcNum base, exp, two, temp;
2843 	BcDig two_digs[2];
2844 
2845 	assert(a != NULL && b != NULL && c != NULL && d != NULL);
2846 	assert(a != d && b != d && c != d);
2847 
2848 	if (BC_ERR(BC_NUM_ZERO(c))) bc_vm_err(BC_ERR_MATH_DIVIDE_BY_ZERO);
2849 	if (BC_ERR(BC_NUM_NEG(b))) bc_vm_err(BC_ERR_MATH_NEGATIVE);
2850 	if (BC_ERR(BC_NUM_RDX_VAL(a) || BC_NUM_RDX_VAL(b) || BC_NUM_RDX_VAL(c)))
2851 		bc_vm_err(BC_ERR_MATH_NON_INTEGER);
2852 
2853 	bc_num_expand(d, c->len);
2854 
2855 	BC_SIG_LOCK;
2856 
2857 	bc_num_init(&base, c->len);
2858 	bc_num_setup(&two, two_digs, sizeof(two_digs) / sizeof(BcDig));
2859 	bc_num_init(&temp, b->len + 1);
2860 	bc_num_createCopy(&exp, b);
2861 
2862 	BC_SETJMP_LOCKED(err);
2863 
2864 	BC_SIG_UNLOCK;
2865 
2866 	bc_num_one(&two);
2867 	two.num[0] = 2;
2868 	bc_num_one(d);
2869 
2870 	// We already checked for 0.
2871 	bc_num_rem(a, c, &base, 0);
2872 
2873 	while (BC_NUM_NONZERO(&exp)) {
2874 
2875 		// Num two cannot be 0, so no errors.
2876 		bc_num_divmod(&exp, &two, &exp, &temp, 0);
2877 
2878 		if (BC_NUM_ONE(&temp) && !BC_NUM_NEG_NP(temp)) {
2879 
2880 			assert(BC_NUM_RDX_VALID(d));
2881 			assert(BC_NUM_RDX_VALID_NP(base));
2882 
2883 			bc_num_mul(d, &base, &temp, 0);
2884 
2885 			// We already checked for 0.
2886 			bc_num_rem(&temp, c, d, 0);
2887 		}
2888 
2889 		assert(BC_NUM_RDX_VALID_NP(base));
2890 
2891 		bc_num_mul(&base, &base, &temp, 0);
2892 
2893 		// We already checked for 0.
2894 		bc_num_rem(&temp, c, &base, 0);
2895 	}
2896 
2897 err:
2898 	BC_SIG_MAYLOCK;
2899 	bc_num_free(&exp);
2900 	bc_num_free(&temp);
2901 	bc_num_free(&base);
2902 	BC_LONGJMP_CONT;
2903 	assert(!BC_NUM_NEG(d) || d->len);
2904 	assert(BC_NUM_RDX_VALID(d));
2905 	assert(!d->len || d->num[d->len - 1] || BC_NUM_RDX_VAL(d) == d->len);
2906 }
2907 #endif // DC_ENABLED
2908 
2909 #if BC_DEBUG_CODE
2910 void bc_num_printDebug(const BcNum *n, const char *name, bool emptyline) {
2911 	bc_file_puts(&vm.fout, name);
2912 	bc_file_puts(&vm.fout, ": ");
2913 	bc_num_printDecimal(n);
2914 	bc_file_putchar(&vm.fout, '\n');
2915 	if (emptyline) bc_file_putchar(&vm.fout, '\n');
2916 	vm.nchars = 0;
2917 }
2918 
2919 void bc_num_printDigs(const BcDig *n, size_t len, bool emptyline) {
2920 
2921 	size_t i;
2922 
2923 	for (i = len - 1; i < len; --i)
2924 		bc_file_printf(&vm.fout, " %lu", (unsigned long) n[i]);
2925 
2926 	bc_file_putchar(&vm.fout, '\n');
2927 	if (emptyline) bc_file_putchar(&vm.fout, '\n');
2928 	vm.nchars = 0;
2929 }
2930 
2931 void bc_num_printWithDigs(const BcNum *n, const char *name, bool emptyline) {
2932 	bc_file_puts(&vm.fout, name);
2933 	bc_file_printf(&vm.fout, " len: %zu, rdx: %zu, scale: %zu\n",
2934 	               name, n->len, BC_NUM_RDX_VAL(n), n->scale);
2935 	bc_num_printDigs(n->num, n->len, emptyline);
2936 }
2937 
2938 void bc_num_dump(const char *varname, const BcNum *n) {
2939 
2940 	ulong i, scale = n->scale;
2941 
2942 	bc_file_printf(&vm.ferr, "\n%s = %s", varname,
2943 	               n->len ? (BC_NUM_NEG(n) ? "-" : "+") : "0 ");
2944 
2945 	for (i = n->len - 1; i < n->len; --i) {
2946 
2947 		if (i + 1 == BC_NUM_RDX_VAL(n)) bc_file_puts(&vm.ferr, ". ");
2948 
2949 		if (scale / BC_BASE_DIGS != BC_NUM_RDX_VAL(n) - i - 1)
2950 			bc_file_printf(&vm.ferr, "%lu ", (unsigned long) n->num[i]);
2951 		else {
2952 
2953 			int mod = scale % BC_BASE_DIGS;
2954 			int d = BC_BASE_DIGS - mod;
2955 			BcDig div;
2956 
2957 			if (mod != 0) {
2958 				div = n->num[i] / ((BcDig) bc_num_pow10[(ulong) d]);
2959 				bc_file_printf(&vm.ferr, "%lu", (unsigned long) div);
2960 			}
2961 
2962 			div = n->num[i] % ((BcDig) bc_num_pow10[(ulong) d]);
2963 			bc_file_printf(&vm.ferr, " ' %lu ", (unsigned long) div);
2964 		}
2965 	}
2966 
2967 	bc_file_printf(&vm.ferr, "(%zu | %zu.%zu / %zu) %lu\n",
2968 	               n->scale, n->len, BC_NUM_RDX_VAL(n), n->cap,
2969 	               (unsigned long) (void*) n->num);
2970 }
2971 #endif // BC_DEBUG_CODE
2972