src/num.c

/*
 * *****************************************************************************
 *
 * SPDX-License-Identifier: BSD-2-Clause
 *
 * Copyright (c) 2018-2021 Gavin D. Howard and contributors.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * * Redistributions of source code must retain the above copyright notice, this
 *   list of conditions and the following disclaimer.
 *
 * * Redistributions in binary form must reproduce the above copyright notice,
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 * *****************************************************************************
 *
 * Code for the number type.
 *
 */

#include <assert.h>
#include <ctype.h>
#include <stdbool.h>
#include <stdlib.h>
#include <string.h>
#include <setjmp.h>
#include <limits.h>

#include <num.h>
#include <rand.h>
#include <vm.h>

// Before you try to understand this code, see the development manual
// (manuals/development.md#numbers).

static void bc_num_m(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale);

/**
 * Multiply two numbers and throw a math error if they overflow.
 * @param a  The first operand.
 * @param b  The second operand.
 * @return   The product of the two operands.
 */
static inline size_t bc_num_mulOverflow(size_t a, size_t b) {
	size_t res = a * b;
	if (BC_ERR(BC_VM_MUL_OVERFLOW(a, b, res))) bc_err(BC_ERR_MATH_OVERFLOW);
	return res;
}

/**
 * Conditionally negate @a n based on @a neg. Algorithm taken from
 * https://graphics.stanford.edu/~seander/bithacks.html#ConditionalNegate .
 * @param n    The value to turn into a signed value and negate.
 * @param neg  The condition to negate or not.
 */
static inline ssize_t bc_num_neg(size_t n, bool neg) {
	return (((ssize_t) n) ^ -((ssize_t) neg)) + neg;
}

/**
 * Compare a BcNum against zero.
 * @param n  The number to compare.
 * @return   -1 if the number is less than 0, 1 if greater, and 0 if equal.
 */
ssize_t bc_num_cmpZero(const BcNum *n) {
	return bc_num_neg((n)->len != 0, BC_NUM_NEG(n));
}

/**
 * Return the number of integer limbs in a BcNum. This is the opposite of rdx.
 * @param n  The number to return the amount of integer limbs for.
 * @return   The amount of integer limbs in @a n.
 */
static inline size_t bc_num_int(const BcNum *n) {
	return n->len ? n->len - BC_NUM_RDX_VAL(n) : 0;
}

/**
 * Expand a number's allocation capacity to at least req limbs.
 * @param n    The number to expand.
 * @param req  The number limbs to expand the allocation capacity to.
 */
static void bc_num_expand(BcNum *restrict n, size_t req) {

	assert(n != NULL);

	req = req >= BC_NUM_DEF_SIZE ? req : BC_NUM_DEF_SIZE;

	if (req > n->cap) {

		BC_SIG_LOCK;

		n->num = bc_vm_realloc(n->num, BC_NUM_SIZE(req));
		n->cap = req;

		BC_SIG_UNLOCK;
	}
}

/**
 * Set a number to 0 with the specified scale.
 * @param n      The number to set to zero.
 * @param scale  The scale to set the number to.
 */
static void bc_num_setToZero(BcNum *restrict n, size_t scale) {
	assert(n != NULL);
	n->scale = scale;
	n->len = n->rdx = 0;
}

void bc_num_zero(BcNum *restrict n) {
	bc_num_setToZero(n, 0);
}

void bc_num_one(BcNum *restrict n) {
	bc_num_zero(n);
	n->len = 1;
	n->num[0] = 1;
}

/**
 * "Cleans" a number, which means reducing the length if the most significant
 * limbs are zero.
 * @param n  The number to clean.
 */
static void bc_num_clean(BcNum *restrict n) {

	// Reduce the length.
	while (BC_NUM_NONZERO(n) && !n->num[n->len - 1]) n->len -= 1;

	// Special cases.
	if (BC_NUM_ZERO(n)) n->rdx = 0;
	else {

		// len must be at least as much as rdx.
		size_t rdx = BC_NUM_RDX_VAL(n);
		if (n->len < rdx) n->len = rdx;
	}
}

/**
 * Returns the log base 10 of @a i. I could have done this with floating-point
 * math, and in fact, I originally did. However, that was the only
 * floating-point code in the entire codebase, and I decided I didn't want any.
 * This is fast enough. Also, it might handle larger numbers better.
 * @param i  The number to return the log base 10 of.
 * @return   The log base 10 of @a i.
 */
static size_t bc_num_log10(size_t i) {
	size_t len;
	for (len = 1; i; i /= BC_BASE, ++len);
	assert(len - 1 <= BC_BASE_DIGS + 1);
	return len - 1;
}

/**
 * Returns the number of decimal digits in a limb that are zero starting at the
 * most significant digits. This basically returns how much of the limb is used.
 * @param n  The number.
 * @return   The number of decimal digits that are 0 starting at the most
 *           significant digits.
 */
static inline size_t bc_num_zeroDigits(const BcDig *n) {
	assert(*n >= 0);
	assert(((size_t) *n) < BC_BASE_POW);
	return BC_BASE_DIGS - bc_num_log10((size_t) *n);
}

/**
 * Return the total number of integer digits in a number. This is the opposite
 * of scale, like bc_num_int() is the opposite of rdx.
 * @param n  The number.
 * @return   The number of integer digits in @a n.
 */
static size_t bc_num_intDigits(const BcNum *n) {
	size_t digits = bc_num_int(n) * BC_BASE_DIGS;
	if (digits > 0) digits -= bc_num_zeroDigits(n->num + n->len - 1);
	return digits;
}

/**
 * Returns the number of limbs of a number that are non-zero starting at the
 * most significant limbs. This expects that there are *no* integer limbs in the
 * number because it is specifically to figure out how many zero limbs after the
 * decimal place to ignore. If there are zero limbs after non-zero limbs, they
 * are counted as non-zero limbs.
 * @param n  The number.
 * @return   The number of non-zero limbs after the decimal point.
 */
static size_t bc_num_nonZeroLen(const BcNum *restrict n) {
	size_t i, len = n->len;
	assert(len == BC_NUM_RDX_VAL(n));
	for (i = len - 1; i < len && !n->num[i]; --i);
	assert(i + 1 > 0);
	return i + 1;
}

/**
 * Performs a one-limb add with a carry.
 * @param a      The first limb.
 * @param b      The second limb.
 * @param carry  An in/out parameter; the carry in from the previous add and the
 *               carry out from this add.
 * @return       The resulting limb sum.
 */
static BcDig bc_num_addDigits(BcDig a, BcDig b, bool *carry) {

	assert(((BcBigDig) BC_BASE_POW) * 2 == ((BcDig) BC_BASE_POW) * 2);
	assert(a < BC_BASE_POW);
	assert(b < BC_BASE_POW);

	a += b + *carry;
	*carry = (a >= BC_BASE_POW);
	if (*carry) a -= BC_BASE_POW;

	assert(a >= 0);
	assert(a < BC_BASE_POW);

	return a;
}

/**
 * Performs a one-limb subtract with a carry.
 * @param a      The first limb.
 * @param b      The second limb.
 * @param carry  An in/out parameter; the carry in from the previous subtract
 *               and the carry out from this subtract.
 * @return       The resulting limb difference.
 */
static BcDig bc_num_subDigits(BcDig a, BcDig b, bool *carry) {

	assert(a < BC_BASE_POW);
	assert(b < BC_BASE_POW);

	b += *carry;
	*carry = (a < b);
	if (*carry) a += BC_BASE_POW;

	assert(a - b >= 0);
	assert(a - b < BC_BASE_POW);

	return a - b;
}

/**
 * Add two BcDig arrays and store the result in the first array.
 * @param a    The first operand and out array.
 * @param b    The second operand.
 * @param len  The length of @a b.
 */
static void bc_num_addArrays(BcDig *restrict a, const BcDig *restrict b,
                             size_t len)
{
	size_t i;
	bool carry = false;

	for (i = 0; i < len; ++i) a[i] = bc_num_addDigits(a[i], b[i], &carry);

	// Take care of the extra limbs in the bigger array.
	for (; carry; ++i) a[i] = bc_num_addDigits(a[i], 0, &carry);
}

/**
 * Subtract two BcDig arrays and store the result in the first array.
 * @param a    The first operand and out array.
 * @param b    The second operand.
 * @param len  The length of @a b.
 */
static void bc_num_subArrays(BcDig *restrict a, const BcDig *restrict b,
                             size_t len)
{
	size_t i;
	bool carry = false;

	for (i = 0; i < len; ++i) a[i] = bc_num_subDigits(a[i], b[i], &carry);

	// Take care of the extra limbs in the bigger array.
	for (; carry; ++i) a[i] = bc_num_subDigits(a[i], 0, &carry);
}

/**
 * Multiply a BcNum array by a one-limb number. This is a faster version of
 * multiplication for when we can use it.
 * @param a  The BcNum to multiply by the one-limb number.
 * @param b  The one limb of the one-limb number.
 * @param c  The return parameter.
 */
static void bc_num_mulArray(const BcNum *restrict a, BcBigDig b,
                            BcNum *restrict c)
{
	size_t i;
	BcBigDig carry = 0;

	assert(b <= BC_BASE_POW);

	// Make sure the return parameter is big enough.
	if (a->len + 1 > c->cap) bc_num_expand(c, a->len + 1);

	// We want the entire return parameter to be zero for cleaning later.
	memset(c->num, 0, BC_NUM_SIZE(c->cap));

	// Actual multiplication loop.
	for (i = 0; i < a->len; ++i) {
		BcBigDig in = ((BcBigDig) a->num[i]) * b + carry;
		c->num[i] = in % BC_BASE_POW;
		carry = in / BC_BASE_POW;
	}

	assert(carry < BC_BASE_POW);

	// Finishing touches.
	c->num[i] = (BcDig) carry;
	c->len = a->len;
	c->len += (carry != 0);

	bc_num_clean(c);

	// Postconditions.
	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
}

/**
 * Divide a BcNum array by a one-limb number. This is a faster version of divide
 * for when we can use it.
 * @param a    The BcNum to multiply by the one-limb number.
 * @param b    The one limb of the one-limb number.
 * @param c    The return parameter for the quotient.
 * @param rem  The return parameter for the remainder.
 */
static void bc_num_divArray(const BcNum *restrict a, BcBigDig b,
                            BcNum *restrict c, BcBigDig *rem)
{
	size_t i;
	BcBigDig carry = 0;

	assert(c->cap >= a->len);

	// Actual division loop.
	for (i = a->len - 1; i < a->len; --i) {
		BcBigDig in = ((BcBigDig) a->num[i]) + carry * BC_BASE_POW;
		assert(in / b < BC_BASE_POW);
		c->num[i] = (BcDig) (in / b);
		carry = in % b;
	}

	// Finishing touches.
	c->len = a->len;
	bc_num_clean(c);
	*rem = carry;

	// Postconditions.
	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
}

/**
 * Compare two BcDig arrays and return >0 if @a b is greater, <0 if @a b is
 * less, and 0 if equal. Both @a a and @a b must have the same length.
 * @param a    The first array.
 * @param b    The second array.
 * @param len  The minimum length of the arrays.
 */
static ssize_t bc_num_compare(const BcDig *restrict a, const BcDig *restrict b,
                              size_t len)
{
	size_t i;
	BcDig c = 0;
	for (i = len - 1; i < len && !(c = a[i] - b[i]); --i);
	return bc_num_neg(i + 1, c < 0);
}

ssize_t bc_num_cmp(const BcNum *a, const BcNum *b) {

	size_t i, min, a_int, b_int, diff, ardx, brdx;
	BcDig *max_num, *min_num;
	bool a_max, neg = false;
	ssize_t cmp;

	assert(a != NULL && b != NULL);

	// Same num? Equal.
	if (a == b) return 0;

	// Easy cases.
	if (BC_NUM_ZERO(a)) return bc_num_neg(b->len != 0, !BC_NUM_NEG(b));
	if (BC_NUM_ZERO(b)) return bc_num_cmpZero(a);
	if (BC_NUM_NEG(a)) {
		if (BC_NUM_NEG(b)) neg = true;
		else return -1;
	}
	else if (BC_NUM_NEG(b)) return 1;

	// Get the number of int limbs in each number and get the difference.
	a_int = bc_num_int(a);
	b_int = bc_num_int(b);
	a_int -= b_int;

	// If there's a difference, then just return the comparison.
	if (a_int) return neg ? -((ssize_t) a_int) : (ssize_t) a_int;

	// Get the rdx's and figure out the max.
	ardx = BC_NUM_RDX_VAL(a);
	brdx = BC_NUM_RDX_VAL(b);
	a_max = (ardx > brdx);

	// Set variables based on the above.
	if (a_max) {
		min = brdx;
		diff = ardx - brdx;
		max_num = a->num + diff;
		min_num = b->num;
	}
	else {
		min = ardx;
		diff = brdx - ardx;
		max_num = b->num + diff;
		min_num = a->num;
	}

	// Do a full limb-by-limb comparison.
	cmp = bc_num_compare(max_num, min_num, b_int + min);

	// If we found a difference, return it based on state.
	if (cmp) return bc_num_neg((size_t) cmp, !a_max == !neg);

	// If there was no difference, then the final step is to check which number
	// has greater or lesser limbs beyond the other's.
	for (max_num -= diff, i = diff - 1; i < diff; --i) {
		if (max_num[i]) return bc_num_neg(1, !a_max == !neg);
	}

	return 0;
}

void bc_num_truncate(BcNum *restrict n, size_t places) {

	size_t nrdx, places_rdx;

	if (!places) return;

	// Grab these needed values; places_rdx is the rdx equivalent to places like
	// rdx is to scale.
	nrdx = BC_NUM_RDX_VAL(n);
	places_rdx = nrdx ? nrdx - BC_NUM_RDX(n->scale - places) : 0;

	// We cannot truncate more places than we have.
	assert(places <= n->scale && (BC_NUM_ZERO(n) || places_rdx <= n->len));

	n->scale -= places;
	BC_NUM_RDX_SET(n, nrdx - places_rdx);

	// Only when the number is nonzero do we need to do the hard stuff.
	if (BC_NUM_NONZERO(n)) {

		size_t pow;

		// This calculates how many decimal digits are in the least significant
		// limb.
		pow = n->scale % BC_BASE_DIGS;
		pow = pow ? BC_BASE_DIGS - pow : 0;
		pow = bc_num_pow10[pow];

		n->len -= places_rdx;

		// We have to move limbs to maintain invariants. The limbs must begin at
		// the beginning of the BcNum array.
		memmove(n->num, n->num + places_rdx, BC_NUM_SIZE(n->len));

		// Clear the lower part of the last digit.
		if (BC_NUM_NONZERO(n)) n->num[0] -= n->num[0] % (BcDig) pow;

		bc_num_clean(n);
	}
}

void bc_num_extend(BcNum *restrict n, size_t places) {

	size_t nrdx, places_rdx;

	if (!places) return;

	// Easy case with zero; set the scale.
	if (BC_NUM_ZERO(n)) {
		n->scale += places;
		return;
	}

	// Grab these needed values; places_rdx is the rdx equivalent to places like
	// rdx is to scale.
	nrdx = BC_NUM_RDX_VAL(n);
	places_rdx = BC_NUM_RDX(places + n->scale) - nrdx;

	// This is the hard case. We need to expand the number, move the limbs, and
	// set the limbs that were just cleared.
	if (places_rdx) {
		bc_num_expand(n, bc_vm_growSize(n->len, places_rdx));
		memmove(n->num + places_rdx, n->num, BC_NUM_SIZE(n->len));
		memset(n->num, 0, BC_NUM_SIZE(places_rdx));
	}

	// Finally, set scale and rdx.
	BC_NUM_RDX_SET(n, nrdx + places_rdx);
	n->scale += places;
	n->len += places_rdx;

	assert(BC_NUM_RDX_VAL(n) == BC_NUM_RDX(n->scale));
}

