biguintcore.d

/** Fundamental operations for arbitrary-precision arithmetic
 *
 * These functions are for internal use only.
 */
/*          Copyright Don Clugston 2008 - 2010.
 * Distributed under the Boost Software License, Version 1.0.
 *    (See accompanying file LICENSE_1_0.txt or copy at
 *          http://www.boost.org/LICENSE_1_0.txt)
 */
/* References:
   "Modern Computer Arithmetic" (MCA) is the primary reference for all
    algorithms used in this library.
  - R.P. Brent and P. Zimmermann, "Modern Computer Arithmetic",
    Version 0.5.9, (Oct 2010).
  - C. Burkinel and J. Ziegler, "Fast Recursive Division", MPI-I-98-1-022,
    Max-Planck Institute fuer Informatik, (Oct 1998).
  - G. Hanrot, M. Quercia, and P. Zimmermann, "The Middle Product Algorithm, I.",
    INRIA 4664, (Dec 2002).
  - M. Bodrato and A. Zanoni, "What about Toom-Cook Matrices Optimality?",
    http://bodrato.it/papers (2006).
  - A. Fog, "Optimizing subroutines in assembly language",
    www.agner.org/optimize (2008).
  - A. Fog, "The microarchitecture of Intel and AMD CPU's",
    www.agner.org/optimize (2008).
  - A. Fog, "Instruction tables: Lists of instruction latencies, throughputs
    and micro-operation breakdowns for Intel and AMD CPU's.", www.agner.org/optimize (2008).

Idioms:
  Many functions in this module use
  'func(Tulong)(Tulong x) if (is(Tulong == ulong))' rather than 'func(ulong x)'
  in order to disable implicit conversion.

*/
module std.internal.math.biguintcore;

version(D_InlineAsm_X86)
{
    import std.internal.math.biguintx86;
}
else
{
    import std.internal.math.biguintnoasm;
}

alias multibyteAdd = multibyteAddSub!('+');
alias multibyteSub = multibyteAddSub!('-');


private import core.cpuid;
private import std.traits : Unqual;

shared static this()
{
    CACHELIMIT = core.cpuid.datacache[0].size*1024/2;
    FASTDIVLIMIT = 100;
}

private:
// Limits for when to switch between algorithms.
immutable size_t CACHELIMIT;   // Half the size of the data cache.
immutable size_t FASTDIVLIMIT; // crossover to recursive division


// These constants are used by shift operations
static if (BigDigit.sizeof == int.sizeof)
{
    enum { LG2BIGDIGITBITS = 5, BIGDIGITSHIFTMASK = 31 };
    alias BIGHALFDIGIT = ushort;
}
else static if (BigDigit.sizeof == long.sizeof)
{
    alias BIGHALFDIGIT = uint;
    enum { LG2BIGDIGITBITS = 6, BIGDIGITSHIFTMASK = 63 };
}
else static assert(0, "Unsupported BigDigit size");

private import std.exception : assumeUnique;
private import std.traits:isIntegral;
enum BigDigitBits = BigDigit.sizeof*8;
template maxBigDigits(T) if (isIntegral!T)
{
    enum maxBigDigits = (T.sizeof+BigDigit.sizeof-1)/BigDigit.sizeof;
}

enum immutable(BigDigit) [] ZERO = [0];
enum immutable(BigDigit) [] ONE = [1];
enum immutable(BigDigit) [] TWO = [2];
enum immutable(BigDigit) [] TEN = [10];


public:

/// BigUint performs memory management and wraps the low-level calls.
struct BigUint
{
private:
    pure invariant()
    {
        assert( data.length >= 1 && (data.length == 1 || data[$-1] != 0 ));
    }
    immutable(BigDigit) [] data = ZERO;
    this(immutable(BigDigit) [] x) pure
    {
       data = x;
    }
    this(T)(T x) pure if (isIntegral!T)
    {
        opAssign(x);
    }
public:
    // Length in uints
    size_t uintLength() pure const
    {
        static if (BigDigit.sizeof == uint.sizeof)
        {
            return data.length;
        }
        else static if (BigDigit.sizeof == ulong.sizeof)
        {
            return data.length * 2 -
            ((data[$-1] & 0xFFFF_FFFF_0000_0000L) ? 1 : 0);
        }
    }
    size_t ulongLength() pure const
    {
        static if (BigDigit.sizeof == uint.sizeof)
        {
            return (data.length + 1) >> 1;
        }
        else static if (BigDigit.sizeof == ulong.sizeof)
        {
            return data.length;
        }
    }

    // The value at (cast(ulong[])data)[n]
    ulong peekUlong(int n) pure const
    {
        static if (BigDigit.sizeof == int.sizeof)
        {
            if (data.length == n*2 + 1) return data[n*2];
            return data[n*2] + ((cast(ulong)data[n*2 + 1]) << 32 );
        }
        else static if (BigDigit.sizeof == long.sizeof)
        {
            return data[n];
        }
    }
    uint peekUint(int n) pure const
    {
        static if (BigDigit.sizeof == int.sizeof)
        {
            return data[n];
        }
        else
        {
            ulong x = data[n >> 1];
            return (n & 1) ? cast(uint)(x >> 32) : cast(uint)x;
        }
    }
public:
    ///
    void opAssign(Tulong)(Tulong u) pure if (is (Tulong == ulong))
    {
        if (u == 0) data = ZERO;
        else if (u == 1) data = ONE;
        else if (u == 2) data = TWO;
        else if (u == 10) data = TEN;
        else
        {
            static if (BigDigit.sizeof == int.sizeof)
            {
                uint ulo = cast(uint)(u & 0xFFFF_FFFF);
                uint uhi = cast(uint)(u >> 32);
                if (uhi == 0)
                {
                    data = [ulo];
                }
                else
                {
                    data = [ulo, uhi];
                }
            }
            else static if (BigDigit.sizeof == long.sizeof)
            {
                data = [u];
            }
        }
    }
    void opAssign(Tdummy = void)(BigUint y) pure
    {
        this.data = y.data;
    }

    ///
    int opCmp(Tdummy = void)(const BigUint y) pure const
    {
        if (data.length != y.data.length)
            return (data.length > y.data.length) ?  1 : -1;
        size_t k = highestDifferentDigit(data, y.data);
        if (data[k] == y.data[k])
            return 0;
        return data[k] > y.data[k] ? 1 : -1;
    }

    ///
    int opCmp(Tulong)(Tulong y) pure const if(is (Tulong == ulong))
    {
        if (data.length > maxBigDigits!Tulong)
            return 1;

        foreach_reverse (i; 0 .. maxBigDigits!Tulong)
        {
            BigDigit tmp = cast(BigDigit)(y>>(i*BigDigitBits));
            if (tmp == 0)
                if (data.length >= i+1)
                {
                    // Since ZERO is [0], so we cannot simply return 1 here, as
                    // data[i] would be 0 for i==0 in that case.
                    return (data[i] > 0) ? 1 : 0;
                }
                else
                    continue;
            else
                if (i+1 > data.length)
                    return -1;
                else if (tmp != data[i])
                    return data[i] > tmp ? 1 : -1;
        }
        return 0;
    }

    bool opEquals(Tdummy = void)(ref const BigUint y) pure const
    {
           return y.data[] == data[];
    }

    bool opEquals(Tdummy = void)(ulong y) pure const
    {
        if (data.length > 2)
            return false;
        uint ylo = cast(uint)(y & 0xFFFF_FFFF);
        uint yhi = cast(uint)(y >> 32);
        if (data.length==2 && data[1]!=yhi)
            return false;
        if (data.length==1 && yhi!=0)
            return false;
        return (data[0] == ylo);
    }

    bool isZero() pure const nothrow @safe
    {
        return data.length == 1 && data[0] == 0;
    }

    size_t numBytes() pure const
    {
        return data.length * BigDigit.sizeof;
    }

    // the extra bytes are added to the start of the string
    char [] toDecimalString(int frontExtraBytes) const pure
    {
        auto predictlength = 20+20*(data.length/2); // just over 19
        char [] buff = new char[frontExtraBytes + predictlength];
        ptrdiff_t sofar = biguintToDecimal(buff, data.dup);
        return buff[sofar-frontExtraBytes..$];
    }

    /** Convert to a hex string, printing a minimum number of digits 'minPadding',
     *  allocating an additional 'frontExtraBytes' at the start of the string.
     *  Padding is done with padChar, which may be '0' or ' '.
     *  'separator' is a digit separation character. If non-zero, it is inserted
     *  between every 8 digits.
     *  Separator characters do not contribute to the minPadding.
     */
    char [] toHexString(int frontExtraBytes, char separator = 0,
            int minPadding=0, char padChar = '0') const pure
    {
        // Calculate number of extra padding bytes
        size_t extraPad = (minPadding > data.length * 2 * BigDigit.sizeof)
            ? minPadding - data.length * 2 * BigDigit.sizeof : 0;

        // Length not including separator bytes
        size_t lenBytes = data.length * 2 * BigDigit.sizeof;

        // Calculate number of separator bytes
        size_t mainSeparatorBytes = separator ? (lenBytes  / 8) - 1 : 0;
        size_t totalSeparatorBytes = separator ? ((extraPad + lenBytes + 7) / 8) - 1: 0;

        char [] buff = new char[lenBytes + extraPad + totalSeparatorBytes + frontExtraBytes];
        biguintToHex(buff[$ - lenBytes - mainSeparatorBytes .. $], data, separator);
        if (extraPad > 0)
        {
            if (separator)
            {
                size_t start = frontExtraBytes; // first index to pad
                if (extraPad &7)
                {
                    // Do 1 to 7 extra zeros.
                    buff[frontExtraBytes .. frontExtraBytes + (extraPad & 7)] = padChar;
                    buff[frontExtraBytes + (extraPad & 7)] = (padChar == ' ' ? ' ' : separator);
                    start += (extraPad & 7) + 1;
                }
                for (int i=0; i< (extraPad >> 3); ++i)
                {
                    buff[start .. start + 8] = padChar;
                    buff[start + 8] = (padChar == ' ' ? ' ' : separator);
                    start += 9;
                }
            }
            else
            {
                buff[frontExtraBytes .. frontExtraBytes + extraPad]=padChar;
            }
        }
        int z = frontExtraBytes;
        if (lenBytes > minPadding)
        {
            // Strip leading zeros.
            ptrdiff_t maxStrip = lenBytes - minPadding;
            while (z< buff.length-1 && (buff[z]=='0' || buff[z]==padChar) && maxStrip>0)
            {
                ++z;
                --maxStrip;
            }
        }
        if (padChar!='0')
        {
            // Convert leading zeros into padChars.
            for (size_t k= z; k< buff.length-1 && (buff[k]=='0' || buff[k]==padChar); ++k)
            {
                if (buff[k]=='0') buff[k]=padChar;
            }
        }
        return buff[z-frontExtraBytes..$];
    }

