Wolfgang Thallergcc16 merge934c3885 months agoHistory
2958 lines · 97.8 KB · D
RawDownload
/** 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){    static import std.internal.math.biguintx86;}static import std.internal.math.biguintnoasm; import std.internal.math.biguintnoasm : BigDigit, KARATSUBALIMIT,    KARATSUBASQUARELIMIT; alias multibyteAdd = multibyteAddSub!('+');alias multibyteSub = multibyteAddSub!('-'); private import std.traits;private import std.range.primitives;public import std.ascii : LetterCase;import std.range.primitives;import std.traits; private: // dipatchers to the right low-level primitives. Added to allow BigInt CTFE for// 32 bit systems (https://issues.dlang.org/show_bug.cgi?id=14767) although it's// used by the other architectures too.// See comments below in case it has to be refactored.version (X86)uint multibyteAddSub(char op)(uint[] dest, const(uint)[] src1, const (uint)[] src2, uint carry){    // must be checked before, otherwise D_InlineAsm_X86 is true.    if (__ctfe)        return std.internal.math.biguintnoasm.multibyteAddSub!op(dest, src1, src2, carry);    // Runtime.    else version (D_InlineAsm_X86)        return std.internal.math.biguintx86.multibyteAddSub!op(dest, src1, src2, carry);    // Runtime if no asm available.    else        return std.internal.math.biguintnoasm.multibyteAddSub!op(dest, src1, src2, carry);}// Any other architectureelse alias multibyteAddSub = std.internal.math.biguintnoasm.multibyteAddSub; version (X86)uint multibyteIncrementAssign(char op)(uint[] dest, uint carry){    if (__ctfe)        return std.internal.math.biguintnoasm.multibyteIncrementAssign!op(dest, carry);    else version (D_InlineAsm_X86)        return std.internal.math.biguintx86.multibyteIncrementAssign!op(dest, carry);    else        return std.internal.math.biguintnoasm.multibyteIncrementAssign!op(dest, carry);}else alias multibyteIncrementAssign = std.internal.math.biguintnoasm.multibyteIncrementAssign; version (X86)uint multibyteShl()(uint[] dest, const(uint)[] src, uint numbits){    if (__ctfe)        return std.internal.math.biguintnoasm.multibyteShl(dest, src, numbits);    else version (D_InlineAsm_X86)        return std.internal.math.biguintx86.multibyteShl(dest, src, numbits);    else        return std.internal.math.biguintnoasm.multibyteShl(dest, src, numbits);}else alias multibyteShl = std.internal.math.biguintnoasm.multibyteShl; version (X86)void multibyteShr()(uint[] dest, const(uint)[] src, uint numbits){    if (__ctfe)        std.internal.math.biguintnoasm.multibyteShr(dest, src, numbits);    else version (D_InlineAsm_X86)        std.internal.math.biguintx86.multibyteShr(dest, src, numbits);    else        std.internal.math.biguintnoasm.multibyteShr(dest, src, numbits);}else alias multibyteShr = std.internal.math.biguintnoasm.multibyteShr; version (X86)uint multibyteMul()(uint[] dest, const(uint)[] src, uint multiplier, uint carry){    if (__ctfe)        return std.internal.math.biguintnoasm.multibyteMul(dest, src, multiplier, carry);    else version (D_InlineAsm_X86)        return std.internal.math.biguintx86.multibyteMul(dest, src, multiplier, carry);    else        return std.internal.math.biguintnoasm.multibyteMul(dest, src, multiplier, carry);}else alias multibyteMul = std.internal.math.biguintnoasm.multibyteMul; version (X86)uint multibyteMulAdd(char op)(uint[] dest, const(uint)[] src, uint multiplier, uint carry){    if (__ctfe)        return std.internal.math.biguintnoasm.multibyteMulAdd!op(dest, src, multiplier, carry);    else version (D_InlineAsm_X86)        return std.internal.math.biguintx86.multibyteMulAdd!op(dest, src, multiplier, carry);    else        return std.internal.math.biguintnoasm.multibyteMulAdd!op(dest, src, multiplier, carry);}else alias multibyteMulAdd = std.internal.math.biguintnoasm.multibyteMulAdd; version (X86)void multibyteMultiplyAccumulate()(uint[] dest, const(uint)[] left, const(uint)[] right){    if (__ctfe)        std.internal.math.biguintnoasm.multibyteMultiplyAccumulate(dest, left, right);    else version (D_InlineAsm_X86)        std.internal.math.biguintx86.multibyteMultiplyAccumulate(dest, left, right);    else        std.internal.math.biguintnoasm.multibyteMultiplyAccumulate(dest, left, right);}else alias multibyteMultiplyAccumulate = std.internal.math.biguintnoasm.multibyteMultiplyAccumulate; version (X86)uint multibyteDivAssign()(uint[] dest, uint divisor, uint overflow){    if (__ctfe)        return std.internal.math.biguintnoasm.multibyteDivAssign(dest, divisor, overflow);    else version (D_InlineAsm_X86)        return std.internal.math.biguintx86.multibyteDivAssign(dest, divisor, overflow);    else        return std.internal.math.biguintnoasm.multibyteDivAssign(dest, divisor, overflow);}else alias multibyteDivAssign = std.internal.math.biguintnoasm.multibyteDivAssign; version (X86)void multibyteAddDiagonalSquares()(uint[] dest, const(uint)[] src){    if (__ctfe)        std.internal.math.biguintnoasm.multibyteAddDiagonalSquares(dest, src);    else version (D_InlineAsm_X86)        std.internal.math.biguintx86.multibyteAddDiagonalSquares(dest, src);    else        std.internal.math.biguintnoasm.multibyteAddDiagonalSquares(dest, src);}else alias multibyteAddDiagonalSquares = std.internal.math.biguintnoasm.multibyteAddDiagonalSquares; version (X86)void multibyteTriangleAccumulate()(uint[] dest, const(uint)[] x){    if (__ctfe)        std.internal.math.biguintnoasm.multibyteTriangleAccumulate(dest, x);    else version (D_InlineAsm_X86)        std.internal.math.biguintx86.multibyteTriangleAccumulate(dest, x);    else        std.internal.math.biguintnoasm.multibyteTriangleAccumulate(dest, x);}else alias multibyteTriangleAccumulate = std.internal.math.biguintnoasm.multibyteTriangleAccumulate; version (X86)void multibyteSquare()(BigDigit[] result, const(BigDigit)[] x){    if (__ctfe)        std.internal.math.biguintnoasm.multibyteSquare(result, x);    else version (D_InlineAsm_X86)        std.internal.math.biguintx86.multibyteSquare(result, x);    else        std.internal.math.biguintnoasm.multibyteSquare(result, x);}else alias multibyteSquare = std.internal.math.biguintnoasm.multibyteSquare; // Limits for when to switch between algorithms.// Half the size of the data cache.@nogc nothrow pure @safe size_t getCacheLimit(){    import core.cpuid : dataCaches;    return dataCaches[0].size * 1024 / 2;}enum size_t FASTDIVLIMIT = 100; // crossover to recursive division  // These constants are used by shift operationsstatic 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"); import std.exception : assumeUnique;import std.traits : isIntegral;enum BigDigitBits = BigDigit.sizeof*8;template maxBigDigits(T)if (isIntegral!T){    enum maxBigDigits = (T.sizeof+BigDigit.sizeof-1)/BigDigit.sizeof;} static immutable BigDigit[] ZERO = [0];static immutable BigDigit[] ONE = [1];static immutable BigDigit[] TWO = [2];static 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 ),                "Invariant requires data to not empty or zero");    }     immutable(BigDigit) [] data = ZERO;     this(return scope immutable(BigDigit) [] x) pure nothrow @nogc @safe    {       data = x;    }  package(std)  // used from: std.bigint    this(T)(T x) pure nothrow @safe scope    if (isIntegral!T)    {        opAssign(x);    }     enum trustedAssumeUnique = function(BigDigit[] input) pure @trusted @nogc {        return assumeUnique(input);    };public:    // Length in uints    @property size_t uintLength() pure nothrow const @safe @nogc scope    {        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);        }    }    @property size_t ulongLength() pure nothrow const @safe @nogc scope    {        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(size_t n) pure nothrow const @safe @nogc scope    {        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(size_t n) pure nothrow const @safe @nogc scope    {        static if (BigDigit.sizeof == int.sizeof)        {            return data[n];        }        else        {            immutable x = data[n >> 1];            return (n & 1) ? cast(uint)(x >> 32) : cast(uint) x;        }    }     ///    void opAssign(Tulong)(Tulong u) pure nothrow @safe scope    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 nothrow @nogc @safe scope    {        this.data = y.data;    }     ///    int opCmp(Tdummy = void)(const BigUint y) pure nothrow @nogc const @safe scope    {        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 nothrow @nogc const @safe scope    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 nothrow @nogc const @safe scope    {           return y.data[] == data[];    }     bool opEquals(Tdummy = void)(ulong y) pure nothrow @nogc const @safe scope    {        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 @nogc scope    {        return data.length == 1 && data[0] == 0;    }     size_t numBytes() pure nothrow const @safe @nogc scope    {        return data.length * BigDigit.sizeof;    }     // the extra bytes are added to the start of the string    char [] toDecimalString(int frontExtraBytes) const pure nothrow @safe scope    {        immutable 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',            LetterCase letterCase = LetterCase.upper) const pure nothrow @safe scope    {        // 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;        immutable totalSeparatorBytes = separator ? ((extraPad + lenBytes + 7) / 8) - 1: 0;         char [] buff = new char[lenBytes + extraPad + totalSeparatorBytes + frontExtraBytes];        biguintToHex(buff[$ - lenBytes - mainSeparatorBytes .. $], data, separator, letterCase);        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..$];    }     /**     * Convert to an octal string.     */    char[] toOctalString() pure nothrow @safe const scope    {        auto predictLength = 1 + data.length*BigDigitBits / 3;

Lines 1–500 of 2958.

Keyboard shortcuts

General

Show keyboard shortcuts
?
Focus the current search
/

Issues and pull requests

Open filters
f

Issues, pull requests and Inbox

Select next item
j↓
Select previous item
k↑
Open selected item
Enter
Dismiss filters or clear selection
Esc

Repository code

Go to file
t