【笔记】Java BigInteger - Multiplication

BigInteger to/from uint[]

用uint[]来表示非负大数，其中数组开头是大数的最高32位，数组结尾是大数最低32位。

/// <summary>
/// <see cref="uint"/>数组转为非负大整数
/// </summary>
private static BigInteger ValueOf(uint[] value)
{
    var result = BigInteger.Zero;
    var array = value.SkipWhile(num => num == 0u).ToArray();
    foreach (var num in array)
    {
        result <<= 32;
        result |= (num & 0xFFFF_FFFF);
    }

    return result;
}

/// <summary>
/// 非负大整数转为<see cref="uint"/>数组
/// </summary>
private static uint[] ToIntArray(BigInteger value)
{
    var byteCount = value.GetByteCount();
    var len = (int)Math.Ceiling(byteCount / 4d);
    var result = new uint[len];
    for (var i = len - 1; i >= 0; --i)
    {
        result[i] = (uint)(value & 0xFFFF_FFFF);
        value >>= 32;
    }

    return result;
}

UnSigned Classical Multiplication

private static readonly long LONG_MASK = 0xFFFF_FFFFL;

/// <summary>
/// 非负大数乘法，数组第一个<see cref="uint"/>存放最高32位，最后一个<see cref="uint"/>存放最低32位。
/// </summary>
public static uint[] MultiplyNonegative(uint[] left, uint[] right)
{
    var xstart = left.Length - 1;
    var ystart = right.Length - 1;
    var result = new uint[left.Length + right.Length];

    var carry = 0L;
    for (int j = ystart, k = ystart + 1 + xstart; j >= 0; j--, k--)
    {
        var product = (right[j] & LONG_MASK) *
                        (left[xstart] & LONG_MASK) + carry;
        result[k] = (uint)product;
        carry = (long)((ulong)product >> 32);
    }
    result[xstart] = (uint)carry;

    for (var i = xstart - 1; i >= 0; i--)
    {
        carry = 0;
        for (int j = ystart, k = ystart + 1 + i; j >= 0; j--, k--)
        {
            var product = (right[j] & LONG_MASK) *
                            (left[i] & LONG_MASK) +
                            (result[k] & LONG_MASK) + carry;
            result[k] = (uint)product;
            carry = (long)((ulong)product >> 32);
        }
        result[i] = (uint)carry;
    }

    return result;
}

BigInteger to/from (uint[], bool)

用(uint[], bool)来表示有符号大数，其中uint[]是大数的绝对值，bool为false时是负数。

/// <summary>
/// (<see cref="uint"/>[], <see cref="bool"/>) to <see cref="BigInteger"/>
/// </summary>
private BigInteger ValueOf((uint[], bool) value)
{
    var result = BigInteger.Zero;
    var array = value.Item1.SkipWhile(num => num == 0u).ToArray();
    foreach (var num in array)
    {
        result <<= 32;
        result |= (num & 0xFFFF_FFFF);
    }

    return value.Item2 ? result : -result;
}


/// <summary>
/// <see cref="BigInteger"/> to (<see cref="uint"/>[], <see cref="bool"/>)
/// </summary>
private (uint[], bool) ToTuple(BigInteger value)
{
    var positive = BigInteger.Abs(value);

    var byteCount = positive.GetByteCount();
    var len = (int)Math.Ceiling(byteCount / 4d);
    var result = new uint[len];
    for (var i = len - 1; i >= 0; --i)
    {
        result[i] = (uint)(positive & 0xFFFF_FFFF);
        positive >>= 32;
    }

    return (result, value >= 0);
}

Signed Karatsuba Multiplication

/// <summary>
/// 非负大数Karatsuba乘法，数组第一个<see cref="uint"/>存放最高32位，最后一个<see cref="uint"/>存放最低32位。
/// </summary>
public static (uint[], bool) MultiplyKaratsuba((uint[], bool) left, (uint[], bool) right)
{
    var xlen = left.Item1.Length;
    var ylen = right.Item1.Length;

    // The number of ints in each half of the number.
    var half = (Math.Max(xlen, ylen) + 1) >> 1;

    // xl and yl are the lower halves of x and y respectively,
    // xh and yh are the upper halves.
    var xl = GetLower(left, half);
    var xh = GetUpper(left, half);
    var yl = GetLower(right, half);
    var yh = GetUpper(right, half);

    var p1 = Multiply(xh, yh);  // p1 = xh*yh
    var p2 = Multiply(xl, yl);  // p2 = xl*yl

    // p3 = (xl + xh)*(y1 + yh) - p1 - p2 = xh*yl + xl*yh
    var t1 = Multiply(xh, yl);
    var t2 = Multiply(xl, yh);
    var p3 = Add(t1, t2);

    // result = p1 * 2^(32*2*half) + p3 * 2^(32*half) + p2
    var result = ShiftLeft(p1, half << 6);
    result = Add(result, ShiftLeft(p3, half << 5));
    result = Add(result, p2);

    return (result.Item1, left.Item2 == right.Item2);
}