/**
 * Retires (finishes) a multiplication or division operation.
 */
static void bc_num_retireMul(BcNum *restrict n, size_t scale,
                             bool neg1, bool neg2)
{
	// Make sure scale is correct.
	if (n->scale < scale) bc_num_extend(n, scale - n->scale);
	else bc_num_truncate(n, n->scale - scale);

	bc_num_clean(n);

	// We need to ensure rdx is correct.
	if (BC_NUM_NONZERO(n)) n->rdx = BC_NUM_NEG_VAL(n, !neg1 != !neg2);
}

/**
 * Splits a number into two BcNum's. This is used in Karatsuba.
 * @param n    The number to split.
 * @param idx  The index to split at.
 * @param a    An out parameter; the low part of @a n.
 * @param b    An out parameter; the high part of @a n.
 */
static void bc_num_split(const BcNum *restrict n, size_t idx,
                         BcNum *restrict a, BcNum *restrict b)
{
	// We want a and b to be clear.
	assert(BC_NUM_ZERO(a));
	assert(BC_NUM_ZERO(b));

	// The usual case.
	if (idx < n->len) {

		// Set the fields first.
		b->len = n->len - idx;
		a->len = idx;
		a->scale = b->scale = 0;
		BC_NUM_RDX_SET(a, 0);
		BC_NUM_RDX_SET(b, 0);

		assert(a->cap >= a->len);
		assert(b->cap >= b->len);

		// Copy the arrays. This is not necessary for safety, but it is faster,
		// for some reason.
		memcpy(b->num, n->num + idx, BC_NUM_SIZE(b->len));
		memcpy(a->num, n->num, BC_NUM_SIZE(idx));

		bc_num_clean(b);
	}
	// If the index is weird, just skip the split.
	else bc_num_copy(a, n);

	bc_num_clean(a);
}

/**
 * Copies a number into another, but shifts the rdx so that the result number
 * only sees the integer part of the source number.
 * @param n  The number to copy.
 * @param r  The result number with a shifted rdx, len, and num.
 */
static void bc_num_shiftRdx(const BcNum *restrict n, BcNum *restrict r) {

	size_t rdx = BC_NUM_RDX_VAL(n);

	r->len = n->len - rdx;
	r->cap = n->cap - rdx;
	r->num = n->num + rdx;

	BC_NUM_RDX_SET_NEG(r, 0, BC_NUM_NEG(n));
	r->scale = 0;
}

/**
 * Shifts a number so that all of the least significant limbs of the number are
 * skipped. This must be undone by bc_num_unshiftZero().
 * @param n  The number to shift.
 */
static size_t bc_num_shiftZero(BcNum *restrict n) {

	size_t i;

	// If we don't have an integer, that is a problem, but it's also a bug
	// because the caller should have set everything up right.
	assert(!BC_NUM_RDX_VAL(n) || BC_NUM_ZERO(n));

	for (i = 0; i < n->len && !n->num[i]; ++i);

	n->len -= i;
	n->num += i;

	return i;
}

/**
 * Undo the damage done by bc_num_unshiftZero(). This must be called like all
 * cleanup functions: after a label used by setjmp() and longjmp().
 * @param n           The number to unshift.
 * @param places_rdx  The amount the number was originally shift.
 */
static void bc_num_unshiftZero(BcNum *restrict n, size_t places_rdx) {
	n->len += places_rdx;
	n->num -= places_rdx;
}

/**
 * Shifts the digits right within a number by no more than BC_BASE_DIGS - 1.
 * This is the final step on shifting numbers with the shift operators. It
 * depends on the caller to shift the limbs properly because all it does is
 * ensure that the rdx point is realigned to be between limbs.
 * @param n    The number to shift digits in.
 * @param dig  The number of places to shift right.
 */
static void bc_num_shift(BcNum *restrict n, BcBigDig dig) {

	size_t i, len = n->len;
	BcBigDig carry = 0, pow;
	BcDig *ptr = n->num;

	assert(dig < BC_BASE_DIGS);

	// Figure out the parameters for division.
	pow = bc_num_pow10[dig];
	dig = bc_num_pow10[BC_BASE_DIGS - dig];

	// Run a series of divisions and mods with carries across the entire number
	// array. This effectively shifts everything over.
	for (i = len - 1; i < len; --i) {
		BcBigDig in, temp;
		in = ((BcBigDig) ptr[i]);
		temp = carry * dig;
		carry = in % pow;
		ptr[i] = ((BcDig) (in / pow)) + (BcDig) temp;
	}

	assert(!carry);
}

/**
 * Shift a number left by a certain number of places. This is the workhorse of
 * the left shift operator.
 * @param n       The number to shift left.
 * @param places  The amount of places to shift @a n left by.
 */
static void bc_num_shiftLeft(BcNum *restrict n, size_t places) {

	BcBigDig dig;
	size_t places_rdx;
	bool shift;

	if (!places) return;

	// Make sure to grow the number if necessary.
	if (places > n->scale) {
		size_t size = bc_vm_growSize(BC_NUM_RDX(places - n->scale), n->len);
		if (size > SIZE_MAX - 1) bc_err(BC_ERR_MATH_OVERFLOW);
	}

	// If zero, we can just set the scale and bail.
	if (BC_NUM_ZERO(n)) {
		if (n->scale >= places) n->scale -= places;
		else n->scale = 0;
		return;
	}

	// When I changed bc to have multiple digits per limb, this was the hardest
	// code to change. This and shift right. Make sure you understand this
	// before attempting anything.

	// This is how many limbs we will shift.
	dig = (BcBigDig) (places % BC_BASE_DIGS);
	shift = (dig != 0);

	// Convert places to a rdx value.
	places_rdx = BC_NUM_RDX(places);

	// If the number is not an integer, we need special care. The reason an
	// integer doesn't is because left shift would only extend the integer,
	// whereas a non-integer might have its fractional part eliminated or only
	// partially eliminated.
	if (n->scale) {

		size_t nrdx = BC_NUM_RDX_VAL(n);

		// If the number's rdx is bigger, that's the hard case.
		if (nrdx >= places_rdx) {

			size_t mod = n->scale % BC_BASE_DIGS, revdig;

			// We want mod to be in the range [1, BC_BASE_DIGS], not
			// [0, BC_BASE_DIGS).
			mod = mod ? mod : BC_BASE_DIGS;

			// We need to reverse dig to get the actual number of digits.
			revdig = dig ? BC_BASE_DIGS - dig : 0;

			// If the two overflow BC_BASE_DIGS, we need to move an extra place.
			if (mod + revdig > BC_BASE_DIGS) places_rdx = 1;
			else places_rdx = 0;
		}
		else places_rdx -= nrdx;
	}

	// If this is non-zero, we need an extra place, so expand, move, and set.
	if (places_rdx) {
		bc_num_expand(n, bc_vm_growSize(n->len, places_rdx));
		memmove(n->num + places_rdx, n->num, BC_NUM_SIZE(n->len));
		memset(n->num, 0, BC_NUM_SIZE(places_rdx));
		n->len += places_rdx;
	}

	// Set the scale appropriately.
	if (places > n->scale) {
		n->scale = 0;
		BC_NUM_RDX_SET(n, 0);
	}
	else {
		n->scale -= places;
		BC_NUM_RDX_SET(n, BC_NUM_RDX(n->scale));
	}

	// Finally, shift within limbs.
	if (shift) bc_num_shift(n, BC_BASE_DIGS - dig);

	bc_num_clean(n);
}

void bc_num_shiftRight(BcNum *restrict n, size_t places) {

	BcBigDig dig;
	size_t places_rdx, scale, scale_mod, int_len, expand;
	bool shift;

	if (!places) return;

	// If zero, we can just set the scale and bail.
	if (BC_NUM_ZERO(n)) {
		n->scale += places;
		bc_num_expand(n, BC_NUM_RDX(n->scale));
		return;
	}

	// Amount within a limb we have to shift by.
	dig = (BcBigDig) (places % BC_BASE_DIGS);
	shift = (dig != 0);

	scale = n->scale;

	// Figure out how the scale is affected.
	scale_mod = scale % BC_BASE_DIGS;
	scale_mod = scale_mod ? scale_mod : BC_BASE_DIGS;

	// We need to know the int length and rdx for places.
	int_len = bc_num_int(n);
	places_rdx = BC_NUM_RDX(places);

	// If we are going to shift past a limb boundary or not, set accordingly.
	if (scale_mod + dig > BC_BASE_DIGS) {
		expand = places_rdx - 1;
		places_rdx = 1;
	}
	else {
		expand = places_rdx;
		places_rdx = 0;
	}

	// Clamp expanding.
	if (expand > int_len) expand -= int_len;
	else expand = 0;

	// Extend, expand, and zero.
	bc_num_extend(n, places_rdx * BC_BASE_DIGS);
	bc_num_expand(n, bc_vm_growSize(expand, n->len));
	memset(n->num + n->len, 0, BC_NUM_SIZE(expand));

	// Set the fields.
	n->len += expand;
	n->scale = 0;
	BC_NUM_RDX_SET(n, 0);

	// Finally, shift within limbs.
	if (shift) bc_num_shift(n, dig);

	n->scale = scale + places;
	BC_NUM_RDX_SET(n, BC_NUM_RDX(n->scale));

	bc_num_clean(n);

	assert(BC_NUM_RDX_VAL(n) <= n->len && n->len <= n->cap);
	assert(BC_NUM_RDX_VAL(n) == BC_NUM_RDX(n->scale));
}

/**
 * Invert @a into @a b at the current scale.
 * @param a      The number to invert.
 * @param b      The return parameter. This must be preallocated.
 * @param scale  The current scale.
 */
static inline void bc_num_inv(BcNum *a, BcNum *b, size_t scale) {
	assert(BC_NUM_NONZERO(a));
	bc_num_div(&vm.one, a, b, scale);
}

/**
 * Tests if a number is a integer with scale or not. Returns true if the number
 * is not an integer. If it is, its integer shifted form is copied into the
 * result parameter for use where only integers are allowed.
 * @param n  The integer to test and shift.
 * @param r  The number to store the shifted result into. This number should
 *           *not* be allocated.
 * @return   True if the number is a non-integer, false otherwise.
 */
static bool bc_num_nonInt(const BcNum *restrict n, BcNum *restrict r) {

	bool zero;
	size_t i, rdx = BC_NUM_RDX_VAL(n);

	if (!rdx) {
		memcpy(r, n, sizeof(BcNum));
		return false;
	}

	zero = true;

	for (i = 0; zero && i < rdx; ++i) zero = (n->num[i] == 0);

	if (BC_ERR(!zero)) return true;

	bc_num_shiftRdx(n, r);

	return false;
}

#if BC_ENABLE_EXTRA_MATH

/**
 * Execute common code for an operater that needs an integer for the second
 * operand and return the integer operand as a BcBigDig.
 * @param a  The first operand.
 * @param b  The second operand.
 * @param c  The result operand.
 * @return   The second operand as a hardware integer.
 */
static BcBigDig bc_num_intop(const BcNum *a, const BcNum *b, BcNum *restrict c)
{
	BcNum temp;

	if (BC_ERR(bc_num_nonInt(b, &temp))) bc_err(BC_ERR_MATH_NON_INTEGER);

	bc_num_copy(c, a);

	return bc_num_bigdig(&temp);
}
#endif // BC_ENABLE_EXTRA_MATH

/**
 * This is the actual implementation of add *and* subtract. Since this function
 * doesn't need to use scale (per the bc spec), I am hijacking it to say whether
 * it's doing an add or a subtract. And then I convert substraction to addition
 * of negative second operand. This is a BcNumBinOp function.
 * @param a    The first operand.
 * @param b    The second operand.
 * @param c    The return parameter.
 * @param sub  Non-zero for a subtract, zero for an add.
 */
static void bc_num_as(BcNum *a, BcNum *b, BcNum *restrict c, size_t sub) {

	BcDig *ptr_c, *ptr_l, *ptr_r;
	size_t i, min_rdx, max_rdx, diff, a_int, b_int, min_len, max_len, max_int;
	size_t len_l, len_r, ardx, brdx;
	bool b_neg, do_sub, do_rev_sub, carry, c_neg;

	if (BC_NUM_ZERO(b)) {
		bc_num_copy(c, a);
		return;
	}

	if (BC_NUM_ZERO(a)) {
		bc_num_copy(c, b);
		c->rdx = BC_NUM_NEG_VAL(c, BC_NUM_NEG(b) != sub);
		return;
	}

	// Invert sign of b if it is to be subtracted. This operation must
	// precede the tests for any of the operands being zero.
	b_neg = (BC_NUM_NEG(b) != sub);

	// Figure out if we will actually add the numbers if their signs are equal
	// or subtract.
	do_sub = (BC_NUM_NEG(a) != b_neg);

	a_int = bc_num_int(a);
	b_int = bc_num_int(b);
	max_int = BC_MAX(a_int, b_int);

	// Figure out which number will have its last limbs copied (for addition) or
	// subtracted (for subtraction).
	ardx = BC_NUM_RDX_VAL(a);
	brdx = BC_NUM_RDX_VAL(b);
	min_rdx = BC_MIN(ardx, brdx);
	max_rdx = BC_MAX(ardx, brdx);
	diff = max_rdx - min_rdx;

	max_len = max_int + max_rdx;

	if (do_sub) {

		// Check whether b has to be subtracted from a or a from b.
		if (a_int != b_int) do_rev_sub = (a_int < b_int);
		else if (ardx > brdx)
			do_rev_sub = (bc_num_compare(a->num + diff, b->num, b->len) < 0);
		else
			do_rev_sub = (bc_num_compare(a->num, b->num + diff, a->len) <= 0);
	}
	else {

		// The result array of the addition might come out one element
		// longer than the bigger of the operand arrays.
		max_len += 1;
		do_rev_sub = (a_int < b_int);
	}

	assert(max_len <= c->cap);

	// Cache values for simple code later.
	if (do_rev_sub) {
		ptr_l = b->num;
		ptr_r = a->num;
		len_l = b->len;
		len_r = a->len;
	}
	else {
		ptr_l = a->num;
		ptr_r = b->num;
		len_l = a->len;
		len_r = b->len;
	}

	ptr_c = c->num;
	carry = false;

	// This is true if the numbers have a different number of limbs after the
	// decimal point.
	if (diff) {

		// If the rdx values of the operands do not match, the result will
		// have low end elements that are the positive or negative trailing
		// elements of the operand with higher rdx value.
		if ((ardx > brdx) != do_rev_sub) {

			// !do_rev_sub && ardx > brdx || do_rev_sub && brdx > ardx
			// The left operand has BcDig values that need to be copied,
			// either from a or from b (in case of a reversed subtraction).
			memcpy(ptr_c, ptr_l, BC_NUM_SIZE(diff));
			ptr_l += diff;
			len_l -= diff;
		}
		else {

			// The right operand has BcDig values that need to be copied
			// or subtracted from zero (in case of a subtraction).
			if (do_sub) {

				// do_sub (do_rev_sub && ardx > brdx ||
				// !do_rev_sub && brdx > ardx)
				for (i = 0; i < diff; i++)
					ptr_c[i] = bc_num_subDigits(0, ptr_r[i], &carry);
			}
			else {

				// !do_sub && brdx > ardx
				memcpy(ptr_c, ptr_r, BC_NUM_SIZE(diff));
			}

			// Future code needs to ignore the limbs we just did.
			ptr_r += diff;
			len_r -= diff;
		}

		// The return value pointer needs to ignore what we just did.
		ptr_c += diff;
	}

	// This is the length that can be directly added/subtracted.
	min_len = BC_MIN(len_l, len_r);

	// After dealing with possible low array elements that depend on only one
	// operand above, the actual add or subtract can be performed as if the rdx
	// of both operands was the same.
	//
	// Inlining takes care of eliminating constant zero arguments to
	// addDigit/subDigit (checked in disassembly of resulting bc binary
	// compiled with gcc and clang).
	if (do_sub) {

		// Actual subtraction.
		for (i = 0; i < min_len; ++i)
			ptr_c[i] = bc_num_subDigits(ptr_l[i], ptr_r[i], &carry);

		// Finishing the limbs beyond the direct subtraction.
		for (; i < len_l; ++i) ptr_c[i] = bc_num_subDigits(ptr_l[i], 0, &carry);
	}
	else {

		// Actual addition.
		for (i = 0; i < min_len; ++i)
			ptr_c[i] = bc_num_addDigits(ptr_l[i], ptr_r[i], &carry);

		// Finishing the limbs beyond the direct addition.
		for (; i < len_l; ++i) ptr_c[i] = bc_num_addDigits(ptr_l[i], 0, &carry);

		// Addition can create an extra limb. We take care of that here.
		ptr_c[i] = bc_num_addDigits(0, 0, &carry);
	}

	assert(carry == false);

	// The result has the same sign as a, unless the operation was a
	// reverse subtraction (b - a).
	c_neg = BC_NUM_NEG(a) != (do_sub && do_rev_sub);
	BC_NUM_RDX_SET_NEG(c, max_rdx, c_neg);
	c->len = max_len;
	c->scale = BC_MAX(a->scale, b->scale);

	bc_num_clean(c);
}

