summaryrefslogtreecommitdiff
path: root/source
diff options
context:
space:
mode:
Diffstat (limited to 'source')
-rw-r--r--source/tanya/math/mp.d2300
-rw-r--r--source/tanya/math/package.d226
2 files changed, 1367 insertions, 1159 deletions
diff --git a/source/tanya/math/mp.d b/source/tanya/math/mp.d
index d36ec18..50efb64 100644
--- a/source/tanya/math/mp.d
+++ b/source/tanya/math/mp.d
@@ -1,1060 +1,1240 @@
-/* This Source Code Form is subject to the terms of the Mozilla Public
- * License, v. 2.0. If a copy of the MPL was not distributed with this
- * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
-
-/**
- * Copyright: Eugene Wissner 2016.
- * License: $(LINK2 https://www.mozilla.org/en-US/MPL/2.0/,
- * Mozilla Public License, v. 2.0).
- * Authors: $(LINK2 mailto:belka@caraus.de, Eugene Wissner)
- */
-module tanya.math.mp;
-
-import std.algorithm.iteration;
-import std.algorithm.searching;
-import std.algorithm.mutation;
-import std.math;
-import std.range;
-import std.traits;
-import tanya.memory;
-
-/**
- * Mutliple precision integer.
- */
-struct Integer
-{
- package ubyte[] rep;
- private bool sign;
-
- /**
- * Creates a multiple precision integer.
- *
- * Params:
- * T = Value type.
- * value = Initial value.
- * allocator = Allocator.
- *
- * Precondition: $(D_INLINECODE allocator !is null)
- */
- this(T)(in auto ref T value, shared Allocator allocator = defaultAllocator)
- nothrow @safe @nogc
- if (isIntegral!T || is(T == Integer))
- {
- this(allocator);
- static if (isIntegral!T)
- {
- assignInt(value);
- }
- else if (value.length > 0)
- {
- allocator.resize!(ubyte, false)(rep, value.length);
- rep[] = value.rep[];
- sign = value.sign;
- }
- }
-
- /// Ditto.
- this(shared Allocator allocator) pure nothrow @safe @nogc
- in
- {
- assert(allocator !is null);
- }
- body
- {
- allocator_ = allocator;
- }
-
- private @nogc unittest
- {
- auto h1 = Integer(79);
- assert(h1.length == 1);
- assert(h1.rep[0] == 79);
- assert(!h1.sign);
-
- auto h2 = Integer(-2);
- assert(h2.length == 1);
- assert(h2.rep[0] == 2);
- assert(h2.sign);
- }
-
- @disable this(this);
-
- /**
- * Destroys the internal representation.
- */
- ~this() nothrow @safe @nogc
- {
- allocator.dispose(rep);
- }
-
- /*
- * Figures out the minimum amount of space this value will take
- * up in bytes and resizes the internal storage. Sets the sign.
- */
- private void assignInt(T)(in ref T value)
- nothrow @safe @nogc
- in
- {
- static assert(isIntegral!T);
- }
- body
- {
- ubyte size = ulong.sizeof;
- ulong mask;
-
- static if (isSigned!T)
- {
- sign = value < 0 ? true : false;
- immutable absolute = value.abs;
- }
- else
- {
- sign = false;
- alias absolute = value;
- }
- for (mask = 0xff00000000000000; mask >= 0xff; mask >>= 8)
- {
- if (absolute & mask)
- {
- break;
- }
- --size;
- }
- allocator.resize!(ubyte, false)(rep, size);
-
- /* 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;
- ushort shift;
- for (auto i = rep.length; i; --i, mask <<= 8, shift += 8)
- {
- rep[i - 1] = cast(ubyte) ((absolute & mask) >> shift);
- }
- }
-
- /**
- * Assigns a new value.
- *
- * Params:
- * T = Value type.
- * value = Value.
- *
- * Returns: $(D_KEYWORD this).
- */
- ref Integer opAssign(T)(in auto ref T value) nothrow @safe @nogc
- if (isIntegral!T || is(T == Integer))
- {
- static if (isIntegral!T)
- {
- assignInt(value);
- }
- else
- {
- allocator.resize!(ubyte, false)(rep, value.length);
- rep[0 .. $] = value.rep[0 .. $];
- sign = value.sign;
- }
- return this;
- }
-
- private @nogc unittest
- {
- auto h = Integer(1019);
- assert(h.length == 2);
- assert(h.rep[0] == 0b00000011 && h.rep[1] == 0b11111011);
-
- h = 3337;
- assert(h.length == 2);
- assert(h.rep[0] == 0b00001101 && h.rep[1] == 0b00001001);
-
- h = 688;
- assert(h.length == 2);
- assert(h.rep[0] == 0b00000010 && h.rep[1] == 0b10110000);
-
- h = 0;
- assert(h.length == 0);
- }
-
- /**
- * Returns: Integer size.
- */
- @property size_t length() const pure nothrow @safe @nogc
- {
- return rep.length;
- }
-
- /**
- * Params:
- * h = The second integer.
- *
- * Returns: Whether the two integers are equal.
- */
- bool opEquals(in Integer h) const nothrow @safe @nogc
- {
- return rep == h.rep;
- }
-
- /// Ditto.
- bool opEquals(in ref Integer h) const nothrow @safe @nogc
- {
- return rep == h.rep;
- }
-
- ///
- unittest
- {
- auto h1 = Integer(1019);
-
- assert(h1 == Integer(1019));
- assert(h1 != Integer(109));
- }
-
- /**
- * Params:
- * h = The second integer.
- *
- * Returns: A positive number if $(D_INLINECODE this > h), a negative
- * number if $(D_INLINECODE this > h), `0` otherwise.
- */
- int opCmp(in ref Integer h) const nothrow @safe @nogc
- {
- if (length > h.length)
- {
- return 1;
- }
- if (length < h.length)
- {
- return -1;
- }
- // Otherwise, keep searching through the representational integers
- // until one is bigger than another - once we've found one, it's
- // safe to stop, since the lower order bytes can't affect the
- // comparison
- for (size_t i, j; i < length && j < h.length; ++i, ++j)
- {
- if (rep[i] < h.rep[j])
- {
- return -1;
- }
- else if (rep[i] > h.rep[j])
- {
- return 1;
- }
- }
- // if we got all the way to the end without a comparison, the
- // two are equal
- return 0;
- }
-
- /// Ditto.
- int opCmp(in Integer h) const nothrow @safe @nogc
- {
- return opCmp(h);
- }
-
- ///
- unittest
- {
- auto h1 = Integer(1019);
- auto h2 = Integer(1019);
- assert(h1 == h2);
-
- h2 = 3337;
- assert(h1 < h2);
-
- h2 = 688;
- assert(h1 > h2);
- }
-
- private void add(in ref ubyte[] h) nothrow @safe @nogc
- {
- uint sum;
- uint carry = 0;
- ubyte[] tmp;
-
- if (h.length > length)
- {
- allocator.resize!(ubyte, false)(tmp, h.length);
- tmp[0 .. h.length] = 0;
- tmp[h.length - length .. $] = rep[0 .. length];
- swap(rep, tmp);
- }
-
- auto i = length;
- auto j = h.length;
-
- do
- {
- --i;
- if (j)
- {
- --j;
- sum = rep[i] + h[j] + carry;
- }
- else
- {
- sum = rep[i] + carry;
- }
- carry = sum > 0xff;
- rep[i] = cast(ubyte) sum;
- }
- while (i);
-
- if (carry)
- {
- // Still overflowed; allocate more space
- allocator.resize!(ubyte, false)(tmp, length + 1);
- tmp[1 .. $] = rep[0 .. length];
- tmp[0] = 0x01;
- swap(rep, tmp);
- }
- allocator.dispose(tmp);
- }
-
- private void subtract(in ref ubyte[] h) nothrow @trusted @nogc
- {
- auto i = rep.length;
- auto j = h.length;
- uint borrow = 0;
-
- do
- {
- int difference;
- --i;
-
- if (j)
- {
- --j;
- difference = rep[i] - h[j] - borrow;
- }
- else
- {
- difference = rep[i] - borrow;
- }
- borrow = difference < 0;
- rep[i] = cast(ubyte) difference;
- }
- while (i);
-
- if (borrow && i && rep[i - 1])
- {
- --rep[i - 1];
- }
- // Go through the representation array and see how many of the
- // left-most bytes are unused. Remove them and resize the array.
- immutable offset = rep.countUntil!((const ref a) => a != 0);
- if (offset > 0)
- {
- ubyte[] tmp = cast(ubyte[]) allocator.allocate(rep.length - offset);
- tmp[0 .. $] = rep[offset .. $];
- allocator.deallocate(rep);
- rep = tmp;
- }
- else if (offset == -1)
- {
- allocator.dispose(rep);
- }
- }
-
- /**
- * Assignment operators with another $(D_PSYMBOL Integer).
- *
- * Params:
- * op = Operation.
- * h = The second integer.
- *
- * Returns: $(D_KEYWORD this).
- */
- ref Integer opOpAssign(string op)(in auto ref Integer h) nothrow @safe @nogc
- if ((op == "+") || (op == "-"))
- out
- {
- assert(rep.length || !sign, "0 should be positive.");
- assert(!rep.length || rep[0]);
- }
- body
- {
- static if (op == "+")
- {
- if (h.sign == sign)
- {
- add(h.rep);
- }
- else
- {
- if (this >= h)
- {
- subtract(h.rep);
- }
- else
- {
- auto tmp = Integer(h);
- tmp.subtract(rep);
- swap(rep, tmp.rep);
- sign = length == 0 ? false : h.sign;
- }
- }
- }
- else
- {
- if (h.sign == sign)
- {
- if (this >= h)
- {
- subtract(h.rep);
- }
- else
- {
- auto tmp = Integer(h);
- tmp.subtract(rep);
- swap(rep, tmp.rep);
- sign = length == 0 ? false : !sign;
- }
- }
- else
- {
- add(h.rep);
- }
- }
- return this;
- }
-
- private @nogc unittest
- {
- {
- auto h1 = Integer(1019);
- auto h2 = Integer(3337);
- h1 += h2;
- assert(h1.length == 2);
- assert(h1.rep[0] == 0x11 && h1.rep[1] == 0x04);
- }
- {
- auto h1 = Integer(4356);
- auto h2 = Integer(2_147_483_647);
- ubyte[4] expected = [0x80, 0x00, 0x11, 0x03];
- h1 += h2;
- assert(h1.rep == expected);
- }
- {
- auto h1 = Integer(2147488003L);
- auto h2 = Integer(2_147_483_647);
- ubyte[5] expected = [0x01, 0x00, 0x00, 0x11, 0x02];
- h1 += h2;
- assert(h1.rep == expected);
- }
- {
- auto h1 = Integer(3);
- auto h2 = Integer(4);
- h1 -= h2;
- assert(h1.length == 1);
- assert(h1.rep[0] == 0x01);
- assert(h1.sign);
- }
- }
-
- private @nogc unittest
- {
- {
- auto h1 = Integer(8589934590L);
- auto h2 = Integer(2147483647);
- ubyte[5] expected = [0x01, 0x7f, 0xff, 0xff, 0xff];
-
- h1 -= h2;
- assert(h1.rep == expected);
- }
- {
- auto h1 = Integer(6442450943);
- auto h2 = Integer(4294967294);
- ubyte[4] expected = [0x80, 0x00, 0x00, 0x01];
- h1 -= h2;
- assert(h1.rep == expected);
- }
- {
- auto h1 = Integer(2147483649);
- auto h2 = Integer(h1);
- h1 -= h2;
- assert(h1.length == 0);
- }
- }
-
- /// Ditto.
- ref Integer opOpAssign(string op)(in auto ref Integer h) nothrow @safe @nogc
- if (op == "*")
- out
- {
- assert(rep.length || !sign, "0 should be positive.");
- assert(!rep.length || rep[0]);
- }
- body
- {
- auto i = h.rep.length;
- if (length == 0)
- {
- return this;
- }
- else if (i == 0)
- {
- opAssign(0);
- return this;
- }
- auto temp = Integer(this, allocator);
- immutable sign = sign == h.sign ? false : true;
-
- opAssign(0);
- do
- {
- --i;
- for (ubyte mask = 0x01; mask; mask <<= 1)
- {
- if (mask & h.rep[i])
- {
- opOpAssign!"+"(temp);
- }
- temp <<= 1;
- }
- }
- while (i);
- this.sign = sign;
-
- return this;
- }
-
- ///
- unittest
- {
- auto h1 = Integer(123);
- auto h2 = Integer(456);
- h1 *= h2;
- assert(cast(long) h1 == 56088);
- }
-
- private unittest
- {
- assert((Integer(1) * Integer()).length == 0);
- }
-
- /// Ditto.
- ref Integer opOpAssign(string op)(in auto ref Integer h) nothrow @safe @nogc
- if ((op == "/") || (op == "%"))
- in
- {
- assert(h.length > 0, "Division by zero.");
- }
- body
- {
- auto divisor = Integer(h, allocator);
- size_t bitSize;
-
- // First, left-shift divisor until it's >= than the dividend
- for (; opCmp(divisor) > 0; ++bitSize)
- {
- divisor <<= 1;
- }
- static if (op == "/")
- {
- ubyte[] quotient;
- allocator.resize!(ubyte, false)(quotient, bitSize / 8 + 1);
- }
-
- // "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 (opCmp(divisor) >= 0)
- {
- opOpAssign!"-"(divisor);
- static if (op == "/")
- {
- quotient[bitPosition / 8] |= (0x80 >> (bitPosition % 8));
- }
- }
- if (bitSize)
- {
- divisor >>= 1;
- }
- else
- {
- break;
- }
- ++bitPosition, --bitSize;
- }
- while (true);
-
- static if (op == "/")
- {
- allocator.dispose(rep);
- rep = quotient;
- sign = sign == h.sign ? false : true;
- }
- return this;
- }
-
- private @nogc unittest
- {
- auto h1 = Integer(18);
- auto h2 = Integer(4);
- h1 %= h2;
- assert(h1.length == 1);
- assert(h1.rep[0] == 0x02);
-
- h1 = 8;
- h1 %= h2;
- assert(h1.length == 0);
-
- h1 = 7;
- h1 %= h2;
- assert(h1.length == 1);
- assert(h1.rep[0] == 0x03);
-
- h1 = 56088;
- h2 = 456;
- h1 /= h2;
- assert(h1.length == 1);
- assert(h1.rep[0] == 0x7b);
- }
-
- /// Ditto.
- ref Integer opOpAssign(string op)(in auto ref Integer exp) nothrow @safe @nogc
- if (op == "^^")
- out
- {
- assert(rep.length || !sign, "0 should be positive.");
- assert(!rep.length || rep[0]);
- }
- body
- {
- auto i = exp.rep.length;
- auto tmp1 = Integer(this, allocator);
- auto tmp2 = Integer(allocator);
-
- opAssign(1);
-
- do
- {
- --i;
- for (ubyte mask = 0x01; mask; mask <<= 1)
- {
- if (exp.rep[i] & mask)
- {
- opOpAssign!"*"(tmp1);
- }
- // Square tmp1
- tmp2 = tmp1;
- tmp1 *= tmp2;
- }
- }
- while (i);
-
- 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 = [0x1b, 0x5c, 0xab, 0x9c, 0x31, 0x10];
- assert(h1.rep == expected);
- }
-
- /**
- * Unary operators.
- *
- * Params:
- * op = Operation.
- *
- * Returns: New $(D_PSYMBOL Integer).
- */
- Integer opUnary(string op)() nothrow @safe @nogc
- if ((op == "+") || (op == "-") || (op == "~"))
- {
- auto h = Integer(this, allocator);
- static if (op == "-")
- {
- h.sign = !h.sign;
- }
- else static if (op == "~")
- {
- h.rep.each!((ref a) => a = ~a);
- }
- return h;
- }
-
- private @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);
-
- h2 = -h1;
- assert(h2.length == 1);
- assert(h2.rep[0] == 79);
- assert(h2.sign);
-
- 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);
-
- h2 = ~h1;
- assert(h2.rep[0] == ~cast(ubyte) 79);
- }
-
- private void decrement() nothrow @safe @nogc
- {
- immutable size = rep.retro.countUntil!((const ref a) => a != 0);
- if (rep[0] == 1)
- {
- allocator.resize!(ubyte, false)(rep, rep.length - 1);
- rep[0 .. $] = typeof(rep[0]).max;
- }
- else
- {
- --rep[$ - size - 1];
- rep[$ - size .. $] = typeof(rep[0]).max;
- }
- }
-
- private void increment() nothrow @safe @nogc
- {
- auto size = rep.retro.countUntil!((const ref a) => a != typeof(rep[0]).max);
- if (size == -1)
- {
- size = length;
- allocator.resize!(ubyte, false)(rep, rep.length + 1);
- rep[0] = 1;
- }
- else
- {
- ++rep[$ - size - 1];
- }
- rep[$ - size .. $] = 0;
- }
-
- /**
- * Increment/decrement.
- *
- * Params:
- * op = Operation.
- *
- * Returns: $(D_KEYWORD this).
- */
- ref Integer opUnary(string op)() nothrow @safe @nogc
- if ((op == "++") || (op == "--"))
- out
- {
- assert(rep.length || !sign, "0 should be positive.");
- assert(!rep.length || rep[0]);
- }
- body
- {
- static if (op == "++")
- {
- if (sign)
- {
- decrement();
- if (length == 0)
- {
- sign = false;
- }
- }
- else
- {
- increment();
- }
- }
- else if (sign)
- {
- increment();
- }
- else
- {
- decrement();
- }
- return this;
- }
-
- private @nogc unittest
- {
- Integer h;
-
- ++h;
- assert(h.length == 1);
- assert(h.rep[0] == 0x01);
-
- --h;
- assert(h.length == 0);
-
- h = 511;
- ++h;
- assert(h.length == 2);
- assert(h.rep[0] == 0x02 && h.rep[1] == 0x00);
-
- --h;
- assert(h.length == 2);
- assert(h.rep[0] == 0x01 && h.rep[1] == 0xff);
-
- h = 79;
- ++h;
- assert(h.length == 1);
- assert(h.rep[0] == 0x50);
-
- --h;
- assert(h.length == 1);
- assert(h.rep[0] == 0x4f);
-
- h = 65535;
- ++h;
- assert(h.length == 3);
- assert(h.rep[0] == 0x01 && h.rep[1] == 0x00 && h.rep[2] == 0x00);
-
- --h;
- assert(h.length == 2);
- assert(h.rep[0] == 0xff && h.rep[1] == 0xff);
-
- h = -2;
- ++h;
- assert(h.length == 1);
- assert(h.rep[0] == 0x01);
- }
-
- /**
- * 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 pure nothrow @safe @nogc
- {
- return length == 0 ? false : true;
- }
-
- /// Ditto.
- T opCast(T)() const pure nothrow @safe @nogc
- if (isIntegral!T && isSigned!T)
- {
- ulong ret;
- for (size_t i = length, j; i > 0 && j <= T.sizeof * 4; --i, j += 8)
- {
- ret |= cast(T) (rep[i - 1]) << j;
- }
- return sign ? -ret : ret;
- }
-
- /// Ditto.
- T opCast(T)() const pure nothrow @safe @nogc
- if (isIntegral!T && isUnsigned!T)
- {
- ulong ret;
- for (size_t i = length, j; i > 0 && j <= T.sizeof * 8; --i, j += 8)
- {
- ret |= cast(T) (rep[i - 1]) << j;
- }
- return ret;
- }
-
- ///
- unittest
- {
- auto h = Integer(79);
- assert(cast(long) h == 79);
-
- h = -79;
- assert(cast(long) h == -79);
-
- h = 4294967295;
- assert(cast(long) h == 4294967295);
-
- h = -4294967295;
- assert(cast(long) h == -4294967295);
- }
-
- /**
- * Shift operations.
- *
- * Params:
- * op = Left or right shift.
- * n = Number of bits to shift by.
- *
- * Returns: An $(D_PSYMBOL Integer) shifted by $(D_PARAM n) bits.
- */
- ref Integer opOpAssign(string op)(in auto ref size_t n) nothrow @safe @nogc
- if (op == ">>")
- out
- {
- assert(rep.length || !sign, "0 should be positive.");
- assert(!rep.length || rep[0]);
- }
- body
- {
- immutable step = n / 8;
-
- if (step >= rep.length)
- {
- allocator.resize!(ubyte, false)(rep, 0);
- return this;
- }
-
- size_t i, j;
- ubyte carry;
- immutable bit = n % 8;
- immutable delta = 8 - bit;
-
- carry = cast(ubyte) (rep[0] << delta);
- rep[0] = (rep[0] >> bit);
- if (rep[0])
- {
- ++j;
- }
- for (i = 1; i < rep.length; ++i)
- {
- immutable oldCarry = carry;
- carry = cast(ubyte) (rep[i] << delta);
- rep[j] = (rep[i] >> bit | oldCarry);
- ++j;
- }
- allocator.resize!(ubyte, false)(rep, rep.length - n / 8 - (i == j ? 0 : 1));
-
- return this;
- }
-
- private @nogc unittest
- {
- auto h1 = Integer(4294967294);
- h1 >>= 10;
- assert(h1.length == 3);
- assert(h1.rep[0] == 0x3f && h1.rep[1] == 0xff && h1.rep[2] == 0xff);
-
- h1 = 27336704;
- h1 >>= 1;
- assert(h1.length == 3);
- assert(h1.rep[0] == 0xd0 && h1.rep[1] == 0x90 && h1.rep[2] == 0x00);
-
- h1 = 4294967294;
- h1 >>= 20;
- assert(h1.length == 2);
- assert(h1.rep[0] == 0x0f && h1.rep[1] == 0xff);
-
- h1 >>= 0;
- assert(h1.length == 2);
- assert(h1.rep[0] == 0x0f && h1.rep[1] == 0xff);
-
- h1 >>= 20;
- assert(h1.length == 0);
-
- h1 >>= 2;
- assert(h1.length == 0);
-
- h1 = 1431655765;
- h1 >>= 16;
- assert(h1.length == 2);
- assert(h1.rep[0] == 0x55 && h1.rep[1] == 0x55);
-
- h1 >>= 16;
- assert(h1.length == 0);
- }
-
- /// Ditto.
- ref Integer opOpAssign(string op)(in auto ref size_t n) nothrow @safe @nogc
- if (op == "<<")
- out
- {
- assert(rep.length || !sign, "0 should be positive.");
- assert(!rep.length || rep[0]);
- }
- body
- {
- ubyte carry;
- auto i = rep.length;
- size_t j;
- immutable bit = n % 8;
- immutable delta = 8 - bit;
-
- if (cast(ubyte) (rep[0] >> delta))
- {
- allocator.resize!(ubyte, false)(rep, i + n / 8 + 1);
- j = i + 1;
- }
- else
- {
- allocator.resize!(ubyte, false)(rep, i + n / 8);
- j = i;
- }
- do
- {
- --i, --j;
- immutable oldCarry = carry;
- carry = rep[i] >> delta;
- rep[j] = cast(ubyte) ((rep[i] << bit) | oldCarry);
- }
- while (i);
- if (carry)
- {
- rep[0] = carry;
- }
- return this;
- }
-
- private @nogc unittest
- {
- auto h1 = Integer(4294967295);
- ubyte[5] expected = [0x01, 0xff, 0xff, 0xff, 0xfe];
- h1 <<= 1;
- assert(h1.rep == expected);
- }
-
- /// Ditto.
- Integer opBinary(string op)(in auto ref size_t n) nothrow @safe @nogc
- if (op == "<<" || op == ">>" || op == "+" || op == "-" || op == "/"
- || op == "*" || op == "^^" || op == "%")
- {
- auto ret = Integer(this, allocator);
- mixin("ret " ~ op ~ "= n;");
- return ret;
- }
-
- ///
- unittest
- {
- auto h1 = Integer(425);
- auto h2 = h1 << 1;
- assert(cast(long) h2 == 850);
-
- h2 = h1 >> 1;
- assert(cast(long) h2 == 212);
- }
-
- /// Ditto.
- Integer opBinary(string op)(in auto ref Integer h) nothrow @safe @nogc
- if (op == "+" || op == "-" || op == "/"
- || op == "*" || op == "^^" || op == "%")
- {
- auto ret = Integer(this, allocator);
- mixin("ret " ~ op ~ "= h;");
- return ret;
- }
-
- mixin DefaultAllocator;
-}
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
+
+/**
+ * Arbitrary precision arithmetic.
+ *
+ * Copyright: Eugene Wissner 2016-2017.
+ * License: $(LINK2 https://www.mozilla.org/en-US/MPL/2.0/,
+ * Mozilla Public License, v. 2.0).
+ * Authors: $(LINK2 mailto:info@caraus.de, Eugene Wissner)
+ */
+module tanya.math.mp;
+
+import core.exception;
+import std.algorithm;
+import std.range;
+import std.traits;
+import tanya.math;
+import tanya.memory;
+
+/**
+ * Algebraic sign.
+ */
+enum Sign : bool
+{
+ /// The value is positive or `0`.
+ positive = false,
+
+ /// The value is negative.
+ negative = true,
+}
+
+/**
+ * Mutliple precision integer.
+ */
+struct Integer
+{
+ private size_t size;
+ package ubyte* rep;
+ package Sign sign;
+
+ pure nothrow @safe @nogc invariant
+ {
+ assert(this.size > 0 || !this.sign, "0 should be positive.");
+ }
+
+ /**
+ * Creates a multiple precision integer.
+ *
+ * Params:
+ * T = Value type.
+ * value = Initial value.
+ * allocator = Allocator.
+ *
+ * Precondition: $(D_INLINECODE allocator !is null)
+ */
+ this(T)(const T value, shared Allocator allocator = defaultAllocator)
+ if (isIntegral!T)
+ {
+ this(allocator);
+ assign(value);
+ }
+
+ /// Ditto.
+ this(const ref Integer value, shared Allocator allocator = defaultAllocator)
+ nothrow @safe @nogc
+ {
+ this(allocator);
+ assign(value);
+ }
+
+ /// Ditto.
+ this(Integer value, shared Allocator allocator = defaultAllocator)
+ nothrow @safe @nogc
+ {
+ this(allocator);
+ if (allocator is value.allocator)
+ {
+ this.rep = value.rep;
+ this.size = value.size;
+ this.sign = value.sign;
+ value.rep = null;
+ value.size = 0;
+ value.sign = Sign.positive;
+ }
+ else
+ {
+ assign(value);
+ }
+ }
+
+ /// Ditto.
+ this(shared Allocator allocator) pure nothrow @safe @nogc
+ in
+ {
+ assert(allocator !is null);
+ }
+ body
+ {
+ 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.
+ *
+ * Params:
+ * R = Range type.
+ * sign = Sign.
+ * value = Range.
+ * allocator = Allocator.
+ *
+ * Precondition: $(D_INLINECODE allocator !is null)
+ */
+ this(R)(const Sign sign, R value, shared Allocator allocator = defaultAllocator)
+ if (isInputRange!R && hasLength!R && is(Unqual!(ElementType!R) == ubyte))
+ out
+ {
+ assert(this.size == 0 || this.rep[0] > 0);
+ }
+ body
+ {
+ this(allocator);
+ while (!value.empty && value.front == 0)
+ {
+ value.popFront();
+ }
+ 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]);
+ }
+
+ /**
+ * Copies the integer.
+ */
+ this(this) nothrow @trusted @nogc
+ {
+ auto tmp = allocator.resize!ubyte(null, this.size);
+ this.rep[0 .. this.size].copy(tmp);
+ this.rep = tmp.ptr;
+ }
+
+ /**
+ * Destroys the 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 (auto i = this.size; i; --i, mask <<= 8, shift += 8)
+ {
+ this.rep[i - 1] = 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]);
+ }
+
+ /**
+ * Returns: Integer size.
+ */
+ @property size_t length() const pure nothrow @safe @nogc
+ {
+ return this.size;
+ }
+
+ /**
+ * Assigns a new value.
+ *
+ * Params:
+ * T = Value type.
+ * value = Value.
+ *
+ * Returns: $(D_KEYWORD this).
+ */
+ ref Integer opAssign(T)(const T value)
+ if (isIntegral!T)
+ {
+ assign(value);
+ return this;
+ }
+
+ /// Ditto.
+ ref Integer opAssign(const ref Integer value) nothrow @safe @nogc
+ {
+ assign(value);
+ return this;
+ }
+
+ /// Ditto.
+ ref Integer opAssign(Integer value) nothrow @safe @nogc
+ {
+ swap(this.rep, value.rep);
+ swap(this.sign, value.sign);
+ swap(this.size, value.size);
+ return this;
+ }
+
+ private @nogc unittest
+ {
+ auto integer = Integer(1019);
+ assert(integer.length == 2);
+ assert(integer.rep[0] == 0b00000011 && integer.rep[1] == 0b11111011);
+
+ integer = 3337;
+ assert(integer.length == 2);
+ assert(integer.rep[0] == 0b00001101 && integer.rep[1] == 0b00001001);
+
+ integer = 688;
+ assert(integer.length == 2);
+ assert(integer.rep[0] == 0b00000010 && integer.rep[1] == 0b10110000);
+
+ integer = 0;
+ assert(integer.length == 0);
+ }
+
+ /*
+ * Extend the space for this by 1 byte and set the LSB to 1.
+ */
+ private void expand() nothrow @trusted @nogc
+ {
+ rep = allocator.resize(this.rep[0 .. this.size], this.size + 1).ptr;
+ ++this.size;
+ auto target = this.rep[1 .. this.size].retro;
+ this.rep[0 .. this.size - 1].retro.copy(target);
+ this.rep[0] = 0x01;
+ }
+
+ /*
+ * 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].countUntil!((const ref a) => a != 0);
+
+ if (i > 0)
+ {
+ this.rep[i .. this.size].copy(this.rep[0 .. this.size - i]);
+ this.rep = allocator.resize(this.rep[0 .. this.size], this.size - i).ptr;
+ this.size -= i;
+ }
+ else if (i == -1)
+ {
+ 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
+ {
+ if (summand.length > this.length)
+ {
+ auto tmp = this.rep[0 .. this.size];
+ this.rep = allocator.resize!ubyte(null, summand.size).ptr;
+ tmp.copy(this.rep[summand.size - this.size .. summand.size]);
+ allocator.resize(tmp[0 .. this.size], 0);
+ this.size = summand.size;
+ }
+
+ bool carry;
+ size_t i = this.size;
+ size_t j = summand.size;
+ do
+ {
+ uint sum;
+ --i;
+ if (j)
+ {
+ --j;
+ sum = this.rep[i] + summand.rep[j] + carry;
+ }
+ else
+ {
+ sum = this.rep[i] + carry;
+ }
+
+ carry = sum > 0xff;
+ this.rep[i] = cast(ubyte) sum;
+ }
+ while (i);
+
+ if (carry)
+ {
+ // Still overflowed; allocate more space
+ expand();
+ }
+ }
+
+ private void subtract(const ref Integer subtrahend) nothrow @trusted @nogc
+ {
+ size_t i = this.size;
+ size_t j = subtrahend.size;
+ bool borrow;
+
+ while (i)
+ {
+ int difference;
+ --i;
+
+ if (j)
+ {
+ --j;
+ difference = this.rep[i] - subtrahend.rep[j] - borrow;
+ }
+ else
+ {
+ difference = this.rep[i] - borrow;
+ }
+ borrow = difference < 0;
+ this.rep[i] = cast(ubyte) difference;
+ }
+
+ if (borrow && i > 0)
+ {
+ if (this.rep[i - 1]) // Don't borrow i
+ {
+ --this.rep[i - 1];
+ }
+ }
+ contract();
+ }
+
+ private int compare(const ref Integer that) const nothrow @trusted @nogc
+ {
+ if (this.length > that.length)
+ {
+ return 1;
+ }
+ else if (this.length < that.length)
+ {
+ return -1;
+ }
+ return this.rep[0 .. this.size].cmp(that.rep[0 .. that.size]);
+ }
+
+ /**
+ * Comparison.
+ *
+ * Params:
+ * 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
+ {
+ if (this.sign != that.sign)
+ {
+ return this.sign ? -1 : 1;
+ }
+ return compare(that);
+ }
+
+ ///
+ @safe @nogc unittest
+ {
+ auto integer1 = Integer(1019);
+ auto integer2 = Integer(1019);
+ assert(integer1 == integer2);
+
+ integer2 = 3337;
+ assert(integer1 < integer2);
+
+ integer2 = 688;
+ assert(integer1 > integer2);
+
+ integer2 = -3337;
+ assert(integer1 > integer2);
+ }
+
+ /// Ditto.
+ int opCmp(I)(const auto ref I that) const
+ if (isIntegral!I)
+ {
+ if (that < 0 && !this.sign)
+ {
+ return 1;
+ }
+ else if (that > 0 && this.sign)
+ {
+ return -1;
+ }
+ else if (this.length > I.sizeof)
+ {
+ return this.sign ? -1 : 1;
+ }
+
+ const diff = (cast(I) this) - that;
+ if (diff > 0)
+ {
+ return 1;
+ }
+ else if (diff < 0)
+ {
+ return -1;
+ }
+ return 0;
+ }
+
+ ///
+ @safe @nogc unittest
+ {
+ auto integer = Integer(1019);
+
+ assert(integer == 1019);
+ assert(integer < 3337);
+ assert(integer > 688);
+ assert(integer > -3337);
+ }
+
+ /**
+ * Params:
+ * that = The second integer.
+ *
+ * Returns: Whether the two integers are equal.
+ */
+ bool opEquals(I)(const auto ref I that) const
+ if (is(I : Integer) || isIntegral!I)
+ {
+ return opCmp!I(that) == 0;
+ }
+
+ ///
+ @safe @nogc unittest
+ {
+ auto integer = Integer(1019);
+
+ assert(integer == Integer(1019));
+ assert(integer != Integer(109));
+ }
+
+ /**
+ * Assignment operators with another $(D_PSYMBOL Integer).
+ *
+ * Params:
+ * op = Operation.
+ * operand = The second operand.
+ *
+ * Returns: $(D_KEYWORD this).
+ */
+ ref Integer opOpAssign(string op : "+")(const auto ref Integer operand)
+ out
+ {
+ assert(this.size == 0 || this.rep[0] > 0);
+ }
+ body
+ {
+ if (this.sign == operand.sign)
+ {
+ add(operand);
+ }
+ 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;
+ }
+
+ }
+ return this;
+ }
+
+ private @nogc unittest
+ {
+ {
+ auto h1 = Integer(1019);
+ auto h2 = Integer(3337);
+ h1 += h2;
+ assert(h1.length == 2);
+ assert(h1.rep[0] == 0x11 && h1.rep[1] == 0x04);
+ }
+ {
+ auto h1 = Integer(4356);
+ auto h2 = Integer(2_147_483_647);
+ ubyte[4] expected = [0x80, 0x00, 0x11, 0x03];
+ h1 += h2;
+ assert(h1.rep[0 .. h1.length] == expected);
+ }
+ {
+ auto h1 = Integer(2147488003L);
+ auto h2 = Integer(2_147_483_647);
+ ubyte[5] expected = [0x01, 0x00, 0x00, 0x11, 0x02];
+ h1 += h2;
+ assert(h1.rep[0 .. h1.length] == expected);
+ }
+ }
+
+ /// Ditto.
+ ref Integer opOpAssign(string op : "-")(const auto ref Integer operand)
+ out
+ {
+ assert(this.size == 0 || this.rep[0] > 0);
+ }
+ body
+ {
+ if (operand.sign == this.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;
+ }
+ }
+ }
+ else
+ {
+ add(operand);
+ }
+ return this;
+ }
+
+ private @nogc unittest
+ {
+ {
+ auto h1 = Integer(3);
+ auto h2 = Integer(4);
+ h1 -= h2;
+ assert(h1.length == 1);
+ assert(h1.rep[0] == 0x01);
+ assert(h1.sign);
+ }
+ {
+ auto h1 = Integer(8589934590L);
+ auto h2 = Integer(2147483647);
+ ubyte[5] expected = [0x01, 0x7f, 0xff, 0xff, 0xff];
+
+ h1 -= h2;
+ assert(h1.rep[0 .. h1.size] == expected);
+ }
+ {
+ auto h1 = Integer(6442450943);
+ auto h2 = Integer(4294967294);
+ ubyte[4] expected = [0x80, 0x00, 0x00, 0x01];
+ h1 -= h2;
+ assert(h1.rep[0 .. h1.size] == expected);
+ }
+ {
+ auto h1 = Integer(2147483649);
+ auto h2 = Integer(h1);
+ h1 -= h2;
+ assert(h1.length == 0);
+ }
+ }
+
+ /// Ditto.
+ ref Integer opOpAssign(string op : "*")(const auto ref Integer operand) @trusted
+ out
+ {
+ assert(this.size == 0 || this.rep[0] > 0);
+ }
+ body
+ {
+ size_t i = operand.size;
+ if (this.length == 0)
+ {
+ return this;
+ }
+ else if (i == 0)
+ {
+ this = 0;
+ return this;
+ }
+ auto temp = Integer(this, this.allocator);
+ auto sign = this.sign != operand.sign;
+
+ this = 0;
+ do
+ {
+ --i;
+ for (ubyte mask = 0x01; mask; mask <<= 1)
+ {
+ if (mask & operand.rep[i])
+ {
+ add(temp);
+ }
+ temp <<= 1;
+ }
+ }
+ while (i);
+
+ this.sign = sign ? Sign.negative : Sign.positive;
+
+ return this;
+ }
+
+ ///
+ @safe @nogc unittest
+ {
+ auto h1 = Integer(123);
+ auto h2 = Integer(456);
+ h1 *= h2;
+ assert(h1 == 56088);
+ }
+
+ private @nogc unittest
+ {
+ assert((Integer(1) * Integer()).length == 0);
+ }
+
+ /// Ditto.
+ ref Integer opOpAssign(string op : "^^")(const auto ref Integer operand)
+ @trusted
+ out
+ {
+ assert(this.size == 0 || this.rep[0] > 0);
+ }
+ body
+ {
+ size_t i = operand.size;
+
+ auto tmp1 = Integer(this, this.allocator);
+ this = 1;
+
+ do
+ {
+ --i;
+ for (ubyte mask = 0x01; mask; mask <<= 1)
+ {
+ if (operand.rep[i] & mask)
+ {
+ this *= tmp1;
+ }
+ // Square tmp1
+ auto tmp2 = tmp1;
+ tmp1 *= tmp2;
+ }
+ }
+ while (i);
+
+ 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 = [0x1b, 0x5c, 0xab, 0x9c, 0x31, 0x10];
+ assert(h1.rep[0 .. h1.size] == expected);
+ }
+
+ /// Ditto.
+ ref Integer opOpAssign(string op)(const auto ref Integer operand) @trusted
+ if ((op == "%") || (op == "/"))
+ in
+ {
+ assert(operand.length > 0, "Division by zero.");
+ }
+ out
+ {
+ assert(this.size == 0 || this.rep[0] > 0);
+ }
+ 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);
+ }
+
+ // 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[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();
+ }
+
+ return this;
+ }
+
+ private @nogc unittest
+ {
+ auto h1 = Integer(18);
+ auto h2 = Integer(4);
+ h1 %= h2;
+ assert(h1.length == 1);
+ assert(h1.rep[0] == 0x02);
+
+ h1 = 8;
+ h1 %= h2;
+ assert(h1.length == 0);
+
+ h1 = 7;
+ h1 %= h2;
+ assert(h1.length == 1);
+ assert(h1.rep[0] == 0x03);
+
+ h1 = 56088;
+ h2 = 456;
+ h1 /= h2;
+ assert(h1.length == 1);
+ assert(h1.rep[0] == 0x7b);
+ }
+
+ /// Ditto.
+ ref Integer opOpAssign(string op : ">>")(const size_t operand) @trusted
+ out
+ {
+ assert(this.size == 0 || this.rep[0] > 0);
+ }
+ body
+ {
+ const step = operand / 8;
+
+ if (step >= this.length)
+ {
+ 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;
+
+ carry = cast(ubyte) (this.rep[0] << delta);
+ this.rep[0] = this.rep[0] >> bit;
+ if (rep[0])
+ {
+ ++j;
+ }
+ for (i = 1; i < this.length; ++i)
+ {
+ immutable oldCarry = carry;
+ carry = cast(ubyte) (this.rep[i] << delta);
+ this.rep[j] = this.rep[i] >> bit | oldCarry;
+ ++j;
+ }
+ const newSize = this.length - operand / 8 - (i == j ? 0 : 1);
+ this.rep = allocator.resize(this.rep[0 .. this.size], newSize).ptr;
+ this.size = newSize;
+
+ return this;
+ }
+
+ private @nogc unittest
+ {
+ auto integer = Integer(4294967294);
+ integer >>= 10;
+ assert(integer.length == 3);
+ assert(integer.rep[0] == 0x3f && integer.rep[1] == 0xff && integer.rep[2] == 0xff);
+
+ integer = 27336704;
+ integer >>= 1;
+ assert(integer.length == 3);
+ assert(integer.rep[0] == 0xd0 && integer.rep[1] == 0x90 && integer.rep[2] == 0x00);
+
+ integer = 4294967294;
+ integer >>= 20;
+ assert(integer.length == 2);
+ assert(integer.rep[0] == 0x0f && integer.rep[1] == 0xff);
+
+ integer >>= 0;
+ assert(integer.length == 2);
+ assert(integer.rep[0] == 0x0f && integer.rep[1] == 0xff);
+
+ integer >>= 20;
+ assert(integer.length == 0);
+
+ integer >>= 2;
+ assert(integer.length == 0);
+
+ integer = 1431655765;
+ integer >>= 16;
+ assert(integer.length == 2);
+ assert(integer.rep[0] == 0x55 && integer.rep[1] == 0x55);
+
+ integer >>= 16;
+ assert(integer.length == 0);
+ }
+
+ /// Ditto.
+ ref Integer opOpAssign(string op : "<<")(const size_t operand) @trusted
+ out
+ {
+ assert(this.size == 0 || this.rep[0] > 0);
+ }
+ body
+ {
+ ubyte carry;
+ auto i = this.length;
+ size_t j;
+ const bit = operand % 8;
+ const delta = 8 - bit;
+
+ if (cast(ubyte) (this.rep[0] >> delta))
+ {
+ const newSize = i + operand / 8 + 1;
+ this.rep = allocator.resize(this.rep[0 .. this.size], newSize).ptr;
+ this.size = newSize;
+ j = i + 1;
+ }
+ else
+ {
+ const newSize = i + operand / 8;
+ this.rep = allocator.resize(this.rep[0 .. this.size], newSize).ptr;
+ j = i;
+ }
+ do
+ {
+ --i, --j;
+ const oldCarry = carry;
+ carry = rep[i] >> delta;
+ rep[j] = cast(ubyte) ((this.rep[i] << bit) | oldCarry);
+ }
+ while (i);
+ if (carry)
+ {
+ rep[0] = carry;
+ }
+ return this;
+ }
+
+ private @nogc unittest
+ {
+ auto integer = Integer(4294967295);
+ ubyte[5] expected = [0x01, 0xff, 0xff, 0xff, 0xfe];
+ integer <<= 1;
+ assert(integer.rep[0 .. integer.size] == expected);
+ }
+
+ private void decrement() nothrow @trusted @nogc
+ in
+ {
+ assert(this.length > 0);
+ }
+ body
+ {
+ auto p = this.rep + this.size;
+ while (p != this.rep)
+ {
+ --p, --*p;
+ if (*p != ubyte.max)
+ {
+ break;
+ }
+ }
+ contract();
+ }
+
+ private void increment() nothrow @trusted @nogc
+ {
+ auto p = this.rep + this.size;
+ while (p != this.rep)
+ {
+ --p, ++*p;
+ if (*p > 0)
+ {
+ return;
+ }
+ }
+ if (p == this.rep)
+ {
+ expand();
+ }
+ }
+
+ /**
+ * Unary operators.
+ *
+ * Params:
+ * op = Operation.
+ *
+ * Returns: New $(D_PSYMBOL Integer).
+ */
+ Integer opUnary(string op : "~")() const
+ {
+ auto ret = Integer(this, this.allocator);
+ ret.rep[0 .. ret.size].each!((ref a) => a = ~a);
+ return ret;
+ }
+
+ /// Ditto.
+ Integer opUnary(string op : "-")() const
+ {
+ auto ret = Integer(this, this.allocator);
+ ret.sign = ret.sign ? Sign.positive : Sign.negative;
+ return ret;
+ }
+
+ /**
+ * Unary operators.
+ *
+ * Params:
+ * op = Operation.
+ *
+ * Returns: $(D_KEYWORD this).
+ */
+ ref inout(Integer) opUnary(string op : "+")() inout
+ {
+ return this;
+ }
+
+ private @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);
+
+ h2 = -h1;
+ assert(h2.length == 1);
+ assert(h2.rep[0] == 79);
+ assert(h2.sign);
+
+ 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);
+
+ h2 = ~h1;
+ assert(h2.rep[0] == ~cast(ubyte) 79);
+ }
+
+ /// Ditto.
+ ref Integer opUnary(string op : "++")()
+ out
+ {
+ assert(this.size == 0 || this.rep[0] > 0);
+ }
+ body
+ {
+ if (this.sign)
+ {
+ decrement();
+ if (this.length == 0)
+ {
+ this.sign = Sign.positive;
+ }
+ }
+ else
+ {
+ increment();
+ }
+ return this;
+ }
+
+ /// Ditto.
+ ref Integer opUnary(string op : "--")()
+ out
+ {
+ assert(this.size == 0 || this.rep[0] > 0);
+ }
+ body
+ {
+ if (this.sign)
+ {
+ increment();
+ }
+ else
+ {
+ decrement();
+ }
+ return this;
+ }
+
+ private @nogc unittest
+ {
+ Integer integer;
+
+ ++integer;
+ assert(integer.length == 1);
+ assert(integer.rep[0] == 0x01);
+
+ --integer;
+ assert(integer.length == 0);
+
+ integer = 511;
+ ++integer;
+ assert(integer.length == 2);
+ assert(integer.rep[0] == 0x02 && integer.rep[1] == 0x00);
+
+ --integer;
+ assert(integer.length == 2);
+ assert(integer.rep[0] == 0x01 && integer.rep[1] == 0xff);
+
+ integer = 79;
+ ++integer;
+ assert(integer.length == 1);
+ assert(integer.rep[0] == 0x50);
+
+ --integer;
+ assert(integer.length == 1);
+ assert(integer.rep[0] == 0x4f);
+
+ integer = 65535;
+ ++integer;
+ assert(integer.length == 3);
+ assert(integer.rep[0] == 0x01 && integer.rep[1] == 0x00 && integer.rep[2] == 0x00);
+
+ --integer;
+ assert(integer.length == 2);
+ assert(integer.rep[0] == 0xff && integer.rep[1] == 0xff);
+
+ integer = -2;
+ ++integer;
+ assert(integer.length == 1);
+ assert(integer.rep[0] == 0x01);
+ }
+
+ /**
+ * Implements binary operators.
+ *
+ * Params:
+ * op = Operation.
+ * operand = The second operand.
+ */
+ Integer opBinary(string op)(const auto ref Integer operand) const
+ if (op == "+" || op == "-" || op == "*" || op == "^^")
+ {
+ auto ret = Integer(this, this.allocator);
+ mixin("ret " ~ op ~ "= operand;");
+ return ret;
+ }
+
+ /// Ditto.
+ Integer opBinary(string op)(const auto ref Integer operand) const
+ if (op == "/" || op == "%")
+ in
+ {
+ assert(operand.length > 0, "Division by zero.");
+ }
+ body
+ {
+ auto ret = Integer(this, this.allocator);
+ mixin("ret " ~ op ~ "= operand;");
+ return ret;
+ }
+
+ /// Ditto.
+ Integer opBinary(string op)(const auto ref size_t operand) const
+ if (op == "<<" || op == ">>")
+ {
+ auto ret = Integer(this, this.allocator);
+ mixin("ret " ~ op ~ "= operand;");
+ return ret;
+ }
+
+ ///
+ @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;
+ }
+ T ret;
+ const(ubyte)* src = this.rep + this.size - 1;
+ for (ubyte shift; src >= this.rep && shift <= T.sizeof * 8; --src, shift += 8)
+ {
+ ret |= cast(T) (*src) << shift;
+ }
+ return ret;
+ }
+
+ ///
+ @safe @nogc unittest
+ {
+ auto integer = Integer(79);
+ assert(cast(long) integer == 79);
+
+ integer = -79;
+ assert(cast(long) integer == -79);
+
+ integer = 4294967295;
+ assert(cast(long) integer == 4294967295);
+
+ integer = -4294967295;
+ assert(cast(long) integer == -4294967295);
+ }
+
+ mixin DefaultAllocator;
+}
diff --git a/source/tanya/math/package.d b/source/tanya/math/package.d
index 8449c34..994d5c7 100644
--- a/source/tanya/math/package.d
+++ b/source/tanya/math/package.d
@@ -3,11 +3,11 @@
* file, You can obtain one at http://mozilla.org/MPL/2.0/. */
/**
- * Copyright: Eugene Wissner 2016.
+ * Copyright: Eugene Wissner 2016-2017.
* License: $(LINK2 https://www.mozilla.org/en-US/MPL/2.0/,
* Mozilla Public License, v. 2.0).
- * Authors: $(LINK2 mailto:belka@caraus.de, Eugene Wissner)
- */
+ * Authors: $(LINK2 mailto:info@caraus.de, Eugene Wissner)
+ */
module tanya.math;
import std.traits;
@@ -16,7 +16,7 @@ public import tanya.math.random;
version (unittest)
{
- import std.algorithm.iteration;
+ import std.algorithm.iteration;
}
/**
@@ -26,12 +26,12 @@ version (unittest)
* is used to allocate the result.
*
* Params:
- * I = Base type.
- * G = Exponent type.
- * H = Divisor type:
- * x = Base.
- * y = Exponent.
- * z = Divisor.
+ * I = Base type.
+ * G = Exponent type.
+ * H = Divisor type:
+ * x = Base.
+ * y = Exponent.
+ * z = Divisor.
*
* Returns: Reminder of the division of $(D_PARAM x) to the power $(D_PARAM y)
* by $(D_PARAM z).
@@ -39,134 +39,162 @@ version (unittest)
* Precondition: $(D_INLINECODE z > 0)
*/
H pow(I, G, H)(in auto ref I x, in auto ref G y, in auto ref H z)
- if (isIntegral!I && isIntegral!G && isIntegral!H)
+ if (isIntegral!I && isIntegral!G && isIntegral!H)
in
{
- assert(z > 0, "Division by zero.");
+ assert(z > 0, "Division by zero.");
}
body
{
- G mask = G.max / 2 + 1;
- H result;
-
- if (y == 0)
- {
- return 1 % z;
- }
- else if (y == 1)
- {
- return x % z;
- }
- do
- {
- immutable bit = y & mask;
- if (!result && bit)
- {
- result = x;
- continue;
- }
-
- result *= result;
- if (bit)
- {
- result *= x;
- }
- result %= z;
- }
- while (mask >>= 1);
-
- return result;
+ G mask = G.max / 2 + 1;
+ H result;
+
+ if (y == 0)
+ {
+ return 1 % z;
+ }
+ else if (y == 1)
+ {
+ return x % z;
+ }
+ do
+ {
+ immutable bit = y & mask;
+ if (!result && bit)
+ {
+ result = x;
+ continue;
+ }
+
+ result *= result;
+ if (bit)
+ {
+ result *= x;
+ }
+ result %= z;
+ }
+ while (mask >>= 1);
+
+ return result;
}
/// Ditto.
-I pow(I)(in auto ref I x, in auto ref I y, in auto ref I z)
- if (is(I == Integer))
+I pow(I)(const auto ref I x, const auto ref I y, const auto ref I z)
+ if (is(I == Integer))
in
{
- assert(z.length > 0, "Division by zero.");
+ assert(z.length > 0, "Division by zero.");
}
body
{
- size_t i = y.length;
- auto tmp2 = Integer(x.allocator), tmp1 = Integer(x, x.allocator);
- Integer result = Integer(x.allocator);
-
- if (x.length == 0 && i != 0)
- {
- i = 0;
- }
- else
- {
- result = 1;
- }
- while (i)
- {
- --i;
- for (ubyte mask = 0x01; mask; mask <<= 1)
- {
- if (y.rep[i] & mask)
- {
- result *= tmp1;
- result %= z;
- }
- tmp2 = tmp1;
- tmp1 *= tmp2;
- tmp1 %= z;
- }
- }
- return result;
+ size_t i = y.length;
+ auto tmp1 = Integer(x, x.allocator);
+ auto result = Integer(x.allocator);
+
+ if (x.length == 0 && i != 0)
+ {
+ i = 0;
+ }
+ else
+ {
+ result = 1;
+ }
+ while (i)
+ {
+ --i;
+ for (ubyte mask = 0x01; mask; mask <<= 1)
+ {
+ if (y.rep[i] & mask)
+ {
+ result *= tmp1;
+ result %= z;
+ }
+ auto tmp2 = tmp1;
+ tmp1 *= tmp2;
+ tmp1 %= z;
+ }
+ }
+ return result;
}
///
pure nothrow @safe @nogc unittest
{
- assert(pow(3, 5, 7) == 5);
- assert(pow(2, 2, 1) == 0);
- assert(pow(3, 3, 3) == 0);
- assert(pow(7, 4, 2) == 1);
- assert(pow(53, 0, 2) == 1);
- assert(pow(53, 1, 3) == 2);
- assert(pow(53, 2, 5) == 4);
- assert(pow(0, 0, 5) == 1);
- assert(pow(0, 5, 5) == 0);
+ assert(pow(3, 5, 7) == 5);
+ assert(pow(2, 2, 1) == 0);
+ assert(pow(3, 3, 3) == 0);
+ assert(pow(7, 4, 2) == 1);
+ assert(pow(53, 0, 2) == 1);
+ assert(pow(53, 1, 3) == 2);
+ assert(pow(53, 2, 5) == 4);
+ assert(pow(0, 0, 5) == 1);
+ assert(pow(0, 5, 5) == 0);
}
///
unittest
{
- assert(cast(long) pow(Integer(3), Integer(5), Integer(7)) == 5);
- assert(cast(long) pow(Integer(2), Integer(2), Integer(1)) == 0);
- assert(cast(long) pow(Integer(3), Integer(3), Integer(3)) == 0);
- assert(cast(long) pow(Integer(7), Integer(4), Integer(2)) == 1);
- assert(cast(long) pow(Integer(53), Integer(0), Integer(2)) == 1);
- assert(cast(long) pow(Integer(53), Integer(1), Integer(3)) == 2);
- assert(cast(long) pow(Integer(53), Integer(2), Integer(5)) == 4);
- assert(cast(long) pow(Integer(0), Integer(0), Integer(5)) == 1);
- assert(cast(long) pow(Integer(0), Integer(5), Integer(5)) == 0);
+ assert(pow(Integer(3), Integer(5), Integer(7)) == 5);
+ assert(pow(Integer(2), Integer(2), Integer(1)) == 0);
+ assert(pow(Integer(3), Integer(3), Integer(3)) == 0);
+ assert(pow(Integer(7), Integer(4), Integer(2)) == 1);
+ assert(pow(Integer(53), Integer(0), Integer(2)) == 1);
+ assert(pow(Integer(53), Integer(1), Integer(3)) == 2);
+ assert(pow(Integer(53), Integer(2), Integer(5)) == 4);
+ assert(pow(Integer(0), Integer(0), Integer(5)) == 1);
+ assert(pow(Integer(0), Integer(5), Integer(5)) == 0);
}
/**
* Checks if $(D_PARAM x) is a prime.
*
* Params:
- * x = The number should be checked.
+ * x = The number should be checked.
*
* Returns: $(D_KEYWORD true) if $(D_PARAM x) is a prime number,
* $(D_KEYWORD false) otherwise.
*/
bool isPseudoprime(ulong x) nothrow pure @safe @nogc
{
- return pow(2, x - 1, x) == 1;
+ return pow(2, x - 1, x) == 1;
}
///
unittest
{
- uint[30] known = [74623, 74653, 74687, 74699, 74707, 74713, 74717, 74719,
- 74843, 74747, 74759, 74761, 74771, 74779, 74797, 74821,
- 74827, 9973, 104729, 15485867, 49979693, 104395303,
- 593441861, 104729, 15485867, 49979693, 104395303,
- 593441861, 899809363, 982451653];
+ uint[30] known = [74623, 74653, 74687, 74699, 74707, 74713, 74717, 74719,
+ 74843, 74747, 74759, 74761, 74771, 74779, 74797, 74821,
+ 74827, 9973, 104729, 15485867, 49979693, 104395303,
+ 593441861, 104729, 15485867, 49979693, 104395303,
+ 593441861, 899809363, 982451653];
+
+ known.each!((ref x) => assert(isPseudoprime(x)));
+}
+
+/**
+ * Params:
+ * I = Value type.
+ * x = Value.
+ *
+ * Returns: The absolute value of a number.
+ */
+I abs(I : Integer)(const auto ref I x)
+{
+ auto result = Integer(x, x.allocator);
+ result.sign = Sign.positive;
+ return result;
+}
+
+/// Ditto.
+I abs(I : Integer)(I x)
+{
+ x.sign = Sign.positive;
+ return x;
+}
- known.each!((ref x) => assert(isPseudoprime(x)));
+/// Ditto.
+I abs(I)(const I x)
+ if (isIntegral!I)
+{
+ return x >= 0 ? x : -x;
}
generated by cgit v1.2.3 (git 2.26.2) at 2026-06-27 05:46:28 +0000