private static (uint[], bool) GetLower((uint[], bool) value, int n)
{
    if (value.Item1.Length <= n)
    {
        return Abs(value);
    }

    var result = value.Item1.Skip(value.Item1.Length - n).Take(n).ToArray();
    result = StripLeadingZeroInts(result);

    return (result, value.Item2);
}

private static (uint[], bool) GetUpper((uint[], bool) value, int n)
{
    if (value.Item1.Length <= n)
    {
        return ZERO;
    }

    var result = value.Item1.Take(value.Item1.Length - n).ToArray();
    result = StripLeadingZeroInts(result);

    return (result, value.Item2);
}

Signed Toom-Cook-3 Multiplication

/// <summary>
/// ToomCook3乘法，数组第一个<see cref="uint"/>存放最高32位，最后一个<see cref="uint"/>存放最低32位。
/// </summary>
public static (uint[], bool) MultiplyToomCook3((uint[], bool) left, (uint[], bool) right)
{
    var alen = left.Item1.Length;
    var blen = right.Item1.Length;

    var largest = Math.Max(alen, blen);

    // k is the size (in ints) of the lower-order slices.
    var k = (largest + 2) / 3;   // Equal to ceil(largest/3)

    // r is the size (in ints) of the highest-order slice.
    var r = largest - (k << 1);
            
    // Obtain slices of the numbers. a2 and b2 are the most significant
    // bits of the numbers a and b, and a0 and b0 the least significant.
    var a2 = GetToomSlice(left, k, r, 0, largest);
    var a1 = GetToomSlice(left, k, r, 1, largest);
    var a0 = GetToomSlice(left, k, r, 2, largest);
    var b2 = GetToomSlice(right, k, r, 0, largest);
    var b1 = GetToomSlice(right, k, r, 1, largest);
    var b0 = GetToomSlice(right, k, r, 2, largest);
            
    var v0 = Multiply(a0, b0);
    var da1 = Add(a2, a0);
    var db1 = Add(b2, b0);
    var vm1 = Multiply(Subtract(da1, a1), Subtract(db1, b1));
    da1 = Add(da1, a1);
    db1 = Add(db1, b1);
    var v1 = Multiply(da1, db1);

    var tmp = Add(db1, b2);
    tmp = ShiftLeft(tmp, 1);
    tmp = Subtract(tmp, b0);
    var v2 = Add(da1, a2);
    v2 = ShiftLeft(v2, 1);
    v2 = Subtract(v2, a0);
    v2 = Multiply(v2, tmp);
    var vinf = Multiply(a2, b2);

    // The algorithm requires two divisions by 2 and one by 3.
    // All divisions are known to be exact, that is, they do not produce
    // remainders, and all results are positive.  The divisions by 2 are
    // implemented as right shifts which are relatively efficient, leaving
    // only an exact division by 3, which is done by a specialized
    // linear-time algorithm.
    var t2 = Subtract(v2, vm1); t2 = ExactDivideBy3(t2);
    var tm1 = Subtract(v1, vm1); tm1 = ShiftRight(tm1, 1);
    var t1 = Subtract(v1, v0);
    t2 = Subtract(t2, t1); t2 = ShiftRight(t2, 1);
    t1 = Subtract(t1, tm1); t1 = Subtract(t1, vinf);
    t2 = Subtract(t2, ShiftLeft(vinf, 1));
    tm1 = Subtract(tm1, t2);

    // Number of bits to shift left.
    var ss = k << 5;

    var result = ShiftLeft(vinf, ss); result = Add(result, t2);
    result = ShiftLeft(result, ss); result = Add(result, t1);
    result = ShiftLeft(result, ss); result = Add(result, tm1);
    result = ShiftLeft(result, ss); result = Add(result, v0);

    return (result.Item1, left.Item2 == right.Item2);
}