/**
 * The simple multiplication that karatsuba dishes out to when the length of the
 * numbers gets low enough. This doesn't use scale because it treats the
 * operands as though they are integers.
 * @param a  The first operand.
 * @param b  The second operand.
 * @param c  The return parameter.
 */
static void bc_num_m_simp(const BcNum *a, const BcNum *b, BcNum *restrict c) {

	size_t i, alen = a->len, blen = b->len, clen;
	BcDig *ptr_a = a->num, *ptr_b = b->num, *ptr_c;
	BcBigDig sum = 0, carry = 0;

	assert(sizeof(sum) >= sizeof(BcDig) * 2);
	assert(!BC_NUM_RDX_VAL(a) && !BC_NUM_RDX_VAL(b));

	// Make sure c is big enough.
	clen = bc_vm_growSize(alen, blen);
	bc_num_expand(c, bc_vm_growSize(clen, 1));

	// If we don't memset, then we might have uninitialized data use later.
	ptr_c = c->num;
	memset(ptr_c, 0, BC_NUM_SIZE(c->cap));

	// This is the actual multiplication loop. It uses the lattice form of long
	// multiplication (see the explanation on the web page at
	// https://knilt.arcc.albany.edu/What_is_Lattice_Multiplication or the
	// explanation at Wikipedia).
	for (i = 0; i < clen; ++i) {

		ssize_t sidx = (ssize_t) (i - blen + 1);
		size_t j, k;

		// These are the start indices.
		j = (size_t) BC_MAX(0, sidx);
		k = BC_MIN(i, blen - 1);

		// On every iteration of this loop, a multiplication happens, then the
		// sum is automatically calculated.
		for (; j < alen && k < blen; ++j, --k) {

			sum += ((BcBigDig) ptr_a[j]) * ((BcBigDig) ptr_b[k]);

			if (sum >= ((BcBigDig) BC_BASE_POW) * BC_BASE_POW) {
				carry += sum / BC_BASE_POW;
				sum %= BC_BASE_POW;
			}
		}

		// Calculate the carry.
		if (sum >= BC_BASE_POW) {
			carry += sum / BC_BASE_POW;
			sum %= BC_BASE_POW;
		}

		// Store and set up for next iteration.
		ptr_c[i] = (BcDig) sum;
		assert(ptr_c[i] < BC_BASE_POW);
		sum = carry;
		carry = 0;
	}

	// This should always be true because there should be no carry on the last
	// digit; multiplication never goes above the sum of both lengths.
	assert(!sum);

	c->len = clen;
}

/**
 * Does a shifted add or subtract for Karatsuba below. This calls either
 * bc_num_addArrays() or bc_num_subArrays().
 * @param n      An in/out parameter; the first operand and return parameter.
 * @param a      The second operand.
 * @param shift  The amount to shift @a n by when adding/subtracting.
 * @param op     The function to call, either bc_num_addArrays() or
 *               bc_num_subArrays().
 */
static void bc_num_shiftAddSub(BcNum *restrict n, const BcNum *restrict a,
                               size_t shift, BcNumShiftAddOp op)
{
	assert(n->len >= shift + a->len);
	assert(!BC_NUM_RDX_VAL(n) && !BC_NUM_RDX_VAL(a));
	op(n->num + shift, a->num, a->len);
}

/**
 * Implements the Karatsuba algorithm.
 */
static void bc_num_k(const BcNum *a, const BcNum *b, BcNum *restrict c) {

	size_t max, max2, total;
	BcNum l1, h1, l2, h2, m2, m1, z0, z1, z2, temp;
	BcDig *digs, *dig_ptr;
	BcNumShiftAddOp op;
	bool aone = BC_NUM_ONE(a);

	assert(BC_NUM_ZERO(c));

	if (BC_NUM_ZERO(a) || BC_NUM_ZERO(b)) return;

	if (aone || BC_NUM_ONE(b)) {
		bc_num_copy(c, aone ? b : a);
		if ((aone && BC_NUM_NEG(a)) || BC_NUM_NEG(b)) BC_NUM_NEG_TGL(c);
		return;
	}

	// Shell out to the simple algorithm with certain conditions.
	if (a->len < BC_NUM_KARATSUBA_LEN || b->len < BC_NUM_KARATSUBA_LEN) {
		bc_num_m_simp(a, b, c);
		return;
	}

	// We need to calculate the max size of the numbers that can result from the
	// operations.
	max = BC_MAX(a->len, b->len);
	max = BC_MAX(max, BC_NUM_DEF_SIZE);
	max2 = (max + 1) / 2;

	// Calculate the space needed for all of the temporary allocations. We do
	// this to just allocate once.
	total = bc_vm_arraySize(BC_NUM_KARATSUBA_ALLOCS, max);

	BC_SIG_LOCK;

	// Allocate space for all of the temporaries.
	digs = dig_ptr = bc_vm_malloc(BC_NUM_SIZE(total));

	// Set up the temporaries.
	bc_num_setup(&l1, dig_ptr, max);
	dig_ptr += max;
	bc_num_setup(&h1, dig_ptr, max);
	dig_ptr += max;
	bc_num_setup(&l2, dig_ptr, max);
	dig_ptr += max;
	bc_num_setup(&h2, dig_ptr, max);
	dig_ptr += max;
	bc_num_setup(&m1, dig_ptr, max);
	dig_ptr += max;
	bc_num_setup(&m2, dig_ptr, max);

	// Some temporaries need the ability to grow, so we allocate them
	// separately.
	max = bc_vm_growSize(max, 1);
	bc_num_init(&z0, max);
	bc_num_init(&z1, max);
	bc_num_init(&z2, max);
	max = bc_vm_growSize(max, max) + 1;
	bc_num_init(&temp, max);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	// First, set up c.
	bc_num_expand(c, max);
	c->len = max;
	memset(c->num, 0, BC_NUM_SIZE(c->len));

	// Split the parameters.
	bc_num_split(a, max2, &l1, &h1);
	bc_num_split(b, max2, &l2, &h2);

	// Do the subtraction.
	bc_num_sub(&h1, &l1, &m1, 0);
	bc_num_sub(&l2, &h2, &m2, 0);

	// The if statements below are there for efficiency reasons. The best way to
	// understand them is to understand the Karatsuba algorithm because now that
	// the ollocations and splits are done, the algorithm is pretty
	// straightforward.

	if (BC_NUM_NONZERO(&h1) && BC_NUM_NONZERO(&h2)) {

		assert(BC_NUM_RDX_VALID_NP(h1));
		assert(BC_NUM_RDX_VALID_NP(h2));

		bc_num_m(&h1, &h2, &z2, 0);
		bc_num_clean(&z2);

		bc_num_shiftAddSub(c, &z2, max2 * 2, bc_num_addArrays);
		bc_num_shiftAddSub(c, &z2, max2, bc_num_addArrays);
	}

	if (BC_NUM_NONZERO(&l1) && BC_NUM_NONZERO(&l2)) {

		assert(BC_NUM_RDX_VALID_NP(l1));
		assert(BC_NUM_RDX_VALID_NP(l2));

		bc_num_m(&l1, &l2, &z0, 0);
		bc_num_clean(&z0);

		bc_num_shiftAddSub(c, &z0, max2, bc_num_addArrays);
		bc_num_shiftAddSub(c, &z0, 0, bc_num_addArrays);
	}

	if (BC_NUM_NONZERO(&m1) && BC_NUM_NONZERO(&m2)) {

		assert(BC_NUM_RDX_VALID_NP(m1));
		assert(BC_NUM_RDX_VALID_NP(m1));

		bc_num_m(&m1, &m2, &z1, 0);
		bc_num_clean(&z1);

		op = (BC_NUM_NEG_NP(m1) != BC_NUM_NEG_NP(m2)) ?
		     bc_num_subArrays : bc_num_addArrays;
		bc_num_shiftAddSub(c, &z1, max2, op);
	}

err:
	BC_SIG_MAYLOCK;
	free(digs);
	bc_num_free(&temp);
	bc_num_free(&z2);
	bc_num_free(&z1);
	bc_num_free(&z0);
	BC_LONGJMP_CONT;
}

/**
 * Does checks for Karatsuba. It also changes things to ensure that the
 * Karatsuba and simple multiplication can treat the numbers as integers. This
 * is a BcNumBinOp function.
 * @param a      The first operand.
 * @param b      The second operand.
 * @param c      The return parameter.
 * @param scale  The current scale.
 */
static void bc_num_m(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {

	BcNum cpa, cpb;
	size_t ascale, bscale, ardx, brdx, azero = 0, bzero = 0, zero, len, rscale;

	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));

	bc_num_zero(c);

	ascale = a->scale;
	bscale = b->scale;

	// This sets the final scale according to the bc spec.
	scale = BC_MAX(scale, ascale);
	scale = BC_MAX(scale, bscale);
	rscale = ascale + bscale;
	scale = BC_MIN(rscale, scale);

	// If this condition is true, we can use bc_num_mulArray(), which would be
	// much faster.
	if ((a->len == 1 || b->len == 1) && !a->rdx && !b->rdx) {

		BcNum *operand;
		BcBigDig dig;

		// Set the correct operands.
		if (a->len == 1) {
			dig = (BcBigDig) a->num[0];
			operand = b;
		}
		else {
			dig = (BcBigDig) b->num[0];
			operand = a;
		}

		bc_num_mulArray(operand, dig, c);

		// Need to make sure the sign is correct.
		if (BC_NUM_NONZERO(c))
			c->rdx = BC_NUM_NEG_VAL(c, BC_NUM_NEG(a) != BC_NUM_NEG(b));

		return;
	}

	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));

	BC_SIG_LOCK;

	// We need copies because of all of the mutation needed to make Karatsuba
	// think the numbers are integers.
	bc_num_init(&cpa, a->len + BC_NUM_RDX_VAL(a));
	bc_num_init(&cpb, b->len + BC_NUM_RDX_VAL(b));

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	bc_num_copy(&cpa, a);
	bc_num_copy(&cpb, b);

	assert(BC_NUM_RDX_VALID_NP(cpa));
	assert(BC_NUM_RDX_VALID_NP(cpb));

	BC_NUM_NEG_CLR_NP(cpa);
	BC_NUM_NEG_CLR_NP(cpb);

	assert(BC_NUM_RDX_VALID_NP(cpa));
	assert(BC_NUM_RDX_VALID_NP(cpb));

	// These are what makes them appear like integers.
	ardx = BC_NUM_RDX_VAL_NP(cpa) * BC_BASE_DIGS;
	bc_num_shiftLeft(&cpa, ardx);

	brdx = BC_NUM_RDX_VAL_NP(cpb) * BC_BASE_DIGS;
	bc_num_shiftLeft(&cpb, brdx);

	// We need to reset the jump here because azero and bzero are used in the
	// cleanup, and local variables are not guaranteed to be the same after a
	// jump.
	BC_SIG_LOCK;

	BC_UNSETJMP;

	// We want to ignore zero limbs.
	azero = bc_num_shiftZero(&cpa);
	bzero = bc_num_shiftZero(&cpb);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	bc_num_clean(&cpa);
	bc_num_clean(&cpb);

	bc_num_k(&cpa, &cpb, c);

	// The return parameter needs to have its scale set. This is the start. It
	// also needs to be shifted by the same amount as a and b have limbs after
	// the decimal point.
	zero = bc_vm_growSize(azero, bzero);
	len = bc_vm_growSize(c->len, zero);

	bc_num_expand(c, len);

	// Shift c based on the limbs after the decimal point in a and b.
	bc_num_shiftLeft(c, (len - c->len) * BC_BASE_DIGS);
	bc_num_shiftRight(c, ardx + brdx);

	bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));

err:
	BC_SIG_MAYLOCK;
	bc_num_unshiftZero(&cpb, bzero);
	bc_num_free(&cpb);
	bc_num_unshiftZero(&cpa, azero);
	bc_num_free(&cpa);
	BC_LONGJMP_CONT;
}

/**
 * Returns true if the BcDig array has non-zero limbs, false otherwise.
 * @param a    The array to test.
 * @param len  The length of the array.
 * @return     True if @a has any non-zero limbs, false otherwise.
 */
static bool bc_num_nonZeroDig(BcDig *restrict a, size_t len) {
	size_t i;
	bool nonzero = false;
	for (i = len - 1; !nonzero && i < len; --i) nonzero = (a[i] != 0);
	return nonzero;
}

/**
 * Compares a BcDig array against a BcNum. This is especially suited for
 * division. Returns >0 if @a a is greater than @a b, <0 if it is less, and =0
 * if they are equal.
 * @param a    The array.
 * @param b    The number.
 * @param len  The length to assume the arrays are. This is always less than the
 *             actual length because of how this is implemented.
 */
static ssize_t bc_num_divCmp(const BcDig *a, const BcNum *b, size_t len) {

	ssize_t cmp;

	if (b->len > len && a[len]) cmp = bc_num_compare(a, b->num, len + 1);
	else if (b->len <= len) {
		if (a[len]) cmp = 1;
		else cmp = bc_num_compare(a, b->num, len);
	}
	else cmp = -1;

	return cmp;
}

/**
 * Extends the two operands of a division by BC_BASE_DIGS minus the number of
 * digits in the divisor estimate. In other words, it is shifting the numbers in
 * order to force the divisor estimate to fill the limb.
 * @param a        The first operand.
 * @param b        The second operand.
 * @param divisor  The divisor estimate.
 */
static void bc_num_divExtend(BcNum *restrict a, BcNum *restrict b,
                             BcBigDig divisor)
{
	size_t pow;

	assert(divisor < BC_BASE_POW);

	pow = BC_BASE_DIGS - bc_num_log10((size_t) divisor);

	bc_num_shiftLeft(a, pow);
	bc_num_shiftLeft(b, pow);
}

/**
 * Actually does division. This is a rewrite of my original code by Stefan Esser
 * from FreeBSD.
 * @param a      The first operand.
 * @param b      The second operand.
 * @param c      The return parameter.
 * @param scale  The current scale.
 */
