diff --git a/source/tanya/math/mp.d b/source/tanya/math/mp.d index 024a0e4..e8dcce3 100644 --- a/source/tanya/math/mp.d +++ b/source/tanya/math/mp.d @@ -12,11 +12,11 @@ */ module tanya.math.mp; -import core.exception; import std.algorithm; +import std.ascii; import std.range; import std.traits; -import tanya.math; +import tanya.container.vector; import tanya.memory; /** @@ -36,8 +36,8 @@ enum Sign : bool */ struct Integer { - private size_t size; - package ubyte* rep; + package digit[] rep; + package ptrdiff_t size; package Sign sign; pure nothrow @safe @nogc invariant @@ -45,6 +45,16 @@ struct Integer assert(this.size > 0 || !this.sign, "0 should be positive."); } + private alias digit = uint; + private alias word = ulong; + + // Count of bits per digit. + private enum : digit + { + digitBitCount = 28, + mask = 0xfffffff, + } + /** * Creates a multiple precision integer. * @@ -59,20 +69,21 @@ struct Integer if (isIntegral!T) { this(allocator); - assign(value); + this = value; } /// Ditto. - this(const ref Integer value, shared Allocator allocator = defaultAllocator) - nothrow @safe @nogc + this(T)(ref T value, shared Allocator allocator = defaultAllocator) + if (is(Unqual!T == Integer)) { this(allocator); - assign(value); + this = value; } /// Ditto. - this(Integer value, shared Allocator allocator = defaultAllocator) + this(T)(T value, shared Allocator allocator = defaultAllocator) nothrow @safe @nogc + if (is(T == Integer)) { this(allocator); if (allocator is value.allocator) @@ -86,7 +97,7 @@ struct Integer } else { - assign(value); + this = value; } } @@ -101,22 +112,6 @@ struct Integer this.allocator_ = allocator; } - private @nogc unittest - { - { - auto integer = Integer(79); - assert(integer.length == 1); - assert(integer.rep[0] == 79); - assert(!integer.sign); - } - { - auto integer = Integer(-2); - assert(integer.length == 1); - assert(integer.rep[0] == 2); - assert(integer.sign); - } - } - /** * Constructs the integer from a sign-magnitude $(D_KEYWORD ubyte) range. * @@ -128,26 +123,41 @@ struct Integer * * Precondition: $(D_INLINECODE allocator !is null) */ - this(R)(const Sign sign, R value, shared Allocator allocator = defaultAllocator) - @trusted - if (isInputRange!R && hasLength!R && is(Unqual!(ElementType!R) == ubyte)) + this(R)(const Sign sign, + R value, + shared Allocator allocator = defaultAllocator) + if (isBidirectionalRange!R && hasLength!R + && is(Unqual!(ElementType!R) == ubyte)) { this(allocator); - while (!value.empty && value.front == 0) + grow(value.length / (digitBitCount / 8) + 1); + + int bit, delta; + + for (; !value.empty; ++this.size) { - value.popFront(); + word w; + for (bit = delta; (bit < digitBitCount) && !value.empty; bit += 8) + { + w |= (cast(word) value.back) << bit; + value.popBack(); + } + + delta = bit - digitBitCount; + this.rep[this.size] |= w & mask; + + if (delta > 0) + { + this.rep[this.size + 1] = (w >> digitBitCount) & mask; + } } - this.rep = allocator.resize(this.rep[0 .. this.size], value.length).ptr; - this.size = value.length; - this.sign = sign; - value.copy(this.rep[0 .. this.size].retro); } - private @nogc unittest + nothrow @safe @nogc unittest { - ubyte[5] range = [ 0x02, 0x11, 0x00, 0x00, 0x01 ]; + ubyte[8] range = [ 0x66, 0x77, 0x88, 0x99, 0xaa, 0xbb, 0xdd, 0xee ]; auto integer = Integer(Sign.positive, range[]); - assert(equal(range[].retro, integer.rep[0 .. integer.size])); + assert(integer == 7383520307673030126); } /** @@ -155,9 +165,9 @@ struct Integer */ this(this) nothrow @trusted @nogc { - auto tmp = allocator.resize!ubyte(null, this.size); + auto tmp = allocator.resize!digit(null, this.size); this.rep[0 .. this.size].copy(tmp); - this.rep = tmp.ptr; + this.rep = tmp; } /** @@ -165,69 +175,15 @@ struct Integer */ ~this() nothrow @trusted @nogc { - allocator.deallocate(this.rep[0 .. this.size]); - } - - /* - * Figures out the minimum amount of space this value will take - * up in bytes and resizes the internal storage. Sets the sign. - */ - private void assign(T)(const ref T value) @trusted - if (isIntegral!T) - { - ubyte size = ulong.sizeof; - ulong mask; - - this.sign = Sign.positive; // Reset the sign. - static if (isSigned!T) - { - const absolute = value.abs; - } - else - { - alias absolute = value; - } - for (mask = 0xff00000000000000; mask >= 0xff; mask >>= 8) - { - if (absolute & mask) - { - break; - } - --size; - } - - this.rep = allocator.resize(this.rep[0 .. this.size], size).ptr; - this.size = size; - static if (isSigned!T) - { - this.sign = value < 0 ? Sign.negative : Sign.positive; - } - - /* Work backward through the int, masking off each byte (up to the - first 0 byte) and copy it into the internal representation in - big-endian format. */ - mask = 0xff; - ubyte shift; - for (size_t i; i < this.size; ++i, mask <<= 8, shift += 8) - { - this.rep[i] = cast(ubyte) ((absolute & mask) >> shift); - } - } - - private void assign(const ref Integer value) nothrow @trusted @nogc - { - this.rep = allocator.resize(this.rep[0 .. this.size], value.size).ptr; - this.size = value.size; - this.sign = value.sign; - value.rep[0 .. value.size].copy(this.rep[0 .. this.size]); + allocator.resize(this.rep, 0); } /** - * Returns: Integer size. + * Returns: Integer byte length. */ @property size_t length() const pure nothrow @safe @nogc { - return this.size; + return (countBits() + 7) / 8; // Round up. } /** @@ -242,154 +198,327 @@ struct Integer ref Integer opAssign(T)(const T value) if (isIntegral!T) { - assign(value); + rep[0 .. this.size].fill(digit.init); + grow(digitBitCount / 8 + 1); + + static if (isSigned!T) + { + ulong absolute; + if (value >= 0) + { + absolute = value; + this.sign = Sign.positive; + } + else + { + absolute = -value; + this.sign = Sign.negative; + } + } + else + { + ulong absolute = value; + this.sign = Sign.positive; + } + + for (this.size = 0; absolute; absolute >>= digitBitCount, ++this.size) + { + this.rep[this.size] = absolute & mask; + } + return this; } /// Ditto. - ref Integer opAssign(const ref Integer value) nothrow @safe @nogc + ref Integer opAssign(T)(ref T value) @trusted + if (is(Unqual!T == Integer)) { - assign(value); + this.rep = allocator.resize(this.rep, value.size); + value.rep[0 .. value.size].copy(this.rep[0 .. value.size]); + this.size = value.size; + this.sign = value.sign; + return this; } /// Ditto. - ref Integer opAssign(Integer value) nothrow @safe @nogc + ref Integer opAssign(T)(T value) nothrow @safe @nogc + if (is(T == Integer)) { - swap(this.rep, value.rep); - swap(this.sign, value.sign); - swap(this.size, value.size); + if (this.allocator is value.allocator) + { + swap(this.rep, value.rep); + swap(this.sign, value.sign); + swap(this.size, value.size); + } + else + { + this = value; + } return this; } - private @nogc unittest + /** + * Casting. + * + * Params: + * T = Target type. + * + * Returns: Converted value. + */ + T opCast(T : bool)() const { - auto integer = Integer(1019); - assert(integer.length == 2); - assert(integer.rep[1] == 0b00000011 && integer.rep[0] == 0b11111011); + return this.size > 0; + } - integer = 3337; - assert(integer.length == 2); - assert(integer.rep[1] == 0b00001101 && integer.rep[0] == 0b00001001); + /// Ditto. + T opCast(T)() const + if (isIntegral!T && isUnsigned!T) + { + T ret; + ubyte shift; + for (size_t i; i < this.size && shift <= T.sizeof * 8; ++i) + { + ret |= (cast(T) this.rep[i]) << shift; + shift += digitBitCount; + } + return ret; + } - integer = 688; - assert(integer.length == 2); - assert(integer.rep[1] == 0b00000010 && integer.rep[0] == 0b10110000); + /// Ditto. + T opCast(T)() const + if (isIntegral!T && isSigned!T) + { + return this.sign ? -(cast(Unsigned!T) this) : cast(Unsigned!T) this; + } + + /// + @safe @nogc unittest + { + auto integer = Integer(79); + assert(cast(ushort) integer == 79); + + integer = -79; + assert(cast(short) integer == -79); + + integer = 4294967295; + assert(cast(long) integer == 4294967295); + + integer = -4294967295; + assert(cast(long) integer == -4294967295); + + integer = long.min; + assert(cast(long) integer == long.min); + + integer = long.min + 1; + assert(cast(long) integer == long.min + 1); integer = 0; - assert(integer.length == 0); + assert(cast(long) integer == 0); } - /* - * Extend the space for this by 1 byte and set the LSB to 1. + /* trim unused digits + * + * This is used to ensure that leading zero digits are + * trimed and the leading "size" digit will be non-zero + * Typically very fast. Also fixes the sign if there + * are no more leading digits */ - private void expand() nothrow @trusted @nogc + void contract() nothrow @safe @nogc { - rep = allocator.resize(this.rep[0 .. this.size], this.size + 1).ptr; - this.rep[this.size] = 0x01; - ++this.size; - } - - /* - * Go through this and see how many of the left-most bytes are unused. - * Remove them and resize this appropriately. - */ - private void contract() nothrow @trusted @nogc - { - const i = this.rep[0 .. this.size] - .retro - .countUntil!((const ref a) => a != 0); - - if (i > 0) + /* decrease size while the most significant digit is + * zero. + */ + while ((this.size > 0) && (this.rep[this.size - 1] == 0)) { - this.rep = allocator.resize(this.rep[0 .. this.size], this.size - i).ptr; - this.size -= i; + --this.size; } - else if (i == -1) + + /* reset the sign flag if size == 0 */ + if (this.size == 0) { this.sign = Sign.positive; - this.rep = allocator.resize(this.rep[0 .. this.size], 0).ptr; - this.size = 0; } } - private void add(const ref Integer summand) nothrow @trusted @nogc + private void grow(const size_t size) nothrow @trusted @nogc { - if (summand.length > this.length) + if (this.rep.length >= size) { - this.rep = allocator.resize!ubyte(this.rep[0 .. this.size], summand.size).ptr; - this.rep[this.size .. summand.size].initializeAll(); + return; + } + const oldLength = this.rep.length; + allocator.resize(this.rep, size); + this.rep[oldLength .. $].fill(digit.init); + } - this.size = summand.size; + private size_t countBits() const pure nothrow @safe @nogc + { + if (this.size == 0) + { + return 0; + } + auto r = (this.size - 1) * digitBitCount; + digit q = this.rep[this.size - 1]; + + while (q > 0) + { + ++r; + q >>= (cast(digit) 1); + } + return r; + } + + private void add(ref const Integer summand, ref Integer sum) + const nothrow @safe @nogc + { + const(digit)[] max, min; + + if (this.size > summand.size) + { + min = summand.rep[0 .. summand.size]; + max = this.rep[0 .. this.size]; + } + else + { + min = this.rep[0 .. this.size]; + max = summand.rep[0 .. summand.size]; + } + sum.grow(max.length + 1); + + const oldSize = sum.size; + sum.size = cast(int) (max.length + 1); + + auto result = sum.rep; + digit carry; + foreach (i, ref const d; min) + { + result.front = d + max.front + carry; + carry = result.front >> digitBitCount; + result.front &= mask; + + max.popFront(); + result.popFront(); } - bool carry; + // Copy higher digests if one of the summands is greater than another + // one. + for (; !max.empty; max.popFront(), result.popFront()) + { + result.front = max.front + carry; + carry = result.front >> digitBitCount; + result.front &= mask; + } + result.front = carry; + + // Clear digits above the old size. + for (auto i = sum.size; i < oldSize; ++i) + { + sum.rep[i] = 0; + } + + sum.contract(); + } + + private void add(const digit summand, ref Integer sum) + const nothrow @safe @nogc + { + sum.grow(this.size + 2); + + sum.rep[0] = this.rep[0] + summand; + auto carry = sum.rep[0] >> digitBitCount; + sum.rep[0] &= mask; + size_t i; - size_t j; - do + for (i = 1; i < this.size; ++i) { - uint sum; - if (j < summand.size) - { - sum = this.rep[i] + summand.rep[j] + carry; - ++j; - } - else - { - sum = this.rep[i] + carry; - } - - carry = sum > 0xff; - this.rep[i] = cast(ubyte) sum; + sum.rep[i] = this.rep[i] + carry; + carry = sum.rep[i] >> digitBitCount; + sum.rep[i] &= mask; } - while (++i < this.size); + sum.rep[i++] = carry; - if (carry) + for (; i < sum.size; ++i) { - // Still overflowed; allocate more space - expand(); + sum.rep[i] = 0; } + sum.size = this.size + 1; + sum.contract(); } - private void subtract(const ref Integer subtrahend) nothrow @trusted @nogc + private void subtract(ref const Integer subtrahend, ref Integer difference) + const nothrow @safe @nogc { + difference.grow(this.size); + + const oldSize = difference.size; + difference.size = this.size; + size_t i; - size_t j; - bool borrow; + digit carry; - while (i < this.size) + for (i = 0; i < subtrahend.size; ++i) { - int difference; - - if (j < subtrahend.size) - { - difference = this.rep[i] - subtrahend.rep[j] - borrow; - ++j; - } - else - { - difference = this.rep[i] - borrow; - } - borrow = difference < 0; - this.rep[i] = cast(ubyte) difference; - - ++i; + difference.rep[i] = (this.rep[i] - subtrahend.rep[i]) - carry; + carry = difference.rep[i] >> (cast(digit) (8 * digit.sizeof) - 1); + difference.rep[i] &= mask; } - if (borrow && i < this.size && this.rep[i]) + // Copy higher digests if the minuend has more digits than the + // subtrahend. + for (; i < this.size; ++i) { - --this.rep[i]; + difference.rep[i] = this.rep[i] - carry; + carry = difference.rep[i] >> (cast(digit) ((8 * digit.sizeof) - 1)); + difference.rep[i] &= mask; } - contract(); + + // Clear digits above the size. + for (i = difference.size; i < oldSize; ++i) + { + difference.rep[i] = 0; + } + + difference.contract(); } - private int compare(const ref Integer that) const nothrow @trusted @nogc + private void subtract(const digit subtrahend, ref Integer difference) + const nothrow @safe @nogc { - if (length > that.length) + difference.grow(this.size); + + const oldSize = difference.size; + + difference.sign = this.sign; + difference.size = this.size; + + difference.rep[0] = this.rep[0] - subtrahend; + auto carry = difference.rep[0] >> ((digit.sizeof * 8) - 1); + difference.rep[0] &= mask; + + size_t i; + for (i = 1; i < this.size; ++i) + { + difference.rep[i] = this.rep[i] - carry; + carry = difference.rep[i] >> ((digit.sizeof * 8) - 1); + difference.rep[i] &= mask; + } + + for (; i < oldSize; ++i) + { + difference.rep[i] = 0; + } + difference.contract(); + } + + // Compare the magnitude. + private int compare(ref const Integer that) const pure nothrow @safe @nogc + { + if (this.size > that.size) { return 1; } - else if (length < that.length) + else if (this.size < that.size) { return -1; } @@ -402,18 +531,33 @@ struct Integer * Comparison. * * Params: + * I = Comparand type. * that = The second integer. * * Returns: A positive number if $(D_INLINECODE this > that), a negative * number if $(D_INLINECODE this < that), `0` otherwise. */ - int opCmp(I : Integer)(const auto ref I that) const + int opCmp(I : Integer)(auto ref const I that) const { if (this.sign != that.sign) { - return this.sign ? -1 : 1; + if (this.sign == Sign.negative) + { + return -1; + } + else + { + return 1; + } + } + if (this.sign == Sign.negative) + { + return that.compare(this); + } + else + { + return compare(that); } - return compare(that); } /// @@ -434,7 +578,7 @@ struct Integer } /// Ditto. - int opCmp(I)(const auto ref I that) const + int opCmp(I)(const I that) const if (isIntegral!I) { if (that < 0 && !this.sign) @@ -445,7 +589,7 @@ struct Integer { return -1; } - else if (this.length > I.sizeof) + else if (this.size > I.sizeof) { return this.sign ? -1 : 1; } @@ -475,11 +619,12 @@ struct Integer /** * Params: + * I = Comparand type. * that = The second integer. * * Returns: Whether the two integers are equal. */ - bool opEquals(I)(const auto ref I that) const + bool opEquals(I)(auto ref const I that) const if (is(I : Integer) || isIntegral!I) { return opCmp!I(that) == 0; @@ -503,160 +648,112 @@ struct Integer * * Returns: $(D_KEYWORD this). */ - ref Integer opOpAssign(string op : "+")(const auto ref Integer operand) + ref Integer opOpAssign(string op : "+")(auto ref const Integer operand) { if (this.sign == operand.sign) { - add(operand); + add(operand, this); + } + else if (compare(operand) < 0) + { + this.sign = operand.sign; + operand.subtract(this, this); } else { - if (this >= operand) - { - subtract(operand); - } - else - { - // Swap the operands. - auto tmp = Integer(this, this.allocator); - this = Integer(operand, this.allocator); - subtract(tmp); - this.sign = operand.sign; - } - + subtract(operand, this); } return this; } - private @nogc unittest + /// + unittest { { auto h1 = Integer(1019); auto h2 = Integer(3337); h1 += h2; - assert(h1.length == 2); - assert(h1.rep[1] == 0x11 && h1.rep[0] == 0x04); + assert(h1 == 4356); } { auto h1 = Integer(4356); auto h2 = Integer(2_147_483_647); - ubyte[4] expected = [ 0x03, 0x11, 0x00, 0x80 ]; h1 += h2; - assert(h1.rep[0 .. h1.size] == expected); + assert(h1 == 2147488003); } { auto h1 = Integer(2147488003L); auto h2 = Integer(2_147_483_647); - ubyte[5] expected = [ 0x02, 0x11, 0x00, 0x00, 0x01 ]; h1 += h2; - assert(h1.rep[0 .. h1.size] == expected); + assert(h1 == 4294971650); } } /// Ditto. - ref Integer opOpAssign(string op : "-")(const auto ref Integer operand) + ref Integer opOpAssign(string op : "-")(auto ref const Integer operand) { - if (operand.sign == this.sign) + if (this.sign != operand.sign) { - if (this >= operand) - { - subtract(operand); - } - else - { - // Swap the operands. - auto tmp = Integer(this, this.allocator); - this = Integer(operand, this.allocator); - subtract(tmp); - - if (operand.sign || this.size == 0) - { - this.sign = Sign.positive; - } - else - { - this.sign = Sign.negative; - } - } + add(operand, this); + } + else if (compare(operand) >= 0) + { + subtract(operand, this); } else { - add(operand); + operand.subtract(this, this); + this.sign = this.sign ? Sign.positive : Sign.negative; } return this; } - private @nogc unittest + /// + unittest { { auto h1 = Integer(3); auto h2 = Integer(4); h1 -= h2; - assert(h1.length == 1); - assert(h1.rep[0] == 0x01); - assert(h1.sign == Sign.negative); + assert(h1 == -1); } { auto h1 = Integer(8589934590L); auto h2 = Integer(2147483647); - ubyte[5] expected = [ 0xff, 0xff, 0xff, 0x7f, 0x01 ]; - h1 -= h2; - assert(h1.rep[0 .. h1.size] == expected); + assert(h1 == 6442450943); } { auto h1 = Integer(6442450943); auto h2 = Integer(4294967294); - ubyte[4] expected = [ 0x01, 0x00, 0x00, 0x80 ]; h1 -= h2; - assert(h1.rep[0 .. h1.size] == expected); + assert(h1 == 2147483649); } { auto h1 = Integer(2147483649); auto h2 = Integer(h1); h1 -= h2; - assert(h1.length == 0); + assert(h1 == 0); } } /// Ditto. - ref Integer opOpAssign(string op : "*")(const auto ref Integer operand) @trusted + ref Integer opOpAssign(string op : "*")(auto ref const Integer operand) { - size_t i; - if (this.length == 0) - { - return this; - } - else if (operand.length == 0) - { - this = 0; - return this; - } - auto temp = Integer(this, this.allocator); - auto sign = this.sign != operand.sign; + const digits = this.size + operand.size + 1; - this = 0; - do - { - for (ubyte mask = 0x01; mask; mask <<= 1) - { - if (mask & operand.rep[i]) - { - add(temp); - } - temp <<= 1; - } - ++i; - } - while (i < operand.size); + multiply(operand, this, digits); - this.sign = sign ? Sign.negative : Sign.positive; + if (this.size > 0) + { + this.sign = this.sign == operand.sign ? Sign.positive : Sign.negative; + } return this; } /// - @safe @nogc unittest + nothrow @safe @nogc unittest { auto h1 = Integer(123); auto h2 = Integer(456); @@ -664,284 +761,155 @@ struct Integer assert(h1 == 56088); } - private @nogc unittest - { - assert((Integer(1) * Integer()).length == 0); - } - /// Ditto. - ref Integer opOpAssign(string op : "^^")(const auto ref Integer operand) - @trusted - { - size_t i; - - auto tmp1 = Integer(this, this.allocator); - this = 1; - - do - { - for (ubyte mask = 0x01; mask; mask <<= 1) - { - if (operand.rep[i] & mask) - { - this *= tmp1; - } - // Square tmp1 - auto tmp2 = tmp1; - tmp1 *= tmp2; - } - ++i; - } - while (i < operand.size); - - return this; - } - - private @nogc unittest - { - auto h1 = Integer(2); - auto h2 = Integer(4); - - h1 ^^= h2; - assert(h1.length == 1); - assert(h1.rep[0] == 0x10); - - h1 = Integer(2342); - h1 ^^= h2; - ubyte[6] expected = [ 0x10, 0x31, 0x9c, 0xab, 0x5c, 0x1b ]; - assert(h1.rep[0 .. h1.size] == expected); - } - - /// Ditto. - ref Integer opOpAssign(string op)(const auto ref Integer operand) @trusted - if ((op == "%") || (op == "/")) + ref Integer opOpAssign(string op : "/")(auto ref const Integer operand) in { assert(operand.length > 0, "Division by zero."); } body { - auto divisor = Integer(operand, this.allocator); - size_t bitSize; - - // First, left-shift divisor until it's >= than the dividend - while (compare(divisor) > 0) - { - divisor <<= 1; - ++bitSize; - } - static if (op == "/") - { - auto quotient = allocator.resize!ubyte(null, bitSize / 8 + 1); - quotient.initializeAll(); - } - - // bitPosition keeps track of which bit, of the quotient, - // is being set or cleared on the current operation. - auto bitPosition = 8 - (bitSize % 8) - 1; - do - { - if (compare(divisor) >= 0) - { - subtract(divisor); // dividend -= divisor - static if (op == "/") - { - quotient[$ - 1 - bitPosition / 8] |= 0x80 >> (bitPosition % 8); - } - } - if (bitSize != 0) - { - divisor >>= 1; - } - ++bitPosition; - } - while (bitSize--); - - static if (op == "/") - { - allocator.resize(this.rep[0 .. this.size], 0); - this.rep = quotient.ptr; - this.size = quotient.length; - this.sign = this.sign == divisor.sign ? Sign.positive : Sign.negative; - contract(); - } - + divide(operand, this); return this; } - private @nogc unittest + /// Ditto. + ref Integer opOpAssign(string op : "%")(auto ref const Integer operand) + in + { + assert(operand.length > 0, "Division by zero."); + } + body + { + divide(operand, null, this); + return this; + } + + nothrow @safe @nogc unittest { auto h1 = Integer(18); auto h2 = Integer(4); h1 %= h2; - assert(h1.length == 1); - assert(h1.rep[0] == 0x02); + assert(h1 == 2); h1 = 8; h1 %= h2; - assert(h1.length == 0); + assert(h1 == 0); h1 = 7; h1 %= h2; - assert(h1.length == 1); - assert(h1.rep[0] == 0x03); + assert(h1 == 3); h1 = 56088; h2 = 456; h1 /= h2; - assert(h1.length == 1); - assert(h1.rep[0] == 0x7b); + assert(h1 == 123); } /// Ditto. - ref Integer opOpAssign(string op : ">>")(const size_t operand) @trusted + ref Integer opOpAssign(string op : ">>")(const size_t operand) { - const step = operand / 8; - - if (step >= this.length) + if (operand == 0) { - this.rep = allocator.resize(this.rep[0 .. this.size], 0).ptr; - this.size = 0; return this; } - - size_t i, j; - ubyte carry; - const bit = operand % 8; - const delta = 8 - bit; - - for (j = step; j < length - 1; ++i, ++j) + if (operand >= digitBitCount) { - carry = cast(ubyte) (this.rep[i + 1] << delta); - this.rep[i] = (this.rep[i] >> bit) | carry; + shiftRight(operand / digitBitCount); } - this.rep[i] = this.rep[j] >> bit; - size_t newSize = length - step; - if (this.rep[i] == 0) + const bit = cast(digit) (operand % digitBitCount); + if (bit != 0) { - --newSize; - } - this.rep = allocator.resize(this.rep[0 .. this.size], newSize).ptr; - this.size = newSize; + const mask = ((cast(digit) 1) << bit) - 1; + const shift = digitBitCount - bit; + digit carry; + foreach (ref d; this.rep[0 .. this.size].retro) + { + const newCarry = d & mask; + d = (d >> bit) | (carry << shift); + carry = newCarry; + } + } + this.contract(); return this; } - private @nogc unittest + /// + nothrow @safe @nogc unittest { auto integer = Integer(4294967294); integer >>= 10; - assert(integer.length == 3); - assert(integer.rep[2] == 0x3f && integer.rep[1] == 0xff && integer.rep[0] == 0xff); + assert(integer == 4194303); integer = 27336704; integer >>= 1; - assert(integer.length == 3); - assert(integer.rep[2] == 0xd0 && integer.rep[1] == 0x90 && integer.rep[0] == 0x00); + assert(integer == 13668352); integer = 4294967294; integer >>= 20; - assert(integer.length == 2); - assert(integer.rep[1] == 0x0f && integer.rep[0] == 0xff); + assert(integer == 4095); integer >>= 0; - assert(integer.length == 2); - assert(integer.rep[1] == 0x0f && integer.rep[0] == 0xff); + assert(integer == 4095); integer >>= 20; - assert(integer.length == 0); + assert(integer == 0); integer >>= 2; - assert(integer.length == 0); + assert(integer == 0); integer = 1431655765; integer >>= 16; - assert(integer.length == 2); - assert(integer.rep[1] == 0x55 && integer.rep[0] == 0x55); + assert(integer == 21845); integer >>= 16; - assert(integer.length == 0); + assert(integer == 0); } /// Ditto. - ref Integer opOpAssign(string op : "<<")(const size_t operand) @trusted + ref Integer opOpAssign(string op : "<<")(const size_t operand) { - auto i = this.length; - size_t j; - size_t newSize; - const bit = operand % 8; - const delta = 8 - bit; - const step = operand / 8; - - auto carry = cast(ubyte) (this.rep[this.size - 1] >> delta); - if (carry != 0) + const step = operand / digitBitCount; + if (this.rep.length < this.size + step + 1) { - newSize = length + step + 1; - this.rep = allocator.resize(this.rep[0 .. this.size], newSize).ptr; - this.size = newSize; - j = newSize - 1; - this.rep[j] = carry; + grow(this.size + step + 1); } - else + if (operand >= digitBitCount) { - newSize = length + step; - this.rep = allocator.resize(this.rep[0 .. this.size], newSize).ptr; - this.size = j = newSize; + shiftLeft(step); } - --i, --j; - for (; i > 0; --i, --j) + const bit = cast(digit) (operand % digitBitCount); + if (bit != 0) { - this.rep[i] = cast(ubyte) (this.rep[j] << bit) | (this.rep[j - 1] >> delta); - } - this.rep[i] = cast(ubyte) (this.rep[j] << bit); - this.rep[0 .. step].fill(cast(ubyte) 0); + const mask = ((cast(digit) 1) << bit) - 1; + const shift = digitBitCount - bit; + digit carry; + foreach (ref d; this.rep[0 .. this.size]) + { + const newCarry = (d >> shift) & mask; + d = ((d << bit) | carry) & this.mask; + carry = newCarry; + } + + if (carry != 0) + { + this.rep[this.size++] = carry; + } + } + this.contract(); return this; } - private @nogc unittest + /// + nothrow @safe @nogc unittest { auto integer = Integer(4294967295); - ubyte[5] expected = [ 0xfe, 0xff, 0xff, 0xff, 0x01 ]; integer <<= 1; - assert(integer.rep[0 .. integer.size] == expected); - } - - private void decrement() nothrow @trusted @nogc - in - { - assert(this.length > 0); - } - body - { - for (ubyte* p = this.rep; p < this.rep + this.size; ++p) - { - --*p; - if (*p != ubyte.max) - { - break; - } - } - contract(); - } - - private void increment() nothrow @trusted @nogc - { - ubyte* p; - for (p = this.rep; p < this.rep + this.size; ++p) - { - ++*p; - if (*p > 0) - { - return; - } - } - if (p == this.rep + this.size) - { - expand(); - } + assert(integer == 8589934590); } /** @@ -954,15 +922,15 @@ struct Integer */ Integer opUnary(string op : "~")() const { - auto ret = Integer(this, this.allocator); - ret.rep[0 .. ret.size].each!((ref a) => a = ~a); + auto ret = Integer(this, allocator); + ret.rep[0 .. ret.size].each!((ref a) => a = ~a & mask); return ret; } /// Ditto. Integer opUnary(string op : "-")() const { - auto ret = Integer(this, this.allocator); + auto ret = Integer(this, allocator); ret.sign = ret.sign ? Sign.positive : Sign.negative; return ret; } @@ -980,35 +948,24 @@ struct Integer return this; } - private @nogc unittest + // + nothrow @safe @nogc unittest { auto h1 = Integer(79); Integer h2; h2 = +h1; - assert(h2.length == 1); - assert(h2.rep[0] == 79); - assert(!h2.sign); - assert(h2 !is h1); + assert(h2 == 79); h2 = -h1; - assert(h2.length == 1); - assert(h2.rep[0] == 79); - assert(h2.sign); + assert(h2 == -79); + assert(h1 == 79); h1 = -h2; - assert(h2.length == 1); - assert(h2.rep[0] == 79); - assert(h2.sign); - assert(h2 !is h1); - - h1 = -h2; - assert(h1.length == 1); - assert(h1.rep[0] == 79); - assert(!h1.sign); + assert(h1 == 79); h2 = ~h1; - assert(h2.rep[0] == ~cast(ubyte) 79); + assert(h2 == ~cast(ubyte) 79); } /// Ditto. @@ -1016,15 +973,11 @@ struct Integer { if (this.sign) { - decrement(); - if (this.length == 0) - { - this.sign = Sign.positive; - } + subtract(1, this); } else { - increment(); + add(1, this); } return this; } @@ -1032,59 +985,56 @@ struct Integer /// Ditto. ref Integer opUnary(string op : "--")() { - if (this.sign) + if (this.size == 0) { - increment(); + add(1, this); + this.sign = Sign.negative; + } + else if (this.sign) + { + add(1, this); } else { - decrement(); + subtract(1, this); } return this; } - private @nogc unittest + /// + nothrow @safe @nogc unittest { Integer integer; ++integer; - assert(integer.length == 1); - assert(integer.rep[0] == 0x01); + assert(integer == 1); --integer; - assert(integer.length == 0); + assert(integer == 0); integer = 511; ++integer; - assert(integer.length == 2); - assert(integer.rep[1] == 0x02 && integer.rep[0] == 0x00); + assert(integer == 512); --integer; - assert(integer.length == 2); - assert(integer.rep[1] == 0x01 && integer.rep[0] == 0xff); + assert(integer == 511); integer = 79; ++integer; - assert(integer.length == 1); - assert(integer.rep[0] == 0x50); + assert(integer == 80); --integer; - assert(integer.length == 1); - assert(integer.rep[0] == 0x4f); - - integer = 65535; - ++integer; - assert(integer.length == 3); - assert(integer.rep[2] == 0x01 && integer.rep[1] == 0x00 && integer.rep[0] == 0x00); - - --integer; - assert(integer.length == 2); - assert(integer.rep[1] == 0xff && integer.rep[0] == 0xff); + assert(integer == 79); integer = -2; ++integer; - assert(integer.length == 1); - assert(integer.rep[0] == 0x01); + assert(integer == -1); + + ++integer; + assert(integer == 0); + + --integer; + assert(integer == -1); } /** @@ -1093,13 +1043,13 @@ struct Integer * Params: * op = Operation. * operand = The second operand. + * + * Returns: Result. */ - Integer opBinary(string op)(const auto ref Integer operand) const - if (op == "+" || op == "-" || op == "*" || op == "^^") + Integer opBinary(string op)(auto ref const Integer operand) const + if ((op == "+" || op == "-") || (op == "*")) { - auto ret = Integer(this, this.allocator); - mixin("ret " ~ op ~ "= operand;"); - return ret; + mixin("return Integer(this, allocator) " ~ op ~ "= operand;"); } /// Ditto. @@ -1111,84 +1061,406 @@ struct Integer } body { - auto ret = Integer(this, this.allocator); - mixin("ret " ~ op ~ "= operand;"); - return ret; + mixin("return Integer(this, allocator) " ~ op ~ "= operand;"); } /// Ditto. - Integer opBinary(string op)(const auto ref size_t operand) const + Integer opBinary(string op)(const size_t operand) const if (op == "<<" || op == ">>") { - auto ret = Integer(this, this.allocator); - mixin("ret " ~ op ~ "= operand;"); - return ret; + mixin("return Integer(this, allocator) " ~ op ~ "= operand;"); + } + + // Shift right a certain amount of digits. + private void shiftRight(const size_t operand) nothrow @safe @nogc + { + if (operand == 0) + { + return; + } + if (this.size <= operand) + { + this = 0; + return; + } + const reducedSize = this.size - operand; + + this.rep[operand .. this.size].copy(this.rep[0 .. reducedSize]); + this.rep[reducedSize .. this.size].fill(digit.init); + this.size = reducedSize; + } + + // Shift left a certain amount of digits. + private void shiftLeft(const size_t operand) nothrow @safe @nogc + { + if (operand == 0) + { + return; + } + const increasedSize = this.size + operand; + grow(increasedSize); + + this.size = increasedSize; + + auto top = this.size - 1; + auto bottom = this.size - 1 - operand; + + for (; top >= operand; --bottom, --top) + { + this.rep[top] = this.rep[bottom]; + } + this.rep[0 .. operand].fill(digit.init); + } + + private void multiply(const digit factor, ref Integer product) + const nothrow @safe @nogc + { + product.grow(this.size + 1); + product.sign = this.sign; + + word carry; + size_t i; + + for (i = 0; i < this.size; ++i) + { + auto newCarry = carry + (cast(word) this.rep[i] * factor); + product.rep[i] = newCarry & mask; + carry = newCarry >> digitBitCount; + } + product.rep[i++] = carry & mask; + + for (; i < this.size; ++i) + { + product.rep[i] = 0; + } + product.size = this.size + 1; + product.contract(); + } + + private void multiply(ref const Integer factor, + ref Integer product, + const size_t digits) const nothrow @safe @nogc + { + Integer intermediate; + intermediate.grow(digits); + intermediate.size = digits; + + for (size_t i; i < this.size; ++i) + { + const limit = min(factor.size, digits - i); + word carry; + auto k = i; + + for (size_t j; j < limit; ++j, ++k) + { + const result = cast(word) intermediate.rep[k] + + (cast(word) this.rep[i] * factor.rep[j]) + + carry; + intermediate.rep[k] = result & mask; + carry = result >> digitBitCount; + } + if (k < digits) + { + intermediate.rep[k] = carry & mask; + } + } + intermediate.contract(); + swap(product, intermediate); + } + + private void divide(Q, ARGS...)(ref const Integer divisor, + auto ref Q quotient, + ref ARGS args) + const nothrow @safe @nogc + if ((is(Q : typeof(null)) + || (is(Q : Integer) && __traits(isRef, quotient))) + && (ARGS.length == 0 || (ARGS.length == 1 && is(ARGS[0] : Integer)))) + in + { + assert(divisor != 0, "Division by zero."); + } + body + { + if (compare(divisor) < 0) + { + static if (ARGS.length == 1) + { + args[0] = this; + } + static if (!is(Q == typeof(null))) + { + quotient = 0; + } + return; + } + + Integer q, t1, t2; + q.grow(this.size + 2); + q.size = this.size + 2; + + t1.grow(2); + t2.grow(3); + + auto x = Integer(this); + auto y = Integer(divisor); + + const sign = this.sign == divisor.sign ? Sign.positive : Sign.negative; + x.sign = y.sign = Sign.positive; + + auto norm = y.countBits() % digitBitCount; + if (norm < digitBitCount - 1) + { + norm = digitBitCount - 1 - norm; + x <<= norm; + y <<= norm; + } + else + { + norm = 0; + } + + auto n = x.size - 1; + auto t = y.size - 1; + + y.shiftLeft(n - t); + + while (x >= y) + { + ++q.rep[n - t]; + x -= y; + } + + y.shiftRight(n - t); + + for (auto i = n; i >= (t + 1); --i) + { + if (i > x.size) + { + continue; + } + if (x.rep[i] == y.rep[t]) + { + q.rep[(i - t) - 1] = (((cast(digit) 1) << digitBitCount) - 1); + } + else + { + word tmp = (cast(word) x.rep[i]) << digitBitCount; + tmp |= x.rep[i - 1]; + tmp /= y.rep[t]; + if (tmp > mask) + { + tmp = mask; + } + q.rep[i - t - 1] = tmp & mask; + } + + q.rep[i - t - 1] = (q.rep[i - t - 1] + 1) & mask; + do + { + q.rep[i - t - 1] = (q.rep[i - t - 1] - 1) & mask; + + // Left hand. + t1 = 0; + t1.rep[0] = ((t - 1) < 0) ? 0 : y.rep[t - 1]; + t1.rep[1] = y.rep[t]; + t1.size = 2; + t1.multiply(q.rep[i - t - 1], t1); + + // Right hand. + t2.rep[0] = ((i - 2) < 0) ? 0 : x.rep[i - 2]; + t2.rep[1] = ((i - 1) < 0) ? 0 : x.rep[i - 1]; + t2.rep[2] = x.rep[i]; + t2.size = 3; + } + while (t1.compare(t2) > 0); + + y.multiply(q.rep[i - t - 1], t1); + + t1.shiftLeft(i - t - 1); + + x -= t1; + + if (x.sign == Sign.negative) + { + t1 = y; + t1.shiftLeft(i - t - 1); + x += t1; + + q.rep[i - t - 1] = (q.rep[i - t - 1] - 1) & mask; + } + } + + x.sign = (x.size == 0) ? Sign.positive : this.sign; + static if (!is(Q == typeof(null))) + { + q.contract(); + swap(q, quotient); + quotient.sign = sign; + } + static if (ARGS.length == 1) + { + x >>= norm; + swap(x, args[0]); + } + } + + private Integer square() nothrow @safe @nogc + { + Integer result; + const resultSize = 2 * this.size + 1; + + result.grow(resultSize); + result.size = resultSize; + + for (size_t i; i < this.size; ++i) + { + const doubleI = 2 * i; + word product = cast(word) result.rep[doubleI] + + (cast(word) this.rep[i] * this.rep[i]); + + result.rep[doubleI] = product & mask; + + word carry = product >> digitBitCount; + size_t k = doubleI + 1; + + for (auto j = i + 1; j < this.size; ++j, ++k) + { + product = (cast(word) this.rep[i]) * (cast(word) this.rep[j]); + product = (cast(word) result.rep[k]) + (2 * product) + carry; + + result.rep[k] = product & mask; + carry = product >> digitBitCount; + } + for (; carry != 0; ++k) + { + product = (cast(word) result.rep[k]) + carry; + result.rep[k] = product & mask; + carry = product >> digitBitCount; + } + } + result.contract(); + + return result; + } + + Vector!ubyte toVector() const nothrow @safe @nogc + { + Vector!ubyte vector; + bool firstBit; + ubyte carry; + + if (this.size == 0) + { + return vector; + } + vector.reserve(this.size * digit.sizeof); + + // The first digit needs extra handling since it can have leading + // non significant zeros. + int digitCount = digitBitCount - 8; + const first = this.rep[this.size - 1]; + const prevBitCount = ((this.size - 1) * digitBitCount); + const fullBytesBitCount = ((prevBitCount - 1) / 8 + 1) * 8; + + // Find out the right alignment of the first byte. + if ((fullBytesBitCount - prevBitCount) == 0) + { + digitCount -= digit.sizeof * 8 - digitBitCount; + } + for (; digitCount >= 0; digitCount -= 8) + { + if (firstBit || ((first >> digitCount) != 0)) + { + firstBit = true; + vector.insertBack(cast(ubyte) (first >> digitCount)); + } + } + if (digitCount >= -8) + { + carry = (first << -digitCount) & 0xff; + digitCount += digitBitCount; + } + else + { + carry = 0; + digitCount = digitBitCount - 8; + } + + foreach_reverse (d; this.rep[0 .. this.size - 1]) + { + if (carry != 0) // Check the carry from the previous digit. + { + vector.insertBack(cast(ubyte) (carry | (d >> digitCount))); + digitCount -= 8; + } + // Write the digit by bytes. + for (; digitCount >= 0; digitCount -= 8) + { + vector.insertBack(cast(ubyte) (d >> digitCount)); + } + // Check for an incomplete byte. + if (digitCount >= -8) + { + carry = (d << -digitCount) & 0xff; + digitCount += digitBitCount; + } + else + { + carry = 0; + digitCount = digitBitCount - 8; + } + } + if (carry != 0) + { + vector.insertBack(cast(ubyte) (carry >> (digitBitCount - digitCount))); + } + + return vector; } /// - @safe @nogc unittest + nothrow @safe @nogc unittest { - auto integer1 = Integer(425); - auto integer2 = integer1 << 1; - assert(integer2 == 850); - - integer2 = integer1 >> 1; - assert(integer2 == 212); - } - - /** - * Casting. - * - * Params: - * T = Target type. - * - * Returns: $(D_KEYWORD false) if the $(D_PSYMBOL Integer) is 0, - * $(D_KEYWORD true) otherwise. - */ - T opCast(T : bool)() const - { - return this.length > 0; - } - - /// Ditto. - T opCast(T)() const - if (isIntegral!T && isSigned!T) - { - return this.sign ? -(cast(Unsigned!T) this) : cast(Unsigned!T) this; - } - - /// Ditto. - T opCast(T)() const @trusted - if (isIntegral!T && isUnsigned!T) - { - if (this.length == 0) { - return 0; + auto integer = Integer(0x66778899aabbddee); + ubyte[8] expected = [ 0x66, 0x77, 0x88, 0x99, 0xaa, 0xbb, 0xdd, 0xee ]; + + auto vector = integer.toVector(); + assert(equal(vector[], expected[])); } - T ret; - const(ubyte)* src = this.rep; - ubyte shift; - for (; src < this.rep + this.size && shift <= T.sizeof * 8; ++src, shift += 8) { - ret |= (cast(T) *src) << shift; + auto integer = Integer(0x03); + ubyte[1] expected = [ 0x03 ]; + + auto vector = integer.toVector(); + assert(equal(vector[], expected[])); } - return ret; - } + { + ubyte[63] expected = [ + 0x02, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, + 0x09, 0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f, 0x10, + 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18, + 0x19, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f, 0x20, + 0x21, 0x22, 0x23, 0x24, 0x25, 0x26, 0x27, 0x28, + 0x29, 0x2a, 0x2b, 0x2c, 0x2d, 0x2e, 0x2f, 0x30, + 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37, 0x38, + 0x39, 0x3a, 0x3b, 0x00, 0x61, 0x62, 0x63, + ]; + auto integer = Integer(Sign.positive, expected[]); - /// - @safe @nogc unittest - { - auto integer = Integer(79); - assert(cast(long) integer == 79); + auto vector = integer.toVector(); + assert(equal(vector[], expected[])); + } + { + ubyte[14] expected = [ + 0x22, 0x33, 0x44, 0x55, 0x05, 0x06, 0x07, + 0x08, 0x3a, 0x3b, 0x00, 0x61, 0x62, 0x63, + ]; + auto integer = Integer(Sign.positive, expected[]); - integer = -79; - assert(cast(long) integer == -79); - - integer = 4294967295; - assert(cast(long) integer == 4294967295); - - integer = -4294967295; - assert(cast(long) integer == -4294967295); + auto vector = integer.toVector(); + assert(equal(vector[], expected[])); + } } mixin DefaultAllocator; diff --git a/source/tanya/math/package.d b/source/tanya/math/package.d index fd5f2a2..3537bc3 100644 --- a/source/tanya/math/package.d +++ b/source/tanya/math/package.d @@ -90,18 +90,19 @@ body size_t i; auto tmp1 = Integer(x, x.allocator); auto result = Integer(x.allocator); + bool firstBit; - if (x.length == 0 && y.length != 0) + if (x.size == 0 && y.size != 0) { - i = y.length; + i = y.size; } else { result = 1; } - while (i < y.length) + while (i < y.size) { - for (ubyte mask = 0x01; mask; mask <<= 1) + for (uint mask = 0x01; mask != 0x10000000; mask <<= 1) { if (y.rep[i] & mask) {