    // return false if invalid character found
    bool fromHexString(const(char)[] s) pure
    {
        //Strip leading zeros
        int firstNonZero = 0;
        while ((firstNonZero < s.length - 1) &&
            (s[firstNonZero]=='0' || s[firstNonZero]=='_'))
        {
                ++firstNonZero;
        }
        auto len = (s.length - firstNonZero + 15)/4;
        auto tmp = new BigDigit[len+1];
        uint part = 0;
        uint sofar = 0;
        uint partcount = 0;
        assert(s.length>0);
        for (ptrdiff_t i = s.length - 1; i>=firstNonZero; --i)
        {
            assert(i>=0);
            char c = s[i];
            if (s[i]=='_') continue;
            uint x = (c>='0' && c<='9') ? c - '0'
                   : (c>='A' && c<='F') ? c - 'A' + 10
                   : (c>='a' && c<='f') ? c - 'a' + 10
                   : 100;
            if (x==100) return false;
            part >>= 4;
            part |= (x<<(32-4));
            ++partcount;
            if (partcount==8)
            {
                tmp[sofar] = part;
                ++sofar;
                partcount = 0;
                part = 0;
            }
        }
        if (part)
        {
            for ( ; partcount != 8; ++partcount) part >>= 4;
            tmp[sofar] = part;
            ++sofar;
        }
        if (sofar == 0) data = ZERO;
        else data = assumeUnique(tmp[0..sofar]);
        return true;
    }

    // return true if OK; false if erroneous characters found
    bool fromDecimalString(const(char)[] s) pure
    {
        //Strip leading zeros
        int firstNonZero = 0;
        while ((firstNonZero < s.length) &&
            (s[firstNonZero]=='0' || s[firstNonZero]=='_'))
        {
                ++firstNonZero;
        }
        if (firstNonZero == s.length && s.length >= 1)
        {
            data = ZERO;
            return true;
        }
        auto predictlength = (18*2 + 2*(s.length-firstNonZero)) / 19;
        auto tmp = new BigDigit[predictlength];
        uint hi = biguintFromDecimal(tmp, s[firstNonZero..$]);
        tmp.length = hi;
        data = assumeUnique(tmp);
        return true;
    }

    ////////////////////////
    //
    // All of these member functions create a new BigUint.

    // return x >> y
    BigUint opShr(Tulong)(Tulong y) pure const if (is (Tulong == ulong))
    {
        assert(y>0);
        uint bits = cast(uint)y & BIGDIGITSHIFTMASK;
        if ((y>>LG2BIGDIGITBITS) >= data.length) return BigUint(ZERO);
        uint words = cast(uint)(y >> LG2BIGDIGITBITS);
        if (bits==0)
        {
            return BigUint(data[words..$]);
        }
        else
        {
            uint [] result = new BigDigit[data.length - words];
            multibyteShr(result, data[words..$], bits);
            if (result.length>1 && result[$-1]==0) return BigUint(assumeUnique(result[0..$-1]));
            else return BigUint(assumeUnique(result));
        }
    }

    // return x << y
    BigUint opShl(Tulong)(Tulong y) pure const if (is (Tulong == ulong))
    {
        assert(y>0);
        if (isZero()) return this;
        uint bits = cast(uint)y & BIGDIGITSHIFTMASK;
        assert ((y>>LG2BIGDIGITBITS) < cast(ulong)(uint.max));
        uint words = cast(uint)(y >> LG2BIGDIGITBITS);
        BigDigit [] result = new BigDigit[data.length + words+1];
        result[0..words] = 0;
        if (bits==0)
        {
            result[words..words+data.length] = data[];
            return BigUint(assumeUnique(result[0..words+data.length]));
        }
        else
        {
            uint c = multibyteShl(result[words..words+data.length], data, bits);
            if (c==0) return BigUint(assumeUnique(result[0..words+data.length]));
            result[$-1] = c;
            return BigUint(assumeUnique(result));
        }
    }

    // If wantSub is false, return x + y, leaving sign unchanged
    // If wantSub is true, return abs(x - y), negating sign if x < y
    static BigUint addOrSubInt(Tulong)(const BigUint x, Tulong y,
            bool wantSub, ref bool sign) pure if (is(Tulong == ulong))
    {
        BigUint r;
        if (wantSub)
        {   // perform a subtraction
            if (x.data.length > 2)
            {
                r.data = subInt(x.data, y);
            }
            else
            {   // could change sign!
                ulong xx = x.data[0];
                if (x.data.length > 1)
                    xx += (cast(ulong)x.data[1]) << 32;
                ulong d;
                if (xx <= y)
                {
                    d = y - xx;
                    sign = !sign;
                }
                else
                {
                    d = xx - y;
                }
                if (d == 0)
                {
                    r = 0UL;
                    sign = false;
                    return r;
                }
                if (d > uint.max)
                {
                    r.data = [cast(uint)(d & 0xFFFF_FFFF), cast(uint)(d>>32)];
                }
                else
                {
                    r.data = [cast(uint)(d & 0xFFFF_FFFF)];
                }
            }
        }
        else
        {
            r.data = addInt(x.data, y);
        }
        return r;
    }

    // If wantSub is false, return x + y, leaving sign unchanged.
    // If wantSub is true, return abs(x - y), negating sign if x<y
    static BigUint addOrSub(BigUint x, BigUint y, bool wantSub, bool *sign)
        pure
    {
        BigUint r;
        if (wantSub)
        {   // perform a subtraction
            bool negative;
            r.data = sub(x.data, y.data, &negative);
            *sign ^= negative;
            if (r.isZero())
            {
                *sign = false;
            }
        }
        else
        {
            r.data = add(x.data, y.data);
        }
        return r;
    }


    //  return x*y.
    //  y must not be zero.
    static BigUint mulInt(T = ulong)(BigUint x, T y) pure
    {
        if (y==0 || x == 0) return BigUint(ZERO);
        uint hi = cast(uint)(y >>> 32);
        uint lo = cast(uint)(y & 0xFFFF_FFFF);
        uint [] result = new BigDigit[x.data.length+1+(hi!=0)];
        result[x.data.length] = multibyteMul(result[0..x.data.length], x.data, lo, 0);
        if (hi!=0)
        {
            result[x.data.length+1] = multibyteMulAdd!('+')(result[1..x.data.length+1],
                x.data, hi, 0);
        }
        return BigUint(removeLeadingZeros(assumeUnique(result)));
    }

    /*  return x * y.
     */
    static BigUint mul(BigUint x, BigUint y) pure
    {
        if (y==0 || x == 0)
            return BigUint(ZERO);
        auto len = x.data.length + y.data.length;
        BigDigit [] result = new BigDigit[len];
        if (y.data.length > x.data.length)
        {
            mulInternal(result, y.data, x.data);
        }
        else
        {
            if (x.data[]==y.data[]) squareInternal(result, x.data);
            else mulInternal(result, x.data, y.data);
        }
        // the highest element could be zero,
        // in which case we need to reduce the length
        return BigUint(removeLeadingZeros(assumeUnique(result)));
    }

    // return x / y
    static BigUint divInt(T)(BigUint x, T y_) pure if ( is(Unqual!T == uint) )
    {
        uint y = y_;
        if (y == 1)
            return x;
        uint [] result = new BigDigit[x.data.length];
        if ((y&(-y))==y)
        {
            assert(y!=0, "BigUint division by zero");
            // perfect power of 2
            uint b = 0;
            for (;y!=1; y>>=1)
            {
                ++b;
            }
            multibyteShr(result, x.data, b);
        }
        else
        {
            result[] = x.data[];
            uint rem = multibyteDivAssign(result, y, 0);
        }
        return BigUint(removeLeadingZeros(assumeUnique(result)));
    }

    // return x % y
    static uint modInt(T)(BigUint x, T y_) pure if ( is(Unqual!T == uint) )
    {
        uint y = y_;
        assert(y!=0);
        if ((y&(-y)) == y)
        {   // perfect power of 2
            return x.data[0] & (y-1);
        }
        else
        {
            // horribly inefficient - malloc, copy, & store are unnecessary.
            uint [] wasteful = new BigDigit[x.data.length];
            wasteful[] = x.data[];
            uint rem = multibyteDivAssign(wasteful, y, 0);
            delete wasteful;
            return rem;
        }
    }

    // return x / y
    static BigUint div(BigUint x, BigUint y) pure
    {
        if (y.data.length > x.data.length)
            return BigUint(ZERO);
        if (y.data.length == 1)
            return divInt(x, y.data[0]);
        BigDigit [] result = new BigDigit[x.data.length - y.data.length + 1];
        divModInternal(result, null, x.data, y.data);
        return BigUint(removeLeadingZeros(assumeUnique(result)));
    }