static void bc_num_d_long(BcNum *restrict a, BcNum *restrict b,
                          BcNum *restrict c, size_t scale)
{
	BcBigDig divisor;
	size_t len, end, i, rdx;
	BcNum cpb;
	bool nonzero = false;

	assert(b->len < a->len);

	len = b->len;
	end = a->len - len;

	assert(len >= 1);

	// This is a final time to make sure c is big enough and that its array is
	// properly zeroed.
	bc_num_expand(c, a->len);
	memset(c->num, 0, c->cap * sizeof(BcDig));

	// Setup.
	BC_NUM_RDX_SET(c, BC_NUM_RDX_VAL(a));
	c->scale = a->scale;
	c->len = a->len;

	// This is pulling the most significant limb of b in order to establish a
	// good "estimate" for the actual divisor.
	divisor = (BcBigDig) b->num[len - 1];

	// The entire bit of code in this if statement is to tighten the estimate of
	// the divisor. The condition asks if b has any other non-zero limbs.
	if (len > 1 && bc_num_nonZeroDig(b->num, len - 1)) {

		// This takes a little bit of understanding. The "10*BC_BASE_DIGS/6+1"
		// results in either 16 for 64-bit 9-digit limbs or 7 for 32-bit 4-digit
		// limbs. Then it shifts a 1 by that many, which in both cases, puts the
		// result above *half* of the max value a limb can store. Basically,
		// this quickly calculates if the divisor is greater than half the max
		// of a limb.
		nonzero = (divisor > 1 << ((10 * BC_BASE_DIGS) / 6 + 1));

		// If the divisor is *not* greater than half the limb...
		if (!nonzero) {

			// Extend the parameters by the number of missing digits in the
			// divisor.
			bc_num_divExtend(a, b, divisor);

			// Check bc_num_d(). In there, we grow a again and again. We do it
			// again here; we *always* want to be sure it is big enough.
			len = BC_MAX(a->len, b->len);
			bc_num_expand(a, len + 1);

			// Make a have a zero most significant limb to match the len.
			if (len + 1 > a->len) a->len = len + 1;

			// Grab the new divisor estimate, new because the shift has made it
			// different.
			len = b->len;
			end = a->len - len;
			divisor = (BcBigDig) b->num[len - 1];

			nonzero = bc_num_nonZeroDig(b->num, len - 1);
		}
	}

	// If b has other nonzero limbs, we want the divisor to be one higher, so
	// that it is an upper bound.
	divisor += nonzero;

	// Make sure c can fit the new length.
	bc_num_expand(c, a->len);
	memset(c->num, 0, BC_NUM_SIZE(c->cap));

	assert(c->scale >= scale);
	rdx = BC_NUM_RDX_VAL(c) - BC_NUM_RDX(scale);

	BC_SIG_LOCK;

	bc_num_init(&cpb, len + 1);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	// This is the actual division loop.
	for (i = end - 1; i < end && i >= rdx && BC_NUM_NONZERO(a); --i) {

		ssize_t cmp;
		BcDig *n;
		BcBigDig result;

		n = a->num + i;
		assert(n >= a->num);
		result = 0;

		cmp = bc_num_divCmp(n, b, len);

		// This is true if n is greater than b, which means that division can
		// proceed, so this inner loop is the part that implements one instance
		// of the division.
		while (cmp >= 0) {

			BcBigDig n1, dividend, quotient;

			// These should be named obviously enough. Just imagine that it's a
			// division of one limb. Because that's what it is.
			n1 = (BcBigDig) n[len];
			dividend = n1 * BC_BASE_POW + (BcBigDig) n[len - 1];
			quotient = (dividend / divisor);

			// If this is true, then we can just subtract. Remember: setting
			// quotient to 1 is not bad because we already know that n is
			// greater than b.
			if (quotient <= 1) {
				quotient = 1;
				bc_num_subArrays(n, b->num, len);
			}
			else {

				assert(quotient <= BC_BASE_POW);

				// We need to multiply and subtract for a quotient above 1.
				bc_num_mulArray(b, (BcBigDig) quotient, &cpb);
				bc_num_subArrays(n, cpb.num, cpb.len);
			}

			// The result is the *real* quotient, by the way, but it might take
			// multiple trips around this loop to get it.
			result += quotient;
			assert(result <= BC_BASE_POW);

			// And here's why it might take multiple trips: n might *still* be
			// greater than b. So we have to loop again. That's what this is
			// setting up for: the condition of the while loop.
			if (nonzero) cmp = bc_num_divCmp(n, b, len);
			else cmp = -1;
		}

		assert(result < BC_BASE_POW);

		// Store the actual limb quotient.
		c->num[i] = (BcDig) result;
	}

err:
	BC_SIG_MAYLOCK;
	bc_num_free(&cpb);
	BC_LONGJMP_CONT;
}

/**
 * Implements division. This is a BcNumBinOp function.
 * @param a      The first operand.
 * @param b      The second operand.
 * @param c      The return parameter.
 * @param scale  The current scale.
 */
static void bc_num_d(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {

	size_t len, cpardx;
	BcNum cpa, cpb;

	if (BC_NUM_ZERO(b)) bc_err(BC_ERR_MATH_DIVIDE_BY_ZERO);

	if (BC_NUM_ZERO(a)) {
		bc_num_setToZero(c, scale);
		return;
	}

	if (BC_NUM_ONE(b)) {
		bc_num_copy(c, a);
		bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));
		return;
	}

	// If this is true, we can use bc_num_divArray(), which would be faster.
	if (!BC_NUM_RDX_VAL(a) && !BC_NUM_RDX_VAL(b) && b->len == 1 && !scale) {
		BcBigDig rem;
		bc_num_divArray(a, (BcBigDig) b->num[0], c, &rem);
		bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));
		return;
	}

	len = bc_num_divReq(a, b, scale);

	BC_SIG_LOCK;

	// Initialize copies of the parameters. We want the length of the first
	// operand copy to be as big as the result because of the way the division
	// is implemented.
	bc_num_init(&cpa, len);
	bc_num_copy(&cpa, a);
	bc_num_createCopy(&cpb, b);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	len = b->len;

	// Like the above comment, we want the copy of the first parameter to be
	// larger than the second parameter.
	if (len > cpa.len) {
		bc_num_expand(&cpa, bc_vm_growSize(len, 2));
		bc_num_extend(&cpa, (len - cpa.len) * BC_BASE_DIGS);
	}

	cpardx = BC_NUM_RDX_VAL_NP(cpa);
	cpa.scale = cpardx * BC_BASE_DIGS;

	// This is just setting up the scale in preparation for the division.
	bc_num_extend(&cpa, b->scale);
	cpardx = BC_NUM_RDX_VAL_NP(cpa) - BC_NUM_RDX(b->scale);
	BC_NUM_RDX_SET_NP(cpa, cpardx);
	cpa.scale = cpardx * BC_BASE_DIGS;

	// Once again, just setting things up, this time to match scale.
	if (scale > cpa.scale) {
		bc_num_extend(&cpa, scale);
		cpardx = BC_NUM_RDX_VAL_NP(cpa);
		cpa.scale = cpardx * BC_BASE_DIGS;
	}

	// Grow if necessary.
	if (cpa.cap == cpa.len) bc_num_expand(&cpa, bc_vm_growSize(cpa.len, 1));

	// We want an extra zero in front to make things simpler.
	cpa.num[cpa.len++] = 0;

	// Still setting things up. Why all of these things are needed is not
	// something that can be easily explained, but it has to do with making the
	// actual algorithm easier to understand because it can assume a lot of
	// things. Thus, you should view all of this setup code as establishing
	// assumptions for bc_num_d_long(), where the actual division happens.
	if (cpardx == cpa.len) cpa.len = bc_num_nonZeroLen(&cpa);
	if (BC_NUM_RDX_VAL_NP(cpb) == cpb.len) cpb.len = bc_num_nonZeroLen(&cpb);
	cpb.scale = 0;
	BC_NUM_RDX_SET_NP(cpb, 0);

	bc_num_d_long(&cpa, &cpb, c, scale);

	bc_num_retireMul(c, scale, BC_NUM_NEG(a), BC_NUM_NEG(b));

err:
	BC_SIG_MAYLOCK;
	bc_num_free(&cpb);
	bc_num_free(&cpa);
	BC_LONGJMP_CONT;
}

/**
 * Implements divmod. This is the actual modulus function; since modulus
 * requires a division anyway, this returns the quotient and modulus. Either can
 * be thrown out as desired.
 * @param a      The first operand.
 * @param b      The second operand.
 * @param c      The return parameter for the quotient.
 * @param d      The return parameter for the modulus.
 * @param scale  The current scale.
 * @param ts     The scale that the operation should be done to. Yes, it's not
 *               necessarily the same as scale, per the bc spec.
 */
static void bc_num_r(BcNum *a, BcNum *b, BcNum *restrict c,
                     BcNum *restrict d, size_t scale, size_t ts)
{
	BcNum temp;
	bool neg;

	if (BC_NUM_ZERO(b)) bc_err(BC_ERR_MATH_DIVIDE_BY_ZERO);

	if (BC_NUM_ZERO(a)) {
		bc_num_setToZero(c, ts);
		bc_num_setToZero(d, ts);
		return;
	}

	BC_SIG_LOCK;

	bc_num_init(&temp, d->cap);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	// Division.
	bc_num_d(a, b, c, scale);

	// We want an extra digit so we can safely truncate.
	if (scale) scale = ts + 1;

	assert(BC_NUM_RDX_VALID(c));
	assert(BC_NUM_RDX_VALID(b));

	// Implement the rest of the (a - (a / b) * b) formula.
	bc_num_m(c, b, &temp, scale);
	bc_num_sub(a, &temp, d, scale);

	// Extend if necessary.
	if (ts > d->scale && BC_NUM_NONZERO(d)) bc_num_extend(d, ts - d->scale);

	neg = BC_NUM_NEG(d);
	bc_num_retireMul(d, ts, BC_NUM_NEG(a), BC_NUM_NEG(b));
	d->rdx = BC_NUM_NEG_VAL(d, BC_NUM_NONZERO(d) ? neg : false);

err:
	BC_SIG_MAYLOCK;
	bc_num_free(&temp);
	BC_LONGJMP_CONT;
}

/**
 * Implements modulus/remainder. (Yes, I know they are different, but not in the
 * context of bc.) This is a BcNumBinOp function.
 * @param a      The first operand.
 * @param b      The second operand.
 * @param c      The return parameter.
 * @param scale  The current scale.
 */
