xref: /freebsd/contrib/bc/src/num.c (revision 78bc019d220e05abb5b12f678f9b4a847019bbcc)
1 /*
2  * *****************************************************************************
3  *
4  * SPDX-License-Identifier: BSD-2-Clause
5  *
6  * Copyright (c) 2018-2021 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 // Before you try to understand this code, see the development manual
49 // (manuals/development.md#numbers).
50 
51 static void
52 bc_num_m(BcNum* a, BcNum* b, BcNum* restrict c, size_t scale);
53 
54 /**
55  * Multiply two numbers and throw a math error if they overflow.
56  * @param a  The first operand.
57  * @param b  The second operand.
58  * @return   The product of the two operands.
59  */
60 static inline size_t
61 bc_num_mulOverflow(size_t a, size_t b)
62 {
63 	size_t res = a * b;
64 	if (BC_ERR(BC_VM_MUL_OVERFLOW(a, b, res))) bc_err(BC_ERR_MATH_OVERFLOW);
65 	return res;
66 }
67 
68 /**
69  * Conditionally negate @a n based on @a neg. Algorithm taken from
70  * https://graphics.stanford.edu/~seander/bithacks.html#ConditionalNegate .
71  * @param n    The value to turn into a signed value and negate.
72  * @param neg  The condition to negate or not.
73  */
74 static inline ssize_t
75 bc_num_neg(size_t n, bool neg)
76 {
77 	return (((ssize_t) n) ^ -((ssize_t) neg)) + neg;
78 }
79 
80 /**
81  * Compare a BcNum against zero.
82  * @param n  The number to compare.
83  * @return   -1 if the number is less than 0, 1 if greater, and 0 if equal.
84  */
85 ssize_t
86 bc_num_cmpZero(const BcNum* n)
87 {
88 	return bc_num_neg((n)->len != 0, BC_NUM_NEG(n));
89 }
90 
91 /**
92  * Return the number of integer limbs in a BcNum. This is the opposite of rdx.
93  * @param n  The number to return the amount of integer limbs for.
94  * @return   The amount of integer limbs in @a n.
95  */
96 static inline size_t
97 bc_num_int(const BcNum* n)
98 {
99 	return n->len ? n->len - BC_NUM_RDX_VAL(n) : 0;
100 }
101 
102 /**
103  * Expand a number's allocation capacity to at least req limbs.
104  * @param n    The number to expand.
105  * @param req  The number limbs to expand the allocation capacity to.
106  */
107 static void
108 bc_num_expand(BcNum* restrict n, size_t req)
109 {
110 	assert(n != NULL);
111 
112 	req = req >= BC_NUM_DEF_SIZE ? req : BC_NUM_DEF_SIZE;
113 
114 	if (req > n->cap)
115 	{
116 		BC_SIG_LOCK;
117 
118 		n->num = bc_vm_realloc(n->num, BC_NUM_SIZE(req));
119 		n->cap = req;
120 
121 		BC_SIG_UNLOCK;
122 	}
123 }
124 
125 /**
126  * Set a number to 0 with the specified scale.
127  * @param n      The number to set to zero.
128  * @param scale  The scale to set the number to.
129  */
130 static void
131 bc_num_setToZero(BcNum* restrict n, size_t scale)
132 {
133 	assert(n != NULL);
134 	n->scale = scale;
135 	n->len = n->rdx = 0;
136 }
137 
138 void
139 bc_num_zero(BcNum* restrict n)
140 {
141 	bc_num_setToZero(n, 0);
142 }
143 
144 void
145 bc_num_one(BcNum* restrict n)
146 {
147 	bc_num_zero(n);
148 	n->len = 1;
149 	n->num[0] = 1;
150 }
151 
152 /**
153  * "Cleans" a number, which means reducing the length if the most significant
154  * limbs are zero.
155  * @param n  The number to clean.
156  */
157 static void
158 bc_num_clean(BcNum* restrict n)
159 {
160 	// Reduce the length.
161 	while (BC_NUM_NONZERO(n) && !n->num[n->len - 1])
162 	{
163 		n->len -= 1;
164 	}
165 
166 	// Special cases.
167 	if (BC_NUM_ZERO(n)) n->rdx = 0;
168 	else
169 	{
170 		// len must be at least as much as rdx.
171 		size_t rdx = BC_NUM_RDX_VAL(n);
172 		if (n->len < rdx) n->len = rdx;
173 	}
174 }
175 
176 /**
177  * Returns the log base 10 of @a i. I could have done this with floating-point
178  * math, and in fact, I originally did. However, that was the only
179  * floating-point code in the entire codebase, and I decided I didn't want any.
180  * This is fast enough. Also, it might handle larger numbers better.
181  * @param i  The number to return the log base 10 of.
182  * @return   The log base 10 of @a i.
183  */
184 static size_t
185 bc_num_log10(size_t i)
186 {
187 	size_t len;
188 
189 	for (len = 1; i; i /= BC_BASE, ++len)
190 	{
191 		continue;
192 	}
193 
194 	assert(len - 1 <= BC_BASE_DIGS + 1);
195 
196 	return len - 1;
197 }
198 
199 /**
200  * Returns the number of decimal digits in a limb that are zero starting at the
201  * most significant digits. This basically returns how much of the limb is used.
202  * @param n  The number.
203  * @return   The number of decimal digits that are 0 starting at the most
204  *           significant digits.
205  */
206 static inline size_t
207 bc_num_zeroDigits(const BcDig* n)
208 {
209 	assert(*n >= 0);
210 	assert(((size_t) *n) < BC_BASE_POW);
211 	return BC_BASE_DIGS - bc_num_log10((size_t) *n);
212 }
213 
214 /**
215  * Return the total number of integer digits in a number. This is the opposite
216  * of scale, like bc_num_int() is the opposite of rdx.
217  * @param n  The number.
218  * @return   The number of integer digits in @a n.
219  */
220 static size_t
221 bc_num_intDigits(const BcNum* n)
222 {
223 	size_t digits = bc_num_int(n) * BC_BASE_DIGS;
224 	if (digits > 0) digits -= bc_num_zeroDigits(n->num + n->len - 1);
225 	return digits;
226 }
227 
228 /**
229  * Returns the number of limbs of a number that are non-zero starting at the
230  * most significant limbs. This expects that there are *no* integer limbs in the
231  * number because it is specifically to figure out how many zero limbs after the
232  * decimal place to ignore. If there are zero limbs after non-zero limbs, they
233  * are counted as non-zero limbs.
234  * @param n  The number.
235  * @return   The number of non-zero limbs after the decimal point.
236  */
237 static size_t
238 bc_num_nonZeroLen(const BcNum* restrict n)
239 {
240 	size_t i, len = n->len;
241 
242 	assert(len == BC_NUM_RDX_VAL(n));
243 
244 	for (i = len - 1; i < len && !n->num[i]; --i)
245 	{
246 		continue;
247 	}
248 
249 	assert(i + 1 > 0);
250 
251 	return i + 1;
252 }
253 
254 /**
255  * Performs a one-limb add with a carry.
256  * @param a      The first limb.
257  * @param b      The second limb.
258  * @param carry  An in/out parameter; the carry in from the previous add and the
259  *               carry out from this add.
260  * @return       The resulting limb sum.
261  */
262 static BcDig
263 bc_num_addDigits(BcDig a, BcDig b, bool* carry)
264 {
265 	assert(((BcBigDig) BC_BASE_POW) * 2 == ((BcDig) BC_BASE_POW) * 2);
266 	assert(a < BC_BASE_POW);
267 	assert(b < BC_BASE_POW);
268 
269 	a += b + *carry;
270 	*carry = (a >= BC_BASE_POW);
271 	if (*carry) a -= BC_BASE_POW;
272 
273 	assert(a >= 0);
274 	assert(a < BC_BASE_POW);
275 
276 	return a;
277 }
278 
279 /**
280  * Performs a one-limb subtract with a carry.
281  * @param a      The first limb.
282  * @param b      The second limb.
283  * @param carry  An in/out parameter; the carry in from the previous subtract
284  *               and the carry out from this subtract.
285  * @return       The resulting limb difference.
286  */
287 static BcDig
288 bc_num_subDigits(BcDig a, BcDig b, bool* carry)
289 {
290 	assert(a < BC_BASE_POW);
291 	assert(b < BC_BASE_POW);
292 
293 	b += *carry;
294 	*carry = (a < b);
295 	if (*carry) a += BC_BASE_POW;
296 
297 	assert(a - b >= 0);
298 	assert(a - b < BC_BASE_POW);
299 
300 	return a - b;
301 }
302 
303 /**
304  * Add two BcDig arrays and store the result in the first array.
305  * @param a    The first operand and out array.
306  * @param b    The second operand.
307  * @param len  The length of @a b.
308  */
309 static void
310 bc_num_addArrays(BcDig* restrict a, const BcDig* restrict b, size_t len)
311 {
312 	size_t i;
313 	bool carry = false;
314 
315 	for (i = 0; i < len; ++i)
316 	{
317 		a[i] = bc_num_addDigits(a[i], b[i], &carry);
318 	}
319 
320 	// Take care of the extra limbs in the bigger array.
321 	for (; carry; ++i)
322 	{
323 		a[i] = bc_num_addDigits(a[i], 0, &carry);
324 	}
325 }
326 
327 /**
328  * Subtract two BcDig arrays and store the result in the first array.
329  * @param a    The first operand and out array.
330  * @param b    The second operand.
331  * @param len  The length of @a b.
332  */
333 static void
334 bc_num_subArrays(BcDig* restrict a, const BcDig* restrict b, size_t len)
335 {
336 	size_t i;
337 	bool carry = false;
338 
339 	for (i = 0; i < len; ++i)
340 	{
341 		a[i] = bc_num_subDigits(a[i], b[i], &carry);
342 	}
343 
344 	// Take care of the extra limbs in the bigger array.
345 	for (; carry; ++i)
346 	{
347 		a[i] = bc_num_subDigits(a[i], 0, &carry);
348 	}
349 }
350 
351 /**
352  * Multiply a BcNum array by a one-limb number. This is a faster version of
353  * multiplication for when we can use it.
354  * @param a  The BcNum to multiply by the one-limb number.
355  * @param b  The one limb of the one-limb number.
356  * @param c  The return parameter.
357  */
358 static void
359 bc_num_mulArray(const BcNum* restrict a, BcBigDig b, BcNum* restrict c)
360 {
361 	size_t i;
362 	BcBigDig carry = 0;
363 
364 	assert(b <= BC_BASE_POW);
365 
366 	// Make sure the return parameter is big enough.
367 	if (a->len + 1 > c->cap) bc_num_expand(c, a->len + 1);
368 
369 	// We want the entire return parameter to be zero for cleaning later.
370 	// NOLINTNEXTLINE
371 	memset(c->num, 0, BC_NUM_SIZE(c->cap));
372 
373 	// Actual multiplication loop.
374 	for (i = 0; i < a->len; ++i)
375 	{
376 		BcBigDig in = ((BcBigDig) a->num[i]) * b + carry;
377 		c->num[i] = in % BC_BASE_POW;
378 		carry = in / BC_BASE_POW;
379 	}
380 
381 	assert(carry < BC_BASE_POW);
382 
383 	// Finishing touches.
384 	c->num[i] = (BcDig) carry;
385 	c->len = a->len;
386 	c->len += (carry != 0);
387 
388 	bc_num_clean(c);
389 
390 	// Postconditions.
391 	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
392 	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
393 	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
394 }
395 
396 /**
397  * Divide a BcNum array by a one-limb number. This is a faster version of divide
398  * for when we can use it.
399  * @param a    The BcNum to multiply by the one-limb number.
400  * @param b    The one limb of the one-limb number.
401  * @param c    The return parameter for the quotient.
402  * @param rem  The return parameter for the remainder.
403  */
404 static void
405 bc_num_divArray(const BcNum* restrict a, BcBigDig b, BcNum* restrict c,
406                 BcBigDig* rem)
407 {
408 	size_t i;
409 	BcBigDig carry = 0;
410 
411 	assert(c->cap >= a->len);
412 
413 	// Actual division loop.
414 	for (i = a->len - 1; i < a->len; --i)
415 	{
416 		BcBigDig in = ((BcBigDig) a->num[i]) + carry * BC_BASE_POW;
417 		assert(in / b < BC_BASE_POW);
418 		c->num[i] = (BcDig) (in / b);
419 		carry = in % b;
420 	}
421 
422 	// Finishing touches.
423 	c->len = a->len;
424 	bc_num_clean(c);
425 	*rem = carry;
426 
427 	// Postconditions.
428 	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
429 	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
430 	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
431 }
432 
433 /**
434  * Compare two BcDig arrays and return >0 if @a b is greater, <0 if @a b is
435  * less, and 0 if equal. Both @a a and @a b must have the same length.
436  * @param a    The first array.
437  * @param b    The second array.
438  * @param len  The minimum length of the arrays.
439  */
440 static ssize_t
441 bc_num_compare(const BcDig* restrict a, const BcDig* restrict b, size_t len)
442 {
443 	size_t i;
444 	BcDig c = 0;
445 	for (i = len - 1; i < len && !(c = a[i] - b[i]); --i)
446 	{
447 		continue;
448 	}
449 	return bc_num_neg(i + 1, c < 0);
450 }
451 
452 ssize_t
453 bc_num_cmp(const BcNum* a, const BcNum* b)
454 {
455 	size_t i, min, a_int, b_int, diff, ardx, brdx;
456 	BcDig* max_num;
457 	BcDig* min_num;
458 	bool a_max, neg = false;
459 	ssize_t cmp;
460 
461 	assert(a != NULL && b != NULL);
462 
463 	// Same num? Equal.
464 	if (a == b) return 0;
465 
466 	// Easy cases.
467 	if (BC_NUM_ZERO(a)) return bc_num_neg(b->len != 0, !BC_NUM_NEG(b));
468 	if (BC_NUM_ZERO(b)) return bc_num_cmpZero(a);
469 	if (BC_NUM_NEG(a))
470 	{
471 		if (BC_NUM_NEG(b)) neg = true;
472 		else return -1;
473 	}
474 	else if (BC_NUM_NEG(b)) return 1;
475 
476 	// Get the number of int limbs in each number and get the difference.
477 	a_int = bc_num_int(a);
478 	b_int = bc_num_int(b);
479 	a_int -= b_int;
480 
481 	// If there's a difference, then just return the comparison.
482 	if (a_int) return neg ? -((ssize_t) a_int) : (ssize_t) a_int;
483 
484 	// Get the rdx's and figure out the max.
485 	ardx = BC_NUM_RDX_VAL(a);
486 	brdx = BC_NUM_RDX_VAL(b);
487 	a_max = (ardx > brdx);
488 
489 	// Set variables based on the above.
490 	if (a_max)
491 	{
492 		min = brdx;
493 		diff = ardx - brdx;
494 		max_num = a->num + diff;
495 		min_num = b->num;
496 	}
497 	else
498 	{
499 		min = ardx;
500 		diff = brdx - ardx;
501 		max_num = b->num + diff;
502 		min_num = a->num;
503 	}
504 
505 	// Do a full limb-by-limb comparison.
506 	cmp = bc_num_compare(max_num, min_num, b_int + min);
507 
508 	// If we found a difference, return it based on state.
509 	if (cmp) return bc_num_neg((size_t) cmp, !a_max == !neg);
510 
511 	// If there was no difference, then the final step is to check which number
512 	// has greater or lesser limbs beyond the other's.
513 	for (max_num -= diff, i = diff - 1; i < diff; --i)
514 	{
515 		if (max_num[i]) return bc_num_neg(1, !a_max == !neg);
516 	}
517 
518 	return 0;
519 }
520 
521 void
522 bc_num_truncate(BcNum* restrict n, size_t places)
523 {
524 	size_t nrdx, places_rdx;
525 
526 	if (!places) return;
527 
528 	// Grab these needed values; places_rdx is the rdx equivalent to places like
529 	// rdx is to scale.
530 	nrdx = BC_NUM_RDX_VAL(n);
531 	places_rdx = nrdx ? nrdx - BC_NUM_RDX(n->scale - places) : 0;
532 
533 	// We cannot truncate more places than we have.
534 	assert(places <= n->scale && (BC_NUM_ZERO(n) || places_rdx <= n->len));
535 
536 	n->scale -= places;
537 	BC_NUM_RDX_SET(n, nrdx - places_rdx);
538 
539 	// Only when the number is nonzero do we need to do the hard stuff.
540 	if (BC_NUM_NONZERO(n))
541 	{
542 		size_t pow;
543 
544 		// This calculates how many decimal digits are in the least significant
545 		// limb.
546 		pow = n->scale % BC_BASE_DIGS;
547 		pow = pow ? BC_BASE_DIGS - pow : 0;
548 		pow = bc_num_pow10[pow];
549 
550 		n->len -= places_rdx;
551 
552 		// We have to move limbs to maintain invariants. The limbs must begin at
553 		// the beginning of the BcNum array.
554 		// NOLINTNEXTLINE
555 		memmove(n->num, n->num + places_rdx, BC_NUM_SIZE(n->len));
556 
557 		// Clear the lower part of the last digit.
558 		if (BC_NUM_NONZERO(n)) n->num[0] -= n->num[0] % (BcDig) pow;
559 
560 		bc_num_clean(n);
561 	}
562 }
563 
564 void
565 bc_num_extend(BcNum* restrict n, size_t places)
566 {
567 	size_t nrdx, places_rdx;
568 
569 	if (!places) return;
570 
571 	// Easy case with zero; set the scale.
572 	if (BC_NUM_ZERO(n))
573 	{
574 		n->scale += places;
575 		return;
576 	}
577 
578 	// Grab these needed values; places_rdx is the rdx equivalent to places like
579 	// rdx is to scale.
580 	nrdx = BC_NUM_RDX_VAL(n);
581 	places_rdx = BC_NUM_RDX(places + n->scale) - nrdx;
582 
583 	// This is the hard case. We need to expand the number, move the limbs, and
584 	// set the limbs that were just cleared.
585 	if (places_rdx)
586 	{
587 		bc_num_expand(n, bc_vm_growSize(n->len, places_rdx));
588 		// NOLINTNEXTLINE
589 		memmove(n->num + places_rdx, n->num, BC_NUM_SIZE(n->len));
590 		// NOLINTNEXTLINE
591 		memset(n->num, 0, BC_NUM_SIZE(places_rdx));
592 	}
593 
594 	// Finally, set scale and rdx.
595 	BC_NUM_RDX_SET(n, nrdx + places_rdx);
596 	n->scale += places;
597 	n->len += places_rdx;
598 
599 	assert(BC_NUM_RDX_VAL(n) == BC_NUM_RDX(n->scale));
600 }
601 
602 /**
603  * Retires (finishes) a multiplication or division operation.
604  */
605 static void
606 bc_num_retireMul(BcNum* restrict n, size_t scale, bool neg1, bool neg2)
607 {
608 	// Make sure scale is correct.
609 	if (n->scale < scale) bc_num_extend(n, scale - n->scale);
610 	else bc_num_truncate(n, n->scale - scale);
611 
612 	bc_num_clean(n);
613 
614 	// We need to ensure rdx is correct.
615 	if (BC_NUM_NONZERO(n)) n->rdx = BC_NUM_NEG_VAL(n, !neg1 != !neg2);
616 }
617 
618 /**
619  * Splits a number into two BcNum's. This is used in Karatsuba.
620  * @param n    The number to split.
621  * @param idx  The index to split at.
622  * @param a    An out parameter; the low part of @a n.
623  * @param b    An out parameter; the high part of @a n.
624  */
625 static void
626 bc_num_split(const BcNum* restrict n, size_t idx, BcNum* restrict a,
627              BcNum* restrict b)
628 {
629 	// We want a and b to be clear.
630 	assert(BC_NUM_ZERO(a));
631 	assert(BC_NUM_ZERO(b));
632 
633 	// The usual case.
634 	if (idx < n->len)
635 	{
636 		// Set the fields first.
637 		b->len = n->len - idx;
638 		a->len = idx;
639 		a->scale = b->scale = 0;
640 		BC_NUM_RDX_SET(a, 0);
641 		BC_NUM_RDX_SET(b, 0);
642 
643 		assert(a->cap >= a->len);
644 		assert(b->cap >= b->len);
645 
646 		// Copy the arrays. This is not necessary for safety, but it is faster,
647 		// for some reason.
648 		// NOLINTNEXTLINE
649 		memcpy(b->num, n->num + idx, BC_NUM_SIZE(b->len));
650 		// NOLINTNEXTLINE
651 		memcpy(a->num, n->num, BC_NUM_SIZE(idx));
652 
653 		bc_num_clean(b);
654 	}
655 	// If the index is weird, just skip the split.
656 	else bc_num_copy(a, n);
657 
658 	bc_num_clean(a);
659 }
660 
661 /**
662  * Copies a number into another, but shifts the rdx so that the result number
663  * only sees the integer part of the source number.
664  * @param n  The number to copy.
665  * @param r  The result number with a shifted rdx, len, and num.
666  */
667 static void
668 bc_num_shiftRdx(const BcNum* restrict n, BcNum* restrict r)
669 {
670 	size_t rdx = BC_NUM_RDX_VAL(n);
671 
672 	r->len = n->len - rdx;
673 	r->cap = n->cap - rdx;
674 	r->num = n->num + rdx;
675 
676 	BC_NUM_RDX_SET_NEG(r, 0, BC_NUM_NEG(n));
677 	r->scale = 0;
678 }
679 
680 /**
681  * Shifts a number so that all of the least significant limbs of the number are
682  * skipped. This must be undone by bc_num_unshiftZero().
683  * @param n  The number to shift.
684  */
685 static size_t
686 bc_num_shiftZero(BcNum* restrict n)
687 {
688 	size_t i;
689 
690 	// If we don't have an integer, that is a problem, but it's also a bug
691 	// because the caller should have set everything up right.
692 	assert(!BC_NUM_RDX_VAL(n) || BC_NUM_ZERO(n));
693 
694 	for (i = 0; i < n->len && !n->num[i]; ++i)
695 	{
696 		continue;
697 	}
698 
699 	n->len -= i;
700 	n->num += i;
701 
702 	return i;
703 }
704 
705 /**
706  * Undo the damage done by bc_num_unshiftZero(). This must be called like all
707  * cleanup functions: after a label used by setjmp() and longjmp().
708  * @param n           The number to unshift.
709  * @param places_rdx  The amount the number was originally shift.
710  */
711 static void
712 bc_num_unshiftZero(BcNum* restrict n, size_t places_rdx)
713 {
714 	n->len += places_rdx;
715 	n->num -= places_rdx;
716 }
717 
718 /**
719  * Shifts the digits right within a number by no more than BC_BASE_DIGS - 1.
720  * This is the final step on shifting numbers with the shift operators. It
721  * depends on the caller to shift the limbs properly because all it does is
722  * ensure that the rdx point is realigned to be between limbs.
723  * @param n    The number to shift digits in.
724  * @param dig  The number of places to shift right.
725  */
726 static void
727 bc_num_shift(BcNum* restrict n, BcBigDig dig)
728 {
729 	size_t i, len = n->len;
730 	BcBigDig carry = 0, pow;
731 	BcDig* ptr = n->num;
732 
733 	assert(dig < BC_BASE_DIGS);
734 
735 	// Figure out the parameters for division.
736 	pow = bc_num_pow10[dig];
737 	dig = bc_num_pow10[BC_BASE_DIGS - dig];
738 
739 	// Run a series of divisions and mods with carries across the entire number
740 	// array. This effectively shifts everything over.
741 	for (i = len - 1; i < len; --i)
742 	{
743 		BcBigDig in, temp;
744 		in = ((BcBigDig) ptr[i]);
745 		temp = carry * dig;
746 		carry = in % pow;
747 		ptr[i] = ((BcDig) (in / pow)) + (BcDig) temp;
748 	}
749 
750 	assert(!carry);
751 }
752 
753 /**
754  * Shift a number left by a certain number of places. This is the workhorse of
755  * the left shift operator.
756  * @param n       The number to shift left.
757  * @param places  The amount of places to shift @a n left by.
758  */
759 static void
760 bc_num_shiftLeft(BcNum* restrict n, size_t places)
761 {
762 	BcBigDig dig;
763 	size_t places_rdx;
764 	bool shift;
765 
766 	if (!places) return;
767 
768 	// Make sure to grow the number if necessary.
769 	if (places > n->scale)
770 	{
771 		size_t size = bc_vm_growSize(BC_NUM_RDX(places - n->scale), n->len);
772 		if (size > SIZE_MAX - 1) bc_err(BC_ERR_MATH_OVERFLOW);
773 	}
774 
775 	// If zero, we can just set the scale and bail.
776 	if (BC_NUM_ZERO(n))
777 	{
778 		if (n->scale >= places) n->scale -= places;
779 		else n->scale = 0;
780 		return;
781 	}
782 
783 	// When I changed bc to have multiple digits per limb, this was the hardest
784 	// code to change. This and shift right. Make sure you understand this
785 	// before attempting anything.
786 
787 	// This is how many limbs we will shift.
788 	dig = (BcBigDig) (places % BC_BASE_DIGS);
789 	shift = (dig != 0);
790 
791 	// Convert places to a rdx value.
792 	places_rdx = BC_NUM_RDX(places);
793 
794 	// If the number is not an integer, we need special care. The reason an
795 	// integer doesn't is because left shift would only extend the integer,
796 	// whereas a non-integer might have its fractional part eliminated or only
797 	// partially eliminated.
798 	if (n->scale)
799 	{
800 		size_t nrdx = BC_NUM_RDX_VAL(n);
801 
802 		// If the number's rdx is bigger, that's the hard case.
803 		if (nrdx >= places_rdx)
804 		{
805 			size_t mod = n->scale % BC_BASE_DIGS, revdig;
806 
807 			// We want mod to be in the range [1, BC_BASE_DIGS], not
808 			// [0, BC_BASE_DIGS).
809 			mod = mod ? mod : BC_BASE_DIGS;
810 
811 			// We need to reverse dig to get the actual number of digits.
812 			revdig = dig ? BC_BASE_DIGS - dig : 0;
813 
814 			// If the two overflow BC_BASE_DIGS, we need to move an extra place.
815 			if (mod + revdig > BC_BASE_DIGS) places_rdx = 1;
816 			else places_rdx = 0;
817 		}
818 		else places_rdx -= nrdx;
819 	}
820 
821 	// If this is non-zero, we need an extra place, so expand, move, and set.
822 	if (places_rdx)
823 	{
824 		bc_num_expand(n, bc_vm_growSize(n->len, places_rdx));
825 		// NOLINTNEXTLINE
826 		memmove(n->num + places_rdx, n->num, BC_NUM_SIZE(n->len));
827 		// NOLINTNEXTLINE
828 		memset(n->num, 0, BC_NUM_SIZE(places_rdx));
829 		n->len += places_rdx;
830 	}
831 
832 	// Set the scale appropriately.
833 	if (places > n->scale)
834 	{
835 		n->scale = 0;
836 		BC_NUM_RDX_SET(n, 0);
837 	}
838 	else
839 	{
840 		n->scale -= places;
841 		BC_NUM_RDX_SET(n, BC_NUM_RDX(n->scale));
842 	}
843 
844 	// Finally, shift within limbs.
845 	if (shift) bc_num_shift(n, BC_BASE_DIGS - dig);
846 
847 	bc_num_clean(n);
848 }
849 
850 void
851 bc_num_shiftRight(BcNum* restrict n, size_t places)
852 {
853 	BcBigDig dig;
854 	size_t places_rdx, scale, scale_mod, int_len, expand;
855 	bool shift;
856 
857 	if (!places) return;
858 
859 	// If zero, we can just set the scale and bail.
860 	if (BC_NUM_ZERO(n))
861 	{
862 		n->scale += places;
863 		bc_num_expand(n, BC_NUM_RDX(n->scale));
864 		return;
865 	}
866 
867 	// Amount within a limb we have to shift by.
868 	dig = (BcBigDig) (places % BC_BASE_DIGS);
869 	shift = (dig != 0);
870 
871 	scale = n->scale;
872 
873 	// Figure out how the scale is affected.
874 	scale_mod = scale % BC_BASE_DIGS;
875 	scale_mod = scale_mod ? scale_mod : BC_BASE_DIGS;
876 
877 	// We need to know the int length and rdx for places.
878 	int_len = bc_num_int(n);
879 	places_rdx = BC_NUM_RDX(places);
880 
881 	// If we are going to shift past a limb boundary or not, set accordingly.
882 	if (scale_mod + dig > BC_BASE_DIGS)
883 	{
884 		expand = places_rdx - 1;
885 		places_rdx = 1;
886 	}
887 	else
888 	{
889 		expand = places_rdx;
890 		places_rdx = 0;
891 	}
892 
893 	// Clamp expanding.
894 	if (expand > int_len) expand -= int_len;
895 	else expand = 0;
896 
897 	// Extend, expand, and zero.
898 	bc_num_extend(n, places_rdx * BC_BASE_DIGS);
899 	bc_num_expand(n, bc_vm_growSize(expand, n->len));
900 	// NOLINTNEXTLINE
901 	memset(n->num + n->len, 0, BC_NUM_SIZE(expand));
902 
903 	// Set the fields.
904 	n->len += expand;
905 	n->scale = 0;
906 	BC_NUM_RDX_SET(n, 0);
907 
908 	// Finally, shift within limbs.
909 	if (shift) bc_num_shift(n, dig);
910 
911 	n->scale = scale + places;
912 	BC_NUM_RDX_SET(n, BC_NUM_RDX(n->scale));
913 
914 	bc_num_clean(n);
915 
916 	assert(BC_NUM_RDX_VAL(n) <= n->len && n->len <= n->cap);
917 	assert(BC_NUM_RDX_VAL(n) == BC_NUM_RDX(n->scale));
918 }
919 
920 /**
921  * Invert @a into @a b at the current scale.
922  * @param a      The number to invert.
923  * @param b      The return parameter. This must be preallocated.
924  * @param scale  The current scale.
925  */
926 static inline void
927 bc_num_inv(BcNum* a, BcNum* b, size_t scale)
928 {
929 	assert(BC_NUM_NONZERO(a));
930 	bc_num_div(&vm.one, a, b, scale);
931 }
932 
933 /**
934  * Tests if a number is a integer with scale or not. Returns true if the number
935  * is not an integer. If it is, its integer shifted form is copied into the
936  * result parameter for use where only integers are allowed.
937  * @param n  The integer to test and shift.
938  * @param r  The number to store the shifted result into. This number should
939  *           *not* be allocated.
940  * @return   True if the number is a non-integer, false otherwise.
941  */
942 static bool
943 bc_num_nonInt(const BcNum* restrict n, BcNum* restrict r)
944 {
945 	bool zero;
946 	size_t i, rdx = BC_NUM_RDX_VAL(n);
947 
948 	if (!rdx)
949 	{
950 		// NOLINTNEXTLINE
951 		memcpy(r, n, sizeof(BcNum));
952 		return false;
953 	}
954 
955 	zero = true;
956 
957 	for (i = 0; zero && i < rdx; ++i)
958 	{
959 		zero = (n->num[i] == 0);
960 	}
961 
962 	if (BC_ERR(!zero)) return true;
963 
964 	bc_num_shiftRdx(n, r);
965 
966 	return false;
967 }
968 
969 #if BC_ENABLE_EXTRA_MATH
970 
971 /**
972  * Execute common code for an operater that needs an integer for the second
973  * operand and return the integer operand as a BcBigDig.
974  * @param a  The first operand.
975  * @param b  The second operand.
976  * @param c  The result operand.
977  * @return   The second operand as a hardware integer.
978  */
979 static BcBigDig
980 bc_num_intop(const BcNum* a, const BcNum* b, BcNum* restrict c)
981 {
982 	BcNum temp;
983 
984 	if (BC_ERR(bc_num_nonInt(b, &temp))) bc_err(BC_ERR_MATH_NON_INTEGER);
985 
986 	bc_num_copy(c, a);
987 
988 	return bc_num_bigdig(&temp);
989 }
990 #endif // BC_ENABLE_EXTRA_MATH
991 
992 /**
993  * This is the actual implementation of add *and* subtract. Since this function
994  * doesn't need to use scale (per the bc spec), I am hijacking it to say whether
995  * it's doing an add or a subtract. And then I convert substraction to addition
996  * of negative second operand. This is a BcNumBinOp function.
997  * @param a    The first operand.
998  * @param b    The second operand.
999  * @param c    The return parameter.
1000  * @param sub  Non-zero for a subtract, zero for an add.
1001  */
1002 static void
1003 bc_num_as(BcNum* a, BcNum* b, BcNum* restrict c, size_t sub)
1004 {
1005 	BcDig* ptr_c;
1006 	BcDig* ptr_l;
1007 	BcDig* ptr_r;
1008 	size_t i, min_rdx, max_rdx, diff, a_int, b_int, min_len, max_len, max_int;
1009 	size_t len_l, len_r, ardx, brdx;
1010 	bool b_neg, do_sub, do_rev_sub, carry, c_neg;
1011 
1012 	if (BC_NUM_ZERO(b))
1013 	{
1014 		bc_num_copy(c, a);
1015 		return;
1016 	}
1017 
1018 	if (BC_NUM_ZERO(a))
1019 	{
1020 		bc_num_copy(c, b);
1021 		c->rdx = BC_NUM_NEG_VAL(c, BC_NUM_NEG(b) != sub);
1022 		return;
1023 	}
1024 
1025 	// Invert sign of b if it is to be subtracted. This operation must
1026 	// precede the tests for any of the operands being zero.
1027 	b_neg = (BC_NUM_NEG(b) != sub);
1028 
1029 	// Figure out if we will actually add the numbers if their signs are equal
1030 	// or subtract.
1031 	do_sub = (BC_NUM_NEG(a) != b_neg);
1032 
1033 	a_int = bc_num_int(a);
1034 	b_int = bc_num_int(b);
1035 	max_int = BC_MAX(a_int, b_int);
1036 
1037 	// Figure out which number will have its last limbs copied (for addition) or
1038 	// subtracted (for subtraction).
1039 	ardx = BC_NUM_RDX_VAL(a);
1040 	brdx = BC_NUM_RDX_VAL(b);
1041 	min_rdx = BC_MIN(ardx, brdx);
1042 	max_rdx = BC_MAX(ardx, brdx);
1043 	diff = max_rdx - min_rdx;
1044 
1045 	max_len = max_int + max_rdx;
1046 
1047 	if (do_sub)
1048 	{
1049 		// Check whether b has to be subtracted from a or a from b.
1050 		if (a_int != b_int) do_rev_sub = (a_int < b_int);
1051 		else if (ardx > brdx)
1052 		{
1053 			do_rev_sub = (bc_num_compare(a->num + diff, b->num, b->len) < 0);
1054 		}
1055 		else do_rev_sub = (bc_num_compare(a->num, b->num + diff, a->len) <= 0);
1056 	}
1057 	else
1058 	{
1059 		// The result array of the addition might come out one element
1060 		// longer than the bigger of the operand arrays.
1061 		max_len += 1;
1062 		do_rev_sub = (a_int < b_int);
1063 	}
1064 
1065 	assert(max_len <= c->cap);
1066 
1067 	// Cache values for simple code later.
1068 	if (do_rev_sub)
1069 	{
1070 		ptr_l = b->num;
1071 		ptr_r = a->num;
1072 		len_l = b->len;
1073 		len_r = a->len;
1074 	}
1075 	else
1076 	{
1077 		ptr_l = a->num;
1078 		ptr_r = b->num;
1079 		len_l = a->len;
1080 		len_r = b->len;
1081 	}
1082 
1083 	ptr_c = c->num;
1084 	carry = false;
1085 
1086 	// This is true if the numbers have a different number of limbs after the
1087 	// decimal point.
1088 	if (diff)
1089 	{
1090 		// If the rdx values of the operands do not match, the result will
1091 		// have low end elements that are the positive or negative trailing
1092 		// elements of the operand with higher rdx value.
1093 		if ((ardx > brdx) != do_rev_sub)
1094 		{
1095 			// !do_rev_sub && ardx > brdx || do_rev_sub && brdx > ardx
1096 			// The left operand has BcDig values that need to be copied,
1097 			// either from a or from b (in case of a reversed subtraction).
1098 			// NOLINTNEXTLINE
1099 			memcpy(ptr_c, ptr_l, BC_NUM_SIZE(diff));
1100 			ptr_l += diff;
1101 			len_l -= diff;
1102 		}
1103 		else
1104 		{
1105 			// The right operand has BcDig values that need to be copied
1106 			// or subtracted from zero (in case of a subtraction).
1107 			if (do_sub)
1108 			{
1109 				// do_sub (do_rev_sub && ardx > brdx ||
1110 				// !do_rev_sub && brdx > ardx)
1111 				for (i = 0; i < diff; i++)
1112 				{
1113 					ptr_c[i] = bc_num_subDigits(0, ptr_r[i], &carry);
1114 				}
1115 			}
1116 			else
1117 			{
1118 				// !do_sub && brdx > ardx
1119 				// NOLINTNEXTLINE
1120 				memcpy(ptr_c, ptr_r, BC_NUM_SIZE(diff));
1121 			}
1122 
1123 			// Future code needs to ignore the limbs we just did.
1124 			ptr_r += diff;
1125 			len_r -= diff;
1126 		}
1127 
1128 		// The return value pointer needs to ignore what we just did.
1129 		ptr_c += diff;
1130 	}
1131 
1132 	// This is the length that can be directly added/subtracted.
1133 	min_len = BC_MIN(len_l, len_r);
1134 
1135 	// After dealing with possible low array elements that depend on only one
1136 	// operand above, the actual add or subtract can be performed as if the rdx
1137 	// of both operands was the same.
1138 	//
1139 	// Inlining takes care of eliminating constant zero arguments to
1140 	// addDigit/subDigit (checked in disassembly of resulting bc binary
1141 	// compiled with gcc and clang).
1142 	if (do_sub)
1143 	{
1144 		// Actual subtraction.
1145 		for (i = 0; i < min_len; ++i)
1146 		{
1147 			ptr_c[i] = bc_num_subDigits(ptr_l[i], ptr_r[i], &carry);
1148 		}
1149 
1150 		// Finishing the limbs beyond the direct subtraction.
1151 		for (; i < len_l; ++i)
1152 		{
1153 			ptr_c[i] = bc_num_subDigits(ptr_l[i], 0, &carry);
1154 		}
1155 	}
1156 	else
1157 	{
1158 		// Actual addition.
1159 		for (i = 0; i < min_len; ++i)
1160 		{
1161 			ptr_c[i] = bc_num_addDigits(ptr_l[i], ptr_r[i], &carry);
1162 		}
1163 
1164 		// Finishing the limbs beyond the direct addition.
1165 		for (; i < len_l; ++i)
1166 		{
1167 			ptr_c[i] = bc_num_addDigits(ptr_l[i], 0, &carry);
1168 		}
1169 
1170 		// Addition can create an extra limb. We take care of that here.
1171 		ptr_c[i] = bc_num_addDigits(0, 0, &carry);
1172 	}
1173 
1174 	assert(carry == false);
1175 
1176 	// The result has the same sign as a, unless the operation was a
1177 	// reverse subtraction (b - a).
1178 	c_neg = BC_NUM_NEG(a) != (do_sub && do_rev_sub);
1179 	BC_NUM_RDX_SET_NEG(c, max_rdx, c_neg);
1180 	c->len = max_len;
1181 	c->scale = BC_MAX(a->scale, b->scale);
1182 
1183 	bc_num_clean(c);
1184 }
1185 
1186 /**
1187  * The simple multiplication that karatsuba dishes out to when the length of the
1188  * numbers gets low enough. This doesn't use scale because it treats the
1189  * operands as though they are integers.
1190  * @param a  The first operand.
1191  * @param b  The second operand.
1192  * @param c  The return parameter.
1193  */
1194 static void
1195 bc_num_m_simp(const BcNum* a, const BcNum* b, BcNum* restrict c)
1196 {
1197 	size_t i, alen = a->len, blen = b->len, clen;
1198 	BcDig* ptr_a = a->num;
1199 	BcDig* ptr_b = b->num;
1200 	BcDig* ptr_c;
1201 	BcBigDig sum = 0, carry = 0;
1202 
1203 	assert(sizeof(sum) >= sizeof(BcDig) * 2);
1204 	assert(!BC_NUM_RDX_VAL(a) && !BC_NUM_RDX_VAL(b));
1205 
1206 	// Make sure c is big enough.
1207 	clen = bc_vm_growSize(alen, blen);
1208 	bc_num_expand(c, bc_vm_growSize(clen, 1));
1209 
1210 	// If we don't memset, then we might have uninitialized data use later.
1211 	ptr_c = c->num;
1212 	// NOLINTNEXTLINE
1213 	memset(ptr_c, 0, BC_NUM_SIZE(c->cap));
1214 
1215 	// This is the actual multiplication loop. It uses the lattice form of long
1216 	// multiplication (see the explanation on the web page at
1217 	// https://knilt.arcc.albany.edu/What_is_Lattice_Multiplication or the
1218 	// explanation at Wikipedia).
1219 	for (i = 0; i < clen; ++i)
1220 	{
1221 		ssize_t sidx = (ssize_t) (i - blen + 1);
1222 		size_t j, k;
1223 
1224 		// These are the start indices.
1225 		j = (size_t) BC_MAX(0, sidx);
1226 		k = BC_MIN(i, blen - 1);
1227 
1228 		// On every iteration of this loop, a multiplication happens, then the
1229 		// sum is automatically calculated.
1230 		for (; j < alen && k < blen; ++j, --k)
1231 		{
1232 			sum += ((BcBigDig) ptr_a[j]) * ((BcBigDig) ptr_b[k]);
1233 
1234 			if (sum >= ((BcBigDig) BC_BASE_POW) * BC_BASE_POW)
1235 			{
1236 				carry += sum / BC_BASE_POW;
1237 				sum %= BC_BASE_POW;
1238 			}
1239 		}
1240 
1241 		// Calculate the carry.
1242 		if (sum >= BC_BASE_POW)
1243 		{
1244 			carry += sum / BC_BASE_POW;
1245 			sum %= BC_BASE_POW;
1246 		}
1247 
1248 		// Store and set up for next iteration.
1249 		ptr_c[i] = (BcDig) sum;
1250 		assert(ptr_c[i] < BC_BASE_POW);
1251 		sum = carry;
1252 		carry = 0;
1253 	}
1254 
1255 	// This should always be true because there should be no carry on the last
1256 	// digit; multiplication never goes above the sum of both lengths.
1257 	assert(!sum);
1258 
1259 	c->len = clen;
1260 }
1261 
1262 /**
1263  * Does a shifted add or subtract for Karatsuba below. This calls either
1264  * bc_num_addArrays() or bc_num_subArrays().
1265  * @param n      An in/out parameter; the first operand and return parameter.
1266  * @param a      The second operand.
1267  * @param shift  The amount to shift @a n by when adding/subtracting.
1268  * @param op     The function to call, either bc_num_addArrays() or
1269  *               bc_num_subArrays().
1270  */
1271 static void
1272 bc_num_shiftAddSub(BcNum* restrict n, const BcNum* restrict a, size_t shift,
1273                    BcNumShiftAddOp op)
1274 {
1275 	assert(n->len >= shift + a->len);
1276 	assert(!BC_NUM_RDX_VAL(n) && !BC_NUM_RDX_VAL(a));
1277 	op(n->num + shift, a->num, a->len);
1278 }
1279 
1280 /**
1281  * Implements the Karatsuba algorithm.
1282  */
1283 static void
1284 bc_num_k(const BcNum* a, const BcNum* b, BcNum* restrict c)
1285 {
1286 	size_t max, max2, total;
1287 	BcNum l1, h1, l2, h2, m2, m1, z0, z1, z2, temp;
1288 	BcDig* digs;
1289 	BcDig* dig_ptr;
1290 	BcNumShiftAddOp op;
1291 	bool aone = BC_NUM_ONE(a);
1292 
1293 	assert(BC_NUM_ZERO(c));
1294 
1295 	if (BC_NUM_ZERO(a) || BC_NUM_ZERO(b)) return;
1296 
1297 	if (aone || BC_NUM_ONE(b))
1298 	{
1299 		bc_num_copy(c, aone ? b : a);
1300 		if ((aone && BC_NUM_NEG(a)) || BC_NUM_NEG(b)) BC_NUM_NEG_TGL(c);
1301 		return;
1302 	}
1303 
1304 	// Shell out to the simple algorithm with certain conditions.
1305 	if (a->len < BC_NUM_KARATSUBA_LEN || b->len < BC_NUM_KARATSUBA_LEN)
1306 	{
1307 		bc_num_m_simp(a, b, c);
1308 		return;
1309 	}
1310 
1311 	// We need to calculate the max size of the numbers that can result from the
1312 	// operations.
1313 	max = BC_MAX(a->len, b->len);
1314 	max = BC_MAX(max, BC_NUM_DEF_SIZE);
1315 	max2 = (max + 1) / 2;
1316 
1317 	// Calculate the space needed for all of the temporary allocations. We do
1318 	// this to just allocate once.
1319 	total = bc_vm_arraySize(BC_NUM_KARATSUBA_ALLOCS, max);
1320 
1321 	BC_SIG_LOCK;
1322 
1323 	// Allocate space for all of the temporaries.
1324 	digs = dig_ptr = bc_vm_malloc(BC_NUM_SIZE(total));
1325 
1326 	// Set up the temporaries.
1327 	bc_num_setup(&l1, dig_ptr, max);
1328 	dig_ptr += max;
1329 	bc_num_setup(&h1, dig_ptr, max);
1330 	dig_ptr += max;
1331 	bc_num_setup(&l2, dig_ptr, max);
1332 	dig_ptr += max;
1333 	bc_num_setup(&h2, dig_ptr, max);
1334 	dig_ptr += max;
1335 	bc_num_setup(&m1, dig_ptr, max);
1336 	dig_ptr += max;
1337 	bc_num_setup(&m2, dig_ptr, max);
1338 
1339 	// Some temporaries need the ability to grow, so we allocate them
1340 	// separately.
1341 	max = bc_vm_growSize(max, 1);
1342 	bc_num_init(&z0, max);
1343 	bc_num_init(&z1, max);
1344 	bc_num_init(&z2, max);
1345 	max = bc_vm_growSize(max, max) + 1;
1346 	bc_num_init(&temp, max);
1347 
1348 	BC_SETJMP_LOCKED(err);
1349 
1350 	BC_SIG_UNLOCK;
1351 
1352 	// First, set up c.
1353 	bc_num_expand(c, max);
1354 	c->len = max;
1355 	// NOLINTNEXTLINE
1356 	memset(c->num, 0, BC_NUM_SIZE(c->len));
1357 
1358 	// Split the parameters.
1359 	bc_num_split(a, max2, &l1, &h1);
1360 	bc_num_split(b, max2, &l2, &h2);
1361 
1362 	// Do the subtraction.
1363 	bc_num_sub(&h1, &l1, &m1, 0);
1364 	bc_num_sub(&l2, &h2, &m2, 0);
1365 
1366 	// The if statements below are there for efficiency reasons. The best way to
1367 	// understand them is to understand the Karatsuba algorithm because now that
1368 	// the ollocations and splits are done, the algorithm is pretty
1369 	// straightforward.
1370 
1371 	if (BC_NUM_NONZERO(&h1) && BC_NUM_NONZERO(&h2))
1372 	{
1373 		assert(BC_NUM_RDX_VALID_NP(h1));
1374 		assert(BC_NUM_RDX_VALID_NP(h2));
1375 
1376 		bc_num_m(&h1, &h2, &z2, 0);
1377 		bc_num_clean(&z2);
1378 
1379 		bc_num_shiftAddSub(c, &z2, max2 * 2, bc_num_addArrays);
1380 		bc_num_shiftAddSub(c, &z2, max2, bc_num_addArrays);
1381 	}
1382 
1383 	if (BC_NUM_NONZERO(&l1) && BC_NUM_NONZERO(&l2))
1384 	{
1385 		assert(BC_NUM_RDX_VALID_NP(l1));
1386 		assert(BC_NUM_RDX_VALID_NP(l2));
1387 
1388 		bc_num_m(&l1, &l2, &z0, 0);
1389 		bc_num_clean(&z0);
1390 
1391 		bc_num_shiftAddSub(c, &z0, max2, bc_num_addArrays);
1392 		bc_num_shiftAddSub(c, &z0, 0, bc_num_addArrays);
1393 	}
1394 
1395 	if (BC_NUM_NONZERO(&m1) && BC_NUM_NONZERO(&m2))
1396 	{
1397 		assert(BC_NUM_RDX_VALID_NP(m1));
1398 		assert(BC_NUM_RDX_VALID_NP(m1));
1399 
1400 		bc_num_m(&m1, &m2, &z1, 0);
1401 		bc_num_clean(&z1);
1402 
1403 		op = (BC_NUM_NEG_NP(m1) != BC_NUM_NEG_NP(m2)) ?
1404 		         bc_num_subArrays :
1405 		         bc_num_addArrays;
1406 		bc_num_shiftAddSub(c, &z1, max2, op);
1407 	}
1408 
1409 err:
1410 	BC_SIG_MAYLOCK;
1411 	free(digs);
1412 	bc_num_free(&temp);
1413 	bc_num_free(&z2);
1414 	bc_num_free(&z1);
1415 	bc_num_free(&z0);
1416 	BC_LONGJMP_CONT;
1417 }
1418 
1419 /**
1420  * Does checks for Karatsuba. It also changes things to ensure that the
1421  * Karatsuba and simple multiplication can treat the numbers as integers. This
1422  * is a BcNumBinOp function.
1423  * @param a      The first operand.
1424  * @param b      The second operand.
1425  * @param c      The return parameter.
1426  * @param scale  The current scale.
1427  */
1428 static void
1429 bc_num_m(BcNum* a, BcNum* b, BcNum* restrict c, size_t scale)
1430 {
1431 	BcNum cpa, cpb;
1432 	size_t ascale, bscale, ardx, brdx, azero = 0, bzero = 0, zero, len, rscale;
1433 
1434 	assert(BC_NUM_RDX_VALID(a));
1435 	assert(BC_NUM_RDX_VALID(b));
1436 
1437 	bc_num_zero(c);
1438 
1439 	ascale = a->scale;
1440 	bscale = b->scale;
1441 
1442 	// This sets the final scale according to the bc spec.
1443 	scale = BC_MAX(scale, ascale);
1444 	scale = BC_MAX(scale, bscale);
1445 	rscale = ascale + bscale;
1446 	scale = BC_MIN(rscale, scale);
1447 
1448 	// If this condition is true, we can use bc_num_mulArray(), which would be
1449 	// much faster.
1450 	if ((a->len == 1 || b->len == 1) && !a->rdx && !b->rdx)
1451 	{
1452 		BcNum* operand;
1453 		BcBigDig dig;
1454 
1455 		// Set the correct operands.
1456 		if (a->len == 1)
1457 		{
1458 			dig = (BcBigDig) a->num[0];
1459 			operand = b;
1460 		}
1461 		else
1462 		{
1463 			dig = (BcBigDig) b->num[0];
1464 			operand = a;
1465 		}
1466 
1467 		bc_num_mulArray(operand, dig, c);
1468 
1469 		// Need to make sure the sign is correct.
1470 		if (BC_NUM_NONZERO(c))
1471 		{
1472 			c->rdx = BC_NUM_NEG_VAL(c, BC_NUM_NEG(a) != BC_NUM_NEG(b));
1473 		}
1474 
1475 		return;
1476 	}
1477 
1478 	assert(BC_NUM_RDX_VALID(a));
1479 	assert(BC_NUM_RDX_VALID(b));
1480 
1481 	BC_SIG_LOCK;
1482 
1483 	// We need copies because of all of the mutation needed to make Karatsuba
1484 	// think the numbers are integers.
1485 	bc_num_init(&cpa, a->len + BC_NUM_RDX_VAL(a));
1486 	bc_num_init(&cpb, b->len + BC_NUM_RDX_VAL(b));
1487 
1488 	BC_SETJMP_LOCKED(err);
1489 
1490 	BC_SIG_UNLOCK;
1491 
1492 	bc_num_copy(&cpa, a);
1493 	bc_num_copy(&cpb, b);
1494 
1495 	assert(BC_NUM_RDX_VALID_NP(cpa));
1496 	assert(BC_NUM_RDX_VALID_NP(cpb));
1497 
1498 	BC_NUM_NEG_CLR_NP(cpa);
1499 	BC_NUM_NEG_CLR_NP(cpb);
1500 
1501 	assert(BC_NUM_RDX_VALID_NP(cpa));
1502 	assert(BC_NUM_RDX_VALID_NP(cpb));
1503 
1504 	// These are what makes them appear like integers.
1505 	ardx = BC_NUM_RDX_VAL_NP(cpa) * BC_BASE_DIGS;
1506 	bc_num_shiftLeft(&cpa, ardx);
1507 
1508 	brdx = BC_NUM_RDX_VAL_NP(cpb) * BC_BASE_DIGS;
1509 	bc_num_shiftLeft(&cpb, brdx);
1510 
1511 	// We need to reset the jump here because azero and bzero are used in the
1512 	// cleanup, and local variables are not guaranteed to be the same after a
1513 	// jump.
1514 	BC_SIG_LOCK;
1515 
1516 	BC_UNSETJMP;
1517 
1518 	// We want to ignore zero limbs.
1519 	azero = bc_num_shiftZero(&cpa);
1520 	bzero = bc_num_shiftZero(&cpb);
1521 
1522 	BC_SETJMP_LOCKED(err);
1523 
1524 	BC_SIG_UNLOCK;
1525 
1526 	bc_num_clean(&cpa);
1527 	bc_num_clean(&cpb);
1528 
1529 	bc_num_k(&cpa, &cpb, c);
1530 
1531 	// The return parameter needs to have its scale set. This is the start. It
1532 	// also needs to be shifted by the same amount as a and b have limbs after
1533 	// the decimal point.
1534 	zero = bc_vm_growSize(azero, bzero);
1535 	len = bc_vm_growSize(c->len, zero);
1536 
1537 	bc_num_expand(c, len);
1538 
1539 	// Shift c based on the limbs after the decimal point in a and b.
1540 	bc_num_shiftLeft(c, (len - c->len) * BC_BASE_DIGS);
1541 	bc_num_shiftRight(c, ardx + brdx);
1542 
1543 	bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));
1544 
1545 err:
1546 	BC_SIG_MAYLOCK;
1547 	bc_num_unshiftZero(&cpb, bzero);
1548 	bc_num_free(&cpb);
1549 	bc_num_unshiftZero(&cpa, azero);
1550 	bc_num_free(&cpa);
1551 	BC_LONGJMP_CONT;
1552 }
1553 
1554 /**
1555  * Returns true if the BcDig array has non-zero limbs, false otherwise.
1556  * @param a    The array to test.
1557  * @param len  The length of the array.
1558  * @return     True if @a has any non-zero limbs, false otherwise.
1559  */
1560 static bool
1561 bc_num_nonZeroDig(BcDig* restrict a, size_t len)
1562 {
1563 	size_t i;
1564 
1565 	bool nonzero = false;
1566 
1567 	for (i = len - 1; !nonzero && i < len; --i)
1568 	{
1569 		nonzero = (a[i] != 0);
1570 	}
1571 
1572 	return nonzero;
1573 }
1574 
1575 /**
1576  * Compares a BcDig array against a BcNum. This is especially suited for
1577  * division. Returns >0 if @a a is greater than @a b, <0 if it is less, and =0
1578  * if they are equal.
1579  * @param a    The array.
1580  * @param b    The number.
1581  * @param len  The length to assume the arrays are. This is always less than the
1582  *             actual length because of how this is implemented.
1583  */
1584 static ssize_t
1585 bc_num_divCmp(const BcDig* a, const BcNum* b, size_t len)
1586 {
1587 	ssize_t cmp;
1588 
1589 	if (b->len > len && a[len]) cmp = bc_num_compare(a, b->num, len + 1);
1590 	else if (b->len <= len)
1591 	{
1592 		if (a[len]) cmp = 1;
1593 		else cmp = bc_num_compare(a, b->num, len);
1594 	}
1595 	else cmp = -1;
1596 
1597 	return cmp;
1598 }
1599 
1600 /**
1601  * Extends the two operands of a division by BC_BASE_DIGS minus the number of
1602  * digits in the divisor estimate. In other words, it is shifting the numbers in
1603  * order to force the divisor estimate to fill the limb.
1604  * @param a        The first operand.
1605  * @param b        The second operand.
1606  * @param divisor  The divisor estimate.
1607  */
1608 static void
1609 bc_num_divExtend(BcNum* restrict a, BcNum* restrict b, BcBigDig divisor)
1610 {
1611 	size_t pow;
1612 
1613 	assert(divisor < BC_BASE_POW);
1614 
1615 	pow = BC_BASE_DIGS - bc_num_log10((size_t) divisor);
1616 
1617 	bc_num_shiftLeft(a, pow);
1618 	bc_num_shiftLeft(b, pow);
1619 }
1620 
1621 /**
1622  * Actually does division. This is a rewrite of my original code by Stefan Esser
1623  * from FreeBSD.
1624  * @param a      The first operand.
1625  * @param b      The second operand.
1626  * @param c      The return parameter.
1627  * @param scale  The current scale.
1628  */
1629 static void
1630 bc_num_d_long(BcNum* restrict a, BcNum* restrict b, BcNum* restrict c,
1631               size_t scale)
1632 {
1633 	BcBigDig divisor;
1634 	size_t len, end, i, rdx;
1635 	BcNum cpb;
1636 	bool nonzero = false;
1637 
1638 	assert(b->len < a->len);
1639 
1640 	len = b->len;
1641 	end = a->len - len;
1642 
1643 	assert(len >= 1);
1644 
1645 	// This is a final time to make sure c is big enough and that its array is
1646 	// properly zeroed.
1647 	bc_num_expand(c, a->len);
1648 	// NOLINTNEXTLINE
1649 	memset(c->num, 0, c->cap * sizeof(BcDig));
1650 
1651 	// Setup.
1652 	BC_NUM_RDX_SET(c, BC_NUM_RDX_VAL(a));
1653 	c->scale = a->scale;
1654 	c->len = a->len;
1655 
1656 	// This is pulling the most significant limb of b in order to establish a
1657 	// good "estimate" for the actual divisor.
1658 	divisor = (BcBigDig) b->num[len - 1];
1659 
1660 	// The entire bit of code in this if statement is to tighten the estimate of
1661 	// the divisor. The condition asks if b has any other non-zero limbs.
1662 	if (len > 1 && bc_num_nonZeroDig(b->num, len - 1))
1663 	{
1664 		// This takes a little bit of understanding. The "10*BC_BASE_DIGS/6+1"
1665 		// results in either 16 for 64-bit 9-digit limbs or 7 for 32-bit 4-digit
1666 		// limbs. Then it shifts a 1 by that many, which in both cases, puts the
1667 		// result above *half* of the max value a limb can store. Basically,
1668 		// this quickly calculates if the divisor is greater than half the max
1669 		// of a limb.
1670 		nonzero = (divisor > 1 << ((10 * BC_BASE_DIGS) / 6 + 1));
1671 
1672 		// If the divisor is *not* greater than half the limb...
1673 		if (!nonzero)
1674 		{
1675 			// Extend the parameters by the number of missing digits in the
1676 			// divisor.
1677 			bc_num_divExtend(a, b, divisor);
1678 
1679 			// Check bc_num_d(). In there, we grow a again and again. We do it
1680 			// again here; we *always* want to be sure it is big enough.
1681 			len = BC_MAX(a->len, b->len);
1682 			bc_num_expand(a, len + 1);
1683 
1684 			// Make a have a zero most significant limb to match the len.
1685 			if (len + 1 > a->len) a->len = len + 1;
1686 
1687 			// Grab the new divisor estimate, new because the shift has made it
1688 			// different.
1689 			len = b->len;
1690 			end = a->len - len;
1691 			divisor = (BcBigDig) b->num[len - 1];
1692 
1693 			nonzero = bc_num_nonZeroDig(b->num, len - 1);
1694 		}
1695 	}
1696 
1697 	// If b has other nonzero limbs, we want the divisor to be one higher, so
1698 	// that it is an upper bound.
1699 	divisor += nonzero;
1700 
1701 	// Make sure c can fit the new length.
1702 	bc_num_expand(c, a->len);
1703 	// NOLINTNEXTLINE
1704 	memset(c->num, 0, BC_NUM_SIZE(c->cap));
1705 
1706 	assert(c->scale >= scale);
1707 	rdx = BC_NUM_RDX_VAL(c) - BC_NUM_RDX(scale);
1708 
1709 	BC_SIG_LOCK;
1710 
1711 	bc_num_init(&cpb, len + 1);
1712 
1713 	BC_SETJMP_LOCKED(err);
1714 
1715 	BC_SIG_UNLOCK;
1716 
1717 	// This is the actual division loop.
1718 	for (i = end - 1; i < end && i >= rdx && BC_NUM_NONZERO(a); --i)
1719 	{
1720 		ssize_t cmp;
1721 		BcDig* n;
1722 		BcBigDig result;
1723 
1724 		n = a->num + i;
1725 		assert(n >= a->num);
1726 		result = 0;
1727 
1728 		cmp = bc_num_divCmp(n, b, len);
1729 
1730 		// This is true if n is greater than b, which means that division can
1731 		// proceed, so this inner loop is the part that implements one instance
1732 		// of the division.
1733 		while (cmp >= 0)
1734 		{
1735 			BcBigDig n1, dividend, quotient;
1736 
1737 			// These should be named obviously enough. Just imagine that it's a
1738 			// division of one limb. Because that's what it is.
1739 			n1 = (BcBigDig) n[len];
1740 			dividend = n1 * BC_BASE_POW + (BcBigDig) n[len - 1];
1741 			quotient = (dividend / divisor);
1742 
1743 			// If this is true, then we can just subtract. Remember: setting
1744 			// quotient to 1 is not bad because we already know that n is
1745 			// greater than b.
1746 			if (quotient <= 1)
1747 			{
1748 				quotient = 1;
1749 				bc_num_subArrays(n, b->num, len);
1750 			}
1751 			else
1752 			{
1753 				assert(quotient <= BC_BASE_POW);
1754 
1755 				// We need to multiply and subtract for a quotient above 1.
1756 				bc_num_mulArray(b, (BcBigDig) quotient, &cpb);
1757 				bc_num_subArrays(n, cpb.num, cpb.len);
1758 			}
1759 
1760 			// The result is the *real* quotient, by the way, but it might take
1761 			// multiple trips around this loop to get it.
1762 			result += quotient;
1763 			assert(result <= BC_BASE_POW);
1764 
1765 			// And here's why it might take multiple trips: n might *still* be
1766 			// greater than b. So we have to loop again. That's what this is
1767 			// setting up for: the condition of the while loop.
1768 			if (nonzero) cmp = bc_num_divCmp(n, b, len);
1769 			else cmp = -1;
1770 		}
1771 
1772 		assert(result < BC_BASE_POW);
1773 
1774 		// Store the actual limb quotient.
1775 		c->num[i] = (BcDig) result;
1776 	}
1777 
1778 err:
1779 	BC_SIG_MAYLOCK;
1780 	bc_num_free(&cpb);
1781 	BC_LONGJMP_CONT;
1782 }
1783 
1784 /**
1785  * Implements division. This is a BcNumBinOp function.
1786  * @param a      The first operand.
1787  * @param b      The second operand.
1788  * @param c      The return parameter.
1789  * @param scale  The current scale.
1790  */
1791 static void
1792 bc_num_d(BcNum* a, BcNum* b, BcNum* restrict c, size_t scale)
1793 {
1794 	size_t len, cpardx;
1795 	BcNum cpa, cpb;
1796 
1797 	if (BC_NUM_ZERO(b)) bc_err(BC_ERR_MATH_DIVIDE_BY_ZERO);
1798 
1799 	if (BC_NUM_ZERO(a))
1800 	{
1801 		bc_num_setToZero(c, scale);
1802 		return;
1803 	}
1804 
1805 	if (BC_NUM_ONE(b))
1806 	{
1807 		bc_num_copy(c, a);
1808 		bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));
1809 		return;
1810 	}
1811 
1812 	// If this is true, we can use bc_num_divArray(), which would be faster.
1813 	if (!BC_NUM_RDX_VAL(a) && !BC_NUM_RDX_VAL(b) && b->len == 1 && !scale)
1814 	{
1815 		BcBigDig rem;
1816 		bc_num_divArray(a, (BcBigDig) b->num[0], c, &rem);
1817 		bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));
1818 		return;
1819 	}
1820 
1821 	len = bc_num_divReq(a, b, scale);
1822 
1823 	BC_SIG_LOCK;
1824 
1825 	// Initialize copies of the parameters. We want the length of the first
1826 	// operand copy to be as big as the result because of the way the division
1827 	// is implemented.
1828 	bc_num_init(&cpa, len);
1829 	bc_num_copy(&cpa, a);
1830 	bc_num_createCopy(&cpb, b);
1831 
1832 	BC_SETJMP_LOCKED(err);
1833 
1834 	BC_SIG_UNLOCK;
1835 
1836 	len = b->len;
1837 
1838 	// Like the above comment, we want the copy of the first parameter to be
1839 	// larger than the second parameter.
1840 	if (len > cpa.len)
1841 	{
1842 		bc_num_expand(&cpa, bc_vm_growSize(len, 2));
1843 		bc_num_extend(&cpa, (len - cpa.len) * BC_BASE_DIGS);
1844 	}
1845 
1846 	cpardx = BC_NUM_RDX_VAL_NP(cpa);
1847 	cpa.scale = cpardx * BC_BASE_DIGS;
1848 
1849 	// This is just setting up the scale in preparation for the division.
1850 	bc_num_extend(&cpa, b->scale);
1851 	cpardx = BC_NUM_RDX_VAL_NP(cpa) - BC_NUM_RDX(b->scale);
1852 	BC_NUM_RDX_SET_NP(cpa, cpardx);
1853 	cpa.scale = cpardx * BC_BASE_DIGS;
1854 
1855 	// Once again, just setting things up, this time to match scale.
1856 	if (scale > cpa.scale)
1857 	{
1858 		bc_num_extend(&cpa, scale);
1859 		cpardx = BC_NUM_RDX_VAL_NP(cpa);
1860 		cpa.scale = cpardx * BC_BASE_DIGS;
1861 	}
1862 
1863 	// Grow if necessary.
1864 	if (cpa.cap == cpa.len) bc_num_expand(&cpa, bc_vm_growSize(cpa.len, 1));
1865 
1866 	// We want an extra zero in front to make things simpler.
1867 	cpa.num[cpa.len++] = 0;
1868 
1869 	// Still setting things up. Why all of these things are needed is not
1870 	// something that can be easily explained, but it has to do with making the
1871 	// actual algorithm easier to understand because it can assume a lot of
1872 	// things. Thus, you should view all of this setup code as establishing
1873 	// assumptions for bc_num_d_long(), where the actual division happens.
1874 	if (cpardx == cpa.len) cpa.len = bc_num_nonZeroLen(&cpa);
1875 	if (BC_NUM_RDX_VAL_NP(cpb) == cpb.len) cpb.len = bc_num_nonZeroLen(&cpb);
1876 	cpb.scale = 0;
1877 	BC_NUM_RDX_SET_NP(cpb, 0);
1878 
1879 	bc_num_d_long(&cpa, &cpb, c, scale);
1880 
1881 	bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));
1882 
1883 err:
1884 	BC_SIG_MAYLOCK;
1885 	bc_num_free(&cpb);
1886 	bc_num_free(&cpa);
1887 	BC_LONGJMP_CONT;
1888 }
1889 
1890 /**
1891  * Implements divmod. This is the actual modulus function; since modulus
1892  * requires a division anyway, this returns the quotient and modulus. Either can
1893  * be thrown out as desired.
1894  * @param a      The first operand.
1895  * @param b      The second operand.
1896  * @param c      The return parameter for the quotient.
1897  * @param d      The return parameter for the modulus.
1898  * @param scale  The current scale.
1899  * @param ts     The scale that the operation should be done to. Yes, it's not
1900  *               necessarily the same as scale, per the bc spec.
1901  */
1902 static void
1903 bc_num_r(BcNum* a, BcNum* b, BcNum* restrict c, BcNum* restrict d, size_t scale,
1904          size_t ts)
1905 {
1906 	BcNum temp;
1907 	bool neg;
1908 
1909 	if (BC_NUM_ZERO(b)) bc_err(BC_ERR_MATH_DIVIDE_BY_ZERO);
1910 
1911 	if (BC_NUM_ZERO(a))
1912 	{
1913 		bc_num_setToZero(c, ts);
1914 		bc_num_setToZero(d, ts);
1915 		return;
1916 	}
1917 
1918 	BC_SIG_LOCK;
1919 
1920 	bc_num_init(&temp, d->cap);
1921 
1922 	BC_SETJMP_LOCKED(err);
1923 
1924 	BC_SIG_UNLOCK;
1925 
1926 	// Division.
1927 	bc_num_d(a, b, c, scale);
1928 
1929 	// We want an extra digit so we can safely truncate.
1930 	if (scale) scale = ts + 1;
1931 
1932 	assert(BC_NUM_RDX_VALID(c));
1933 	assert(BC_NUM_RDX_VALID(b));
1934 
1935 	// Implement the rest of the (a - (a / b) * b) formula.
1936 	bc_num_m(c, b, &temp, scale);
1937 	bc_num_sub(a, &temp, d, scale);
1938 
1939 	// Extend if necessary.
1940 	if (ts > d->scale && BC_NUM_NONZERO(d)) bc_num_extend(d, ts - d->scale);
1941 
1942 	neg = BC_NUM_NEG(d);
1943 	bc_num_retireMul(d, ts, BC_NUM_NEG(a), BC_NUM_NEG(b));
1944 	d->rdx = BC_NUM_NEG_VAL(d, BC_NUM_NONZERO(d) ? neg : false);
1945 
1946 err:
1947 	BC_SIG_MAYLOCK;
1948 	bc_num_free(&temp);
1949 	BC_LONGJMP_CONT;
1950 }
1951 
1952 /**
1953  * Implements modulus/remainder. (Yes, I know they are different, but not in the
1954  * context of bc.) This is a BcNumBinOp function.
1955  * @param a      The first operand.
1956  * @param b      The second operand.
1957  * @param c      The return parameter.
1958  * @param scale  The current scale.
1959  */
1960 static void
1961 bc_num_rem(BcNum* a, BcNum* b, BcNum* restrict c, size_t scale)
1962 {
1963 	BcNum c1;
1964 	size_t ts;
1965 
1966 	ts = bc_vm_growSize(scale, b->scale);
1967 	ts = BC_MAX(ts, a->scale);
1968 
1969 	BC_SIG_LOCK;
1970 
1971 	// Need a temp for the quotient.
1972 	bc_num_init(&c1, bc_num_mulReq(a, b, ts));
1973 
1974 	BC_SETJMP_LOCKED(err);
1975 
1976 	BC_SIG_UNLOCK;
1977 
1978 	bc_num_r(a, b, &c1, c, scale, ts);
1979 
1980 err:
1981 	BC_SIG_MAYLOCK;
1982 	bc_num_free(&c1);
1983 	BC_LONGJMP_CONT;
1984 }
1985 
1986 /**
1987  * Implements power (exponentiation). This is a BcNumBinOp function.
1988  * @param a      The first operand.
1989  * @param b      The second operand.
1990  * @param c      The return parameter.
1991  * @param scale  The current scale.
1992  */
1993 static void
1994 bc_num_p(BcNum* a, BcNum* b, BcNum* restrict c, size_t scale)
1995 {
1996 	BcNum copy, btemp;
1997 	BcBigDig exp;
1998 	size_t powrdx, resrdx;
1999 	bool neg;
2000 
2001 	if (BC_ERR(bc_num_nonInt(b, &btemp))) bc_err(BC_ERR_MATH_NON_INTEGER);
2002 
2003 	if (BC_NUM_ZERO(&btemp))
2004 	{
2005 		bc_num_one(c);
2006 		return;
2007 	}
2008 
2009 	if (BC_NUM_ZERO(a))
2010 	{
2011 		if (BC_NUM_NEG_NP(btemp)) bc_err(BC_ERR_MATH_DIVIDE_BY_ZERO);
2012 		bc_num_setToZero(c, scale);
2013 		return;
2014 	}
2015 
2016 	if (BC_NUM_ONE(&btemp))
2017 	{
2018 		if (!BC_NUM_NEG_NP(btemp)) bc_num_copy(c, a);
2019 		else bc_num_inv(a, c, scale);
2020 		return;
2021 	}
2022 
2023 	neg = BC_NUM_NEG_NP(btemp);
2024 	BC_NUM_NEG_CLR_NP(btemp);
2025 
2026 	exp = bc_num_bigdig(&btemp);
2027 
2028 	BC_SIG_LOCK;
2029 
2030 	bc_num_createCopy(&copy, a);
2031 
2032 	BC_SETJMP_LOCKED(err);
2033 
2034 	BC_SIG_UNLOCK;
2035 
2036 	// If this is true, then we do not have to do a division, and we need to
2037 	// set scale accordingly.
2038 	if (!neg)
2039 	{
2040 		size_t max = BC_MAX(scale, a->scale), scalepow;
2041 		scalepow = bc_num_mulOverflow(a->scale, exp);
2042 		scale = BC_MIN(scalepow, max);
2043 	}
2044 
2045 	// This is only implementing the first exponentiation by squaring, until it
2046 	// reaches the first time where the square is actually used.
2047 	for (powrdx = a->scale; !(exp & 1); exp >>= 1)
2048 	{
2049 		powrdx <<= 1;
2050 		assert(BC_NUM_RDX_VALID_NP(copy));
2051 		bc_num_mul(&copy, &copy, &copy, powrdx);
2052 	}
2053 
2054 	// Make c a copy of copy for the purpose of saving the squares that should
2055 	// be saved.
2056 	bc_num_copy(c, &copy);
2057 	resrdx = powrdx;
2058 
2059 	// Now finish the exponentiation by squaring, this time saving the squares
2060 	// as necessary.
2061 	while (exp >>= 1)
2062 	{
2063 		powrdx <<= 1;
2064 		assert(BC_NUM_RDX_VALID_NP(copy));
2065 		bc_num_mul(&copy, &copy, &copy, powrdx);
2066 
2067 		// If this is true, we want to save that particular square. This does
2068 		// that by multiplying c with copy.
2069 		if (exp & 1)
2070 		{
2071 			resrdx += powrdx;
2072 			assert(BC_NUM_RDX_VALID(c));
2073 			assert(BC_NUM_RDX_VALID_NP(copy));
2074 			bc_num_mul(c, &copy, c, resrdx);
2075 		}
2076 	}
2077 
2078 	// Invert if necessary.
2079 	if (neg) bc_num_inv(c, c, scale);
2080 
2081 	// Truncate if necessary.
2082 	if (c->scale > scale) bc_num_truncate(c, c->scale - scale);
2083 
2084 	bc_num_clean(c);
2085 
2086 err:
2087 	BC_SIG_MAYLOCK;
2088 	bc_num_free(&copy);
2089 	BC_LONGJMP_CONT;
2090 }
2091 
2092 #if BC_ENABLE_EXTRA_MATH
2093 /**
2094  * Implements the places operator. This is a BcNumBinOp function.
2095  * @param a      The first operand.
2096  * @param b      The second operand.
2097  * @param c      The return parameter.
2098  * @param scale  The current scale.
2099  */
2100 static void
2101 bc_num_place(BcNum* a, BcNum* b, BcNum* restrict c, size_t scale)
2102 {
2103 	BcBigDig val;
2104 
2105 	BC_UNUSED(scale);
2106 
2107 	val = bc_num_intop(a, b, c);
2108 
2109 	// Just truncate or extend as appropriate.
2110 	if (val < c->scale) bc_num_truncate(c, c->scale - val);
2111 	else if (val > c->scale) bc_num_extend(c, val - c->scale);
2112 }
2113 
2114 /**
2115  * Implements the left shift operator. This is a BcNumBinOp function.
2116  */
2117 static void
2118 bc_num_left(BcNum* a, BcNum* b, BcNum* restrict c, size_t scale)
2119 {
2120 	BcBigDig val;
2121 
2122 	BC_UNUSED(scale);
2123 
2124 	val = bc_num_intop(a, b, c);
2125 
2126 	bc_num_shiftLeft(c, (size_t) val);
2127 }
2128 
2129 /**
2130  * Implements the right shift operator. This is a BcNumBinOp function.
2131  */
2132 static void
2133 bc_num_right(BcNum* a, BcNum* b, BcNum* restrict c, size_t scale)
2134 {
2135 	BcBigDig val;
2136 
2137 	BC_UNUSED(scale);
2138 
2139 	val = bc_num_intop(a, b, c);
2140 
2141 	if (BC_NUM_ZERO(c)) return;
2142 
2143 	bc_num_shiftRight(c, (size_t) val);
2144 }
2145 #endif // BC_ENABLE_EXTRA_MATH
2146 
2147 /**
2148  * Prepares for, and calls, a binary operator function. This is probably the
2149  * most important function in the entire file because it establishes assumptions
2150  * that make the rest of the code so easy. Those assumptions include:
2151  *
2152  * - a is not the same pointer as c.
2153  * - b is not the same pointer as c.
2154  * - there is enough room in c for the result.
2155  *
2156  * Without these, this whole function would basically have to be duplicated for
2157  * *all* binary operators.
2158  *
2159  * @param a      The first operand.
2160  * @param b      The second operand.
2161  * @param c      The return parameter.
2162  * @param scale  The current scale.
2163  * @param req    The number of limbs needed to fit the result.
2164  */
2165 static void
2166 bc_num_binary(BcNum* a, BcNum* b, BcNum* c, size_t scale, BcNumBinOp op,
2167               size_t req)
2168 {
2169 	BcNum* ptr_a;
2170 	BcNum* ptr_b;
2171 	BcNum num2;
2172 	bool init = false;
2173 
2174 	assert(a != NULL && b != NULL && c != NULL && op != NULL);
2175 
2176 	assert(BC_NUM_RDX_VALID(a));
2177 	assert(BC_NUM_RDX_VALID(b));
2178 
2179 	BC_SIG_LOCK;
2180 
2181 	// Reallocate if c == a.
2182 	if (c == a)
2183 	{
2184 		ptr_a = &num2;
2185 
2186 		// NOLINTNEXTLINE
2187 		memcpy(ptr_a, c, sizeof(BcNum));
2188 		init = true;
2189 	}
2190 	else
2191 	{
2192 		ptr_a = a;
2193 	}
2194 
2195 	// Also reallocate if c == b.
2196 	if (c == b)
2197 	{
2198 		ptr_b = &num2;
2199 
2200 		if (c != a)
2201 		{
2202 			// NOLINTNEXTLINE
2203 			memcpy(ptr_b, c, sizeof(BcNum));
2204 			init = true;
2205 		}
2206 	}
2207 	else
2208 	{
2209 		ptr_b = b;
2210 	}
2211 
2212 	// Actually reallocate. If we don't reallocate, we want to expand at the
2213 	// very least.
2214 	if (init)
2215 	{
2216 		bc_num_init(c, req);
2217 
2218 		// Must prepare for cleanup. We want this here so that locals that got
2219 		// set stay set since a longjmp() is not guaranteed to preserve locals.
2220 		BC_SETJMP_LOCKED(err);
2221 		BC_SIG_UNLOCK;
2222 	}
2223 	else
2224 	{
2225 		BC_SIG_UNLOCK;
2226 		bc_num_expand(c, req);
2227 	}
2228 
2229 	// It is okay for a and b to be the same. If a binary operator function does
2230 	// need them to be different, the binary operator function is responsible
2231 	// for that.
2232 
2233 	// Call the actual binary operator function.
2234 	op(ptr_a, ptr_b, c, scale);
2235 
2236 	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
2237 	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
2238 	assert(BC_NUM_RDX_VALID(c));
2239 	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
2240 
2241 err:
2242 	// Cleanup only needed if we initialized c to a new number.
2243 	if (init)
2244 	{
2245 		BC_SIG_MAYLOCK;
2246 		bc_num_free(&num2);
2247 		BC_LONGJMP_CONT;
2248 	}
2249 }
2250 
2251 /**
2252  * Tests a number string for validity. This function has a history; I originally
2253  * wrote it because I did not trust my parser. Over time, however, I came to
2254  * trust it, so I was able to relegate this function to debug builds only, and I
2255  * used it in assert()'s. But then I created the library, and well, I can't
2256  * trust users, so I reused this for yelling at users.
2257  * @param val  The string to check to see if it's a valid number string.
2258  * @return     True if the string is a valid number string, false otherwise.
2259  */
2260 bool
2261 bc_num_strValid(const char* restrict val)
2262 {
2263 	bool radix = false;
2264 	size_t i, len = strlen(val);
2265 
2266 	// Notice that I don't check if there is a negative sign. That is not part
2267 	// of a valid number, except in the library. The library-specific code takes
2268 	// care of that part.
2269 
2270 	// Nothing in the string is okay.
2271 	if (!len) return true;
2272 
2273 	// Loop through the characters.
2274 	for (i = 0; i < len; ++i)
2275 	{
2276 		BcDig c = val[i];
2277 
2278 		// If we have found a radix point...
2279 		if (c == '.')
2280 		{
2281 			// We don't allow two radices.
2282 			if (radix) return false;
2283 
2284 			radix = true;
2285 			continue;
2286 		}
2287 
2288 		// We only allow digits and uppercase letters.
2289 		if (!(isdigit(c) || isupper(c))) return false;
2290 	}
2291 
2292 	return true;
2293 }
2294 
2295 /**
2296  * Parses one character and returns the digit that corresponds to that
2297  * character according to the base.
2298  * @param c     The character to parse.
2299  * @param base  The base.
2300  * @return      The character as a digit.
2301  */
2302 static BcBigDig
2303 bc_num_parseChar(char c, size_t base)
2304 {
2305 	assert(isupper(c) || isdigit(c));
2306 
2307 	// If a letter...
2308 	if (isupper(c))
2309 	{
2310 		// This returns the digit that directly corresponds with the letter.
2311 		c = BC_NUM_NUM_LETTER(c);
2312 
2313 		// If the digit is greater than the base, we clamp.
2314 		c = ((size_t) c) >= base ? (char) base - 1 : c;
2315 	}
2316 	// Straight convert the digit to a number.
2317 	else c -= '0';
2318 
2319 	return (BcBigDig) (uchar) c;
2320 }
2321 
2322 /**
2323  * Parses a string as a decimal number. This is separate because it's going to
2324  * be the most used, and it can be heavily optimized for decimal only.
2325  * @param n    The number to parse into and return. Must be preallocated.
2326  * @param val  The string to parse.
2327  */
2328 static void
2329 bc_num_parseDecimal(BcNum* restrict n, const char* restrict val)
2330 {
2331 	size_t len, i, temp, mod;
2332 	const char* ptr;
2333 	bool zero = true, rdx;
2334 
2335 	// Eat leading zeroes.
2336 	for (i = 0; val[i] == '0'; ++i)
2337 	{
2338 		continue;
2339 	}
2340 
2341 	val += i;
2342 	assert(!val[0] || isalnum(val[0]) || val[0] == '.');
2343 
2344 	// All 0's. We can just return, since this procedure expects a virgin
2345 	// (already 0) BcNum.
2346 	if (!val[0]) return;
2347 
2348 	// The length of the string is the length of the number, except it might be
2349 	// one bigger because of a decimal point.
2350 	len = strlen(val);
2351 
2352 	// Find the location of the decimal point.
2353 	ptr = strchr(val, '.');
2354 	rdx = (ptr != NULL);
2355 
2356 	// We eat leading zeroes again. These leading zeroes are different because
2357 	// they will come after the decimal point if they exist, and since that's
2358 	// the case, they must be preserved.
2359 	for (i = 0; i < len && (zero = (val[i] == '0' || val[i] == '.')); ++i)
2360 	{
2361 		continue;
2362 	}
2363 
2364 	// Set the scale of the number based on the location of the decimal point.
2365 	// The casts to uintptr_t is to ensure that bc does not hit undefined
2366 	// behavior when doing math on the values.
2367 	n->scale = (size_t) (rdx *
2368 	                     (((uintptr_t) (val + len)) - (((uintptr_t) ptr) + 1)));
2369 
2370 	// Set rdx.
2371 	BC_NUM_RDX_SET(n, BC_NUM_RDX(n->scale));
2372 
2373 	// Calculate length. First, the length of the integer, then the number of
2374 	// digits in the last limb, then the length.
2375 	i = len - (ptr == val ? 0 : i) - rdx;
2376 	temp = BC_NUM_ROUND_POW(i);
2377 	mod = n->scale % BC_BASE_DIGS;
2378 	i = mod ? BC_BASE_DIGS - mod : 0;
2379 	n->len = ((temp + i) / BC_BASE_DIGS);
2380 
2381 	// Expand and zero.
2382 	bc_num_expand(n, n->len);
2383 	// NOLINTNEXTLINE
2384 	memset(n->num, 0, BC_NUM_SIZE(n->len));
2385 
2386 	if (zero)
2387 	{
2388 		// I think I can set rdx directly to zero here because n should be a
2389 		// new number with sign set to false.
2390 		n->len = n->rdx = 0;
2391 	}
2392 	else
2393 	{
2394 		// There is actually stuff to parse if we make it here. Yay...
2395 		BcBigDig exp, pow;
2396 
2397 		assert(i <= BC_NUM_BIGDIG_MAX);
2398 
2399 		// The exponent and power.
2400 		exp = (BcBigDig) i;
2401 		pow = bc_num_pow10[exp];
2402 
2403 		// Parse loop. We parse backwards because numbers are stored little
2404 		// endian.
2405 		for (i = len - 1; i < len; --i, ++exp)
2406 		{
2407 			char c = val[i];
2408 
2409 			// Skip the decimal point.
2410 			if (c == '.') exp -= 1;
2411 			else
2412 			{
2413 				// The index of the limb.
2414 				size_t idx = exp / BC_BASE_DIGS;
2415 
2416 				// Clamp for the base.
2417 				if (isupper(c)) c = '9';
2418 
2419 				// Add the digit to the limb.
2420 				n->num[idx] += (((BcBigDig) c) - '0') * pow;
2421 
2422 				// Adjust the power and exponent.
2423 				if ((exp + 1) % BC_BASE_DIGS == 0) pow = 1;
2424 				else pow *= BC_BASE;
2425 			}
2426 		}
2427 	}
2428 }
2429 
2430 /**
2431  * Parse a number in any base (besides decimal).
2432  * @param n     The number to parse into and return. Must be preallocated.
2433  * @param val   The string to parse.
2434  * @param base  The base to parse as.
2435  */
2436 static void
2437 bc_num_parseBase(BcNum* restrict n, const char* restrict val, BcBigDig base)
2438 {
2439 	BcNum temp, mult1, mult2, result1, result2;
2440 	BcNum* m1;
2441 	BcNum* m2;
2442 	BcNum* ptr;
2443 	char c = 0;
2444 	bool zero = true;
2445 	BcBigDig v;
2446 	size_t i, digs, len = strlen(val);
2447 
2448 	// If zero, just return because the number should be virgin (already 0).
2449 	for (i = 0; zero && i < len; ++i)
2450 	{
2451 		zero = (val[i] == '.' || val[i] == '0');
2452 	}
2453 	if (zero) return;
2454 
2455 	BC_SIG_LOCK;
2456 
2457 	bc_num_init(&temp, BC_NUM_BIGDIG_LOG10);
2458 	bc_num_init(&mult1, BC_NUM_BIGDIG_LOG10);
2459 
2460 	BC_SETJMP_LOCKED(int_err);
2461 
2462 	BC_SIG_UNLOCK;
2463 
2464 	// We split parsing into parsing the integer and parsing the fractional
2465 	// part.
2466 
2467 	// Parse the integer part. This is the easy part because we just multiply
2468 	// the number by the base, then add the digit.
2469 	for (i = 0; i < len && (c = val[i]) && c != '.'; ++i)
2470 	{
2471 		// Convert the character to a digit.
2472 		v = bc_num_parseChar(c, base);
2473 
2474 		// Multiply the number.
2475 		bc_num_mulArray(n, base, &mult1);
2476 
2477 		// Convert the digit to a number and add.
2478 		bc_num_bigdig2num(&temp, v);
2479 		bc_num_add(&mult1, &temp, n, 0);
2480 	}
2481 
2482 	// If this condition is true, then we are done. We still need to do cleanup
2483 	// though.
2484 	if (i == len && !val[i]) goto int_err;
2485 
2486 	// If we get here, we *must* be at the radix point.
2487 	assert(val[i] == '.');
2488 
2489 	BC_SIG_LOCK;
2490 
2491 	// Unset the jump to reset in for these new initializations.
2492 	BC_UNSETJMP;
2493 
2494 	bc_num_init(&mult2, BC_NUM_BIGDIG_LOG10);
2495 	bc_num_init(&result1, BC_NUM_DEF_SIZE);
2496 	bc_num_init(&result2, BC_NUM_DEF_SIZE);
2497 	bc_num_one(&mult1);
2498 
2499 	BC_SETJMP_LOCKED(err);
2500 
2501 	BC_SIG_UNLOCK;
2502 
2503 	// Pointers for easy switching.
2504 	m1 = &mult1;
2505 	m2 = &mult2;
2506 
2507 	// Parse the fractional part. This is the hard part.
2508 	for (i += 1, digs = 0; i < len && (c = val[i]); ++i, ++digs)
2509 	{
2510 		size_t rdx;
2511 
2512 		// Convert the character to a digit.
2513 		v = bc_num_parseChar(c, base);
2514 
2515 		// We keep growing result2 according to the base because the more digits
2516 		// after the radix, the more significant the digits close to the radix
2517 		// should be.
2518 		bc_num_mulArray(&result1, base, &result2);
2519 
2520 		// Convert the digit to a number.
2521 		bc_num_bigdig2num(&temp, v);
2522 
2523 		// Add the digit into the fraction part.
2524 		bc_num_add(&result2, &temp, &result1, 0);
2525 
2526 		// Keep growing m1 and m2 for use after the loop.
2527 		bc_num_mulArray(m1, base, m2);
2528 
2529 		rdx = BC_NUM_RDX_VAL(m2);
2530 
2531 		if (m2->len < rdx) m2->len = rdx;
2532 
2533 		// Switch.
2534 		ptr = m1;
2535 		m1 = m2;
2536 		m2 = ptr;
2537 	}
2538 
2539 	// This one cannot be a divide by 0 because mult starts out at 1, then is
2540 	// multiplied by base, and base cannot be 0, so mult cannot be 0. And this
2541 	// is the reason we keep growing m1 and m2; this division is what converts
2542 	// the parsed fractional part from an integer to a fractional part.
2543 	bc_num_div(&result1, m1, &result2, digs * 2);
2544 
2545 	// Pretruncate.
2546 	bc_num_truncate(&result2, digs);
2547 
2548 	// The final add of the integer part to the fractional part.
2549 	bc_num_add(n, &result2, n, digs);
2550 
2551 	// Basic cleanup.
2552 	if (BC_NUM_NONZERO(n))
2553 	{
2554 		if (n->scale < digs) bc_num_extend(n, digs - n->scale);
2555 	}
2556 	else bc_num_zero(n);
2557 
2558 err:
2559 	BC_SIG_MAYLOCK;
2560 	bc_num_free(&result2);
2561 	bc_num_free(&result1);
2562 	bc_num_free(&mult2);
2563 int_err:
2564 	BC_SIG_MAYLOCK;
2565 	bc_num_free(&mult1);
2566 	bc_num_free(&temp);
2567 	BC_LONGJMP_CONT;
2568 }
2569 
2570 /**
2571  * Prints a backslash+newline combo if the number of characters needs it. This
2572  * is really a convenience function.
2573  */
2574 static inline void
2575 bc_num_printNewline(void)
2576 {
2577 #if !BC_ENABLE_LIBRARY
2578 	if (vm.nchars >= vm.line_len - 1 && vm.line_len)
2579 	{
2580 		bc_vm_putchar('\\', bc_flush_none);
2581 		bc_vm_putchar('\n', bc_flush_err);
2582 	}
2583 #endif // !BC_ENABLE_LIBRARY
2584 }
2585 
2586 /**
2587  * Prints a character after a backslash+newline, if needed.
2588  * @param c       The character to print.
2589  * @param bslash  Whether to print a backslash+newline.
2590  */
2591 static void
2592 bc_num_putchar(int c, bool bslash)
2593 {
2594 	if (c != '\n' && bslash) bc_num_printNewline();
2595 	bc_vm_putchar(c, bc_flush_save);
2596 }
2597 
2598 #if !BC_ENABLE_LIBRARY
2599 
2600 /**
2601  * Prints a character for a number's digit. This is for printing for dc's P
2602  * command. This function does not need to worry about radix points. This is a
2603  * BcNumDigitOp.
2604  * @param n       The "digit" to print.
2605  * @param len     The "length" of the digit, or number of characters that will
2606  *                need to be printed for the digit.
2607  * @param rdx     True if a decimal (radix) point should be printed.
2608  * @param bslash  True if a backslash+newline should be printed if the character
2609  *                limit for the line is reached, false otherwise.
2610  */
2611 static void
2612 bc_num_printChar(size_t n, size_t len, bool rdx, bool bslash)
2613 {
2614 	BC_UNUSED(rdx);
2615 	BC_UNUSED(len);
2616 	BC_UNUSED(bslash);
2617 	assert(len == 1);
2618 	bc_vm_putchar((uchar) n, bc_flush_save);
2619 }
2620 
2621 #endif // !BC_ENABLE_LIBRARY
2622 
2623 /**
2624  * Prints a series of characters for large bases. This is for printing in bases
2625  * above hexadecimal. This is a BcNumDigitOp.
2626  * @param n       The "digit" to print.
2627  * @param len     The "length" of the digit, or number of characters that will
2628  *                need to be printed for the digit.
2629  * @param rdx     True if a decimal (radix) point should be printed.
2630  * @param bslash  True if a backslash+newline should be printed if the character
2631  *                limit for the line is reached, false otherwise.
2632  */
2633 static void
2634 bc_num_printDigits(size_t n, size_t len, bool rdx, bool bslash)
2635 {
2636 	size_t exp, pow;
2637 
2638 	// If needed, print the radix; otherwise, print a space to separate digits.
2639 	bc_num_putchar(rdx ? '.' : ' ', true);
2640 
2641 	// Calculate the exponent and power.
2642 	for (exp = 0, pow = 1; exp < len - 1; ++exp, pow *= BC_BASE)
2643 	{
2644 		continue;
2645 	}
2646 
2647 	// Print each character individually.
2648 	for (exp = 0; exp < len; pow /= BC_BASE, ++exp)
2649 	{
2650 		// The individual subdigit.
2651 		size_t dig = n / pow;
2652 
2653 		// Take the subdigit away.
2654 		n -= dig * pow;
2655 
2656 		// Print the subdigit.
2657 		bc_num_putchar(((uchar) dig) + '0', bslash || exp != len - 1);
2658 	}
2659 }
2660 
2661 /**
2662  * Prints a character for a number's digit. This is for printing in bases for
2663  * hexadecimal and below because they always print only one character at a time.
2664  * This is a BcNumDigitOp.
2665  * @param n       The "digit" to print.
2666  * @param len     The "length" of the digit, or number of characters that will
2667  *                need to be printed for the digit.
2668  * @param rdx     True if a decimal (radix) point should be printed.
2669  * @param bslash  True if a backslash+newline should be printed if the character
2670  *                limit for the line is reached, false otherwise.
2671  */
2672 static void
2673 bc_num_printHex(size_t n, size_t len, bool rdx, bool bslash)
2674 {
2675 	BC_UNUSED(len);
2676 	BC_UNUSED(bslash);
2677 
2678 	assert(len == 1);
2679 
2680 	if (rdx) bc_num_putchar('.', true);
2681 
2682 	bc_num_putchar(bc_num_hex_digits[n], bslash);
2683 }
2684 
2685 /**
2686  * Prints a decimal number. This is specially written for optimization since
2687  * this will be used the most and because bc's numbers are already in decimal.
2688  * @param n        The number to print.
2689  * @param newline  Whether to print backslash+newlines on long enough lines.
2690  */
2691 static void
2692 bc_num_printDecimal(const BcNum* restrict n, bool newline)
2693 {
2694 	size_t i, j, rdx = BC_NUM_RDX_VAL(n);
2695 	bool zero = true;
2696 	size_t buffer[BC_BASE_DIGS];
2697 
2698 	// Print loop.
2699 	for (i = n->len - 1; i < n->len; --i)
2700 	{
2701 		BcDig n9 = n->num[i];
2702 		size_t temp;
2703 		bool irdx = (i == rdx - 1);
2704 
2705 		// Calculate the number of digits in the limb.
2706 		zero = (zero & !irdx);
2707 		temp = n->scale % BC_BASE_DIGS;
2708 		temp = i || !temp ? 0 : BC_BASE_DIGS - temp;
2709 
2710 		// NOLINTNEXTLINE
2711 		memset(buffer, 0, BC_BASE_DIGS * sizeof(size_t));
2712 
2713 		// Fill the buffer with individual digits.
2714 		for (j = 0; n9 && j < BC_BASE_DIGS; ++j)
2715 		{
2716 			buffer[j] = ((size_t) n9) % BC_BASE;
2717 			n9 /= BC_BASE;
2718 		}
2719 
2720 		// Print the digits in the buffer.
2721 		for (j = BC_BASE_DIGS - 1; j < BC_BASE_DIGS && j >= temp; --j)
2722 		{
2723 			// Figure out whether to print the decimal point.
2724 			bool print_rdx = (irdx & (j == BC_BASE_DIGS - 1));
2725 
2726 			// The zero variable helps us skip leading zero digits in the limb.
2727 			zero = (zero && buffer[j] == 0);
2728 
2729 			if (!zero)
2730 			{
2731 				// While the first three arguments should be self-explanatory,
2732 				// the last needs explaining. I don't want to print a newline
2733 				// when the last digit to be printed could take the place of the
2734 				// backslash rather than being pushed, as a single character, to
2735 				// the next line. That's what that last argument does for bc.
2736 				bc_num_printHex(buffer[j], 1, print_rdx,
2737 				                !newline || (j > temp || i != 0));
2738 			}
2739 		}
2740 	}
2741 }
2742 
2743 #if BC_ENABLE_EXTRA_MATH
2744 
2745 /**
2746  * Prints a number in scientific or engineering format. When doing this, we are
2747  * always printing in decimal.
2748  * @param n        The number to print.
2749  * @param eng      True if we are in engineering mode.
2750  * @param newline  Whether to print backslash+newlines on long enough lines.
2751  */
2752 static void
2753 bc_num_printExponent(const BcNum* restrict n, bool eng, bool newline)
2754 {
2755 	size_t places, mod, nrdx = BC_NUM_RDX_VAL(n);
2756 	bool neg = (n->len <= nrdx);
2757 	BcNum temp, exp;
2758 	BcDig digs[BC_NUM_BIGDIG_LOG10];
2759 
2760 	BC_SIG_LOCK;
2761 
2762 	bc_num_createCopy(&temp, n);
2763 
2764 	BC_SETJMP_LOCKED(exit);
2765 
2766 	BC_SIG_UNLOCK;
2767 
2768 	// We need to calculate the exponents, and they change based on whether the
2769 	// number is all fractional or not, obviously.
2770 	if (neg)
2771 	{
2772 		// Figure out how many limbs after the decimal point is zero.
2773 		size_t i, idx = bc_num_nonZeroLen(n) - 1;
2774 
2775 		places = 1;
2776 
2777 		// Figure out how much in the last limb is zero.
2778 		for (i = BC_BASE_DIGS - 1; i < BC_BASE_DIGS; --i)
2779 		{
2780 			if (bc_num_pow10[i] > (BcBigDig) n->num[idx]) places += 1;
2781 			else break;
2782 		}
2783 
2784 		// Calculate the combination of zero limbs and zero digits in the last
2785 		// limb.
2786 		places += (nrdx - (idx + 1)) * BC_BASE_DIGS;
2787 		mod = places % 3;
2788 
2789 		// Calculate places if we are in engineering mode.
2790 		if (eng && mod != 0) places += 3 - mod;
2791 
2792 		// Shift the temp to the right place.
2793 		bc_num_shiftLeft(&temp, places);
2794 	}
2795 	else
2796 	{
2797 		// This is the number of digits that we are supposed to put behind the
2798 		// decimal point.
2799 		places = bc_num_intDigits(n) - 1;
2800 
2801 		// Calculate the true number based on whether engineering mode is
2802 		// activated.
2803 		mod = places % 3;
2804 		if (eng && mod != 0) places -= 3 - (3 - mod);
2805 
2806 		// Shift the temp to the right place.
2807 		bc_num_shiftRight(&temp, places);
2808 	}
2809 
2810 	// Print the shifted number.
2811 	bc_num_printDecimal(&temp, newline);
2812 
2813 	// Print the e.
2814 	bc_num_putchar('e', !newline);
2815 
2816 	// Need to explicitly print a zero exponent.
2817 	if (!places)
2818 	{
2819 		bc_num_printHex(0, 1, false, !newline);
2820 		goto exit;
2821 	}
2822 
2823 	// Need to print sign for the exponent.
2824 	if (neg) bc_num_putchar('-', true);
2825 
2826 	// Create a temporary for the exponent...
2827 	bc_num_setup(&exp, digs, BC_NUM_BIGDIG_LOG10);
2828 	bc_num_bigdig2num(&exp, (BcBigDig) places);
2829 
2830 	/// ..and print it.
2831 	bc_num_printDecimal(&exp, newline);
2832 
2833 exit:
2834 	BC_SIG_MAYLOCK;
2835 	bc_num_free(&temp);
2836 	BC_LONGJMP_CONT;
2837 }
2838 #endif // BC_ENABLE_EXTRA_MATH
2839 
2840 /**
2841  * Converts a number from limbs with base BC_BASE_POW to base @a pow, where
2842  * @a pow is obase^N.
2843  * @param n    The number to convert.
2844  * @param rem  BC_BASE_POW - @a pow.
2845  * @param pow  The power of obase we will convert the number to.
2846  * @param idx  The index of the number to start converting at. Doing the
2847  *             conversion is O(n^2); we have to sweep through starting at the
2848  *             least significant limb
2849  */
2850 static void
2851 bc_num_printFixup(BcNum* restrict n, BcBigDig rem, BcBigDig pow, size_t idx)
2852 {
2853 	size_t i, len = n->len - idx;
2854 	BcBigDig acc;
2855 	BcDig* a = n->num + idx;
2856 
2857 	// Ignore if there's just one limb left. This is the part that requires the
2858 	// extra loop after the one calling this function in bc_num_printPrepare().
2859 	if (len < 2) return;
2860 
2861 	// Loop through the remaining limbs and convert. We start at the second limb
2862 	// because we pull the value from the previous one as well.
2863 	for (i = len - 1; i > 0; --i)
2864 	{
2865 		// Get the limb and add it to the previous, along with multiplying by
2866 		// the remainder because that's the proper overflow. "acc" means
2867 		// "accumulator," by the way.
2868 		acc = ((BcBigDig) a[i]) * rem + ((BcBigDig) a[i - 1]);
2869 
2870 		// Store a value in base pow in the previous limb.
2871 		a[i - 1] = (BcDig) (acc % pow);
2872 
2873 		// Divide by the base and accumulate the remaining value in the limb.
2874 		acc /= pow;
2875 		acc += (BcBigDig) a[i];
2876 
2877 		// If the accumulator is greater than the base...
2878 		if (acc >= BC_BASE_POW)
2879 		{
2880 			// Do we need to grow?
2881 			if (i == len - 1)
2882 			{
2883 				// Grow.
2884 				len = bc_vm_growSize(len, 1);
2885 				bc_num_expand(n, bc_vm_growSize(len, idx));
2886 
2887 				// Update the pointer because it may have moved.
2888 				a = n->num + idx;
2889 
2890 				// Zero out the last limb.
2891 				a[len - 1] = 0;
2892 			}
2893 
2894 			// Overflow into the next limb since we are over the base.
2895 			a[i + 1] += acc / BC_BASE_POW;
2896 			acc %= BC_BASE_POW;
2897 		}
2898 
2899 		assert(acc < BC_BASE_POW);
2900 
2901 		// Set the limb.
2902 		a[i] = (BcDig) acc;
2903 	}
2904 
2905 	// We may have grown the number, so adjust the length.
2906 	n->len = len + idx;
2907 }
2908 
2909 /**
2910  * Prepares a number for printing in a base that is not a divisor of
2911  * BC_BASE_POW. This basically converts the number from having limbs of base
2912  * BC_BASE_POW to limbs of pow, where pow is obase^N.
2913  * @param n    The number to prepare for printing.
2914  * @param rem  The remainder of BC_BASE_POW when divided by a power of the base.
2915  * @param pow  The power of the base.
2916  */
2917 static void
2918 bc_num_printPrepare(BcNum* restrict n, BcBigDig rem, BcBigDig pow)
2919 {
2920 	size_t i;
2921 
2922 	// Loop from the least significant limb to the most significant limb and
2923 	// convert limbs in each pass.
2924 	for (i = 0; i < n->len; ++i)
2925 	{
2926 		bc_num_printFixup(n, rem, pow, i);
2927 	}
2928 
2929 	// bc_num_printFixup() does not do everything it is supposed to, so we do
2930 	// the last bit of cleanup here. That cleanup is to ensure that each limb
2931 	// is less than pow and to expand the number to fit new limbs as necessary.
2932 	for (i = 0; i < n->len; ++i)
2933 	{
2934 		assert(pow == ((BcBigDig) ((BcDig) pow)));
2935 
2936 		// If the limb needs fixing...
2937 		if (n->num[i] >= (BcDig) pow)
2938 		{
2939 			// Do we need to grow?
2940 			if (i + 1 == n->len)
2941 			{
2942 				// Grow the number.
2943 				n->len = bc_vm_growSize(n->len, 1);
2944 				bc_num_expand(n, n->len);
2945 
2946 				// Without this, we might use uninitialized data.
2947 				n->num[i + 1] = 0;
2948 			}
2949 
2950 			assert(pow < BC_BASE_POW);
2951 
2952 			// Overflow into the next limb.
2953 			n->num[i + 1] += n->num[i] / ((BcDig) pow);
2954 			n->num[i] %= (BcDig) pow;
2955 		}
2956 	}
2957 }
2958 
2959 static void
2960 bc_num_printNum(BcNum* restrict n, BcBigDig base, size_t len,
2961                 BcNumDigitOp print, bool newline)
2962 {
2963 	BcVec stack;
2964 	BcNum intp, fracp1, fracp2, digit, flen1, flen2;
2965 	BcNum* n1;
2966 	BcNum* n2;
2967 	BcNum* temp;
2968 	BcBigDig dig = 0, acc, exp;
2969 	BcBigDig* ptr;
2970 	size_t i, j, nrdx, idigits;
2971 	bool radix;
2972 	BcDig digit_digs[BC_NUM_BIGDIG_LOG10 + 1];
2973 
2974 	assert(base > 1);
2975 
2976 	// Easy case. Even with scale, we just print this.
2977 	if (BC_NUM_ZERO(n))
2978 	{
2979 		print(0, len, false, !newline);
2980 		return;
2981 	}
2982 
2983 	// This function uses an algorithm that Stefan Esser <se@freebsd.org> came
2984 	// up with to print the integer part of a number. What it does is convert
2985 	// intp into a number of the specified base, but it does it directly,
2986 	// instead of just doing a series of divisions and printing the remainders
2987 	// in reverse order.
2988 	//
2989 	// Let me explain in a bit more detail:
2990 	//
2991 	// The algorithm takes the current least significant limb (after intp has
2992 	// been converted to an integer) and the next to least significant limb, and
2993 	// it converts the least significant limb into one of the specified base,
2994 	// putting any overflow into the next to least significant limb. It iterates
2995 	// through the whole number, from least significant to most significant,
2996 	// doing this conversion. At the end of that iteration, the least
2997 	// significant limb is converted, but the others are not, so it iterates
2998 	// again, starting at the next to least significant limb. It keeps doing
2999 	// that conversion, skipping one more limb than the last time, until all
3000 	// limbs have been converted. Then it prints them in reverse order.
3001 	//
3002 	// That is the gist of the algorithm. It leaves out several things, such as
3003 	// the fact that limbs are not always converted into the specified base, but
3004 	// into something close, basically a power of the specified base. In
3005 	// Stefan's words, "You could consider BcDigs to be of base 10^BC_BASE_DIGS
3006 	// in the normal case and obase^N for the largest value of N that satisfies
3007 	// obase^N <= 10^BC_BASE_DIGS. [This means that] the result is not in base
3008 	// "obase", but in base "obase^N", which happens to be printable as a number
3009 	// of base "obase" without consideration for neighbouring BcDigs." This fact
3010 	// is what necessitates the existence of the loop later in this function.
3011 	//
3012 	// The conversion happens in bc_num_printPrepare() where the outer loop
3013 	// happens and bc_num_printFixup() where the inner loop, or actual
3014 	// conversion, happens. In other words, bc_num_printPrepare() is where the
3015 	// loop that starts at the least significant limb and goes to the most
3016 	// significant limb. Then, on every iteration of its loop, it calls
3017 	// bc_num_printFixup(), which has the inner loop of actually converting
3018 	// the limbs it passes into limbs of base obase^N rather than base
3019 	// BC_BASE_POW.
3020 
3021 	nrdx = BC_NUM_RDX_VAL(n);
3022 
3023 	BC_SIG_LOCK;
3024 
3025 	// The stack is what allows us to reverse the digits for printing.
3026 	bc_vec_init(&stack, sizeof(BcBigDig), BC_DTOR_NONE);
3027 	bc_num_init(&fracp1, nrdx);
3028 
3029 	// intp will be the "integer part" of the number, so copy it.
3030 	bc_num_createCopy(&intp, n);
3031 
3032 	BC_SETJMP_LOCKED(err);
3033 
3034 	BC_SIG_UNLOCK;
3035 
3036 	// Make intp an integer.
3037 	bc_num_truncate(&intp, intp.scale);
3038 
3039 	// Get the fractional part out.
3040 	bc_num_sub(n, &intp, &fracp1, 0);
3041 
3042 	// If the base is not the same as the last base used for printing, we need
3043 	// to update the cached exponent and power. Yes, we cache the values of the
3044 	// exponent and power. That is to prevent us from calculating them every
3045 	// time because printing will probably happen multiple times on the same
3046 	// base.
3047 	if (base != vm.last_base)
3048 	{
3049 		vm.last_pow = 1;
3050 		vm.last_exp = 0;
3051 
3052 		// Calculate the exponent and power.
3053 		while (vm.last_pow * base <= BC_BASE_POW)
3054 		{
3055 			vm.last_pow *= base;
3056 			vm.last_exp += 1;
3057 		}
3058 
3059 		// Also, the remainder and base itself.
3060 		vm.last_rem = BC_BASE_POW - vm.last_pow;
3061 		vm.last_base = base;
3062 	}
3063 
3064 	exp = vm.last_exp;
3065 
3066 	// If vm.last_rem is 0, then the base we are printing in is a divisor of
3067 	// BC_BASE_POW, which is the easy case because it means that BC_BASE_POW is
3068 	// a power of obase, and no conversion is needed. If it *is* 0, then we have
3069 	// the hard case, and we have to prepare the number for the base.
3070 	if (vm.last_rem != 0) bc_num_printPrepare(&intp, vm.last_rem, vm.last_pow);
3071 
3072 	// After the conversion comes the surprisingly easy part. From here on out,
3073 	// this is basically naive code that I wrote, adjusted for the larger bases.
3074 
3075 	// Fill the stack of digits for the integer part.
3076 	for (i = 0; i < intp.len; ++i)
3077 	{
3078 		// Get the limb.
3079 		acc = (BcBigDig) intp.num[i];
3080 
3081 		// Turn the limb into digits of base obase.
3082 		for (j = 0; j < exp && (i < intp.len - 1 || acc != 0); ++j)
3083 		{
3084 			// This condition is true if we are not at the last digit.
3085 			if (j != exp - 1)
3086 			{
3087 				dig = acc % base;
3088 				acc /= base;
3089 			}
3090 			else
3091 			{
3092 				dig = acc;
3093 				acc = 0;
3094 			}
3095 
3096 			assert(dig < base);
3097 
3098 			// Push the digit onto the stack.
3099 			bc_vec_push(&stack, &dig);
3100 		}
3101 
3102 		assert(acc == 0);
3103 	}
3104 
3105 	// Go through the stack backwards and print each digit.
3106 	for (i = 0; i < stack.len; ++i)
3107 	{
3108 		ptr = bc_vec_item_rev(&stack, i);
3109 
3110 		assert(ptr != NULL);
3111 
3112 		// While the first three arguments should be self-explanatory, the last
3113 		// needs explaining. I don't want to print a newline when the last digit
3114 		// to be printed could take the place of the backslash rather than being
3115 		// pushed, as a single character, to the next line. That's what that
3116 		// last argument does for bc.
3117 		print(*ptr, len, false,
3118 		      !newline || (n->scale != 0 || i == stack.len - 1));
3119 	}
3120 
3121 	// We are done if there is no fractional part.
3122 	if (!n->scale) goto err;
3123 
3124 	BC_SIG_LOCK;
3125 
3126 	// Reset the jump because some locals are changing.
3127 	BC_UNSETJMP;
3128 
3129 	bc_num_init(&fracp2, nrdx);
3130 	bc_num_setup(&digit, digit_digs, sizeof(digit_digs) / sizeof(BcDig));
3131 	bc_num_init(&flen1, BC_NUM_BIGDIG_LOG10);
3132 	bc_num_init(&flen2, BC_NUM_BIGDIG_LOG10);
3133 
3134 	BC_SETJMP_LOCKED(frac_err);
3135 
3136 	BC_SIG_UNLOCK;
3137 
3138 	bc_num_one(&flen1);
3139 
3140 	radix = true;
3141 
3142 	// Pointers for easy switching.
3143 	n1 = &flen1;
3144 	n2 = &flen2;
3145 
3146 	fracp2.scale = n->scale;
3147 	BC_NUM_RDX_SET_NP(fracp2, BC_NUM_RDX(fracp2.scale));
3148 
3149 	// As long as we have not reached the scale of the number, keep printing.
3150 	while ((idigits = bc_num_intDigits(n1)) <= n->scale)
3151 	{
3152 		// These numbers will keep growing.
3153 		bc_num_expand(&fracp2, fracp1.len + 1);
3154 		bc_num_mulArray(&fracp1, base, &fracp2);
3155 
3156 		nrdx = BC_NUM_RDX_VAL_NP(fracp2);
3157 
3158 		// Ensure an invariant.
3159 		if (fracp2.len < nrdx) fracp2.len = nrdx;
3160 
3161 		// fracp is guaranteed to be non-negative and small enough.
3162 		dig = bc_num_bigdig2(&fracp2);
3163 
3164 		// Convert the digit to a number and subtract it from the number.
3165 		bc_num_bigdig2num(&digit, dig);
3166 		bc_num_sub(&fracp2, &digit, &fracp1, 0);
3167 
3168 		// While the first three arguments should be self-explanatory, the last
3169 		// needs explaining. I don't want to print a newline when the last digit
3170 		// to be printed could take the place of the backslash rather than being
3171 		// pushed, as a single character, to the next line. That's what that
3172 		// last argument does for bc.
3173 		print(dig, len, radix, !newline || idigits != n->scale);
3174 
3175 		// Update the multipliers.
3176 		bc_num_mulArray(n1, base, n2);
3177 
3178 		radix = false;
3179 
3180 		// Switch.
3181 		temp = n1;
3182 		n1 = n2;
3183 		n2 = temp;
3184 	}
3185 
3186 frac_err:
3187 	BC_SIG_MAYLOCK;
3188 	bc_num_free(&flen2);
3189 	bc_num_free(&flen1);
3190 	bc_num_free(&fracp2);
3191 err:
3192 	BC_SIG_MAYLOCK;
3193 	bc_num_free(&fracp1);
3194 	bc_num_free(&intp);
3195 	bc_vec_free(&stack);
3196 	BC_LONGJMP_CONT;
3197 }
3198 
3199 /**
3200  * Prints a number in the specified base, or rather, figures out which function
3201  * to call to print the number in the specified base and calls it.
3202  * @param n        The number to print.
3203  * @param base     The base to print in.
3204  * @param newline  Whether to print backslash+newlines on long enough lines.
3205  */
3206 static void
3207 bc_num_printBase(BcNum* restrict n, BcBigDig base, bool newline)
3208 {
3209 	size_t width;
3210 	BcNumDigitOp print;
3211 	bool neg = BC_NUM_NEG(n);
3212 
3213 	// Clear the sign because it makes the actual printing easier when we have
3214 	// to do math.
3215 	BC_NUM_NEG_CLR(n);
3216 
3217 	// Bases at hexadecimal and below are printed as one character, larger bases
3218 	// are printed as a series of digits separated by spaces.
3219 	if (base <= BC_NUM_MAX_POSIX_IBASE)
3220 	{
3221 		width = 1;
3222 		print = bc_num_printHex;
3223 	}
3224 	else
3225 	{
3226 		assert(base <= BC_BASE_POW);
3227 		width = bc_num_log10(base - 1);
3228 		print = bc_num_printDigits;
3229 	}
3230 
3231 	// Print.
3232 	bc_num_printNum(n, base, width, print, newline);
3233 
3234 	// Reset the sign.
3235 	n->rdx = BC_NUM_NEG_VAL(n, neg);
3236 }
3237 
3238 #if !BC_ENABLE_LIBRARY
3239 
3240 void
3241 bc_num_stream(BcNum* restrict n)
3242 {
3243 	bc_num_printNum(n, BC_NUM_STREAM_BASE, 1, bc_num_printChar, false);
3244 }
3245 
3246 #endif // !BC_ENABLE_LIBRARY
3247 
3248 void
3249 bc_num_setup(BcNum* restrict n, BcDig* restrict num, size_t cap)
3250 {
3251 	assert(n != NULL);
3252 	n->num = num;
3253 	n->cap = cap;
3254 	bc_num_zero(n);
3255 }
3256 
3257 void
3258 bc_num_init(BcNum* restrict n, size_t req)
3259 {
3260 	BcDig* num;
3261 
3262 	BC_SIG_ASSERT_LOCKED;
3263 
3264 	assert(n != NULL);
3265 
3266 	// BC_NUM_DEF_SIZE is set to be about the smallest allocation size that
3267 	// malloc() returns in practice, so just use it.
3268 	req = req >= BC_NUM_DEF_SIZE ? req : BC_NUM_DEF_SIZE;
3269 
3270 	// If we can't use a temp, allocate.
3271 	if (req != BC_NUM_DEF_SIZE || (num = bc_vm_takeTemp()) == NULL)
3272 	{
3273 		num = bc_vm_malloc(BC_NUM_SIZE(req));
3274 	}
3275 
3276 	bc_num_setup(n, num, req);
3277 }
3278 
3279 void
3280 bc_num_clear(BcNum* restrict n)
3281 {
3282 	n->num = NULL;
3283 	n->cap = 0;
3284 }
3285 
3286 void
3287 bc_num_free(void* num)
3288 {
3289 	BcNum* n = (BcNum*) num;
3290 
3291 	BC_SIG_ASSERT_LOCKED;
3292 
3293 	assert(n != NULL);
3294 
3295 	if (n->cap == BC_NUM_DEF_SIZE) bc_vm_addTemp(n->num);
3296 	else free(n->num);
3297 }
3298 
3299 void
3300 bc_num_copy(BcNum* d, const BcNum* s)
3301 {
3302 	assert(d != NULL && s != NULL);
3303 
3304 	if (d == s) return;
3305 
3306 	bc_num_expand(d, s->len);
3307 	d->len = s->len;
3308 
3309 	// I can just copy directly here because the sign *and* rdx will be
3310 	// properly preserved.
3311 	d->rdx = s->rdx;
3312 	d->scale = s->scale;
3313 	// NOLINTNEXTLINE
3314 	memcpy(d->num, s->num, BC_NUM_SIZE(d->len));
3315 }
3316 
3317 void
3318 bc_num_createCopy(BcNum* d, const BcNum* s)
3319 {
3320 	BC_SIG_ASSERT_LOCKED;
3321 	bc_num_init(d, s->len);
3322 	bc_num_copy(d, s);
3323 }
3324 
3325 void
3326 bc_num_createFromBigdig(BcNum* restrict n, BcBigDig val)
3327 {
3328 	BC_SIG_ASSERT_LOCKED;
3329 	bc_num_init(n, BC_NUM_BIGDIG_LOG10);
3330 	bc_num_bigdig2num(n, val);
3331 }
3332 
3333 size_t
3334 bc_num_scale(const BcNum* restrict n)
3335 {
3336 	return n->scale;
3337 }
3338 
3339 size_t
3340 bc_num_len(const BcNum* restrict n)
3341 {
3342 	size_t len = n->len;
3343 
3344 	// Always return at least 1.
3345 	if (BC_NUM_ZERO(n)) return n->scale ? n->scale : 1;
3346 
3347 	// If this is true, there is no integer portion of the number.
3348 	if (BC_NUM_RDX_VAL(n) == len)
3349 	{
3350 		// We have to take into account the fact that some of the digits right
3351 		// after the decimal could be zero. If that is the case, we need to
3352 		// ignore them until we hit the first non-zero digit.
3353 
3354 		size_t zero, scale;
3355 
3356 		// The number of limbs with non-zero digits.
3357 		len = bc_num_nonZeroLen(n);
3358 
3359 		// Get the number of digits in the last limb.
3360 		scale = n->scale % BC_BASE_DIGS;
3361 		scale = scale ? scale : BC_BASE_DIGS;
3362 
3363 		// Get the number of zero digits.
3364 		zero = bc_num_zeroDigits(n->num + len - 1);
3365 
3366 		// Calculate the true length.
3367 		len = len * BC_BASE_DIGS - zero - (BC_BASE_DIGS - scale);
3368 	}
3369 	// Otherwise, count the number of int digits and return that plus the scale.
3370 	else len = bc_num_intDigits(n) + n->scale;
3371 
3372 	return len;
3373 }
3374 
3375 void
3376 bc_num_parse(BcNum* restrict n, const char* restrict val, BcBigDig base)
3377 {
3378 	assert(n != NULL && val != NULL && base);
3379 	assert(base >= BC_NUM_MIN_BASE && base <= vm.maxes[BC_PROG_GLOBALS_IBASE]);
3380 	assert(bc_num_strValid(val));
3381 
3382 	// A one character number is *always* parsed as though the base was the
3383 	// maximum allowed ibase, per the bc spec.
3384 	if (!val[1])
3385 	{
3386 		BcBigDig dig = bc_num_parseChar(val[0], BC_NUM_MAX_LBASE);
3387 		bc_num_bigdig2num(n, dig);
3388 	}
3389 	else if (base == BC_BASE) bc_num_parseDecimal(n, val);
3390 	else bc_num_parseBase(n, val, base);
3391 
3392 	assert(BC_NUM_RDX_VALID(n));
3393 }
3394 
3395 void
3396 bc_num_print(BcNum* restrict n, BcBigDig base, bool newline)
3397 {
3398 	assert(n != NULL);
3399 	assert(BC_ENABLE_EXTRA_MATH || base >= BC_NUM_MIN_BASE);
3400 
3401 	// We may need a newline, just to start.
3402 	bc_num_printNewline();
3403 
3404 	if (BC_NUM_NONZERO(n))
3405 	{
3406 		// Print the sign.
3407 		if (BC_NUM_NEG(n)) bc_num_putchar('-', true);
3408 
3409 		// Print the leading zero if necessary.
3410 		if (BC_Z && BC_NUM_RDX_VAL(n) == n->len)
3411 		{
3412 			bc_num_printHex(0, 1, false, !newline);
3413 		}
3414 	}
3415 
3416 	// Short-circuit 0.
3417 	if (BC_NUM_ZERO(n)) bc_num_printHex(0, 1, false, !newline);
3418 	else if (base == BC_BASE) bc_num_printDecimal(n, newline);
3419 #if BC_ENABLE_EXTRA_MATH
3420 	else if (base == 0 || base == 1)
3421 	{
3422 		bc_num_printExponent(n, base != 0, newline);
3423 	}
3424 #endif // BC_ENABLE_EXTRA_MATH
3425 	else bc_num_printBase(n, base, newline);
3426 
3427 	if (newline) bc_num_putchar('\n', false);
3428 }
3429 
3430 BcBigDig
3431 bc_num_bigdig2(const BcNum* restrict n)
3432 {
3433 	// This function returns no errors because it's guaranteed to succeed if
3434 	// its preconditions are met. Those preconditions include both n needs to
3435 	// be non-NULL, n being non-negative, and n being less than vm.max. If all
3436 	// of that is true, then we can just convert without worrying about negative
3437 	// errors or overflow.
3438 
3439 	BcBigDig r = 0;
3440 	size_t nrdx = BC_NUM_RDX_VAL(n);
3441 
3442 	assert(n != NULL);
3443 	assert(!BC_NUM_NEG(n));
3444 	assert(bc_num_cmp(n, &vm.max) < 0);
3445 	assert(n->len - nrdx <= 3);
3446 
3447 	// There is a small speed win from unrolling the loop here, and since it
3448 	// only adds 53 bytes, I decided that it was worth it.
3449 	switch (n->len - nrdx)
3450 	{
3451 		case 3:
3452 		{
3453 			r = (BcBigDig) n->num[nrdx + 2];
3454 
3455 			// Fallthrough.
3456 			BC_FALLTHROUGH
3457 		}
3458 
3459 		case 2:
3460 		{
3461 			r = r * BC_BASE_POW + (BcBigDig) n->num[nrdx + 1];
3462 
3463 			// Fallthrough.
3464 			BC_FALLTHROUGH
3465 		}
3466 
3467 		case 1:
3468 		{
3469 			r = r * BC_BASE_POW + (BcBigDig) n->num[nrdx];
3470 		}
3471 	}
3472 
3473 	return r;
3474 }
3475 
3476 BcBigDig
3477 bc_num_bigdig(const BcNum* restrict n)
3478 {
3479 	assert(n != NULL);
3480 
3481 	// This error checking is extremely important, and if you do not have a
3482 	// guarantee that converting a number will always succeed in a particular
3483 	// case, you *must* call this function to get these error checks. This
3484 	// includes all instances of numbers inputted by the user or calculated by
3485 	// the user. Otherwise, you can call the faster bc_num_bigdig2().
3486 	if (BC_ERR(BC_NUM_NEG(n))) bc_err(BC_ERR_MATH_NEGATIVE);
3487 	if (BC_ERR(bc_num_cmp(n, &vm.max) >= 0)) bc_err(BC_ERR_MATH_OVERFLOW);
3488 
3489 	return bc_num_bigdig2(n);
3490 }
3491 
3492 void
3493 bc_num_bigdig2num(BcNum* restrict n, BcBigDig val)
3494 {
3495 	BcDig* ptr;
3496 	size_t i;
3497 
3498 	assert(n != NULL);
3499 
3500 	bc_num_zero(n);
3501 
3502 	// Already 0.
3503 	if (!val) return;
3504 
3505 	// Expand first. This is the only way this function can fail, and it's a
3506 	// fatal error.
3507 	bc_num_expand(n, BC_NUM_BIGDIG_LOG10);
3508 
3509 	// The conversion is easy because numbers are laid out in little-endian
3510 	// order.
3511 	for (ptr = n->num, i = 0; val; ++i, val /= BC_BASE_POW)
3512 	{
3513 		ptr[i] = val % BC_BASE_POW;
3514 	}
3515 
3516 	n->len = i;
3517 }
3518 
3519 #if BC_ENABLE_EXTRA_MATH
3520 
3521 void
3522 bc_num_rng(const BcNum* restrict n, BcRNG* rng)
3523 {
3524 	BcNum temp, temp2, intn, frac;
3525 	BcRand state1, state2, inc1, inc2;
3526 	size_t nrdx = BC_NUM_RDX_VAL(n);
3527 
3528 	// This function holds the secret of how I interpret a seed number for the
3529 	// PRNG. Well, it's actually in the development manual
3530 	// (manuals/development.md#pseudo-random-number-generator), so look there
3531 	// before you try to understand this.
3532 
3533 	BC_SIG_LOCK;
3534 
3535 	bc_num_init(&temp, n->len);
3536 	bc_num_init(&temp2, n->len);
3537 	bc_num_init(&frac, nrdx);
3538 	bc_num_init(&intn, bc_num_int(n));
3539 
3540 	BC_SETJMP_LOCKED(err);
3541 
3542 	BC_SIG_UNLOCK;
3543 
3544 	assert(BC_NUM_RDX_VALID_NP(vm.max));
3545 
3546 	// NOLINTNEXTLINE
3547 	memcpy(frac.num, n->num, BC_NUM_SIZE(nrdx));
3548 	frac.len = nrdx;
3549 	BC_NUM_RDX_SET_NP(frac, nrdx);
3550 	frac.scale = n->scale;
3551 
3552 	assert(BC_NUM_RDX_VALID_NP(frac));
3553 	assert(BC_NUM_RDX_VALID_NP(vm.max2));
3554 
3555 	// Multiply the fraction and truncate so that it's an integer. The
3556 	// truncation is what clamps it, by the way.
3557 	bc_num_mul(&frac, &vm.max2, &temp, 0);
3558 	bc_num_truncate(&temp, temp.scale);
3559 	bc_num_copy(&frac, &temp);
3560 
3561 	// Get the integer.
3562 	// NOLINTNEXTLINE
3563 	memcpy(intn.num, n->num + nrdx, BC_NUM_SIZE(bc_num_int(n)));
3564 	intn.len = bc_num_int(n);
3565 
3566 	// This assert is here because it has to be true. It is also here to justify
3567 	// some optimizations.
3568 	assert(BC_NUM_NONZERO(&vm.max));
3569 
3570 	// If there *was* a fractional part...
3571 	if (BC_NUM_NONZERO(&frac))
3572 	{
3573 		// This divmod splits frac into the two state parts.
3574 		bc_num_divmod(&frac, &vm.max, &temp, &temp2, 0);
3575 
3576 		// frac is guaranteed to be smaller than vm.max * vm.max (pow).
3577 		// This means that when dividing frac by vm.max, as above, the
3578 		// quotient and remainder are both guaranteed to be less than vm.max,
3579 		// which means we can use bc_num_bigdig2() here and not worry about
3580 		// overflow.
3581 		state1 = (BcRand) bc_num_bigdig2(&temp2);
3582 		state2 = (BcRand) bc_num_bigdig2(&temp);
3583 	}
3584 	else state1 = state2 = 0;
3585 
3586 	// If there *was* an integer part...
3587 	if (BC_NUM_NONZERO(&intn))
3588 	{
3589 		// This divmod splits intn into the two inc parts.
3590 		bc_num_divmod(&intn, &vm.max, &temp, &temp2, 0);
3591 
3592 		// Because temp2 is the mod of vm.max, from above, it is guaranteed
3593 		// to be small enough to use bc_num_bigdig2().
3594 		inc1 = (BcRand) bc_num_bigdig2(&temp2);
3595 
3596 		// Clamp the second inc part.
3597 		if (bc_num_cmp(&temp, &vm.max) >= 0)
3598 		{
3599 			bc_num_copy(&temp2, &temp);
3600 			bc_num_mod(&temp2, &vm.max, &temp, 0);
3601 		}
3602 
3603 		// The if statement above ensures that temp is less than vm.max, which
3604 		// means that we can use bc_num_bigdig2() here.
3605 		inc2 = (BcRand) bc_num_bigdig2(&temp);
3606 	}
3607 	else inc1 = inc2 = 0;
3608 
3609 	bc_rand_seed(rng, state1, state2, inc1, inc2);
3610 
3611 err:
3612 	BC_SIG_MAYLOCK;
3613 	bc_num_free(&intn);
3614 	bc_num_free(&frac);
3615 	bc_num_free(&temp2);
3616 	bc_num_free(&temp);
3617 	BC_LONGJMP_CONT;
3618 }
3619 
3620 void
3621 bc_num_createFromRNG(BcNum* restrict n, BcRNG* rng)
3622 {
3623 	BcRand s1, s2, i1, i2;
3624 	BcNum conv, temp1, temp2, temp3;
3625 	BcDig temp1_num[BC_RAND_NUM_SIZE], temp2_num[BC_RAND_NUM_SIZE];
3626 	BcDig conv_num[BC_NUM_BIGDIG_LOG10];
3627 
3628 	BC_SIG_LOCK;
3629 
3630 	bc_num_init(&temp3, 2 * BC_RAND_NUM_SIZE);
3631 
3632 	BC_SETJMP_LOCKED(err);
3633 
3634 	BC_SIG_UNLOCK;
3635 
3636 	bc_num_setup(&temp1, temp1_num, sizeof(temp1_num) / sizeof(BcDig));
3637 	bc_num_setup(&temp2, temp2_num, sizeof(temp2_num) / sizeof(BcDig));
3638 	bc_num_setup(&conv, conv_num, sizeof(conv_num) / sizeof(BcDig));
3639 
3640 	// This assert is here because it has to be true. It is also here to justify
3641 	// the assumption that vm.max is not zero.
3642 	assert(BC_NUM_NONZERO(&vm.max));
3643 
3644 	// Because this is true, we can just ignore math errors that would happen
3645 	// otherwise.
3646 	assert(BC_NUM_NONZERO(&vm.max2));
3647 
3648 	bc_rand_getRands(rng, &s1, &s2, &i1, &i2);
3649 
3650 	// Put the second piece of state into a number.
3651 	bc_num_bigdig2num(&conv, (BcBigDig) s2);
3652 
3653 	assert(BC_NUM_RDX_VALID_NP(conv));
3654 
3655 	// Multiply by max to make room for the first piece of state.
3656 	bc_num_mul(&conv, &vm.max, &temp1, 0);
3657 
3658 	// Add in the first piece of state.
3659 	bc_num_bigdig2num(&conv, (BcBigDig) s1);
3660 	bc_num_add(&conv, &temp1, &temp2, 0);
3661 
3662 	// Divide to make it an entirely fractional part.
3663 	bc_num_div(&temp2, &vm.max2, &temp3, BC_RAND_STATE_BITS);
3664 
3665 	// Now start on the increment parts. It's the same process without the
3666 	// divide, so put the second piece of increment into a number.
3667 	bc_num_bigdig2num(&conv, (BcBigDig) i2);
3668 
3669 	assert(BC_NUM_RDX_VALID_NP(conv));
3670 
3671 	// Multiply by max to make room for the first piece of increment.
3672 	bc_num_mul(&conv, &vm.max, &temp1, 0);
3673 
3674 	// Add in the first piece of increment.
3675 	bc_num_bigdig2num(&conv, (BcBigDig) i1);
3676 	bc_num_add(&conv, &temp1, &temp2, 0);
3677 
3678 	// Now add the two together.
3679 	bc_num_add(&temp2, &temp3, n, 0);
3680 
3681 	assert(BC_NUM_RDX_VALID(n));
3682 
3683 err:
3684 	BC_SIG_MAYLOCK;
3685 	bc_num_free(&temp3);
3686 	BC_LONGJMP_CONT;
3687 }
3688 
3689 void
3690 bc_num_irand(BcNum* restrict a, BcNum* restrict b, BcRNG* restrict rng)
3691 {
3692 	BcNum atemp;
3693 	size_t i, len;
3694 
3695 	assert(a != b);
3696 
3697 	if (BC_ERR(BC_NUM_NEG(a))) bc_err(BC_ERR_MATH_NEGATIVE);
3698 
3699 	// If either of these are true, then the numbers are integers.
3700 	if (BC_NUM_ZERO(a) || BC_NUM_ONE(a)) return;
3701 
3702 	if (BC_ERR(bc_num_nonInt(a, &atemp))) bc_err(BC_ERR_MATH_NON_INTEGER);
3703 
3704 	assert(atemp.len);
3705 
3706 	len = atemp.len - 1;
3707 
3708 	// Just generate a random number for each limb.
3709 	for (i = 0; i < len; ++i)
3710 	{
3711 		b->num[i] = (BcDig) bc_rand_bounded(rng, BC_BASE_POW);
3712 	}
3713 
3714 	// Do the last digit explicitly because the bound must be right. But only
3715 	// do it if the limb does not equal 1. If it does, we have already hit the
3716 	// limit.
3717 	if (atemp.num[i] != 1)
3718 	{
3719 		b->num[i] = (BcDig) bc_rand_bounded(rng, (BcRand) atemp.num[i]);
3720 		b->len = atemp.len;
3721 	}
3722 	// We want 1 less len in the case where we skip the last limb.
3723 	else b->len = len;
3724 
3725 	bc_num_clean(b);
3726 
3727 	assert(BC_NUM_RDX_VALID(b));
3728 }
3729 #endif // BC_ENABLE_EXTRA_MATH
3730 
3731 size_t
3732 bc_num_addReq(const BcNum* a, const BcNum* b, size_t scale)
3733 {
3734 	size_t aint, bint, ardx, brdx;
3735 
3736 	// Addition and subtraction require the max of the length of the two numbers
3737 	// plus 1.
3738 
3739 	BC_UNUSED(scale);
3740 
3741 	ardx = BC_NUM_RDX_VAL(a);
3742 	aint = bc_num_int(a);
3743 	assert(aint <= a->len && ardx <= a->len);
3744 
3745 	brdx = BC_NUM_RDX_VAL(b);
3746 	bint = bc_num_int(b);
3747 	assert(bint <= b->len && brdx <= b->len);
3748 
3749 	ardx = BC_MAX(ardx, brdx);
3750 	aint = BC_MAX(aint, bint);
3751 
3752 	return bc_vm_growSize(bc_vm_growSize(ardx, aint), 1);
3753 }
3754 
3755 size_t
3756 bc_num_mulReq(const BcNum* a, const BcNum* b, size_t scale)
3757 {
3758 	size_t max, rdx;
3759 
3760 	// Multiplication requires the sum of the lengths of the numbers.
3761 
3762 	rdx = bc_vm_growSize(BC_NUM_RDX_VAL(a), BC_NUM_RDX_VAL(b));
3763 
3764 	max = BC_NUM_RDX(scale);
3765 
3766 	max = bc_vm_growSize(BC_MAX(max, rdx), 1);
3767 	rdx = bc_vm_growSize(bc_vm_growSize(bc_num_int(a), bc_num_int(b)), max);
3768 
3769 	return rdx;
3770 }
3771 
3772 size_t
3773 bc_num_divReq(const BcNum* a, const BcNum* b, size_t scale)
3774 {
3775 	size_t max, rdx;
3776 
3777 	// Division requires the length of the dividend plus the scale.
3778 
3779 	rdx = bc_vm_growSize(BC_NUM_RDX_VAL(a), BC_NUM_RDX_VAL(b));
3780 
3781 	max = BC_NUM_RDX(scale);
3782 
3783 	max = bc_vm_growSize(BC_MAX(max, rdx), 1);
3784 	rdx = bc_vm_growSize(bc_num_int(a), max);
3785 
3786 	return rdx;
3787 }
3788 
3789 size_t
3790 bc_num_powReq(const BcNum* a, const BcNum* b, size_t scale)
3791 {
3792 	BC_UNUSED(scale);
3793 	return bc_vm_growSize(bc_vm_growSize(a->len, b->len), 1);
3794 }
3795 
3796 #if BC_ENABLE_EXTRA_MATH
3797 size_t
3798 bc_num_placesReq(const BcNum* a, const BcNum* b, size_t scale)
3799 {
3800 	BC_UNUSED(scale);
3801 	return a->len + b->len - BC_NUM_RDX_VAL(a) - BC_NUM_RDX_VAL(b);
3802 }
3803 #endif // BC_ENABLE_EXTRA_MATH
3804 
3805 void
3806 bc_num_add(BcNum* a, BcNum* b, BcNum* c, size_t scale)
3807 {
3808 	assert(BC_NUM_RDX_VALID(a));
3809 	assert(BC_NUM_RDX_VALID(b));
3810 	bc_num_binary(a, b, c, false, bc_num_as, bc_num_addReq(a, b, scale));
3811 }
3812 
3813 void
3814 bc_num_sub(BcNum* a, BcNum* b, BcNum* c, size_t scale)
3815 {
3816 	assert(BC_NUM_RDX_VALID(a));
3817 	assert(BC_NUM_RDX_VALID(b));
3818 	bc_num_binary(a, b, c, true, bc_num_as, bc_num_addReq(a, b, scale));
3819 }
3820 
3821 void
3822 bc_num_mul(BcNum* a, BcNum* b, BcNum* c, size_t scale)
3823 {
3824 	assert(BC_NUM_RDX_VALID(a));
3825 	assert(BC_NUM_RDX_VALID(b));
3826 	bc_num_binary(a, b, c, scale, bc_num_m, bc_num_mulReq(a, b, scale));
3827 }
3828 
3829 void
3830 bc_num_div(BcNum* a, BcNum* b, BcNum* c, size_t scale)
3831 {
3832 	assert(BC_NUM_RDX_VALID(a));
3833 	assert(BC_NUM_RDX_VALID(b));
3834 	bc_num_binary(a, b, c, scale, bc_num_d, bc_num_divReq(a, b, scale));
3835 }
3836 
3837 void
3838 bc_num_mod(BcNum* a, BcNum* b, BcNum* c, size_t scale)
3839 {
3840 	assert(BC_NUM_RDX_VALID(a));
3841 	assert(BC_NUM_RDX_VALID(b));
3842 	bc_num_binary(a, b, c, scale, bc_num_rem, bc_num_divReq(a, b, scale));
3843 }
3844 
3845 void
3846 bc_num_pow(BcNum* a, BcNum* b, BcNum* c, size_t scale)
3847 {
3848 	assert(BC_NUM_RDX_VALID(a));
3849 	assert(BC_NUM_RDX_VALID(b));
3850 	bc_num_binary(a, b, c, scale, bc_num_p, bc_num_powReq(a, b, scale));
3851 }
3852 
3853 #if BC_ENABLE_EXTRA_MATH
3854 void
3855 bc_num_places(BcNum* a, BcNum* b, BcNum* c, size_t scale)
3856 {
3857 	assert(BC_NUM_RDX_VALID(a));
3858 	assert(BC_NUM_RDX_VALID(b));
3859 	bc_num_binary(a, b, c, scale, bc_num_place, bc_num_placesReq(a, b, scale));
3860 }
3861 
3862 void
3863 bc_num_lshift(BcNum* a, BcNum* b, BcNum* c, size_t scale)
3864 {
3865 	assert(BC_NUM_RDX_VALID(a));
3866 	assert(BC_NUM_RDX_VALID(b));
3867 	bc_num_binary(a, b, c, scale, bc_num_left, bc_num_placesReq(a, b, scale));
3868 }
3869 
3870 void
3871 bc_num_rshift(BcNum* a, BcNum* b, BcNum* c, size_t scale)
3872 {
3873 	assert(BC_NUM_RDX_VALID(a));
3874 	assert(BC_NUM_RDX_VALID(b));
3875 	bc_num_binary(a, b, c, scale, bc_num_right, bc_num_placesReq(a, b, scale));
3876 }
3877 #endif // BC_ENABLE_EXTRA_MATH
3878 
3879 void
3880 bc_num_sqrt(BcNum* restrict a, BcNum* restrict b, size_t scale)
3881 {
3882 	BcNum num1, num2, half, f, fprime;
3883 	BcNum* x0;
3884 	BcNum* x1;
3885 	BcNum* temp;
3886 	size_t pow, len, rdx, req, resscale;
3887 	BcDig half_digs[1];
3888 
3889 	assert(a != NULL && b != NULL && a != b);
3890 
3891 	if (BC_ERR(BC_NUM_NEG(a))) bc_err(BC_ERR_MATH_NEGATIVE);
3892 
3893 	// We want to calculate to a's scale if it is bigger so that the result will
3894 	// truncate properly.
3895 	if (a->scale > scale) scale = a->scale;
3896 
3897 	// Set parameters for the result.
3898 	len = bc_vm_growSize(bc_num_intDigits(a), 1);
3899 	rdx = BC_NUM_RDX(scale);
3900 
3901 	// Square root needs half of the length of the parameter.
3902 	req = bc_vm_growSize(BC_MAX(rdx, BC_NUM_RDX_VAL(a)), len >> 1);
3903 
3904 	BC_SIG_LOCK;
3905 
3906 	// Unlike the binary operators, this function is the only single parameter
3907 	// function and is expected to initialize the result. This means that it
3908 	// expects that b is *NOT* preallocated. We allocate it here.
3909 	bc_num_init(b, bc_vm_growSize(req, 1));
3910 
3911 	BC_SIG_UNLOCK;
3912 
3913 	assert(a != NULL && b != NULL && a != b);
3914 	assert(a->num != NULL && b->num != NULL);
3915 
3916 	// Easy case.
3917 	if (BC_NUM_ZERO(a))
3918 	{
3919 		bc_num_setToZero(b, scale);
3920 		return;
3921 	}
3922 
3923 	// Another easy case.
3924 	if (BC_NUM_ONE(a))
3925 	{
3926 		bc_num_one(b);
3927 		bc_num_extend(b, scale);
3928 		return;
3929 	}
3930 
3931 	// Set the parameters again.
3932 	rdx = BC_NUM_RDX(scale);
3933 	rdx = BC_MAX(rdx, BC_NUM_RDX_VAL(a));
3934 	len = bc_vm_growSize(a->len, rdx);
3935 
3936 	BC_SIG_LOCK;
3937 
3938 	bc_num_init(&num1, len);
3939 	bc_num_init(&num2, len);
3940 	bc_num_setup(&half, half_digs, sizeof(half_digs) / sizeof(BcDig));
3941 
3942 	// There is a division by two in the formula. We setup a number that's 1/2
3943 	// so that we can use multiplication instead of heavy division.
3944 	bc_num_one(&half);
3945 	half.num[0] = BC_BASE_POW / 2;
3946 	half.len = 1;
3947 	BC_NUM_RDX_SET_NP(half, 1);
3948 	half.scale = 1;
3949 
3950 	bc_num_init(&f, len);
3951 	bc_num_init(&fprime, len);
3952 
3953 	BC_SETJMP_LOCKED(err);
3954 
3955 	BC_SIG_UNLOCK;
3956 
3957 	// Pointers for easy switching.
3958 	x0 = &num1;
3959 	x1 = &num2;
3960 
3961 	// Start with 1.
3962 	bc_num_one(x0);
3963 
3964 	// The power of the operand is needed for the estimate.
3965 	pow = bc_num_intDigits(a);
3966 
3967 	// The code in this if statement calculates the initial estimate. First, if
3968 	// a is less than 0, then 0 is a good estimate. Otherwise, we want something
3969 	// in the same ballpark. That ballpark is pow.
3970 	if (pow)
3971 	{
3972 		// An odd number is served by starting with 2^((pow-1)/2), and an even
3973 		// number is served by starting with 6^((pow-2)/2). Why? Because math.
3974 		if (pow & 1) x0->num[0] = 2;
3975 		else x0->num[0] = 6;
3976 
3977 		pow -= 2 - (pow & 1);
3978 		bc_num_shiftLeft(x0, pow / 2);
3979 	}
3980 
3981 	// I can set the rdx here directly because neg should be false.
3982 	x0->scale = x0->rdx = 0;
3983 	resscale = (scale + BC_BASE_DIGS) + 2;
3984 
3985 	// This is the calculation loop. This compare goes to 0 eventually as the
3986 	// difference between the two numbers gets smaller than resscale.
3987 	while (bc_num_cmp(x1, x0))
3988 	{
3989 		assert(BC_NUM_NONZERO(x0));
3990 
3991 		// This loop directly corresponds to the iteration in Newton's method.
3992 		// If you know the formula, this loop makes sense. Go study the formula.
3993 
3994 		bc_num_div(a, x0, &f, resscale);
3995 		bc_num_add(x0, &f, &fprime, resscale);
3996 
3997 		assert(BC_NUM_RDX_VALID_NP(fprime));
3998 		assert(BC_NUM_RDX_VALID_NP(half));
3999 
4000 		bc_num_mul(&fprime, &half, x1, resscale);
4001 
4002 		// Switch.
4003 		temp = x0;
4004 		x0 = x1;
4005 		x1 = temp;
4006 	}
4007 
4008 	// Copy to the result and truncate.
4009 	bc_num_copy(b, x0);
4010 	if (b->scale > scale) bc_num_truncate(b, b->scale - scale);
4011 
4012 	assert(!BC_NUM_NEG(b) || BC_NUM_NONZERO(b));
4013 	assert(BC_NUM_RDX_VALID(b));
4014 	assert(BC_NUM_RDX_VAL(b) <= b->len || !b->len);
4015 	assert(!b->len || b->num[b->len - 1] || BC_NUM_RDX_VAL(b) == b->len);
4016 
4017 err:
4018 	BC_SIG_MAYLOCK;
4019 	bc_num_free(&fprime);
4020 	bc_num_free(&f);
4021 	bc_num_free(&num2);
4022 	bc_num_free(&num1);
4023 	BC_LONGJMP_CONT;
4024 }
4025 
4026 void
4027 bc_num_divmod(BcNum* a, BcNum* b, BcNum* c, BcNum* d, size_t scale)
4028 {
4029 	size_t ts, len;
4030 	BcNum *ptr_a, num2;
4031 	bool init = false;
4032 
4033 	// The bulk of this function is just doing what bc_num_binary() does for the
4034 	// binary operators. However, it assumes that only c and a can be equal.
4035 
4036 	// Set up the parameters.
4037 	ts = BC_MAX(scale + b->scale, a->scale);
4038 	len = bc_num_mulReq(a, b, ts);
4039 
4040 	assert(a != NULL && b != NULL && c != NULL && d != NULL);
4041 	assert(c != d && a != d && b != d && b != c);
4042 
4043 	// Initialize or expand as necessary.
4044 	if (c == a)
4045 	{
4046 		// NOLINTNEXTLINE
4047 		memcpy(&num2, c, sizeof(BcNum));
4048 		ptr_a = &num2;
4049 
4050 		BC_SIG_LOCK;
4051 
4052 		bc_num_init(c, len);
4053 
4054 		init = true;
4055 
4056 		BC_SETJMP_LOCKED(err);
4057 
4058 		BC_SIG_UNLOCK;
4059 	}
4060 	else
4061 	{
4062 		ptr_a = a;
4063 		bc_num_expand(c, len);
4064 	}
4065 
4066 	// Do the quick version if possible.
4067 	if (BC_NUM_NONZERO(a) && !BC_NUM_RDX_VAL(a) && !BC_NUM_RDX_VAL(b) &&
4068 	    b->len == 1 && !scale)
4069 	{
4070 		BcBigDig rem;
4071 
4072 		bc_num_divArray(ptr_a, (BcBigDig) b->num[0], c, &rem);
4073 
4074 		assert(rem < BC_BASE_POW);
4075 
4076 		d->num[0] = (BcDig) rem;
4077 		d->len = (rem != 0);
4078 	}
4079 	// Do the slow method.
4080 	else bc_num_r(ptr_a, b, c, d, scale, ts);
4081 
4082 	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
4083 	assert(BC_NUM_RDX_VALID(c));
4084 	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
4085 	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
4086 	assert(!BC_NUM_NEG(d) || BC_NUM_NONZERO(d));
4087 	assert(BC_NUM_RDX_VALID(d));
4088 	assert(BC_NUM_RDX_VAL(d) <= d->len || !d->len);
4089 	assert(!d->len || d->num[d->len - 1] || BC_NUM_RDX_VAL(d) == d->len);
4090 
4091 err:
4092 	// Only cleanup if we initialized.
4093 	if (init)
4094 	{
4095 		BC_SIG_MAYLOCK;
4096 		bc_num_free(&num2);
4097 		BC_LONGJMP_CONT;
4098 	}
4099 }
4100 
4101 void
4102 bc_num_modexp(BcNum* a, BcNum* b, BcNum* c, BcNum* restrict d)
4103 {
4104 	BcNum base, exp, two, temp, atemp, btemp, ctemp;
4105 	BcDig two_digs[2];
4106 
4107 	assert(a != NULL && b != NULL && c != NULL && d != NULL);
4108 	assert(a != d && b != d && c != d);
4109 
4110 	if (BC_ERR(BC_NUM_ZERO(c))) bc_err(BC_ERR_MATH_DIVIDE_BY_ZERO);
4111 
4112 	if (BC_ERR(BC_NUM_NEG(b))) bc_err(BC_ERR_MATH_NEGATIVE);
4113 
4114 #ifndef NDEBUG
4115 	// This is entirely for quieting a useless scan-build error.
4116 	btemp.len = 0;
4117 	ctemp.len = 0;
4118 #endif // NDEBUG
4119 
4120 	// Eliminate fractional parts that are zero or error if they are not zero.
4121 	if (BC_ERR(bc_num_nonInt(a, &atemp) || bc_num_nonInt(b, &btemp) ||
4122 	           bc_num_nonInt(c, &ctemp)))
4123 	{
4124 		bc_err(BC_ERR_MATH_NON_INTEGER);
4125 	}
4126 
4127 	bc_num_expand(d, ctemp.len);
4128 
4129 	BC_SIG_LOCK;
4130 
4131 	bc_num_init(&base, ctemp.len);
4132 	bc_num_setup(&two, two_digs, sizeof(two_digs) / sizeof(BcDig));
4133 	bc_num_init(&temp, btemp.len + 1);
4134 	bc_num_createCopy(&exp, &btemp);
4135 
4136 	BC_SETJMP_LOCKED(err);
4137 
4138 	BC_SIG_UNLOCK;
4139 
4140 	bc_num_one(&two);
4141 	two.num[0] = 2;
4142 	bc_num_one(d);
4143 
4144 	// We already checked for 0.
4145 	bc_num_rem(&atemp, &ctemp, &base, 0);
4146 
4147 	// If you know the algorithm I used, the memory-efficient method, then this
4148 	// loop should be self-explanatory because it is the calculation loop.
4149 	while (BC_NUM_NONZERO(&exp))
4150 	{
4151 		// Num two cannot be 0, so no errors.
4152 		bc_num_divmod(&exp, &two, &exp, &temp, 0);
4153 
4154 		if (BC_NUM_ONE(&temp) && !BC_NUM_NEG_NP(temp))
4155 		{
4156 			assert(BC_NUM_RDX_VALID(d));
4157 			assert(BC_NUM_RDX_VALID_NP(base));
4158 
4159 			bc_num_mul(d, &base, &temp, 0);
4160 
4161 			// We already checked for 0.
4162 			bc_num_rem(&temp, &ctemp, d, 0);
4163 		}
4164 
4165 		assert(BC_NUM_RDX_VALID_NP(base));
4166 
4167 		bc_num_mul(&base, &base, &temp, 0);
4168 
4169 		// We already checked for 0.
4170 		bc_num_rem(&temp, &ctemp, &base, 0);
4171 	}
4172 
4173 err:
4174 	BC_SIG_MAYLOCK;
4175 	bc_num_free(&exp);
4176 	bc_num_free(&temp);
4177 	bc_num_free(&base);
4178 	BC_LONGJMP_CONT;
4179 	assert(!BC_NUM_NEG(d) || d->len);
4180 	assert(BC_NUM_RDX_VALID(d));
4181 	assert(!d->len || d->num[d->len - 1] || BC_NUM_RDX_VAL(d) == d->len);
4182 }
4183 
4184 #if BC_DEBUG_CODE
4185 void
4186 bc_num_printDebug(const BcNum* n, const char* name, bool emptyline)
4187 {
4188 	bc_file_puts(&vm.fout, bc_flush_none, name);
4189 	bc_file_puts(&vm.fout, bc_flush_none, ": ");
4190 	bc_num_printDecimal(n, true);
4191 	bc_file_putchar(&vm.fout, bc_flush_err, '\n');
4192 	if (emptyline) bc_file_putchar(&vm.fout, bc_flush_err, '\n');
4193 	vm.nchars = 0;
4194 }
4195 
4196 void
4197 bc_num_printDigs(const BcDig* n, size_t len, bool emptyline)
4198 {
4199 	size_t i;
4200 
4201 	for (i = len - 1; i < len; --i)
4202 	{
4203 		bc_file_printf(&vm.fout, " %lu", (unsigned long) n[i]);
4204 	}
4205 
4206 	bc_file_putchar(&vm.fout, bc_flush_err, '\n');
4207 	if (emptyline) bc_file_putchar(&vm.fout, bc_flush_err, '\n');
4208 	vm.nchars = 0;
4209 }
4210 
4211 void
4212 bc_num_printWithDigs(const BcNum* n, const char* name, bool emptyline)
4213 {
4214 	bc_file_puts(&vm.fout, bc_flush_none, name);
4215 	bc_file_printf(&vm.fout, " len: %zu, rdx: %zu, scale: %zu\n", name, n->len,
4216 	               BC_NUM_RDX_VAL(n), n->scale);
4217 	bc_num_printDigs(n->num, n->len, emptyline);
4218 }
4219 
4220 void
4221 bc_num_dump(const char* varname, const BcNum* n)
4222 {
4223 	ulong i, scale = n->scale;
4224 
4225 	bc_file_printf(&vm.ferr, "\n%s = %s", varname,
4226 	               n->len ? (BC_NUM_NEG(n) ? "-" : "+") : "0 ");
4227 
4228 	for (i = n->len - 1; i < n->len; --i)
4229 	{
4230 		if (i + 1 == BC_NUM_RDX_VAL(n))
4231 		{
4232 			bc_file_puts(&vm.ferr, bc_flush_none, ". ");
4233 		}
4234 
4235 		if (scale / BC_BASE_DIGS != BC_NUM_RDX_VAL(n) - i - 1)
4236 		{
4237 			bc_file_printf(&vm.ferr, "%lu ", (unsigned long) n->num[i]);
4238 		}
4239 		else
4240 		{
4241 			int mod = scale % BC_BASE_DIGS;
4242 			int d = BC_BASE_DIGS - mod;
4243 			BcDig div;
4244 
4245 			if (mod != 0)
4246 			{
4247 				div = n->num[i] / ((BcDig) bc_num_pow10[(ulong) d]);
4248 				bc_file_printf(&vm.ferr, "%lu", (unsigned long) div);
4249 			}
4250 
4251 			div = n->num[i] % ((BcDig) bc_num_pow10[(ulong) d]);
4252 			bc_file_printf(&vm.ferr, " ' %lu ", (unsigned long) div);
4253 		}
4254 	}
4255 
4256 	bc_file_printf(&vm.ferr, "(%zu | %zu.%zu / %zu) %lu\n", n->scale, n->len,
4257 	               BC_NUM_RDX_VAL(n), n->cap, (unsigned long) (void*) n->num);
4258 
4259 	bc_file_flush(&vm.ferr, bc_flush_err);
4260 }
4261 #endif // BC_DEBUG_CODE
4262