private static (uint[], bool) GetToomSlice((uint[], bool) value, int lowerSize, int upperSize, int slice, int fullsize)
{
    var start = 0; var end = 0;
    var len = value.Item1.Length;
    var offset = fullsize - len;

    if (slice == 0)
    {
        start = 0 - offset;
        end = upperSize - 1 - offset;
    }
    else
    {
        start = upperSize + (slice - 1) * lowerSize - offset;
        end = start + lowerSize - 1;
    }

    if (start < 0)
    {
        start = 0;
    }
    if (end < 0)
    {
        return (new uint[] { 0 }, true);
    }

    var sliceSize = (end - start) + 1;

    if (sliceSize <= 0)
    {
        return (new uint[] { 0 }, true);
    }

    // While performing Toom-Cook, all slices are positive and
    // the sign is adjusted when the final number is composed.
    if (start == 0 && sliceSize >= len)
    {
        return (value.Item1, true);
    }

    var result = value.Item1.Skip(start).Take(sliceSize).ToArray();
    result = StripLeadingZeroInts(result);

    return (result, true);
}
        
private static (uint[], bool) ExactDivideBy3((uint[], bool) value)
{
    var len = value.Item1.Length;
    var result = new uint[len];
    long x, w, q;
    var borrow = 0L;
    for (var i = len - 1; i >= 0; i--)
    {
        x = (value.Item1[i] & LONG_MASK);
        w = x - borrow;
        if (borrow > x)
        {      // Did we make the number go negative?
            borrow = 1L;
        }
        else
        {
            borrow = 0L;
        }

        // 0xAAAAAAAB is the modular inverse of 3 (mod 2^32).  Thus,
        // the effect of this is to divide by 3 (mod 2^32).
        // This is much faster than division on most architectures.
        q = (w * 0xAAAA_AAABL) & LONG_MASK;
        result[i] = (uint)q;

        // Now check the borrow. The second check can of course be
        // eliminated if the first fails.
        if (q >= 0x5555_5556L)
        {
            borrow++;
            if (q >= 0xAAAA_AAABL)
                borrow++;
        }
    }
    result = StripLeadingZeroInts(result);

    return (result, value.Item2);
}

/// <summary>
/// <paramref name="value"/> &lt;&lt; <paramref name="k"/> (<paramref name="k"/> &gt; 0)
/// </summary>
public static (uint[], bool) ShiftLeft((uint[], bool) value, int k)
{
    var nInts = k >> 5;
    var nBits = k & 0x1F;
    var len = value.Item1.Length;
    uint[] result;

    if (nBits == 0)//k = 0 (mod 32)
    {
        result = value.Item1.ToArray();
        Array.Resize(ref result, len + nInts);
    }
    else
    {
        var i = 0;
        var nBits2 = 32 - nBits;
        var highBits = value.Item1[0] >> nBits2;
        if (highBits != 0)
        {
            result = new uint[len + nInts + 1];
            result[i++] = highBits;
        }
        else
        {
            result = new uint[len + nInts];
        }
        var j = 0;
        while (j < len - 1)
            result[i++] = (value.Item1[j++] << nBits) | (value.Item1[j] >> nBits2);

        result[i] = value.Item1[j] << nBits;
    }

    return (StripLeadingZeroInts(result),  value.Item2);
}

/// <summary>
/// <paramref name="value"/> &gt;&gt; <paramref name="k"/> (<paramref name="k"/> &gt; 0)
/// </summary>
public static (uint[], bool) ShiftRight((uint[], bool) value, int k)
{
    var nInts = k >> 5;
    var nBits = k & 0x1F;
    var len = value.Item1.Length;
    uint[] result;

    // Special case: entire contents shifted off the end
    if (nInts >= len)
        return (new uint[] { 0 }, true);

    if (nBits == 0)// k = 0 (mod 32)
    {
        var newMagLen = len - nInts;
        result = value.Item1.Take(newMagLen).ToArray();
    }
    else
    {
        var i = 0;
        var highBits = value.Item1[0] >> nBits;
        if (highBits != 0)
        {
            result = new uint[len - nInts];
            result[i++] = highBits;
        }
        else
        {
            result = new uint[len - nInts - 1];
        }

        int nBits2 = 32 - nBits;
        int j = 0;
        while (j < len - nInts - 1)
            result[i++] = (value.Item1[j++] << nBits2) | (value.Item1[j] >> nBits);
    }

    if (!value.Item2)
    {
        // Find out whether any one-bits were shifted off the end.
        var onesLost = false;
        for (int i = len - 1, j = len - nInts; i >= j && !onesLost; i--)
            onesLost = (value.Item1[i] != 0);
        if (!onesLost && nBits != 0)
            onesLost = (value.Item1[len - nInts - 1] << (32 - nBits) != 0);

        if (onesLost)
            result = Increment(result);
    }

    return (StripLeadingZeroInts(result), value.Item2);
}

private static uint[] Increment(uint[] value)
{
    var lastSum = 0u;
    for (var i = value.Length - 1; i >= 0 && lastSum == 0; i--)
        lastSum = (value[i] += 1);
    if (lastSum == 0)
    {
        value = new uint[value.Length + 1];
        value[0] = 1;
    }
    return value;
}