static void bc_num_rem(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {

	BcNum c1;
	size_t ts;

	ts = bc_vm_growSize(scale, b->scale);
	ts = BC_MAX(ts, a->scale);

	BC_SIG_LOCK;

	// Need a temp for the quotient.
	bc_num_init(&c1, bc_num_mulReq(a, b, ts));

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	bc_num_r(a, b, &c1, c, scale, ts);

err:
	BC_SIG_MAYLOCK;
	bc_num_free(&c1);
	BC_LONGJMP_CONT;
}

/**
 * Implements power (exponentiation). This is a BcNumBinOp function.
 * @param a      The first operand.
 * @param b      The second operand.
 * @param c      The return parameter.
 * @param scale  The current scale.
 */
static void bc_num_p(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {

	BcNum copy, btemp;
	BcBigDig exp;
	size_t powrdx, resrdx;
	bool neg;

	if (BC_ERR(bc_num_nonInt(b, &btemp))) bc_err(BC_ERR_MATH_NON_INTEGER);

	if (BC_NUM_ZERO(&btemp)) {
		bc_num_one(c);
		return;
	}

	if (BC_NUM_ZERO(a)) {
		if (BC_NUM_NEG_NP(btemp)) bc_err(BC_ERR_MATH_DIVIDE_BY_ZERO);
		bc_num_setToZero(c, scale);
		return;
	}

	if (BC_NUM_ONE(&btemp)) {
		if (!BC_NUM_NEG_NP(btemp)) bc_num_copy(c, a);
		else bc_num_inv(a, c, scale);
		return;
	}

	neg = BC_NUM_NEG_NP(btemp);
	BC_NUM_NEG_CLR_NP(btemp);

	exp = bc_num_bigdig(&btemp);

	BC_SIG_LOCK;

	bc_num_createCopy(&copy, a);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	// If this is true, then we do not have to do a division, and we need to
	// set scale accordingly.
	if (!neg) {
		size_t max = BC_MAX(scale, a->scale), scalepow;
		scalepow = bc_num_mulOverflow(a->scale, exp);
		scale = BC_MIN(scalepow, max);
	}

	// This is only implementing the first exponentiation by squaring, until it
	// reaches the first time where the square is actually used.
	for (powrdx = a->scale; !(exp & 1); exp >>= 1) {
		powrdx <<= 1;
		assert(BC_NUM_RDX_VALID_NP(copy));
		bc_num_mul(&copy, &copy, &copy, powrdx);
	}

	// Make c a copy of copy for the purpose of saving the squares that should
	// be saved.
	bc_num_copy(c, &copy);
	resrdx = powrdx;

	// Now finish the exponentiation by squaring, this time saving the squares
	// as necessary.
	while (exp >>= 1) {

		powrdx <<= 1;
		assert(BC_NUM_RDX_VALID_NP(copy));
		bc_num_mul(&copy, &copy, &copy, powrdx);

		// If this is true, we want to save that particular square. This does
		// that by multiplying c with copy.
		if (exp & 1) {
			resrdx += powrdx;
			assert(BC_NUM_RDX_VALID(c));
			assert(BC_NUM_RDX_VALID_NP(copy));
			bc_num_mul(c, &copy, c, resrdx);
		}
	}

	// Invert if necessary.
	if (neg) bc_num_inv(c, c, scale);

	// Truncate if necessary.
	if (c->scale > scale) bc_num_truncate(c, c->scale - scale);

	bc_num_clean(c);

err:
	BC_SIG_MAYLOCK;
	bc_num_free(&copy);
	BC_LONGJMP_CONT;
}

#if BC_ENABLE_EXTRA_MATH
/**
 * Implements the places operator. This is a BcNumBinOp function.
 * @param a      The first operand.
 * @param b      The second operand.
 * @param c      The return parameter.
 * @param scale  The current scale.
 */
static void bc_num_place(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {

	BcBigDig val;

	BC_UNUSED(scale);

	val = bc_num_intop(a, b, c);

	// Just truncate or extend as appropriate.
	if (val < c->scale) bc_num_truncate(c, c->scale - val);
	else if (val > c->scale) bc_num_extend(c, val - c->scale);
}

/**
 * Implements the left shift operator. This is a BcNumBinOp function.
 */
static void bc_num_left(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {

	BcBigDig val;

	BC_UNUSED(scale);

	val = bc_num_intop(a, b, c);

	bc_num_shiftLeft(c, (size_t) val);
}

/**
 * Implements the right shift operator. This is a BcNumBinOp function.
 */
static void bc_num_right(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {

	BcBigDig val;

	BC_UNUSED(scale);

	val = bc_num_intop(a, b, c);

	if (BC_NUM_ZERO(c)) return;

	bc_num_shiftRight(c, (size_t) val);
}
#endif // BC_ENABLE_EXTRA_MATH

/**
 * Prepares for, and calls, a binary operator function. This is probably the
 * most important function in the entire file because it establishes assumptions
 * that make the rest of the code so easy. Those assumptions include:
 *
 * - a is not the same pointer as c.
 * - b is not the same pointer as c.
 * - there is enough room in c for the result.
 *
 * Without these, this whole function would basically have to be duplicated for
 * *all* binary operators.
 *
 * @param a      The first operand.
 * @param b      The second operand.
 * @param c      The return parameter.
 * @param scale  The current scale.
 * @param req    The number of limbs needed to fit the result.
 */
static void bc_num_binary(BcNum *a, BcNum *b, BcNum *c, size_t scale,
                          BcNumBinOp op, size_t req)
{
	BcNum *ptr_a, *ptr_b, num2;
	bool init = false;

	assert(a != NULL && b != NULL && c != NULL && op != NULL);

	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));

	BC_SIG_LOCK;

	// Reallocate if c == a.
	if (c == a) {

		ptr_a = &num2;

		memcpy(ptr_a, c, sizeof(BcNum));
		init = true;
	}
	else {
		ptr_a = a;
	}

	// Also reallocate if c == b.
	if (c == b) {

		ptr_b = &num2;

		if (c != a) {
			memcpy(ptr_b, c, sizeof(BcNum));
			init = true;
		}
	}
	else {
		ptr_b = b;
	}

	// Actually reallocate. If we don't reallocate, we want to expand at the
	// very least.
	if (init) {

		bc_num_init(c, req);

		// Must prepare for cleanup. We want this here so that locals that got
		// set stay set since a longjmp() is not guaranteed to preserve locals.
		BC_SETJMP_LOCKED(err);
		BC_SIG_UNLOCK;
	}
	else {
		BC_SIG_UNLOCK;
		bc_num_expand(c, req);
	}

	// It is okay for a and b to be the same. If a binary operator function does
	// need them to be different, the binary operator function is responsible
	// for that.

	// Call the actual binary operator function.
	op(ptr_a, ptr_b, c, scale);

	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
	assert(BC_NUM_RDX_VALID(c));
	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);

err:
	// Cleanup only needed if we initialized c to a new number.
	if (init) {
		BC_SIG_MAYLOCK;
		bc_num_free(&num2);
		BC_LONGJMP_CONT;
	}
}

#if !defined(NDEBUG) || BC_ENABLE_LIBRARY

/**
 * Tests a number string for validity. This function has a history; I originally
 * wrote it because I did not trust my parser. Over time, however, I came to
 * trust it, so I was able to relegate this function to debug builds only, and I
 * used it in assert()'s. But then I created the library, and well, I can't
 * trust users, so I reused this for yelling at users.
 * @param val  The string to check to see if it's a valid number string.
 * @return     True if the string is a valid number string, false otherwise.
 */
bool bc_num_strValid(const char *restrict val) {

	bool radix = false;
	size_t i, len = strlen(val);

	// Notice that I don't check if there is a negative sign. That is not part
	// of a valid number, except in the library. The library-specific code takes
	// care of that part.

	// Nothing in the string is okay.
	if (!len) return true;

	// Loop through the characters.
	for (i = 0; i < len; ++i) {

		BcDig c = val[i];

		// If we have found a radix point...
		if (c == '.') {

			// We don't allow two radices.
			if (radix) return false;

			radix = true;
			continue;
		}

		// We only allow digits and uppercase letters.
		if (!(isdigit(c) || isupper(c))) return false;
	}

	return true;
}
#endif // !defined(NDEBUG) || BC_ENABLE_LIBRARY

/**
 * Parses one character and returns the digit that corresponds to that
 * character according to the base.
 * @param c     The character to parse.
 * @param base  The base.
 * @return      The character as a digit.
 */
static BcBigDig bc_num_parseChar(char c, size_t base) {

	assert(isupper(c) || isdigit(c));

	// If a letter...
	if (isupper(c)) {

		// This returns the digit that directly corresponds with the letter.
		c = BC_NUM_NUM_LETTER(c);

		// If the digit is greater than the base, we clamp.
		c = ((size_t) c) >= base ? (char) base - 1 : c;
	}
	// Straight convert the digit to a number.
	else c -= '0';

	return (BcBigDig) (uchar) c;
}

/**
 * Parses a string as a decimal number. This is separate because it's going to
 * be the most used, and it can be heavily optimized for decimal only.
 * @param n    The number to parse into and return. Must be preallocated.
 * @param val  The string to parse.
 */
static void bc_num_parseDecimal(BcNum *restrict n, const char *restrict val) {

	size_t len, i, temp, mod;
	const char *ptr;
	bool zero = true, rdx;

	// Eat leading zeroes.
	for (i = 0; val[i] == '0'; ++i);

	val += i;
	assert(!val[0] || isalnum(val[0]) || val[0] == '.');

	// All 0's. We can just return, since this procedure expects a virgin
	// (already 0) BcNum.
	if (!val[0]) return;

	// The length of the string is the length of the number, except it might be
	// one bigger because of a decimal point.
	len = strlen(val);

	// Find the location of the decimal point.
	ptr = strchr(val, '.');
	rdx = (ptr != NULL);

	// We eat leading zeroes again. These leading zeroes are different because
	// they will come after the decimal point if they exist, and since that's
	// the case, they must be preserved.
	for (i = 0; i < len && (zero = (val[i] == '0' || val[i] == '.')); ++i);

	// Set the scale of the number based on the location of the decimal point.
	// The casts to uintptr_t is to ensure that bc does not hit undefined
	// behavior when doing math on the values.
	n->scale = (size_t) (rdx * (((uintptr_t) (val + len)) -
	                            (((uintptr_t) ptr) + 1)));

	// Set rdx.
	BC_NUM_RDX_SET(n, BC_NUM_RDX(n->scale));

	// Calculate length. First, the length of the integer, then the number of
	// digits in the last limb, then the length.
	i = len - (ptr == val ? 0 : i) - rdx;
	temp = BC_NUM_ROUND_POW(i);
	mod = n->scale % BC_BASE_DIGS;
	i = mod ? BC_BASE_DIGS - mod : 0;
	n->len = ((temp + i) / BC_BASE_DIGS);

	// Expand and zero.
	bc_num_expand(n, n->len);
	memset(n->num, 0, BC_NUM_SIZE(n->len));

	if (zero) {
		// I think I can set rdx directly to zero here because n should be a
		// new number with sign set to false.
		n->len = n->rdx = 0;
	}
	else {

		// There is actually stuff to parse if we make it here. Yay...
		BcBigDig exp, pow;

		assert(i <= BC_NUM_BIGDIG_MAX);

		// The exponent and power.
		exp = (BcBigDig) i;
		pow = bc_num_pow10[exp];

		// Parse loop. We parse backwards because numbers are stored little
		// endian.
		for (i = len - 1; i < len; --i, ++exp) {

			char c = val[i];

			// Skip the decimal point.
			if (c == '.') exp -= 1;
			else {

				// The index of the limb.
				size_t idx = exp / BC_BASE_DIGS;

				// Clamp for the base.
				if (isupper(c)) c = '9';

				// Add the digit to the limb.
				n->num[idx] += (((BcBigDig) c) - '0') * pow;

				// Adjust the power and exponent.
				if ((exp + 1) % BC_BASE_DIGS == 0) pow = 1;
				else pow *= BC_BASE;
			}
		}
	}
}

/**
 * Parse a number in any base (besides decimal).
 * @param n     The number to parse into and return. Must be preallocated.
 * @param val   The string to parse.
 * @param base  The base to parse as.
 */
static void bc_num_parseBase(BcNum *restrict n, const char *restrict val,
                             BcBigDig base)
{
	BcNum temp, mult1, mult2, result1, result2, *m1, *m2, *ptr;
	char c = 0;
	bool zero = true;
	BcBigDig v;
	size_t i, digs, len = strlen(val);

	// If zero, just return because the number should be virgin (already 0).
	for (i = 0; zero && i < len; ++i) zero = (val[i] == '.' || val[i] == '0');
	if (zero) return;

	BC_SIG_LOCK;

	bc_num_init(&temp, BC_NUM_BIGDIG_LOG10);
	bc_num_init(&mult1, BC_NUM_BIGDIG_LOG10);

	BC_SETJMP_LOCKED(int_err);

	BC_SIG_UNLOCK;

	// We split parsing into parsing the integer and parsing the fractional
	// part.

	// Parse the integer part. This is the easy part because we just multiply
	// the number by the base, then add the digit.
	for (i = 0; i < len && (c = val[i]) && c != '.'; ++i) {

		// Convert the character to a digit.
		v = bc_num_parseChar(c, base);

		// Multiply the number.
		bc_num_mulArray(n, base, &mult1);

		// Convert the digit to a number and add.
		bc_num_bigdig2num(&temp, v);
		bc_num_add(&mult1, &temp, n, 0);
	}

	// If this condition is true, then we are done. We still need to do cleanup
	// though.
	if (i == len && !val[i]) goto int_err;

	// If we get here, we *must* be at the radix point.
	assert(val[i] == '.');

	BC_SIG_LOCK;

	// Unset the jump to reset in for these new initializations.
	BC_UNSETJMP;

	bc_num_init(&mult2, BC_NUM_BIGDIG_LOG10);
	bc_num_init(&result1, BC_NUM_DEF_SIZE);
	bc_num_init(&result2, BC_NUM_DEF_SIZE);
	bc_num_one(&mult1);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	// Pointers for easy switching.
	m1 = &mult1;
	m2 = &mult2;

	// Parse the fractional part. This is the hard part.
	for (i += 1, digs = 0; i < len && (c = val[i]); ++i, ++digs) {

		size_t rdx;

		// Convert the character to a digit.
		v = bc_num_parseChar(c, base);

		// We keep growing result2 according to the base because the more digits
		// after the radix, the more significant the digits close to the radix
		// should be.
		bc_num_mulArray(&result1, base, &result2);

		// Convert the digit to a number.
		bc_num_bigdig2num(&temp, v);

		// Add the digit into the fraction part.
		bc_num_add(&result2, &temp, &result1, 0);

		// Keep growing m1 and m2 for use after the loop.
		bc_num_mulArray(m1, base, m2);

		rdx = BC_NUM_RDX_VAL(m2);

		if (m2->len < rdx) m2->len = rdx;

		// Switch.
		ptr = m1;
		m1 = m2;
		m2 = ptr;
	}

	// This one cannot be a divide by 0 because mult starts out at 1, then is
	// multiplied by base, and base cannot be 0, so mult cannot be 0. And this
	// is the reason we keep growing m1 and m2; this division is what converts
	// the parsed fractional part from an integer to a fractional part.
	bc_num_div(&result1, m1, &result2, digs * 2);

	// Pretruncate.
	bc_num_truncate(&result2, digs);

	// The final add of the integer part to the fractional part.
	bc_num_add(n, &result2, n, digs);

	// Basic cleanup.
	if (BC_NUM_NONZERO(n)) {
		if (n->scale < digs) bc_num_extend(n, digs - n->scale);
	}
	else bc_num_zero(n);

err:
	BC_SIG_MAYLOCK;
	bc_num_free(&result2);
	bc_num_free(&result1);
	bc_num_free(&mult2);
int_err:
	BC_SIG_MAYLOCK;
	bc_num_free(&mult1);
	bc_num_free(&temp);
	BC_LONGJMP_CONT;
}

/**
 * Prints a backslash+newline combo if the number of characters needs it. This
 * is really a convenience function.
 */
static inline void bc_num_printNewline(void) {
#if !BC_ENABLE_LIBRARY
	if (vm.nchars >= vm.line_len - 1) {
		bc_vm_putchar('\\', bc_flush_none);
		bc_vm_putchar('\n', bc_flush_err);
	}
#endif // !BC_ENABLE_LIBRARY
}

/**
 * Prints a character after a backslash+newline, if needed.
 * @param c       The character to print.
 * @param bslash  Whether to print a backslash+newline.
 */
static void bc_num_putchar(int c, bool bslash) {
	if (c != '\n' && bslash) bc_num_printNewline();
	bc_vm_putchar(c, bc_flush_save);
}

#if !BC_ENABLE_LIBRARY

/**
 * Prints a character for a number's digit. This is for printing for dc's P
 * command. This function does not need to worry about radix points. This is a
 * BcNumDigitOp.
 * @param n       The "digit" to print.
 * @param len     The "length" of the digit, or number of characters that will
 *                need to be printed for the digit.
 * @param rdx     True if a decimal (radix) point should be printed.
 * @param bslash  True if a backslash+newline should be printed if the character
 *                limit for the line is reached, false otherwise.
 */
static void bc_num_printChar(size_t n, size_t len, bool rdx, bool bslash) {
	BC_UNUSED(rdx);
	BC_UNUSED(len);
	BC_UNUSED(bslash);
	assert(len == 1);
	bc_vm_putchar((uchar) n, bc_flush_save);
}

#endif // !BC_ENABLE_LIBRARY

/**
 * Prints a series of characters for large bases. This is for printing in bases
 * above hexadecimal. This is a BcNumDigitOp.
 * @param n       The "digit" to print.
 * @param len     The "length" of the digit, or number of characters that will
 *                need to be printed for the digit.
 * @param rdx     True if a decimal (radix) point should be printed.
 * @param bslash  True if a backslash+newline should be printed if the character
 *                limit for the line is reached, false otherwise.
 */
static void bc_num_printDigits(size_t n, size_t len, bool rdx, bool bslash) {

	size_t exp, pow;

	// If needed, print the radix; otherwise, print a space to separate digits.
	bc_num_putchar(rdx ? '.' : ' ', true);

	// Calculate the exponent and power.
	for (exp = 0, pow = 1; exp < len - 1; ++exp, pow *= BC_BASE);

	// Print each character individually.
	for (exp = 0; exp < len; pow /= BC_BASE, ++exp) {

		// The individual subdigit.
		size_t dig = n / pow;

		// Take the subdigit away.
		n -= dig * pow;

		// Print the subdigit.
		bc_num_putchar(((uchar) dig) + '0', bslash || exp != len - 1);
	}
}

/**
 * Prints a character for a number's digit. This is for printing in bases for
 * hexadecimal and below because they always print only one character at a time.
 * This is a BcNumDigitOp.
 * @param n       The "digit" to print.
 * @param len     The "length" of the digit, or number of characters that will
 *                need to be printed for the digit.
 * @param rdx     True if a decimal (radix) point should be printed.
 * @param bslash  True if a backslash+newline should be printed if the character
 *                limit for the line is reached, false otherwise.
 */
static void bc_num_printHex(size_t n, size_t len, bool rdx, bool bslash) {

	BC_UNUSED(len);
	BC_UNUSED(bslash);

	assert(len == 1);

	if (rdx) bc_num_putchar('.', true);

	bc_num_putchar(bc_num_hex_digits[n], bslash);
}

/**
 * Prints a decimal number. This is specially written for optimization since
 * this will be used the most and because bc's numbers are already in decimal.
 * @param n        The number to print.
 * @param newline  Whether to print backslash+newlines on long enough lines.
 */
static void bc_num_printDecimal(const BcNum *restrict n, bool newline) {

	size_t i, j, rdx = BC_NUM_RDX_VAL(n);
	bool zero = true;
	size_t buffer[BC_BASE_DIGS];

	// Print the sign.
	if (BC_NUM_NEG(n)) bc_num_putchar('-', true);

	// Print loop.
	for (i = n->len - 1; i < n->len; --i) {

		BcDig n9 = n->num[i];
		size_t temp;
		bool irdx = (i == rdx - 1);

		// Calculate the number of digits in the limb.
		zero = (zero & !irdx);
		temp = n->scale % BC_BASE_DIGS;
		temp = i || !temp ? 0 : BC_BASE_DIGS - temp;

		memset(buffer, 0, BC_BASE_DIGS * sizeof(size_t));

		// Fill the buffer with individual digits.
		for (j = 0; n9 && j < BC_BASE_DIGS; ++j) {
			buffer[j] = ((size_t) n9) % BC_BASE;
			n9 /= BC_BASE;
		}

		// Print the digits in the buffer.
		for (j = BC_BASE_DIGS - 1; j < BC_BASE_DIGS && j >= temp; --j) {

			// Figure out whether to print the decimal point.
			bool print_rdx = (irdx & (j == BC_BASE_DIGS - 1));

			// The zero variable helps us skip leading zero digits in the limb.
			zero = (zero && buffer[j] == 0);

			if (!zero) {

				// While the first three arguments should be self-explanatory,
				// the last needs explaining. I don't want to print a newline
				// when the last digit to be printed could take the place of the
				// backslash rather than being pushed, as a single character, to
				// the next line. That's what that last argument does for bc.
				bc_num_printHex(buffer[j], 1, print_rdx,
				                !newline || (j > temp || i != 0));
			}
		}
	}
}

#if BC_ENABLE_EXTRA_MATH

/**
 * Prints a number in scientific or engineering format. When doing this, we are
 * always printing in decimal.
 * @param n        The number to print.
 * @param eng      True if we are in engineering mode.
 * @param newline  Whether to print backslash+newlines on long enough lines.
 */
static void bc_num_printExponent(const BcNum *restrict n,
                                 bool eng, bool newline)
{
	size_t places, mod, nrdx = BC_NUM_RDX_VAL(n);
	bool neg = (n->len <= nrdx);
	BcNum temp, exp;
	BcDig digs[BC_NUM_BIGDIG_LOG10];

	BC_SIG_LOCK;

	bc_num_createCopy(&temp, n);

	BC_SETJMP_LOCKED(exit);

	BC_SIG_UNLOCK;

	// We need to calculate the exponents, and they change based on whether the
	// number is all fractional or not, obviously.
	if (neg) {

		// Figure out how many limbs after the decimal point is zero.
		size_t i, idx = bc_num_nonZeroLen(n) - 1;

		places = 1;

		// Figure out how much in the last limb is zero.
		for (i = BC_BASE_DIGS - 1; i < BC_BASE_DIGS; --i) {
			if (bc_num_pow10[i] > (BcBigDig) n->num[idx]) places += 1;
			else break;
		}

		// Calculate the combination of zero limbs and zero digits in the last
		// limb.
		places += (nrdx - (idx + 1)) * BC_BASE_DIGS;
		mod = places % 3;

		// Calculate places if we are in engineering mode.
		if (eng && mod != 0) places += 3 - mod;

		// Shift the temp to the right place.
		bc_num_shiftLeft(&temp, places);
	}
	else {

		// This is the number of digits that we are supposed to put behind the
		// decimal point.
		places = bc_num_intDigits(n) - 1;

		// Calculate the true number based on whether engineering mode is
		// activated.
		mod = places % 3;
		if (eng && mod != 0) places -= 3 - (3 - mod);

		// Shift the temp to the right place.
		bc_num_shiftRight(&temp, places);
	}

	// Print the shifted number.
	bc_num_printDecimal(&temp, newline);

	// Print the e.
	bc_num_putchar('e', !newline);

	// Need to explicitly print a zero exponent.
	if (!places) {
		bc_num_printHex(0, 1, false, !newline);
		goto exit;
	}

	// Need to print sign for the exponent.
	if (neg) bc_num_putchar('-', true);

	// Create a temporary for the exponent...
	bc_num_setup(&exp, digs, BC_NUM_BIGDIG_LOG10);
	bc_num_bigdig2num(&exp, (BcBigDig) places);

	/// ..and print it.
	bc_num_printDecimal(&exp, newline);

exit:
	BC_SIG_MAYLOCK;
	bc_num_free(&temp);
	BC_LONGJMP_CONT;
}
#endif // BC_ENABLE_EXTRA_MATH

/**
 * Converts a number from limbs with base BC_BASE_POW to base @a pow, where
 * @a pow is obase^N.
 * @param n    The number to convert.
 * @param rem  BC_BASE_POW - @a pow.
 * @param pow  The power of obase we will convert the number to.
 * @param idx  The index of the number to start converting at. Doing the
 *             conversion is O(n^2); we have to sweep through starting at the
 *             least significant limb
 */
static void bc_num_printFixup(BcNum *restrict n, BcBigDig rem,
                              BcBigDig pow, size_t idx)
{
	size_t i, len = n->len - idx;
	BcBigDig acc;
	BcDig *a = n->num + idx;

	// Ignore if there's just one limb left. This is the part that requires the
	// extra loop after the one calling this function in bc_num_printPrepare().
	if (len < 2) return;

	// Loop through the remaining limbs and convert. We start at the second limb
	// because we pull the value from the previous one as well.
	for (i = len - 1; i > 0; --i) {

		// Get the limb and add it to the previous, along with multiplying by
		// the remainder because that's the proper overflow. "acc" means
		// "accumulator," by the way.
		acc = ((BcBigDig) a[i]) * rem + ((BcBigDig) a[i - 1]);

		// Store a value in base pow in the previous limb.
		a[i - 1] = (BcDig) (acc % pow);

		// Divide by the base and accumulate the remaining value in the limb.
		acc /= pow;
		acc += (BcBigDig) a[i];

		// If the accumulator is greater than the base...
		if (acc >= BC_BASE_POW) {

			// Do we need to grow?
			if (i == len - 1) {

				// Grow.
				len = bc_vm_growSize(len, 1);
				bc_num_expand(n, bc_vm_growSize(len, idx));

				// Update the pointer because it may have moved.
				a = n->num + idx;

				// Zero out the last limb.
				a[len - 1] = 0;
			}

			// Overflow into the next limb since we are over the base.
			a[i + 1] += acc / BC_BASE_POW;
			acc %= BC_BASE_POW;
		}

		assert(acc < BC_BASE_POW);

		// Set the limb.
		a[i] = (BcDig) acc;
	}

	// We may have grown the number, so adjust the length.
	n->len = len + idx;
}

/**
 * Prepares a number for printing in a base that is not a divisor of
 * BC_BASE_POW. This basically converts the number from having limbs of base
 * BC_BASE_POW to limbs of pow, where pow is obase^N.
 * @param n    The number to prepare for printing.
 * @param rem  The remainder of BC_BASE_POW when divided by a power of the base.
 * @param pow  The power of the base.
 */
static void bc_num_printPrepare(BcNum *restrict n, BcBigDig rem, BcBigDig pow) {

	size_t i;

	// Loop from the least significant limb to the most significant limb and
	// convert limbs in each pass.
	for (i = 0; i < n->len; ++i) bc_num_printFixup(n, rem, pow, i);

	// bc_num_printFixup() does not do everything it is supposed to, so we do
	// the last bit of cleanup here. That cleanup is to ensure that each limb
	// is less than pow and to expand the number to fit new limbs as necessary.
	for (i = 0; i < n->len; ++i) {

		assert(pow == ((BcBigDig) ((BcDig) pow)));

		// If the limb needs fixing...
		if (n->num[i] >= (BcDig) pow) {

			// Do we need to grow?
			if (i + 1 == n->len) {

				// Grow the number.
				n->len = bc_vm_growSize(n->len, 1);
				bc_num_expand(n, n->len);

				// Without this, we might use uninitialized data.
				n->num[i + 1] = 0;
			}

			assert(pow < BC_BASE_POW);

			// Overflow into the next limb.
			n->num[i + 1] += n->num[i] / ((BcDig) pow);
			n->num[i] %= (BcDig) pow;
		}
	}
}

static void bc_num_printNum(BcNum *restrict n, BcBigDig base, size_t len,
                            BcNumDigitOp print, bool newline)
{
	BcVec stack;
	BcNum intp, fracp1, fracp2, digit, flen1, flen2, *n1, *n2, *temp;
	BcBigDig dig = 0, *ptr, acc, exp;
	size_t i, j, nrdx, idigits;
	bool radix;
	BcDig digit_digs[BC_NUM_BIGDIG_LOG10 + 1];

	assert(base > 1);

	// Easy case. Even with scale, we just print this.
	if (BC_NUM_ZERO(n)) {
		print(0, len, false, !newline);
		return;
	}

	// This function uses an algorithm that Stefan Esser <se@freebsd.org> came
	// up with to print the integer part of a number. What it does is convert
	// intp into a number of the specified base, but it does it directly,
	// instead of just doing a series of divisions and printing the remainders
	// in reverse order.
	//
	// Let me explain in a bit more detail:
	//
	// The algorithm takes the current least significant limb (after intp has
	// been converted to an integer) and the next to least significant limb, and
	// it converts the least significant limb into one of the specified base,
	// putting any overflow into the next to least significant limb. It iterates
	// through the whole number, from least significant to most significant,
	// doing this conversion. At the end of that iteration, the least
	// significant limb is converted, but the others are not, so it iterates
	// again, starting at the next to least significant limb. It keeps doing
	// that conversion, skipping one more limb than the last time, until all
	// limbs have been converted. Then it prints them in reverse order.
	//
	// That is the gist of the algorithm. It leaves out several things, such as
	// the fact that limbs are not always converted into the specified base, but
	// into something close, basically a power of the specified base. In
	// Stefan's words, "You could consider BcDigs to be of base 10^BC_BASE_DIGS
	// in the normal case and obase^N for the largest value of N that satisfies
	// obase^N <= 10^BC_BASE_DIGS. [This means that] the result is not in base
	// "obase", but in base "obase^N", which happens to be printable as a number
	// of base "obase" without consideration for neighbouring BcDigs." This fact
	// is what necessitates the existence of the loop later in this function.
	//
	// The conversion happens in bc_num_printPrepare() where the outer loop
	// happens and bc_num_printFixup() where the inner loop, or actual
	// conversion, happens. In other words, bc_num_printPrepare() is where the
	// loop that starts at the least significant limb and goes to the most
	// significant limb. Then, on every iteration of its loop, it calls
	// bc_num_printFixup(), which has the inner loop of actually converting
	// the limbs it passes into limbs of base obase^N rather than base
	// BC_BASE_POW.

	nrdx = BC_NUM_RDX_VAL(n);

	BC_SIG_LOCK;

	// The stack is what allows us to reverse the digits for printing.
	bc_vec_init(&stack, sizeof(BcBigDig), BC_DTOR_NONE);
	bc_num_init(&fracp1, nrdx);

	// intp will be the "integer part" of the number, so copy it.
	bc_num_createCopy(&intp, n);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	// Make intp an integer.
	bc_num_truncate(&intp, intp.scale);

	// Get the fractional part out.
	bc_num_sub(n, &intp, &fracp1, 0);

	// If the base is not the same as the last base used for printing, we need
	// to update the cached exponent and power. Yes, we cache the values of the
	// exponent and power. That is to prevent us from calculating them every
	// time because printing will probably happen multiple times on the same
	// base.
	if (base != vm.last_base) {

		vm.last_pow = 1;
		vm.last_exp = 0;

		// Calculate the exponent and power.
		while (vm.last_pow * base <= BC_BASE_POW) {
			vm.last_pow *= base;
			vm.last_exp += 1;
		}

		// Also, the remainder and base itself.
		vm.last_rem = BC_BASE_POW - vm.last_pow;
		vm.last_base = base;
	}

	exp = vm.last_exp;

	// If vm.last_rem is 0, then the base we are printing in is a divisor of
	// BC_BASE_POW, which is the easy case because it means that BC_BASE_POW is
	// a power of obase, and no conversion is needed. If it *is* 0, then we have
	// the hard case, and we have to prepare the number for the base.
	if (vm.last_rem != 0) bc_num_printPrepare(&intp, vm.last_rem, vm.last_pow);

	// After the conversion comes the surprisingly easy part. From here on out,
	// this is basically naive code that I wrote, adjusted for the larger bases.

	// Fill the stack of digits for the integer part.
	for (i = 0; i < intp.len; ++i) {

		// Get the limb.
		acc = (BcBigDig) intp.num[i];

		// Turn the limb into digits of base obase.
		for (j = 0; j < exp && (i < intp.len - 1 || acc != 0); ++j)
		{
			// This condition is true if we are not at the last digit.
			if (j != exp - 1) {
				dig = acc % base;
				acc /= base;
			}
			else {
				dig = acc;
				acc = 0;
			}

			assert(dig < base);

			// Push the digit onto the stack.
			bc_vec_push(&stack, &dig);
		}

		assert(acc == 0);
	}

	// Go through the stack backwards and print each digit.
	for (i = 0; i < stack.len; ++i) {

		ptr = bc_vec_item_rev(&stack, i);

		assert(ptr != NULL);

		// While the first three arguments should be self-explanatory, the last
		// needs explaining. I don't want to print a newline when the last digit
		// to be printed could take the place of the backslash rather than being
		// pushed, as a single character, to the next line. That's what that
		// last argument does for bc.
		print(*ptr, len, false, !newline ||
		      (n->scale != 0 || i == stack.len - 1));
	}

	// We are done if there is no fractional part.
	if (!n->scale) goto err;

	BC_SIG_LOCK;

	// Reset the jump because some locals are changing.
	BC_UNSETJMP;

	bc_num_init(&fracp2, nrdx);
	bc_num_setup(&digit, digit_digs, sizeof(digit_digs) / sizeof(BcDig));
	bc_num_init(&flen1, BC_NUM_BIGDIG_LOG10);
	bc_num_init(&flen2, BC_NUM_BIGDIG_LOG10);

	BC_SETJMP_LOCKED(frac_err);

	BC_SIG_UNLOCK;

	bc_num_one(&flen1);

	radix = true;

	// Pointers for easy switching.
	n1 = &flen1;
	n2 = &flen2;

	fracp2.scale = n->scale;
	BC_NUM_RDX_SET_NP(fracp2, BC_NUM_RDX(fracp2.scale));

	// As long as we have not reached the scale of the number, keep printing.
	while ((idigits = bc_num_intDigits(n1)) <= n->scale) {

		// These numbers will keep growing.
		bc_num_expand(&fracp2, fracp1.len + 1);
		bc_num_mulArray(&fracp1, base, &fracp2);

		nrdx = BC_NUM_RDX_VAL_NP(fracp2);

		// Ensure an invariant.
		if (fracp2.len < nrdx) fracp2.len = nrdx;

		// fracp is guaranteed to be non-negative and small enough.
		dig = bc_num_bigdig2(&fracp2);

		// Convert the digit to a number and subtract it from the number.
		bc_num_bigdig2num(&digit, dig);
		bc_num_sub(&fracp2, &digit, &fracp1, 0);

		// While the first three arguments should be self-explanatory, the last
		// needs explaining. I don't want to print a newline when the last digit
		// to be printed could take the place of the backslash rather than being
		// pushed, as a single character, to the next line. That's what that
		// last argument does for bc.
		print(dig, len, radix, !newline || idigits != n->scale);

		// Update the multipliers.
		bc_num_mulArray(n1, base, n2);

		radix = false;

		// Switch.
		temp = n1;
		n1 = n2;
		n2 = temp;
	}

frac_err:
	BC_SIG_MAYLOCK;
	bc_num_free(&flen2);
	bc_num_free(&flen1);
	bc_num_free(&fracp2);
err:
	BC_SIG_MAYLOCK;
	bc_num_free(&fracp1);
	bc_num_free(&intp);
	bc_vec_free(&stack);
	BC_LONGJMP_CONT;
}

/**
 * Prints a number in the specified base, or rather, figures out which function
 * to call to print the number in the specified base and calls it.
 * @param n        The number to print.
 * @param base     The base to print in.
 * @param newline  Whether to print backslash+newlines on long enough lines.
 */
static void bc_num_printBase(BcNum *restrict n, BcBigDig base, bool newline) {

	size_t width;
	BcNumDigitOp print;
	bool neg = BC_NUM_NEG(n);

	// Just take care of the sign right here.
	if (neg) bc_num_putchar('-', true);

	// Clear the sign because it makes the actual printing easier when we have
	// to do math.
	BC_NUM_NEG_CLR(n);

	// Bases at hexadecimal and below are printed as one character, larger bases
	// are printed as a series of digits separated by spaces.
	if (base <= BC_NUM_MAX_POSIX_IBASE) {
		width = 1;
		print = bc_num_printHex;
	}
	else {
		assert(base <= BC_BASE_POW);
		width = bc_num_log10(base - 1);
		print = bc_num_printDigits;
	}

	// Print.
	bc_num_printNum(n, base, width, print, newline);

	// Reset the sign.
	n->rdx = BC_NUM_NEG_VAL(n, neg);
}

#if !BC_ENABLE_LIBRARY

void bc_num_stream(BcNum *restrict n) {
	bc_num_printNum(n, BC_NUM_STREAM_BASE, 1, bc_num_printChar, false);
}

#endif // !BC_ENABLE_LIBRARY

void bc_num_setup(BcNum *restrict n, BcDig *restrict num, size_t cap) {
	assert(n != NULL);
	n->num = num;
	n->cap = cap;
	bc_num_zero(n);
}

void bc_num_init(BcNum *restrict n, size_t req) {

	BcDig *num;

	BC_SIG_ASSERT_LOCKED;

	assert(n != NULL);

	// BC_NUM_DEF_SIZE is set to be about the smallest allocation size that
	// malloc() returns in practice, so just use it.
	req = req >= BC_NUM_DEF_SIZE ? req : BC_NUM_DEF_SIZE;

	// If we can't use a temp, allocate.
	if (req != BC_NUM_DEF_SIZE || (num = bc_vm_takeTemp()) == NULL)
		num = bc_vm_malloc(BC_NUM_SIZE(req));

	bc_num_setup(n, num, req);
}

void bc_num_clear(BcNum *restrict n) {
	n->num = NULL;
	n->cap = 0;
}

void bc_num_free(void *num) {

	BcNum *n = (BcNum*) num;

	BC_SIG_ASSERT_LOCKED;

	assert(n != NULL);

	if (n->cap == BC_NUM_DEF_SIZE) bc_vm_addTemp(n->num);
	else free(n->num);
}

void bc_num_copy(BcNum *d, const BcNum *s) {

	assert(d != NULL && s != NULL);

	if (d == s) return;

	bc_num_expand(d, s->len);
	d->len = s->len;

	// I can just copy directly here because the sign *and* rdx will be
	// properly preserved.
	d->rdx = s->rdx;
	d->scale = s->scale;
	memcpy(d->num, s->num, BC_NUM_SIZE(d->len));
}

void bc_num_createCopy(BcNum *d, const BcNum *s) {
	BC_SIG_ASSERT_LOCKED;
	bc_num_init(d, s->len);
	bc_num_copy(d, s);
}

void bc_num_createFromBigdig(BcNum *restrict n, BcBigDig val) {
	BC_SIG_ASSERT_LOCKED;
	bc_num_init(n, BC_NUM_BIGDIG_LOG10);
	bc_num_bigdig2num(n, val);
}

size_t bc_num_scale(const BcNum *restrict n) {
	return n->scale;
}

size_t bc_num_len(const BcNum *restrict n) {

	size_t len = n->len;

	// Always return at least 1.
	if (BC_NUM_ZERO(n)) return n->scale ? n->scale : 1;

	// If this is true, there is no integer portion of the number.
	if (BC_NUM_RDX_VAL(n) == len) {

		// We have to take into account the fact that some of the digits right
		// after the decimal could be zero. If that is the case, we need to
		// ignore them until we hit the first non-zero digit.

		size_t zero, scale;

		// The number of limbs with non-zero digits.
		len = bc_num_nonZeroLen(n);

		// Get the number of digits in the last limb.
		scale = n->scale % BC_BASE_DIGS;
		scale = scale ? scale : BC_BASE_DIGS;

		// Get the number of zero digits.
		zero = bc_num_zeroDigits(n->num + len - 1);

		// Calculate the true length.
		len = len * BC_BASE_DIGS - zero - (BC_BASE_DIGS - scale);
	}
	// Otherwise, count the number of int digits and return that plus the scale.
	else len = bc_num_intDigits(n) + n->scale;

	return len;
}

void bc_num_parse(BcNum *restrict n, const char *restrict val, BcBigDig base) {

	assert(n != NULL && val != NULL && base);
	assert(base >= BC_NUM_MIN_BASE && base <= vm.maxes[BC_PROG_GLOBALS_IBASE]);
	assert(bc_num_strValid(val));

	// A one character number is *always* parsed as though the base was the
	// maximum allowed ibase, per the bc spec.
	if (!val[1]) {
		BcBigDig dig = bc_num_parseChar(val[0], BC_NUM_MAX_LBASE);
		bc_num_bigdig2num(n, dig);
	}
	else if (base == BC_BASE) bc_num_parseDecimal(n, val);
	else bc_num_parseBase(n, val, base);

	assert(BC_NUM_RDX_VALID(n));
}

void bc_num_print(BcNum *restrict n, BcBigDig base, bool newline) {

	assert(n != NULL);
	assert(BC_ENABLE_EXTRA_MATH || base >= BC_NUM_MIN_BASE);

	// We may need a newline, just to start.
	bc_num_printNewline();

	// Short-circuit 0.
	if (BC_NUM_ZERO(n)) bc_num_printHex(0, 1, false, !newline);
	else if (base == BC_BASE) bc_num_printDecimal(n, newline);
#if BC_ENABLE_EXTRA_MATH
	else if (base == 0 || base == 1)
		bc_num_printExponent(n, base != 0, newline);
#endif // BC_ENABLE_EXTRA_MATH
	else bc_num_printBase(n, base, newline);

	if (newline) bc_num_putchar('\n', false);
}

BcBigDig bc_num_bigdig2(const BcNum *restrict n) {

	// This function returns no errors because it's guaranteed to succeed if
	// its preconditions are met. Those preconditions include both n needs to
	// be non-NULL, n being non-negative, and n being less than vm.max. If all
	// of that is true, then we can just convert without worrying about negative
	// errors or overflow.

	BcBigDig r = 0;
	size_t nrdx = BC_NUM_RDX_VAL(n);

	assert(n != NULL);
	assert(!BC_NUM_NEG(n));
	assert(bc_num_cmp(n, &vm.max) < 0);
	assert(n->len - nrdx <= 3);

	// There is a small speed win from unrolling the loop here, and since it
	// only adds 53 bytes, I decided that it was worth it.
	switch (n->len - nrdx) {

		case 3:
		{
			r = (BcBigDig) n->num[nrdx + 2];
		}
		// Fallthrough.
		BC_FALLTHROUGH

		case 2:
		{
			r = r * BC_BASE_POW + (BcBigDig) n->num[nrdx + 1];
		}
		// Fallthrough.
		BC_FALLTHROUGH

		case 1:
		{
			r = r * BC_BASE_POW + (BcBigDig) n->num[nrdx];
		}
	}

	return r;
}

BcBigDig bc_num_bigdig(const BcNum *restrict n) {

	assert(n != NULL);

	// This error checking is extremely important, and if you do not have a
	// guarantee that converting a number will always succeed in a particular
	// case, you *must* call this function to get these error checks. This
	// includes all instances of numbers inputted by the user or calculated by
	// the user. Otherwise, you can call the faster bc_num_bigdig2().
	if (BC_ERR(BC_NUM_NEG(n))) bc_err(BC_ERR_MATH_NEGATIVE);
	if (BC_ERR(bc_num_cmp(n, &vm.max) >= 0)) bc_err(BC_ERR_MATH_OVERFLOW);

	return bc_num_bigdig2(n);
}

void bc_num_bigdig2num(BcNum *restrict n, BcBigDig val) {

	BcDig *ptr;
	size_t i;

	assert(n != NULL);

	bc_num_zero(n);

	// Already 0.
	if (!val) return;

	// Expand first. This is the only way this function can fail, and it's a
	// fatal error.
	bc_num_expand(n, BC_NUM_BIGDIG_LOG10);

	// The conversion is easy because numbers are laid out in little-endian
	// order.
	for (ptr = n->num, i = 0; val; ++i, val /= BC_BASE_POW)
		ptr[i] = val % BC_BASE_POW;

	n->len = i;
}

#if BC_ENABLE_EXTRA_MATH

void bc_num_rng(const BcNum *restrict n, BcRNG *rng) {

	BcNum temp, temp2, intn, frac;
	BcRand state1, state2, inc1, inc2;
	size_t nrdx = BC_NUM_RDX_VAL(n);

	// This function holds the secret of how I interpret a seed number for the
	// PRNG. Well, it's actually in the development manual
	// (manuals/development.md#pseudo-random-number-generator), so look there
	// before you try to understand this.

	BC_SIG_LOCK;

	bc_num_init(&temp, n->len);
	bc_num_init(&temp2, n->len);
	bc_num_init(&frac, nrdx);
	bc_num_init(&intn, bc_num_int(n));

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	assert(BC_NUM_RDX_VALID_NP(vm.max));

	memcpy(frac.num, n->num, BC_NUM_SIZE(nrdx));
	frac.len = nrdx;
	BC_NUM_RDX_SET_NP(frac, nrdx);
	frac.scale = n->scale;

	assert(BC_NUM_RDX_VALID_NP(frac));
	assert(BC_NUM_RDX_VALID_NP(vm.max2));

	// Multiply the fraction and truncate so that it's an integer. The
	// truncation is what clamps it, by the way.
	bc_num_mul(&frac, &vm.max2, &temp, 0);
	bc_num_truncate(&temp, temp.scale);
	bc_num_copy(&frac, &temp);

	// Get the integer.
	memcpy(intn.num, n->num + nrdx, BC_NUM_SIZE(bc_num_int(n)));
	intn.len = bc_num_int(n);

	// This assert is here because it has to be true. It is also here to justify
	// some optimizations.
	assert(BC_NUM_NONZERO(&vm.max));

	// If there *was* a fractional part...
	if (BC_NUM_NONZERO(&frac)) {

		// This divmod splits frac into the two state parts.
		bc_num_divmod(&frac, &vm.max, &temp, &temp2, 0);

		// frac is guaranteed to be smaller than vm.max * vm.max (pow).
		// This means that when dividing frac by vm.max, as above, the
		// quotient and remainder are both guaranteed to be less than vm.max,
		// which means we can use bc_num_bigdig2() here and not worry about
		// overflow.
		state1 = (BcRand) bc_num_bigdig2(&temp2);
		state2 = (BcRand) bc_num_bigdig2(&temp);
	}
	else state1 = state2 = 0;

	// If there *was* an integer part...
	if (BC_NUM_NONZERO(&intn)) {

		// This divmod splits intn into the two inc parts.
		bc_num_divmod(&intn, &vm.max, &temp, &temp2, 0);

		// Because temp2 is the mod of vm.max, from above, it is guaranteed
		// to be small enough to use bc_num_bigdig2().
		inc1 = (BcRand) bc_num_bigdig2(&temp2);

		// Clamp the second inc part.
		if (bc_num_cmp(&temp, &vm.max) >= 0) {
			bc_num_copy(&temp2, &temp);
			bc_num_mod(&temp2, &vm.max, &temp, 0);
		}

		// The if statement above ensures that temp is less than vm.max, which
		// means that we can use bc_num_bigdig2() here.
		inc2 = (BcRand) bc_num_bigdig2(&temp);
	}
	else inc1 = inc2 = 0;

	bc_rand_seed(rng, state1, state2, inc1, inc2);

err:
	BC_SIG_MAYLOCK;
	bc_num_free(&intn);
	bc_num_free(&frac);
	bc_num_free(&temp2);
	bc_num_free(&temp);
	BC_LONGJMP_CONT;
}

void bc_num_createFromRNG(BcNum *restrict n, BcRNG *rng) {

	BcRand s1, s2, i1, i2;
	BcNum conv, temp1, temp2, temp3;
	BcDig temp1_num[BC_RAND_NUM_SIZE], temp2_num[BC_RAND_NUM_SIZE];
	BcDig conv_num[BC_NUM_BIGDIG_LOG10];

	BC_SIG_LOCK;

	bc_num_init(&temp3, 2 * BC_RAND_NUM_SIZE);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	bc_num_setup(&temp1, temp1_num, sizeof(temp1_num) / sizeof(BcDig));
	bc_num_setup(&temp2, temp2_num, sizeof(temp2_num) / sizeof(BcDig));
	bc_num_setup(&conv, conv_num, sizeof(conv_num) / sizeof(BcDig));

	// This assert is here because it has to be true. It is also here to justify
	// the assumption that vm.max is not zero.
	assert(BC_NUM_NONZERO(&vm.max));

	// Because this is true, we can just ignore math errors that would happen
	// otherwise.
	assert(BC_NUM_NONZERO(&vm.max2));

	bc_rand_getRands(rng, &s1, &s2, &i1, &i2);

	// Put the second piece of state into a number.
	bc_num_bigdig2num(&conv, (BcBigDig) s2);

	assert(BC_NUM_RDX_VALID_NP(conv));

	// Multiply by max to make room for the first piece of state.
	bc_num_mul(&conv, &vm.max, &temp1, 0);

	// Add in the first piece of state.
	bc_num_bigdig2num(&conv, (BcBigDig) s1);
	bc_num_add(&conv, &temp1, &temp2, 0);

	// Divide to make it an entirely fractional part.
	bc_num_div(&temp2, &vm.max2, &temp3, BC_RAND_STATE_BITS);

	// Now start on the increment parts. It's the same process without the
	// divide, so put the second piece of increment into a number.
	bc_num_bigdig2num(&conv, (BcBigDig) i2);

	assert(BC_NUM_RDX_VALID_NP(conv));

	// Multiply by max to make room for the first piece of increment.
	bc_num_mul(&conv, &vm.max, &temp1, 0);

	// Add in the first piece of increment.
	bc_num_bigdig2num(&conv, (BcBigDig) i1);
	bc_num_add(&conv, &temp1, &temp2, 0);

	// Now add the two together.
	bc_num_add(&temp2, &temp3, n, 0);

	assert(BC_NUM_RDX_VALID(n));

err:
	BC_SIG_MAYLOCK;
	bc_num_free(&temp3);
	BC_LONGJMP_CONT;
}

void bc_num_irand(BcNum *restrict a, BcNum *restrict b, BcRNG *restrict rng) {

	BcNum atemp;
	size_t i, len;

	assert(a != b);

	if (BC_ERR(BC_NUM_NEG(a))) bc_err(BC_ERR_MATH_NEGATIVE);

	// If either of these are true, then the numbers are integers.
	if (BC_NUM_ZERO(a) || BC_NUM_ONE(a)) return;

	if (BC_ERR(bc_num_nonInt(a, &atemp))) bc_err(BC_ERR_MATH_NON_INTEGER);

	assert(atemp.len);

	len = atemp.len - 1;

	// Just generate a random number for each limb.
	for (i = 0; i < len; ++i)
		b->num[i] = (BcDig) bc_rand_bounded(rng, BC_BASE_POW);

	// Do the last digit explicitly because the bound must be right. But only
	// do it if the limb does not equal 1. If it does, we have already hit the
	// limit.
	if (atemp.num[i] != 1) {
		b->num[i] = (BcDig) bc_rand_bounded(rng, (BcRand) atemp.num[i]);
		b->len = atemp.len;
	}
	// We want 1 less len in the case where we skip the last limb.
	else b->len = len;

	bc_num_clean(b);

	assert(BC_NUM_RDX_VALID(b));
}
#endif // BC_ENABLE_EXTRA_MATH

size_t bc_num_addReq(const BcNum *a, const BcNum *b, size_t scale) {

	size_t aint, bint, ardx, brdx;

	// Addition and subtraction require the max of the length of the two numbers
	// plus 1.

	BC_UNUSED(scale);

	ardx = BC_NUM_RDX_VAL(a);
	aint = bc_num_int(a);
	assert(aint <= a->len && ardx <= a->len);

	brdx = BC_NUM_RDX_VAL(b);
	bint = bc_num_int(b);
	assert(bint <= b->len && brdx <= b->len);

	ardx = BC_MAX(ardx, brdx);
	aint = BC_MAX(aint, bint);

	return bc_vm_growSize(bc_vm_growSize(ardx, aint), 1);
}

size_t bc_num_mulReq(const BcNum *a, const BcNum *b, size_t scale) {

	size_t max, rdx;

	// Multiplication requires the sum of the lengths of the numbers.

	rdx = bc_vm_growSize(BC_NUM_RDX_VAL(a), BC_NUM_RDX_VAL(b));

	max = BC_NUM_RDX(scale);

	max = bc_vm_growSize(BC_MAX(max, rdx), 1);
	rdx = bc_vm_growSize(bc_vm_growSize(bc_num_int(a), bc_num_int(b)), max);

	return rdx;
}

size_t bc_num_divReq(const BcNum *a, const BcNum *b, size_t scale) {

	size_t max, rdx;

	// Division requires the length of the dividend plus the scale.

	rdx = bc_vm_growSize(BC_NUM_RDX_VAL(a), BC_NUM_RDX_VAL(b));

	max = BC_NUM_RDX(scale);

	max = bc_vm_growSize(BC_MAX(max, rdx), 1);
	rdx = bc_vm_growSize(bc_num_int(a), max);

	return rdx;
}

size_t bc_num_powReq(const BcNum *a, const BcNum *b, size_t scale) {
	BC_UNUSED(scale);
	return bc_vm_growSize(bc_vm_growSize(a->len, b->len), 1);
}

#if BC_ENABLE_EXTRA_MATH
size_t bc_num_placesReq(const BcNum *a, const BcNum *b, size_t scale) {
	BC_UNUSED(scale);
	return a->len + b->len - BC_NUM_RDX_VAL(a) - BC_NUM_RDX_VAL(b);
}
#endif // BC_ENABLE_EXTRA_MATH

void bc_num_add(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));
	bc_num_binary(a, b, c, false, bc_num_as, bc_num_addReq(a, b, scale));
}

