contrib/bc: merge version 5.1.0 from vendor branch

[FreeBSD/FreeBSD.git] / contrib / bc / 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 && vm.line_len) {
                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 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);

        // 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();

        if (BC_NUM_NONZERO(n)) {

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

                // Print the leading zero if necessary.
                if (BC_Z && BC_NUM_RDX_VAL(n) == n->len)
                        bc_num_printHex(0, 1, false, !newline);
        }

        // 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