    // return x % y
    static BigUint mod(BigUint x, BigUint y) pure
    {
        if (y.data.length > x.data.length) return x;
        if (y.data.length == 1)
        {
            return BigUint([modInt(x, y.data[0])]);
        }
        BigDigit [] result = new BigDigit[x.data.length - y.data.length + 1];
        BigDigit [] rem = new BigDigit[y.data.length];
        divModInternal(result, rem, x.data, y.data);
        return BigUint(removeLeadingZeros(assumeUnique(rem)));
    }

    // return x op y
    static BigUint bitwiseOp(string op)(BigUint x, BigUint y, bool xSign, bool ySign, ref bool resultSign) pure if (op == "|" || op == "^" || op == "&")
    {
        auto d1 = includeSign(x.data, y.uintLength, xSign);
        auto d2 = includeSign(y.data, x.uintLength, ySign);

        foreach (i; 0..d1.length)
        {
            mixin("d1[i] " ~ op ~ "= d2[i];");
        }

        mixin("resultSign = xSign " ~ op ~ " ySign;");

        if (resultSign) {
            twosComplement(d1, d1);
        }

        return BigUint(removeLeadingZeros(assumeUnique(d1)));
    }

    /**
     * Return a BigUint which is x raised to the power of y.
     * Method: Powers of 2 are removed from x, then left-to-right binary
     * exponentiation is used.
     * Memory allocation is minimized: at most one temporary BigUint is used.
     */
    static BigUint pow(BigUint x, ulong y) pure
    {
        // Deal with the degenerate cases first.
        if (y==0) return BigUint(ONE);
        if (y==1) return x;
        if (x==0 || x==1) return x;

        BigUint result;

        // Simplify, step 1: Remove all powers of 2.
        uint firstnonzero = firstNonZeroDigit(x.data);
        // Now we know x = x[firstnonzero..$] * (2^^(firstnonzero*BigDigitBits))
        // where BigDigitBits = BigDigit.sizeof * 8

        // See if x[firstnonzero..$] can now fit into a single digit.
        bool singledigit = ((x.data.length - firstnonzero) == 1);
        // If true, then x0 is that digit
        // and the result will be (x0 ^^ y) * (2^^(firstnonzero*y*BigDigitBits))
        BigDigit x0 = x.data[firstnonzero];
        assert(x0 !=0);
        // Length of the non-zero portion
        size_t nonzerolength = x.data.length - firstnonzero;
        ulong y0;
        uint evenbits = 0; // number of even bits in the bottom of x
        while (!(x0 & 1))
        {
            x0 >>= 1;
            ++evenbits;
        }

        if ((x.data.length- firstnonzero == 2))
        {
            // Check for a single digit straddling a digit boundary
            BigDigit x1 = x.data[firstnonzero+1];
            if ((x1 >> evenbits) == 0)
            {
                x0 |= (x1 << (BigDigit.sizeof * 8 - evenbits));
                singledigit = true;
            }
        }
        // Now if (singledigit), x^^y  = (x0 ^^ y) * 2^^(evenbits * y) * 2^^(firstnonzero*y*BigDigitBits))

        uint evenshiftbits = 0; // Total powers of 2 to shift by, at the end

        // Simplify, step 2: For singledigits, see if we can trivially reduce y

        BigDigit finalMultiplier = 1UL;

        if (singledigit)
        {
            // x fits into a single digit. Raise it to the highest power we can
            // that still fits into a single digit, then reduce the exponent accordingly.
            // We're quite likely to have a residual multiply at the end.
            // For example, 10^^100 = (((5^^13)^^7) * 5^^9) * 2^^100.
            // and 5^^13 still fits into a uint.
            evenshiftbits  = cast(uint)( (evenbits * y) & BIGDIGITSHIFTMASK);
            if (x0 == 1)
            {   // Perfect power of 2
                result = 1UL;
                return result << (evenbits + firstnonzero * 8 * BigDigit.sizeof) * y;
            }
            int p = highestPowerBelowUintMax(x0);
            if (y <= p)
            {   // Just do it with pow
                result = cast(ulong)intpow(x0, y);
                if (evenbits + firstnonzero == 0)
                    return result;
                return result << (evenbits + firstnonzero * 8 * BigDigit.sizeof) * y;
            }
            y0 = y / p;
            finalMultiplier = intpow(x0, y - y0*p);
            x0 = intpow(x0, p);
            // Result is x0
            nonzerolength = 1;
        }
        // Now if (singledigit), x^^y  = finalMultiplier * (x0 ^^ y0) * 2^^(evenbits * y) * 2^^(firstnonzero*y*BigDigitBits))

        // Perform a crude check for overflow and allocate result buffer.
        // The length required is y * lg2(x) bits.
        // which will always fit into y*x.length digits. But this is
        // a gross overestimate if x is small (length 1 or 2) and the highest
        // digit is nearly empty.
        // A better estimate is:
        //   y * lg2(x[$-1]/BigDigit.max) + y * (x.length - 1) digits,
        //  and the first term is always between
        //  y * (bsr(x.data[$-1]) + 1) / BIGDIGITBITS and
        //  y * (bsr(x.data[$-1]) + 2) / BIGDIGITBITS
        // For single digit payloads, we already have
        //   x^^y  = finalMultiplier * (x0 ^^ y0) * 2^^(evenbits * y) * 2^^(firstnonzero*y*BigDigitBits))
        // and x0 is almost a full digit, so it's a tight estimate.
        // Number of digits is therefore 1 + x0.length*y0 + (evenbits*y)/BIGDIGIT + firstnonzero*y
        // Note that the divisions must be rounded up.

        // Estimated length in BigDigits
        ulong estimatelength = singledigit
            ? 1 + y0 + ((evenbits*y  + BigDigit.sizeof * 8 - 1) / (BigDigit.sizeof *8)) + firstnonzero*y
            :  x.data.length * y;
        // Imprecise check for overflow. Makes the extreme cases easier to debug
        // (less extreme overflow will result in an out of memory error).
        if (estimatelength > uint.max/(4*BigDigit.sizeof))
            assert(0, "Overflow in BigInt.pow");

        // The result buffer includes space for all the trailing zeros
        BigDigit [] resultBuffer = new BigDigit[cast(size_t)estimatelength];

        // Do all the powers of 2!
        size_t result_start = cast(size_t)( firstnonzero * y
            + (singledigit ? ((evenbits * y) >> LG2BIGDIGITBITS) : 0));

        resultBuffer[0..result_start] = 0;
        BigDigit [] t1 = resultBuffer[result_start..$];
        BigDigit [] r1;

        if (singledigit)
        {
            r1 = t1[0..1];
            r1[0] = x0;
            y = y0;
        }
        else
        {
            // It's not worth right shifting by evenbits unless we also shrink the length after each
            // multiply or squaring operation. That might still be worthwhile for large y.
            r1 = t1[0..x.data.length - firstnonzero];
            r1[0..$] = x.data[firstnonzero..$];
        }

        if (y>1)
        {   // Set r1 = r1 ^^ y.
            // The secondary buffer only needs space for the multiplication results
            BigDigit [] secondaryBuffer = new BigDigit[resultBuffer.length - result_start];
            BigDigit [] t2 = secondaryBuffer;
            BigDigit [] r2;

            int shifts = 63; // num bits in a long
            while(!(y & 0x8000_0000_0000_0000L))
            {
                y <<= 1;
                --shifts;
            }
            y <<=1;

            while(y!=0)
            {
                // For each bit of y: Set r1 =  r1 * r1
                // If the bit is 1, set r1 = r1 * x
                // Eg, if y is 0b101, result = ((x^^2)^^2)*x == x^^5.
                // Optimization opportunity: if more than 2 bits in y are set,
                // it's usually possible to reduce the number of multiplies
                // by caching odd powers of x. eg for y = 54,
                // (0b110110), set u = x^^3, and result is ((u^^8)*u)^^2
                r2 = t2[0 .. r1.length*2];
                squareInternal(r2, r1);
                if (y & 0x8000_0000_0000_0000L)
                {
                    r1 = t1[0 .. r2.length + nonzerolength];
                    if (singledigit)
                    {
                        r1[$-1] = multibyteMul(r1[0 .. $-1], r2, x0, 0);
                    }
                    else
                    {
                        mulInternal(r1, r2, x.data[firstnonzero..$]);
                    }
                }
                else
                {
                    r1 = t1[0 .. r2.length];
                    r1[] = r2[];
                }
                y <<=1;
                shifts--;
            }
            while (shifts>0)
            {
                r2 = t2[0 .. r1.length * 2];
                squareInternal(r2, r1);
                r1 = t1[0 .. r2.length];
                r1[] = r2[];
                --shifts;
            }
        }

        if (finalMultiplier!=1)
        {
            BigDigit carry = multibyteMul(r1, r1, finalMultiplier, 0);
            if (carry)
            {
                r1 = t1[0 .. r1.length + 1];
                r1[$-1] = carry;
            }
        }
        if (evenshiftbits)
        {
            BigDigit carry = multibyteShl(r1, r1, evenshiftbits);
            if (carry!=0)
            {
                r1 = t1[0 .. r1.length + 1];
                r1[$ - 1] = carry;
            }
        }
        while(r1[$ - 1]==0)
        {
            r1=r1[0 .. $ - 1];
        }
        return BigUint(assumeUnique(resultBuffer[0 .. result_start + r1.length]));
    }

    // Implement toHash so that BigUint works properly as an AA key.
    size_t toHash() const @trusted nothrow
    {
        return typeid(data).getHash(&data);
    }

} // end BigUint

unittest
{
    // ulong comparison test
    BigUint a = [1];
    assert(a == 1);
    assert(a < 0x8000_0000_0000_0000UL); // bug 9548

    // bug 12234
    BigUint z = [0];
    assert(z == 0UL);
    assert(!(z > 0UL));
    assert(!(z < 0UL));
}