void bc_num_sub(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));
	bc_num_binary(a, b, c, true, bc_num_as, bc_num_addReq(a, b, scale));
}

void bc_num_mul(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));
	bc_num_binary(a, b, c, scale, bc_num_m, bc_num_mulReq(a, b, scale));
}

void bc_num_div(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));
	bc_num_binary(a, b, c, scale, bc_num_d, bc_num_divReq(a, b, scale));
}

void bc_num_mod(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));
	bc_num_binary(a, b, c, scale, bc_num_rem, bc_num_divReq(a, b, scale));
}

void bc_num_pow(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));
	bc_num_binary(a, b, c, scale, bc_num_p, bc_num_powReq(a, b, scale));
}

#if BC_ENABLE_EXTRA_MATH
void bc_num_places(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));
	bc_num_binary(a, b, c, scale, bc_num_place, bc_num_placesReq(a, b, scale));
}

void bc_num_lshift(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));
	bc_num_binary(a, b, c, scale, bc_num_left, bc_num_placesReq(a, b, scale));
}

void bc_num_rshift(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
	assert(BC_NUM_RDX_VALID(a));
	assert(BC_NUM_RDX_VALID(b));
	bc_num_binary(a, b, c, scale, bc_num_right, bc_num_placesReq(a, b, scale));
}
#endif // BC_ENABLE_EXTRA_MATH