Signed Mulitiplication

private static readonly int KARATSUBA_THRESHOLD = 80;
private static readonly int TOOM_COOK_THRESHOLD = 240;
private static readonly (uint[], bool) ZERO = (new uint[] { 0 }, true);
private static readonly (uint[], bool) ONE = (new uint[] { 1 }, true);
private static readonly (uint[], bool) MINUSONE = (new uint[] { 1 }, false);

/// <summary>
/// 大数乘法，数组第一个<see cref="uint"/>存放最高32位，最后一个<see cref="uint"/>存放最低32位。
/// </summary>
public static (uint[], bool) Multiply((uint[], bool) left, (uint[], bool) right)
{
    left = StripLeadingZeroInts(left);
    right = StripLeadingZeroInts(right);

    if (IsZero(left))
        return ZERO;
    if (IsZero(right))
        return ZERO;
    if (IsOne(left))
        return right;
    if (IsMinusOne(left))
        return Negate(right);
    if (IsOne(right))
        return left;
    if (IsMinusOne(right))
        return Negate(left);

    if (left.Item1.Length < KARATSUBA_THRESHOLD ||
        right.Item1.Length < KARATSUBA_THRESHOLD)
    {
        var result = MultiplyNonegative(left.Item1, right.Item1);
        result = StripLeadingZeroInts(result);

        return (result, left.Item2 == right.Item2);
    }
    else if (left.Item1.Length < TOOM_COOK_THRESHOLD ||
            right.Item1.Length < TOOM_COOK_THRESHOLD)
    {
        return MultiplyKaratsuba(left, right);
    }
    else
    {
        return MultiplyToomCook3(left, right);
    }                       
}

private static bool IsZero((uint[], bool) value)
    => value.Item1.Count() == 0 || value == ZERO;

private static bool IsOne((uint[], bool) value)
    => value == ONE;

private static bool IsMinusOne((uint[], bool) value)
    => value == MINUSONE;

private static uint[] StripLeadingZeroInts(uint[] value)
    => value.SkipWhile(num => num == 0u).ToArray();

private static (uint[], bool) StripLeadingZeroInts((uint[], bool) value)
    => (StripLeadingZeroInts(value.Item1), value.Item2);

private static (uint[], bool) Negate((uint[], bool) value)
    => (value.Item1, !value.Item2);

测试

[TestMethod]
public void MultiplyNonegativeTest()
{
    for (var i = 0; i < 100; ++i)
    {
        var left = BigInteger.Abs(RandomBigInteger());
        var right = BigInteger.Abs(RandomBigInteger());
        var test = MultiplyNonegative(ToIntArray(left), ToIntArray(right));
        var expected = left * right;
        Assert.AreEqual(expected, ValueOf(test));
    }
}

[TestMethod]
public void MultiplyKaratsubaTest()
{
    for (var i = 0; i < 100; ++i)
    {
        var left = RandomBigInteger();
        var right = RandomBigInteger();
        var test = MultiplyKaratsuba(ToTuple(left), ToTuple(right));
        var expected = left * right;
        Assert.AreEqual(expected, ValueOf(test));
    }
}

[TestMethod]
public void MultiplyToomCook3Test()
{
    for (var i = 0; i < 100; ++i)
    {
        var left = RandomBigInteger();
        var right = RandomBigInteger();
        var test = MultiplyToomCook3(ToTuple(left), ToTuple(right));
        var expected = left * right;
        Assert.AreEqual(expected, ValueOf(test));
    }
}

[TestMethod]
public void MultiplyTest()
{
    for (var i = 0; i < 100; ++i)
    {
        var left = RandomBigInteger();
        var right = RandomBigInteger();
        var test = Multiply(ToTuple(left), ToTuple(right));
        var expected = left * right;
        Assert.AreEqual(expected, ValueOf(test));
    }
}

[TestMethod]
public void ConvertTest()
{
    for (var i = 0; i < 100; ++i)
    {
        var value = BigInteger.Abs(RandomBigInteger());
        var test = ToIntArray(value);
        Assert.AreEqual(value, ValueOf(test));

        var value2 = RandomBigInteger();
        var test2 = ToTuple(value2);
        Assert.AreEqual(value2, ValueOf(test2));
    }
}

private BigInteger RandomBigInteger()
{
    var ran = new Random(Guid.NewGuid().GetHashCode());
    var bytes = new byte[ran.Next(100, 3000)];
    ran.NextBytes(bytes);

    return new BigInteger(bytes);
}