// Remove leading zeros from x, to restore the BigUint invariant
inout(BigDigit) [] removeLeadingZeros(inout(BigDigit) [] x) pure
{
    size_t k = x.length;
    while(k>1 && x[k - 1]==0) --k;
    return x[0 .. k];
}

unittest
{
   BigUint r = BigUint([5]);
   BigUint t = BigUint([7]);
   BigUint s = BigUint.mod(r, t);
   assert(s==5);
}


unittest
{
    BigUint r;
    r = 5UL;
    assert(r.peekUlong(0) == 5UL);
    assert(r.peekUint(0) == 5U);
    r = 0x1234_5678_9ABC_DEF0UL;
    assert(r.peekUlong(0) == 0x1234_5678_9ABC_DEF0UL);
    assert(r.peekUint(0) == 0x9ABC_DEF0U);
}


// Pow tests
unittest
{
    BigUint r, s;
    r.fromHexString("80000000_00000001");
    s = BigUint.pow(r, 5);
    r.fromHexString("08000000_00000000_50000000_00000001_40000000_00000002_80000000"
      ~ "_00000002_80000000_00000001");
    assert(s == r);
    s = 10UL;
    s = BigUint.pow(s, 39);
    r.fromDecimalString("1000000000000000000000000000000000000000");
    assert(s == r);
    r.fromHexString("1_E1178E81_00000000");
    s = BigUint.pow(r, 15); // Regression test: this used to overflow array bounds

    r.fromDecimalString("000_000_00");
    assert(r == 0);
    r.fromDecimalString("0007");
    assert(r == 7);
    r.fromDecimalString("0");
    assert(r == 0);
}

// Radix conversion tests
unittest
{
    BigUint r;
    r.fromHexString("1_E1178E81_00000000");
    assert(r.toHexString(0, '_', 0) == "1_E1178E81_00000000");
    assert(r.toHexString(0, '_', 20) == "0001_E1178E81_00000000");
    assert(r.toHexString(0, '_', 16+8) == "00000001_E1178E81_00000000");
    assert(r.toHexString(0, '_', 16+9) == "0_00000001_E1178E81_00000000");
    assert(r.toHexString(0, '_', 16+8+8) ==   "00000000_00000001_E1178E81_00000000");
    assert(r.toHexString(0, '_', 16+8+8+1) ==      "0_00000000_00000001_E1178E81_00000000");
    assert(r.toHexString(0, '_', 16+8+8+1, ' ') == "                  1_E1178E81_00000000");
    assert(r.toHexString(0, 0, 16+8+8+1) == "00000000000000001E1178E8100000000");
    r = 0UL;
    assert(r.toHexString(0, '_', 0) == "0");
    assert(r.toHexString(0, '_', 7) == "0000000");
    assert(r.toHexString(0, '_', 7, ' ') == "      0");
    assert(r.toHexString(0, '#', 9) == "0#00000000");
    assert(r.toHexString(0, 0, 9) == "000000000");

}


private:
void twosComplement(const(BigDigit) [] x, BigDigit[] result) pure
{
    foreach (i; 0..x.length)
    {
        result[i] = ~x[i];
    }
    result[x.length..$] = BigDigit.max;

    bool sgn = false;

    foreach (i; 0..result.length) {
        if (result[i] == BigDigit.max) {
            result[i] = 0;
        } else {
            result[i] += 1;
            break;
        }
    }
}

// Encode BigInt as BigDigit array (sign and 2's complement)
BigDigit[] includeSign(const(BigDigit) [] x, size_t minSize, bool sign) pure
{
    size_t length = (x.length > minSize) ? x.length : minSize;
    BigDigit [] result = new BigDigit[length];
    if (sign)
    {
        twosComplement(x, result);
    }
    else
    {
        result[0..x.length] = x;
    }
    return result;
}

// works for any type
T intpow(T)(T x, ulong n) pure
{
    T p;

    switch (n)
    {
    case 0:
        p = 1;
        break;

    case 1:
        p = x;
        break;

    case 2:
        p = x * x;
        break;

    default:
        p = 1;
        while (1){
            if (n & 1)
                p *= x;
            n >>= 1;
            if (!n)
                break;
            x *= x;
        }
        break;
    }
    return p;
}


//  returns the maximum power of x that will fit in a uint.
int highestPowerBelowUintMax(uint x) pure
{
     assert(x>1);
     static immutable ubyte [22] maxpwr = [ 31, 20, 15, 13, 12, 11, 10, 10, 9, 9,
                                          8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 7];
     if (x<24) return maxpwr[x-2];
     if (x<41) return 6;
     if (x<85) return 5;
     if (x<256) return 4;
     if (x<1626) return 3;
     if (x<65536) return 2;
     return 1;
}

//  returns the maximum power of x that will fit in a ulong.
int highestPowerBelowUlongMax(uint x) pure
{
     assert(x>1);
     static immutable ubyte [39] maxpwr = [ 63, 40, 31, 27, 24, 22, 21, 20, 19, 18,
                                         17, 17, 16, 16, 15, 15, 15, 15, 14, 14,
                                         14, 14, 13, 13, 13, 13, 13, 13, 13, 12,
                                         12, 12, 12, 12, 12, 12, 12, 12, 12];
     if (x<41) return maxpwr[x-2];
     if (x<57) return 11;
     if (x<85) return 10;
     if (x<139) return 9;
     if (x<256) return 8;
     if (x<566) return 7;
     if (x<1626) return 6;
     if (x<7132) return 5;
     if (x<65536) return 4;
     if (x<2642246) return 3;
     return 2;
}

version(unittest) {

int slowHighestPowerBelowUintMax(uint x) pure
{
     int pwr = 1;
     for (ulong q = x;x*q < cast(ulong)uint.max; ) {
         q*=x; ++pwr;
     }
     return pwr;
}

unittest
{
  assert(highestPowerBelowUintMax(10)==9);
  for (int k=82; k<88; ++k) {assert(highestPowerBelowUintMax(k)== slowHighestPowerBelowUintMax(k)); }
}

}


/*  General unsigned subtraction routine for bigints.
 *  Sets result = x - y. If the result is negative, negative will be true.
 */
BigDigit [] sub(const BigDigit [] x, const BigDigit [] y, bool *negative) pure
{
    if (x.length == y.length)
    {
        // There's a possibility of cancellation, if x and y are almost equal.
        ptrdiff_t last = highestDifferentDigit(x, y);
        BigDigit [] result = new BigDigit[last+1];
        if (x[last] < y[last])
        {   // we know result is negative
            multibyteSub(result[0..last+1], y[0..last+1], x[0..last+1], 0);
            *negative = true;
        }
        else
        {   // positive or zero result
            multibyteSub(result[0..last+1], x[0..last+1], y[0..last+1], 0);
            *negative = false;
        }
        while (result.length > 1 && result[$-1] == 0)
        {
            result = result[0..$-1];
        }
//        if (result.length >1 && result[$-1]==0) return result[0..$-1];
        return result;
    }
    // Lengths are different
    const(BigDigit) [] large, small;
    if (x.length < y.length)
    {
        *negative = true;
        large = y; small = x;
    }
    else
    {
        *negative = false;
        large = x; small = y;
    }
    // result.length will be equal to larger length, or could decrease by 1.


    BigDigit [] result = new BigDigit[large.length];
    BigDigit carry = multibyteSub(result[0..small.length], large[0..small.length], small, 0);
    result[small.length..$] = large[small.length..$];
    if (carry)
    {
        multibyteIncrementAssign!('-')(result[small.length..$], carry);
    }
    while (result.length > 1 && result[$-1] == 0)
    {
        result = result[0..$-1];
    }
    return result;
}


// return a + b
BigDigit [] add(const BigDigit [] a, const BigDigit [] b) pure
{
    const(BigDigit) [] x, y;
    if (a.length < b.length)
    {
        x = b; y = a;
    }
    else
    {
        x = a; y = b;
    }
    // now we know x.length > y.length
    // create result. add 1 in case it overflows
    BigDigit [] result = new BigDigit[x.length + 1];

    BigDigit carry = multibyteAdd(result[0..y.length], x[0..y.length], y, 0);
    if (x.length != y.length)
    {
        result[y.length..$-1]= x[y.length..$];
        carry  = multibyteIncrementAssign!('+')(result[y.length..$-1], carry);
    }
    if (carry)
    {
        result[$-1] = carry;
        return result;
    }
    else
        return result[0..$-1];
}

/**  return x + y
 */
BigDigit [] addInt(const BigDigit[] x, ulong y) pure
{
    uint hi = cast(uint)(y >>> 32);
    uint lo = cast(uint)(y& 0xFFFF_FFFF);
    auto len = x.length;
    if (x.length < 2 && hi!=0) ++len;
    BigDigit [] result = new BigDigit[len+1];
    result[0..x.length] = x[];
    if (x.length < 2 && hi!=0) { result[1]=hi; hi=0; }
    uint carry = multibyteIncrementAssign!('+')(result[0..$-1], lo);
    if (hi!=0) carry += multibyteIncrementAssign!('+')(result[1..$-1], hi);
    if (carry)
    {
        result[$-1] = carry;
        return result;
    }
    else
        return result[0..$-1];
}

/** Return x - y.
 *  x must be greater than y.
 */
BigDigit [] subInt(const BigDigit[] x, ulong y) pure
{
    uint hi = cast(uint)(y >>> 32);
    uint lo = cast(uint)(y & 0xFFFF_FFFF);
    BigDigit [] result = new BigDigit[x.length];
    result[] = x[];
    multibyteIncrementAssign!('-')(result[], lo);
    if (hi)
        multibyteIncrementAssign!('-')(result[1..$], hi);
    if (result[$-1] == 0)
        return result[0..$-1];
    else
        return result;
}