void bc_num_sqrt(BcNum *restrict a, BcNum *restrict b, size_t scale) {

	BcNum num1, num2, half, f, fprime, *x0, *x1, *temp;
	size_t pow, len, rdx, req, resscale;
	BcDig half_digs[1];

	assert(a != NULL && b != NULL && a != b);

	if (BC_ERR(BC_NUM_NEG(a))) bc_err(BC_ERR_MATH_NEGATIVE);

	// We want to calculate to a's scale if it is bigger so that the result will
	// truncate properly.
	if (a->scale > scale) scale = a->scale;

	// Set parameters for the result.
	len = bc_vm_growSize(bc_num_intDigits(a), 1);
	rdx = BC_NUM_RDX(scale);

	// Square root needs half of the length of the parameter.
	req = bc_vm_growSize(BC_MAX(rdx, BC_NUM_RDX_VAL(a)), len >> 1);

	BC_SIG_LOCK;

	// Unlike the binary operators, this function is the only single parameter
	// function and is expected to initialize the result. This means that it
	// expects that b is *NOT* preallocated. We allocate it here.
	bc_num_init(b, bc_vm_growSize(req, 1));

	BC_SIG_UNLOCK;

	assert(a != NULL && b != NULL && a != b);
	assert(a->num != NULL && b->num != NULL);

	// Easy case.
	if (BC_NUM_ZERO(a)) {
		bc_num_setToZero(b, scale);
		return;
	}

	// Another easy case.
	if (BC_NUM_ONE(a)) {
		bc_num_one(b);
		bc_num_extend(b, scale);
		return;
	}

	// Set the parameters again.
	rdx = BC_NUM_RDX(scale);
	rdx = BC_MAX(rdx, BC_NUM_RDX_VAL(a));
	len = bc_vm_growSize(a->len, rdx);

	BC_SIG_LOCK;

	bc_num_init(&num1, len);
	bc_num_init(&num2, len);
	bc_num_setup(&half, half_digs, sizeof(half_digs) / sizeof(BcDig));

	// There is a division by two in the formula. We setup a number that's 1/2
	// so that we can use multiplication instead of heavy division.
	bc_num_one(&half);
	half.num[0] = BC_BASE_POW / 2;
	half.len = 1;
	BC_NUM_RDX_SET_NP(half, 1);
	half.scale = 1;

	bc_num_init(&f, len);
	bc_num_init(&fprime, len);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	// Pointers for easy switching.
	x0 = &num1;
	x1 = &num2;

	// Start with 1.
	bc_num_one(x0);

	// The power of the operand is needed for the estimate.
	pow = bc_num_intDigits(a);

	// The code in this if statement calculates the initial estimate. First, if
	// a is less than 0, then 0 is a good estimate. Otherwise, we want something
	// in the same ballpark. That ballpark is pow.
	if (pow) {

		// An odd number is served by starting with 2^((pow-1)/2), and an even
		// number is served by starting with 6^((pow-2)/2). Why? Because math.
		if (pow & 1) x0->num[0] = 2;
		else x0->num[0] = 6;

		pow -= 2 - (pow & 1);
		bc_num_shiftLeft(x0, pow / 2);
	}

	// I can set the rdx here directly because neg should be false.
	x0->scale = x0->rdx = 0;
	resscale = (scale + BC_BASE_DIGS) + 2;

	// This is the calculation loop. This compare goes to 0 eventually as the
	// difference between the two numbers gets smaller than resscale.
	while (bc_num_cmp(x1, x0)) {

		assert(BC_NUM_NONZERO(x0));

		// This loop directly corresponds to the iteration in Newton's method.
		// If you know the formula, this loop makes sense. Go study the formula.

		bc_num_div(a, x0, &f, resscale);
		bc_num_add(x0, &f, &fprime, resscale);

		assert(BC_NUM_RDX_VALID_NP(fprime));
		assert(BC_NUM_RDX_VALID_NP(half));

		bc_num_mul(&fprime, &half, x1, resscale);

		// Switch.
		temp = x0;
		x0 = x1;
		x1 = temp;
	}

	// Copy to the result and truncate.
	bc_num_copy(b, x0);
	if (b->scale > scale) bc_num_truncate(b, b->scale - scale);

	assert(!BC_NUM_NEG(b) || BC_NUM_NONZERO(b));
	assert(BC_NUM_RDX_VALID(b));
	assert(BC_NUM_RDX_VAL(b) <= b->len || !b->len);
	assert(!b->len || b->num[b->len - 1] || BC_NUM_RDX_VAL(b) == b->len);

err:
	BC_SIG_MAYLOCK;
	bc_num_free(&fprime);
	bc_num_free(&f);
	bc_num_free(&num2);
	bc_num_free(&num1);
	BC_LONGJMP_CONT;
}