/**  General unsigned multiply routine for bigints.
 *  Sets result = x * y.
 *
 *  The length of y must not be larger than the length of x.
 *  Different algorithms are used, depending on the lengths of x and y.
 *  TODO: "Modern Computer Arithmetic" suggests the OddEvenKaratsuba algorithm for the
 *  unbalanced case. (But I doubt it would be faster in practice).
 *
 */
void mulInternal(BigDigit[] result, const(BigDigit)[] x, const(BigDigit)[] y)
    pure
{
    assert( result.length == x.length + y.length );
    assert( y.length > 0 );
    assert( x.length >= y.length);
    if (y.length <= KARATSUBALIMIT)
    {
        // Small multiplier, we'll just use the asm classic multiply.
        if (y.length == 1)
        {   // Trivial case, no cache effects to worry about
            result[x.length] = multibyteMul(result[0..x.length], x, y[0], 0);
            return;
        }

        if (x.length + y.length < CACHELIMIT)
            return mulSimple(result, x, y);

        // If x is so big that it won't fit into the cache, we divide it into chunks
        // Every chunk must be greater than y.length.
        // We make the first chunk shorter, if necessary, to ensure this.

        auto chunksize = CACHELIMIT / y.length;
        auto residual  =  x.length % chunksize;
        if (residual < y.length)
        {
            chunksize -= y.length;
        }

        // Use schoolbook multiply.
        mulSimple(result[0 .. chunksize + y.length], x[0..chunksize], y);
        auto done = chunksize;

        while (done < x.length)
        {
            // result[done .. done+ylength] already has a value.
            chunksize = (done + (CACHELIMIT / y.length) < x.length) ? (CACHELIMIT / y.length) :  x.length - done;
            BigDigit [KARATSUBALIMIT] partial;
            partial[0..y.length] = result[done..done+y.length];
            mulSimple(result[done..done+chunksize+y.length], x[done..done+chunksize], y);
            addAssignSimple(result[done..done+chunksize + y.length], partial[0..y.length]);
            done += chunksize;
        }
        return;
    }

    auto half = (x.length >> 1) + (x.length & 1);
    if (2*y.length*y.length <= x.length*x.length)
    {
        // UNBALANCED MULTIPLY
        // Use school multiply to cut into quasi-squares of Karatsuba-size
        // or larger. The ratio of the two sides of the 'square' must be
        // between 1.414:1 and 1:1. Use Karatsuba on each chunk.
        //
        // For maximum performance, we want the ratio to be as close to
        // 1:1 as possible. To achieve this, we can either pad x or y.
        // The best choice depends on the modulus x%y.
        auto numchunks = x.length / y.length;
        auto chunksize = y.length;
        auto extra =  x.length % y.length;
        auto maxchunk = chunksize + extra;
        bool paddingY; // true = we're padding Y, false = we're padding X.
        if (extra * extra * 2 < y.length*y.length)
        {
            // The leftover bit is small enough that it should be incorporated
            // in the existing chunks.
            // Make all the chunks a tiny bit bigger
            // (We're padding y with zeros)
            chunksize += extra / cast(double)numchunks;
            extra = x.length - chunksize*numchunks;
            // there will probably be a few left over.
            // Every chunk will either have size chunksize, or chunksize+1.
            maxchunk = chunksize + 1;
            paddingY = true;
            assert(chunksize + extra + chunksize *(numchunks-1) == x.length );
        }
        else
        {
            // the extra bit is large enough that it's worth making a new chunk.
            // (This means we're padding x with zeros, when doing the first one).
            maxchunk = chunksize;
            ++numchunks;
            paddingY = false;
            assert(extra + chunksize *(numchunks-1) == x.length );
        }
        // We make the buffer a bit bigger so we have space for the partial sums.
        BigDigit [] scratchbuff = new BigDigit[karatsubaRequiredBuffSize(maxchunk) + y.length];
        BigDigit [] partial = scratchbuff[$ - y.length .. $];
        size_t done; // how much of X have we done so far?
        double residual = 0;
        if (paddingY)
        {
            // If the first chunk is bigger, do it first. We're padding y.
            mulKaratsuba(result[0 .. y.length + chunksize + (extra > 0 ? 1 : 0 )],
                x[0 .. chunksize + (extra>0?1:0)], y, scratchbuff);
            done = chunksize + (extra > 0 ? 1 : 0);
            if (extra) --extra;
        }
        else
        {   // We're padding X. Begin with the extra bit.
            mulKaratsuba(result[0 .. y.length + extra], y, x[0..extra], scratchbuff);
            done = extra;
            extra = 0;
        }
        auto basechunksize = chunksize;
        while (done < x.length)
        {
            chunksize = basechunksize + (extra > 0 ? 1 : 0);
            if (extra) --extra;
            partial[] = result[done .. done+y.length];
            mulKaratsuba(result[done .. done + y.length + chunksize],
                       x[done .. done+chunksize], y, scratchbuff);
            addAssignSimple(result[done .. done + y.length + chunksize], partial);
            done += chunksize;
        }
        delete scratchbuff;
    }
    else
    {
        // Balanced. Use Karatsuba directly.
        BigDigit [] scratchbuff = new BigDigit[karatsubaRequiredBuffSize(x.length)];
        mulKaratsuba(result, x, y, scratchbuff);
        delete scratchbuff;
    }
}

/**  General unsigned squaring routine for BigInts.
 *   Sets result = x*x.
 *   NOTE: If the highest half-digit of x is zero, the highest digit of result will
 *   also be zero.
 */
void squareInternal(BigDigit[] result, const BigDigit[] x) pure
{
  // Squaring is potentially half a multiply, plus add the squares of
  // the diagonal elements.
  assert(result.length == 2*x.length);
  if (x.length <= KARATSUBASQUARELIMIT)
  {
      if (x.length==1) {
         result[1] = multibyteMul(result[0..1], x, x[0], 0);
         return;
      }
      return squareSimple(result, x);
  }
  // The nice thing about squaring is that it always stays balanced
  BigDigit [] scratchbuff = new BigDigit[karatsubaRequiredBuffSize(x.length)];
  squareKaratsuba(result, x, scratchbuff);
  delete scratchbuff;
}


import core.bitop : bsr;

/// if remainder is null, only calculate quotient.
void divModInternal(BigDigit [] quotient, BigDigit[] remainder, const BigDigit [] u,
        const BigDigit [] v) pure
{
    assert(quotient.length == u.length - v.length + 1);
    assert(remainder == null || remainder.length == v.length);
    assert(v.length > 1);
    assert(u.length >= v.length);

    // Normalize by shifting v left just enough so that
    // its high-order bit is on, and shift u left the
    // same amount. The highest bit of u will never be set.

    BigDigit [] vn = new BigDigit[v.length];
    BigDigit [] un = new BigDigit[u.length + 1];
    // How much to left shift v, so that its MSB is set.
    uint s = BIGDIGITSHIFTMASK - bsr(v[$-1]);
    if (s!=0)
    {
        multibyteShl(vn, v, s);
        un[$-1] = multibyteShl(un[0..$-1], u, s);
    }
    else
    {
        vn[] = v[];
        un[0..$-1] = u[];
        un[$-1] = 0;
    }
    if (quotient.length<FASTDIVLIMIT)
    {
        schoolbookDivMod(quotient, un, vn);
    }
    else
    {
        blockDivMod(quotient, un, vn);
    }

    // Unnormalize remainder, if required.
    if (remainder != null)
    {
        if (s == 0) remainder[] = un[0..vn.length];
        else multibyteShr(remainder, un[0..vn.length+1], s);
    }
    delete un;
    delete vn;
}

unittest
{
    immutable(uint) [] u = [0, 0xFFFF_FFFE, 0x8000_0000];
    immutable(uint) [] v = [0xFFFF_FFFF, 0x8000_0000];
    uint [] q = new uint[u.length - v.length + 1];
    uint [] r = new uint[2];
    divModInternal(q, r, u, v);
    assert(q[]==[0xFFFF_FFFFu, 0]);
    assert(r[]==[0xFFFF_FFFFu, 0x7FFF_FFFF]);
    u = [0, 0xFFFF_FFFE, 0x8000_0001];
    v = [0xFFFF_FFFF, 0x8000_0000];
    divModInternal(q, r, u, v);
}


private:
// Converts a big uint to a hexadecimal string.
//
// Optionally, a separator character (eg, an underscore) may be added between
// every 8 digits.
// buff.length must be data.length*8 if separator is zero,
// or data.length*9 if separator is non-zero. It will be completely filled.
char [] biguintToHex(char [] buff, const BigDigit [] data, char separator=0)
    pure
{
    int x=0;
    for (ptrdiff_t i=data.length - 1; i>=0; --i)
    {
        toHexZeroPadded(buff[x..x+8], data[i]);
        x+=8;
        if (separator)
        {
            if (i>0) buff[x] = separator;
            ++x;
        }
    }
    return buff;
}

/** Convert a big uint into a decimal string.
 *
 * Params:
 *  data    The biguint to be converted. Will be destroyed.
 *  buff    The destination buffer for the decimal string. Must be
 *          large enough to store the result, including leading zeros.
 *          Will be filled backwards, starting from buff[$-1].
 *
 * buff.length must be >= (data.length*32)/log2(10) = 9.63296 * data.length.
 * Returns:
 *    the lowest index of buff which was used.
 */
size_t biguintToDecimal(char [] buff, BigDigit [] data) pure
{
    ptrdiff_t sofar = buff.length;
    // Might be better to divide by (10^38/2^32) since that gives 38 digits for
    // the price of 3 divisions and a shr; this version only gives 27 digits
    // for 3 divisions.
    while(data.length>1)
    {
        uint rem = multibyteDivAssign(data, 10_0000_0000, 0);
        itoaZeroPadded(buff[sofar-9 .. sofar], rem);
        sofar -= 9;
        if (data[$-1] == 0 && data.length > 1)
        {
            data.length = data.length - 1;
        }
    }
    itoaZeroPadded(buff[sofar-10 .. sofar], data[0]);
    sofar -= 10;
    // and strip off the leading zeros
    while(sofar!= buff.length-1 && buff[sofar] == '0')
        sofar++;
    return sofar;
}

/** Convert a decimal string into a big uint.
 *
 * Params:
 *  data    The biguint to be receive the result. Must be large enough to
 *          store the result.
 *  s       The decimal string. May contain _ or 0..9
 *
 * The required length for the destination buffer is slightly less than
 *  1 + s.length/log2(10) = 1 + s.length/3.3219.
 *
 * Returns:
 *    the highest index of data which was used.
 */
int biguintFromDecimal(BigDigit [] data, const(char)[] s) pure
in
{
    assert((data.length >= 2) || (data.length == 1 && s.length == 1));
}
body
{
    import std.conv : ConvException;

    // Convert to base 1e19 = 10_000_000_000_000_000_000.
    // (this is the largest power of 10 that will fit into a long).
    // The length will be less than 1 + s.length/log2(10) = 1 + s.length/3.3219.
    // 485 bits will only just fit into 146 decimal digits.
    // As we convert the string, we record the number of digits we've seen in base 19:
    // hi is the number of digits/19, lo is the extra digits (0 to 18).
    // TODO: This is inefficient for very large strings (it is O(n^^2)).
    // We should take advantage of fast multiplication once the numbers exceed
    // Karatsuba size.
    uint lo = 0; // number of powers of digits, 0..18
    uint x = 0;
    ulong y = 0;
    uint hi = 0; // number of base 1e19 digits
    data[0] = 0; // initially number is 0.
    if (data.length > 1)
        data[1] = 0;

    for (int i= (s[0]=='-' || s[0]=='+')? 1 : 0; i<s.length; ++i)
    {
        if (s[i] == '_')
            continue;
        if (s[i] < '0' || s[i] > '9')
            throw new ConvException("invalid digit");
        x *= 10;
        x += s[i] - '0';
        ++lo;
        if (lo == 9)
        {
            y = x;
            x = 0;
        }
        if (lo == 18)
        {
            y *= 10_0000_0000;
            y += x;
            x = 0;
        }
        if (lo == 19)
        {
            y *= 10;
            y += x;
            x = 0;
            // Multiply existing number by 10^19, then add y1.
            if (hi>0)
            {
                data[hi] = multibyteMul(data[0..hi], data[0..hi], 1220703125*2u, 0); // 5^13*2 = 0x9184_E72A
                ++hi;
                data[hi] = multibyteMul(data[0..hi], data[0..hi], 15625*262144u, 0); // 5^6*2^18 = 0xF424_0000
                ++hi;
            }
            else
                hi = 2;
            uint c = multibyteIncrementAssign!('+')(data[0..hi], cast(uint)(y&0xFFFF_FFFF));
            c += multibyteIncrementAssign!('+')(data[1..hi], cast(uint)(y>>32));
            if (c!=0)
            {
                data[hi]=c;
                ++hi;
            }
            y = 0;
            lo = 0;
        }
    }
    // Now set y = all remaining digits.
    if (lo>=18)
    {
    }
    else if (lo>=9)
    {
        for (int k=9; k<lo; ++k) y*=10;
        y+=x;
    }
    else
    {
        for (int k=0; k<lo; ++k) y*=10;
        y+=x;
    }
    if (lo != 0)
    {
        if (hi == 0)
        {
            data[0] = cast(uint)y;
            if (data.length == 1)
            {
                hi = 1;
            }
            else
            {
                data[1] = cast(uint)(y >>> 32);
                hi=2;
            }
        }
        else
        {
            while (lo>0)
            {
                uint c = multibyteMul(data[0..hi], data[0..hi], 10, 0);
                if (c!=0)
                {
                    data[hi]=c;
                    ++hi;
                }
                --lo;
            }
            uint c = multibyteIncrementAssign!('+')(data[0..hi], cast(uint)(y&0xFFFF_FFFF));
            if (y > 0xFFFF_FFFFL)
            {
                c += multibyteIncrementAssign!('+')(data[1..hi], cast(uint)(y>>32));
            }
            if (c!=0)
            {
                data[hi]=c;
                ++hi;
            }
        }
    }
    while (hi>1 && data[hi-1]==0)
        --hi;
    return hi;
}


private:
// ------------------------
// These in-place functions are only for internal use; they are incompatible
// with COW.

// Classic 'schoolbook' multiplication.
void mulSimple(BigDigit[] result, const(BigDigit) [] left,
        const(BigDigit)[] right) pure
in
{
    assert(result.length == left.length + right.length);
    assert(right.length>1);
}
body
{
    result[left.length] = multibyteMul(result[0..left.length], left, right[0], 0);
    multibyteMultiplyAccumulate(result[1..$], left, right[1..$]);
}

// Classic 'schoolbook' squaring
void squareSimple(BigDigit[] result, const(BigDigit) [] x) pure
in
{
    assert(result.length == 2*x.length);
    assert(x.length>1);
}
body
{
    multibyteSquare(result, x);
}


// add two uints of possibly different lengths. Result must be as long
// as the larger length.
// Returns carry (0 or 1).
uint addSimple(BigDigit[] result, const BigDigit [] left, const BigDigit [] right) pure
in
{
    assert(result.length == left.length);
    assert(left.length >= right.length);
    assert(right.length>0);
}
body
{
    uint carry = multibyteAdd(result[0..right.length],
            left[0..right.length], right, 0);
    if (right.length < left.length)
    {
        result[right.length..left.length] = left[right.length .. $];
        carry = multibyteIncrementAssign!('+')(result[right.length..$], carry);
    }
    return carry;
}

//  result = left - right
// returns carry (0 or 1)
BigDigit subSimple(BigDigit [] result,const(BigDigit) [] left,
        const(BigDigit) [] right) pure
in
{
    assert(result.length == left.length);
    assert(left.length >= right.length);
    assert(right.length>0);
}
body
{
    BigDigit carry = multibyteSub(result[0..right.length],
            left[0..right.length], right, 0);
    if (right.length < left.length)
    {
        result[right.length..left.length] = left[right.length .. $];
        carry = multibyteIncrementAssign!('-')(result[right.length..$], carry);
    } //else if (result.length==left.length+1) { result[$-1] = carry; carry=0; }
    return carry;
}


/* result = result - right
 * Returns carry = 1 if result was less than right.
*/
BigDigit subAssignSimple(BigDigit [] result, const(BigDigit) [] right) pure
{
    assert(result.length >= right.length);
    uint c = multibyteSub(result[0..right.length], result[0..right.length], right, 0);
    if (c && result.length > right.length)
        c = multibyteIncrementAssign!('-')(result[right.length .. $], c);
    return c;
}

/* result = result + right
*/
BigDigit addAssignSimple(BigDigit [] result, const(BigDigit) [] right) pure
{
    assert(result.length >= right.length);
    uint c = multibyteAdd(result[0..right.length], result[0..right.length], right, 0);
    if (c && result.length > right.length)
       c = multibyteIncrementAssign!('+')(result[right.length .. $], c);
    return c;
}

/* performs result += wantSub? - right : right;
*/
BigDigit addOrSubAssignSimple(BigDigit [] result, const(BigDigit) [] right,
        bool wantSub) pure
{
    if (wantSub)
        return subAssignSimple(result, right);
    else
        return addAssignSimple(result, right);
}


// return true if x<y, considering leading zeros
bool less(const(BigDigit)[] x, const(BigDigit)[] y) pure
{
    assert(x.length >= y.length);
    auto k = x.length-1;
    while(x[k]==0 && k>=y.length)
        --k;
    if (k>=y.length)
        return false;
    while (k>0 && x[k]==y[k])
        --k;
    return x[k] < y[k];
}

// Set result = abs(x-y), return true if result is negative(x<y), false if x<=y.
bool inplaceSub(BigDigit[] result, const(BigDigit)[] x, const(BigDigit)[] y)
    pure
{
    assert(result.length == (x.length >= y.length) ? x.length : y.length);

    size_t minlen;
    bool negative;
    if (x.length >= y.length)
    {
        minlen = y.length;
        negative = less(x, y);
    }
    else
    {
       minlen = x.length;
       negative = !less(y, x);
    }
    const (BigDigit)[] large, small;
    if (negative)
    {
        large = y; small = x;
    }
    else
    {
        large = x; small = y;
    }

    BigDigit carry = multibyteSub(result[0..minlen], large[0..minlen], small[0..minlen], 0);
    if (x.length != y.length)
    {
        result[minlen..large.length]= large[minlen..$];
        result[large.length..$] = 0;
        if (carry)
            multibyteIncrementAssign!('-')(result[minlen..$], carry);
    }
    return negative;
}

/* Determine how much space is required for the temporaries
 * when performing a Karatsuba multiplication.
 */
size_t karatsubaRequiredBuffSize(size_t xlen) pure
{
    return xlen <= KARATSUBALIMIT ? 0 : 2*xlen; // - KARATSUBALIMIT+2;
}