void bc_num_divmod(BcNum *a, BcNum *b, BcNum *c, BcNum *d, size_t scale) {

	size_t ts, len;
	BcNum *ptr_a, num2;
	bool init = false;

	// The bulk of this function is just doing what bc_num_binary() does for the
	// binary operators. However, it assumes that only c and a can be equal.

	// Set up the parameters.
	ts = BC_MAX(scale + b->scale, a->scale);
	len = bc_num_mulReq(a, b, ts);

	assert(a != NULL && b != NULL && c != NULL && d != NULL);
	assert(c != d && a != d && b != d && b != c);

	// Initialize or expand as necessary.
	if (c == a) {

		memcpy(&num2, c, sizeof(BcNum));
		ptr_a = &num2;

		BC_SIG_LOCK;

		bc_num_init(c, len);

		init = true;

		BC_SETJMP_LOCKED(err);

		BC_SIG_UNLOCK;
	}
	else {
		ptr_a = a;
		bc_num_expand(c, len);
	}

	// Do the quick version if possible.
	if (BC_NUM_NONZERO(a) && !BC_NUM_RDX_VAL(a) &&
	    !BC_NUM_RDX_VAL(b) && b->len == 1 && !scale)
	{
		BcBigDig rem;

		bc_num_divArray(ptr_a, (BcBigDig) b->num[0], c, &rem);

		assert(rem < BC_BASE_POW);

		d->num[0] = (BcDig) rem;
		d->len = (rem != 0);
	}
	// Do the slow method.
	else bc_num_r(ptr_a, b, c, d, scale, ts);

	assert(!BC_NUM_NEG(c) || BC_NUM_NONZERO(c));
	assert(BC_NUM_RDX_VALID(c));
	assert(BC_NUM_RDX_VAL(c) <= c->len || !c->len);
	assert(!c->len || c->num[c->len - 1] || BC_NUM_RDX_VAL(c) == c->len);
	assert(!BC_NUM_NEG(d) || BC_NUM_NONZERO(d));
	assert(BC_NUM_RDX_VALID(d));
	assert(BC_NUM_RDX_VAL(d) <= d->len || !d->len);
	assert(!d->len || d->num[d->len - 1] || BC_NUM_RDX_VAL(d) == d->len);

err:
	// Only cleanup if we initialized.
	if (init) {
		BC_SIG_MAYLOCK;
		bc_num_free(&num2);
		BC_LONGJMP_CONT;
	}
}

void bc_num_modexp(BcNum *a, BcNum *b, BcNum *c, BcNum *restrict d) {

	BcNum base, exp, two, temp, atemp, btemp, ctemp;
	BcDig two_digs[2];

	assert(a != NULL && b != NULL && c != NULL && d != NULL);
	assert(a != d && b != d && c != d);

	if (BC_ERR(BC_NUM_ZERO(c))) bc_err(BC_ERR_MATH_DIVIDE_BY_ZERO);

	if (BC_ERR(BC_NUM_NEG(b))) bc_err(BC_ERR_MATH_NEGATIVE);

#ifndef NDEBUG
	// This is entirely for quieting a useless scan-build error.
	btemp.len = 0;
	ctemp.len = 0;
#endif // NDEBUG

	// Eliminate fractional parts that are zero or error if they are not zero.
	if (BC_ERR(bc_num_nonInt(a, &atemp) || bc_num_nonInt(b, &btemp) ||
	           bc_num_nonInt(c, &ctemp)))
	{
		bc_err(BC_ERR_MATH_NON_INTEGER);
	}

	bc_num_expand(d, ctemp.len);

	BC_SIG_LOCK;

	bc_num_init(&base, ctemp.len);
	bc_num_setup(&two, two_digs, sizeof(two_digs) / sizeof(BcDig));
	bc_num_init(&temp, btemp.len + 1);
	bc_num_createCopy(&exp, &btemp);

	BC_SETJMP_LOCKED(err);

	BC_SIG_UNLOCK;

	bc_num_one(&two);
	two.num[0] = 2;
	bc_num_one(d);

	// We already checked for 0.
	bc_num_rem(&atemp, &ctemp, &base, 0);

	// If you know the algorithm I used, the memory-efficient method, then this
	// loop should be self-explanatory because it is the calculation loop.
	while (BC_NUM_NONZERO(&exp)) {

		// Num two cannot be 0, so no errors.
		bc_num_divmod(&exp, &two, &exp, &temp, 0);

		if (BC_NUM_ONE(&temp) && !BC_NUM_NEG_NP(temp)) {

			assert(BC_NUM_RDX_VALID(d));
			assert(BC_NUM_RDX_VALID_NP(base));

			bc_num_mul(d, &base, &temp, 0);

			// We already checked for 0.
			bc_num_rem(&temp, &ctemp, d, 0);
		}

		assert(BC_NUM_RDX_VALID_NP(base));

		bc_num_mul(&base, &base, &temp, 0);

		// We already checked for 0.
		bc_num_rem(&temp, &ctemp, &base, 0);
	}

err:
	BC_SIG_MAYLOCK;
	bc_num_free(&exp);
	bc_num_free(&temp);
	bc_num_free(&base);
	BC_LONGJMP_CONT;
	assert(!BC_NUM_NEG(d) || d->len);
	assert(BC_NUM_RDX_VALID(d));
	assert(!d->len || d->num[d->len - 1] || BC_NUM_RDX_VAL(d) == d->len);
}

#if BC_DEBUG_CODE
void bc_num_printDebug(const BcNum *n, const char *name, bool emptyline) {
	bc_file_puts(&vm.fout, bc_flush_none, name);
	bc_file_puts(&vm.fout, bc_flush_none, ": ");
	bc_num_printDecimal(n, true);
	bc_file_putchar(&vm.fout, bc_flush_err, '\n');
	if (emptyline) bc_file_putchar(&vm.fout, bc_flush_err, '\n');
	vm.nchars = 0;
}

void bc_num_printDigs(const BcDig *n, size_t len, bool emptyline) {

	size_t i;

	for (i = len - 1; i < len; --i)
		bc_file_printf(&vm.fout, " %lu", (unsigned long) n[i]);

	bc_file_putchar(&vm.fout, bc_flush_err, '\n');
	if (emptyline) bc_file_putchar(&vm.fout, bc_flush_err, '\n');
	vm.nchars = 0;
}

void bc_num_printWithDigs(const BcNum *n, const char *name, bool emptyline) {
	bc_file_puts(&vm.fout, bc_flush_none, name);
	bc_file_printf(&vm.fout, " len: %zu, rdx: %zu, scale: %zu\n",
	               name, n->len, BC_NUM_RDX_VAL(n), n->scale);
	bc_num_printDigs(n->num, n->len, emptyline);
}

void bc_num_dump(const char *varname, const BcNum *n) {

	ulong i, scale = n->scale;

	bc_file_printf(&vm.ferr, "\n%s = %s", varname,
	               n->len ? (BC_NUM_NEG(n) ? "-" : "+") : "0 ");

	for (i = n->len - 1; i < n->len; --i) {

		if (i + 1 == BC_NUM_RDX_VAL(n))
			bc_file_puts(&vm.ferr, bc_flush_none, ". ");

		if (scale / BC_BASE_DIGS != BC_NUM_RDX_VAL(n) - i - 1)
			bc_file_printf(&vm.ferr, "%lu ", (unsigned long) n->num[i]);
		else {

			int mod = scale % BC_BASE_DIGS;
			int d = BC_BASE_DIGS - mod;
			BcDig div;

			if (mod != 0) {
				div = n->num[i] / ((BcDig) bc_num_pow10[(ulong) d]);
				bc_file_printf(&vm.ferr, "%lu", (unsigned long) div);
			}

			div = n->num[i] % ((BcDig) bc_num_pow10[(ulong) d]);
			bc_file_printf(&vm.ferr, " ' %lu ", (unsigned long) div);
		}
	}

	bc_file_printf(&vm.ferr, "(%zu | %zu.%zu / %zu) %lu\n",
	               n->scale, n->len, BC_NUM_RDX_VAL(n), n->cap,
	               (unsigned long) (void*) n->num);

	bc_file_flush(&vm.ferr, bc_flush_err);
}
#endif // BC_DEBUG_CODE