/* Sets result = x*y, using Karatsuba multiplication.
* x must be longer or equal to y.
* Valid only for balanced multiplies, where x is not shorter than y.
* It is superior to schoolbook multiplication if and only if
*    sqrt(2)*y.length > x.length > y.length.
* Karatsuba multiplication is O(n^1.59), whereas schoolbook is O(n^2)
* The maximum allowable length of x and y is uint.max; but better algorithms
* should be used far before that length is reached.
* Params:
* scratchbuff      An array long enough to store all the temporaries. Will be destroyed.
*/
void mulKaratsuba(BigDigit [] result, const(BigDigit) [] x,
        const(BigDigit)[] y, BigDigit [] scratchbuff) pure
{
    assert(x.length >= y.length);
          assert(result.length < uint.max, "Operands too large");
    assert(result.length == x.length + y.length);
    if (x.length <= KARATSUBALIMIT)
    {
        return mulSimple(result, x, y);
    }
    // Must be almost square (otherwise, a schoolbook iteration is better)
    assert(2L * y.length * y.length > (x.length-1) * (x.length-1),
        "Bigint Internal Error: Asymmetric Karatsuba");

    // The subtractive version of Karatsuba multiply uses the following result:
    // (Nx1 + x0)*(Ny1 + y0) = (N*N)*x1y1 + x0y0 + N * (x0y0 + x1y1 - mid)
    // where mid = (x0-x1)*(y0-y1)
    // requiring 3 multiplies of length N, instead of 4.
    // The advantage of the subtractive over the additive version is that
    // the mid multiply cannot exceed length N. But there are subtleties:
    // (x0-x1),(y0-y1) may be negative or zero. To keep it simple, we
    // retain all of the leading zeros in the subtractions

    // half length, round up.
    auto half = (x.length >> 1) + (x.length & 1);

    const(BigDigit) [] x0 = x[0 .. half];
    const(BigDigit) [] x1 = x[half .. $];
    const(BigDigit) [] y0 = y[0 .. half];
    const(BigDigit) [] y1 = y[half .. $];
    BigDigit [] mid = scratchbuff[0 .. half*2];
    BigDigit [] newscratchbuff = scratchbuff[half*2 .. $];
    BigDigit [] resultLow = result[0 .. 2*half];
    BigDigit [] resultHigh = result[2*half .. $];
     // initially use result to store temporaries
    BigDigit [] xdiff= result[0 .. half];
    BigDigit [] ydiff = result[half .. half*2];

    // First, we calculate mid, and sign of mid
    bool midNegative = inplaceSub(xdiff, x0, x1)
                      ^ inplaceSub(ydiff, y0, y1);
    mulKaratsuba(mid, xdiff, ydiff, newscratchbuff);

    // Low half of result gets x0 * y0. High half gets x1 * y1

    mulKaratsuba(resultLow, x0, y0, newscratchbuff);

    if (2L * y1.length * y1.length < x1.length * x1.length)
    {
        // an asymmetric situation has been created.
        // Worst case is if x:y = 1.414 : 1, then x1:y1 = 2.41 : 1.
        // Applying one schoolbook multiply gives us two pieces each 1.2:1
        if (y1.length <= KARATSUBALIMIT)
            mulSimple(resultHigh, x1, y1);
        else
        {
            // divide x1 in two, then use schoolbook multiply on the two pieces.
            auto quarter = (x1.length >> 1) + (x1.length & 1);
            bool ysmaller = (quarter >= y1.length);
            mulKaratsuba(resultHigh[0..quarter+y1.length], ysmaller ? x1[0..quarter] : y1,
                ysmaller ? y1 : x1[0..quarter], newscratchbuff);
            // Save the part which will be overwritten.
            bool ysmaller2 = ((x1.length - quarter) >= y1.length);
            newscratchbuff[0..y1.length] = resultHigh[quarter..quarter + y1.length];
            mulKaratsuba(resultHigh[quarter..$], ysmaller2 ? x1[quarter..$] : y1,
                ysmaller2 ? y1 : x1[quarter..$], newscratchbuff[y1.length..$]);

            resultHigh[quarter..$].addAssignSimple(newscratchbuff[0..y1.length]);
        }
    }
    else
        mulKaratsuba(resultHigh, x1, y1, newscratchbuff);

    /* We now have result = x0y0 + (N*N)*x1y1
       Before adding or subtracting mid, we must calculate
       result += N * (x0y0 + x1y1)
       We can do this with three half-length additions. With a = x0y0, b = x1y1:
                      aHI aLO
        +       aHI   aLO
        +       bHI   bLO
        +  bHI  bLO
        =  R3   R2    R1   R0
        R1 = aHI + bLO + aLO
        R2 = aHI + bLO + aHI + carry_from_R1
        R3 = bHi + carry_from_R2

     It might actually be quicker to do it in two full-length additions:
     newscratchbuff[2*half] = addSimple(newscratchbuff[0..2*half], result[0..2*half], result[2*half..$]);
     addAssignSimple(result[half..$], newscratchbuff[0..2*half+1]);
   */
    BigDigit[] R1 = result[half..half*2];
    BigDigit[] R2 = result[half*2..half*3];
    BigDigit[] R3 = result[half*3..$];
    BigDigit c1 = multibyteAdd(R2, R2, R1, 0); // c1:R2 = R2 + R1
    BigDigit c2 = multibyteAdd(R1, R2, result[0..half], 0); // c2:R1 = R2 + R1 + R0
    BigDigit c3 = addAssignSimple(R2, R3); // R2 = R2 + R1 + R3
    if (c1+c2)
        multibyteIncrementAssign!('+')(result[half*2..$], c1+c2);
    if (c1+c3)
        multibyteIncrementAssign!('+')(R3, c1+c3);

    // And finally we subtract mid
    addOrSubAssignSimple(result[half..$], mid, !midNegative);
}

void squareKaratsuba(BigDigit [] result, const BigDigit [] x,
        BigDigit [] scratchbuff) pure
{
    // See mulKaratsuba for implementation comments.
    // Squaring is simpler, since it never gets asymmetric.
    assert(result.length < uint.max, "Operands too large");
    assert(result.length == 2*x.length);
    if (x.length <= KARATSUBASQUARELIMIT)
    {
        return squareSimple(result, x);
    }
    // half length, round up.
    auto half = (x.length >> 1) + (x.length & 1);

    const(BigDigit)[] x0 = x[0 .. half];
    const(BigDigit)[] x1 = x[half .. $];
    BigDigit [] mid = scratchbuff[0 .. half*2];
    BigDigit [] newscratchbuff = scratchbuff[half*2 .. $];
     // initially use result to store temporaries
    BigDigit [] xdiff= result[0 .. half];
    BigDigit [] ydiff = result[half .. half*2];

    // First, we calculate mid. We don't need its sign
    inplaceSub(xdiff, x0, x1);
    squareKaratsuba(mid, xdiff, newscratchbuff);

    // Set result = x0x0 + (N*N)*x1x1
    squareKaratsuba(result[0 .. 2*half], x0, newscratchbuff);
    squareKaratsuba(result[2*half .. $], x1, newscratchbuff);

    /* result += N * (x0x0 + x1x1)
       Do this with three half-length additions. With a = x0x0, b = x1x1:
        R1 = aHI + bLO + aLO
        R2 = aHI + bLO + aHI + carry_from_R1
        R3 = bHi + carry_from_R2
    */
    BigDigit[] R1 = result[half..half*2];
    BigDigit[] R2 = result[half*2..half*3];
    BigDigit[] R3 = result[half*3..$];
    BigDigit c1 = multibyteAdd(R2, R2, R1, 0); // c1:R2 = R2 + R1
    BigDigit c2 = multibyteAdd(R1, R2, result[0..half], 0); // c2:R1 = R2 + R1 + R0
    BigDigit c3 = addAssignSimple(R2, R3); // R2 = R2 + R1 + R3
    if (c1+c2) multibyteIncrementAssign!('+')(result[half*2..$], c1+c2);
    if (c1+c3) multibyteIncrementAssign!('+')(R3, c1+c3);

    // And finally we subtract mid, which is always positive
    subAssignSimple(result[half..$], mid);
}

/* Knuth's Algorithm D, as presented in
 * H.S. Warren, "Hacker's Delight", Addison-Wesley Professional (2002).
 * Also described in "Modern Computer Arithmetic" 0.2, Exercise 1.8.18.
 * Given u and v, calculates  quotient  = u / v, u = u % v.
 * v must be normalized (ie, the MSB of v must be 1).
 * The most significant words of quotient and u may be zero.
 * u[0..v.length] holds the remainder.
 */
void schoolbookDivMod(BigDigit [] quotient, BigDigit [] u, in BigDigit [] v)
    pure
{
    assert(quotient.length == u.length - v.length);
    assert(v.length > 1);
    assert(u.length >= v.length);
    assert((v[$-1]&0x8000_0000)!=0);
    assert(u[$-1] < v[$-1]);
    // BUG: This code only works if BigDigit is uint.
    uint vhi = v[$-1];
    uint vlo = v[$-2];

    for (ptrdiff_t j = u.length - v.length - 1; j >= 0; j--)
    {
        // Compute estimate of quotient[j],
        // qhat = (three most significant words of u)/(two most sig words of v).
        uint qhat;
        if (u[j + v.length] == vhi)
        {
            // uu/vhi could exceed uint.max (it will be 0x8000_0000 or 0x8000_0001)
            qhat = uint.max;
        }
        else
        {
            uint ulo = u[j + v.length - 2];
            version(D_InlineAsm_X86)
            {
                // Note: On DMD, this is only ~10% faster than the non-asm code.
                uint *p = &u[j + v.length - 1];
                asm
                {
                    mov EAX, p;
                    mov EDX, [EAX+4];
                    mov EAX, [EAX];
                    div dword ptr [vhi];
                    mov qhat, EAX;
                    mov ECX, EDX;
div3by2correction:
                    mul dword ptr [vlo]; // EDX:EAX = qhat * vlo
                    sub EAX, ulo;
                    sbb EDX, ECX;
                    jbe div3by2done;
                    mov EAX, qhat;
                    dec EAX;
                    mov qhat, EAX;
                    add ECX, dword ptr [vhi];
                    jnc div3by2correction;
div3by2done:    ;
                }
            }
            else
            { // version(InlineAsm)
                ulong uu = (cast(ulong)(u[j + v.length]) << 32) | u[j + v.length - 1];
                ulong bigqhat = uu / vhi;
                ulong rhat =  uu - bigqhat * vhi;
                qhat = cast(uint)bigqhat;
again:
                if (cast(ulong)qhat * vlo > ((rhat << 32) + ulo))
                {
                    --qhat;
                    rhat += vhi;
                    if (!(rhat & 0xFFFF_FFFF_0000_0000L))
                        goto again;
                }
            } // version(InlineAsm)
        }
        // Multiply and subtract.
        uint carry = multibyteMulAdd!('-')(u[j..j + v.length], v, qhat, 0);

        if (u[j+v.length] < carry)
        {
            // If we subtracted too much, add back
            --qhat;
            carry -= multibyteAdd(u[j..j + v.length],u[j..j + v.length], v, 0);
        }
        quotient[j] = qhat;
        u[j + v.length] = u[j + v.length] - carry;
    }
}

private:

// TODO: Replace with a library call
void itoaZeroPadded(char[] output, uint value, int radix = 10) pure
{
    ptrdiff_t x = output.length - 1;
    for( ; x >= 0; --x)
    {
        output[x]= cast(char)(value % radix + '0');
        value /= radix;
    }
}

void toHexZeroPadded(char[] output, uint value) pure
{
    ptrdiff_t x = output.length - 1;
    static immutable string hexDigits = "0123456789ABCDEF";
    for( ; x>=0; --x)
    {
        output[x] = hexDigits[value & 0xF];
        value >>= 4;
    }
}

private:

// Returns the highest value of i for which left[i]!=right[i],
// or 0 if left[] == right[]
size_t highestDifferentDigit(const BigDigit [] left, const BigDigit [] right) pure
{
    assert(left.length == right.length);
    for (ptrdiff_t i = left.length - 1; i>0; --i)
    {
        if (left[i] != right[i])
            return i;
    }
    return 0;
}

// Returns the lowest value of i for which x[i]!=0.
int firstNonZeroDigit(const BigDigit [] x) pure
{
    int k = 0;
    while (x[k]==0)
    {
        ++k;
        assert(k<x.length);
    }
    return k;
}
import core.stdc.stdio;
/*
    Calculate quotient and remainder of u / v using fast recursive division.
    v must be normalised, and must be at least half as long as u.
    Given u and v, v normalised, calculates  quotient  = u/v, u = u%v.
    scratch is temporary storage space, length must be >= quotient + 1.

Returns:
    u[0..v.length] is the remainder. u[v.length..$] is corrupted.

    Implements algorithm 1.8 from MCA.
    This algorithm has an annoying special case. After the first recursion, the
    highest bit of the quotient may be set. This means that in the second
    recursive call, the 'in' contract would be violated. (This happens only
    when the top quarter of u is equal to the top half of v. A base 10
    equivalent example of this situation is 5517/56; the first step gives
    55/5 = 11). To maintain the in contract, we pad a zero to the top of both
    u and the quotient. 'mayOverflow' indicates that that the special case
    has occurred.
    (In MCA, a different strategy is used: the in contract is weakened, and
    schoolbookDivMod is more general: it allows the high bit of u to be set).
    See also:
    - C. Burkinel and J. Ziegler, "Fast Recursive Division", MPI-I-98-1-022,
      Max-Planck Institute fuer Informatik, (Oct 1998).
*/
void recursiveDivMod(BigDigit[] quotient, BigDigit[] u, const(BigDigit)[] v,
                     BigDigit[] scratch, bool mayOverflow = false) pure
in
{
    // v must be normalized
    assert(v.length > 1);
    assert((v[$ - 1] & 0x8000_0000) != 0);
    assert(!(u[$ - 1] & 0x8000_0000));
    assert(quotient.length == u.length - v.length);
    if (mayOverflow)
    {
        assert(u[$-1] == 0);
        assert(u[$-2] & 0x8000_0000);
    }

    // Must be symmetric. Use block schoolbook division if not.
    assert((mayOverflow ? u.length-1 : u.length) <= 2 * v.length);
    assert((mayOverflow ? u.length-1 : u.length) >= v.length);
    assert(scratch.length >= quotient.length + (mayOverflow ? 0 : 1));
}
body
{
    if (quotient.length < FASTDIVLIMIT)
    {
        return schoolbookDivMod(quotient, u, v);
    }

    // Split quotient into two halves, but keep padding in the top half
    auto k = (mayOverflow ?  quotient.length - 1 : quotient.length) >> 1;

    // RECURSION 1: Calculate the high half of the quotient

    // Note that if u and quotient were padded, they remain padded during
    // this call, so in contract is satisfied.
    recursiveDivMod(quotient[k .. $], u[2 * k .. $], v[k .. $],
        scratch, mayOverflow);

    // quotient[k..$] is our guess at the high quotient.
    // u[2*k.. 2.*k + v.length - k = k + v.length] is the high part of the
    // first remainder. u[0..2*k] is the low part.

    // Calculate the full first remainder to be
    //    remainder - highQuotient * lowDivisor
    // reducing highQuotient until the remainder is positive.
    // The low part of the remainder, u[0..k], cannot be altered by this.

    adjustRemainder(quotient[k .. $], u[k .. k + v.length], v, k,
            scratch[0 .. quotient.length], mayOverflow);

    // RECURSION 2: Calculate the low half of the quotient
    // The full first remainder is now in u[0..k + v.length].

    if (u[k + v.length - 1] & 0x8000_0000)
    {
        // Special case. The high quotient is 0x1_00...000 or 0x1_00...001.
        // This means we need an extra quotient word for the next recursion.
        // We need to restore the invariant for the recursive calls.
        // We do this by padding both u and quotient. Extending u is trivial,
        // because the higher words will not be used again. But for the
        // quotient, we're clobbering the low word of the high quotient,
        // so we need save it, and add it back in after the recursive call.

        auto clobberedQuotient = quotient[k];
        u[k+v.length] = 0;

        recursiveDivMod(quotient[0 .. k+1], u[k .. k + v.length+1],
            v[k .. $], scratch, true);
        adjustRemainder(quotient[0 .. k+1], u[0 .. v.length], v, k,
            scratch[0 .. 2 * k+1], true);

        // Now add the quotient word that got clobbered earlier.
        multibyteIncrementAssign!('+')(quotient[k..$], clobberedQuotient);
    }
    else
    {
        // The special case has NOT happened.
        recursiveDivMod(quotient[0 .. k], u[k .. k + v.length], v[k .. $],
            scratch, false);

        // high remainder is in u[k..k+(v.length-k)] == u[k .. v.length]

        adjustRemainder(quotient[0 .. k], u[0 .. v.length], v, k,
            scratch[0 .. 2 * k]);
    }
}

// rem -= quot * v[0..k].
// If would make rem negative, decrease quot until rem is >=0.
// Needs (quot.length * k) scratch space to store the result of the multiply.
void adjustRemainder(BigDigit[] quot, BigDigit[] rem, const(BigDigit)[] v,
        ptrdiff_t k,
        BigDigit[] scratch, bool mayOverflow = false) pure
{
    assert(rem.length == v.length);
    mulInternal(scratch, quot, v[0 .. k]);
    uint carry = 0;
    if (mayOverflow)
        carry = scratch[$-1] + subAssignSimple(rem, scratch[0..$-1]);
    else
        carry = subAssignSimple(rem, scratch);
    while(carry)
    {
        multibyteIncrementAssign!('-')(quot, 1); // quot--
        carry -= multibyteAdd(rem, rem, v, 0);
    }
}

// Cope with unbalanced division by performing block schoolbook division.
void blockDivMod(BigDigit [] quotient, BigDigit [] u, in BigDigit [] v) pure
{
    assert(quotient.length == u.length - v.length);
    assert(v.length > 1);
    assert(u.length >= v.length);
    assert((v[$-1] & 0x8000_0000)!=0);
    assert((u[$-1] & 0x8000_0000)==0);
    BigDigit [] scratch = new BigDigit[v.length + 1];

    // Perform block schoolbook division, with 'v.length' blocks.
    auto m = u.length - v.length;
    while (m > v.length)
    {
        bool mayOverflow = (u[m + v.length -1 ] & 0x8000_0000)!=0;
        BigDigit saveq;
        if (mayOverflow)
        {
            u[m + v.length] = 0;
            saveq = quotient[m];
        }
        recursiveDivMod(quotient[m-v.length..m + (mayOverflow? 1: 0)],
            u[m - v.length..m + v.length + (mayOverflow? 1: 0)], v, scratch, mayOverflow);
        if (mayOverflow)
        {
            assert(quotient[m] == 0);
            quotient[m] = saveq;
        }
        m -= v.length;
    }
    recursiveDivMod(quotient[0..m], u[0..m + v.length], v, scratch);
    delete scratch;
}

version(unittest)
{
    import std.c.stdio;
}

unittest
{
    void printBiguint(const uint [] data)
    {
        char [] buff = biguintToHex(new char[data.length*9], data, '_');
        printf("%.*s\n", buff.length, buff.ptr);
    }

    void printDecimalBigUint(BigUint data)
    {
        auto str = data.toDecimalString(0);
        printf("%.*s\n", str.length, str.ptr);
    }

    uint [] a, b;
    a = new uint[43];
    b = new uint[179];
    for (int i=0; i<a.length; ++i) a[i] = 0x1234_B6E9 + i;
    for (int i=0; i<b.length; ++i) b[i] = 0x1BCD_8763 - i*546;

    a[$-1] |= 0x8000_0000;
    uint [] r = new uint[a.length];
    uint [] q = new uint[b.length-a.length+1];

    divModInternal(q, r, b, a);
    q = q[0..$-1];
    uint [] r1 = r.dup;
    uint [] q1 = q.dup;
    blockDivMod(q, b, a);
    r = b[0..a.length];
    assert(r[] == r1[]);
    assert(q[] == q1[]);
}