【笔记】Java BigInteger - Square

(uint[], bool)表示大数，uint[]大数的绝对值，bool为false时是负数。

Classical Square

/// <summary>
/// Classical平方，数组第一个<see cref="uint"/>存放最高32位，最后一个<see cref="uint"/>存放最低32位。
/// </summary>
private static uint[] SquareClassical(uint[] value)
{
    var len = value.Length;
    var result = new uint[len << 1];
    // Store the squares, right shifted one bit (i.e., divided by 2)
    var lastProductLowWord = 0u;
    for (int j = 0, i = 0; j < len; j++)
    {
        var piece = (value[j] & LONG_MASK);
        var product = piece * piece;
        result[i++] = (lastProductLowWord << 31) | (uint)((ulong)product >> 33);
        result[i++] = (uint)((ulong)product >> 1);
        lastProductLowWord = (uint)product;
    }

    // Add in off-diagonal sums
    for (int i = len, offset = 1; i > 0; i--, offset += 2)
    {
        var t = (int)value[i - 1];
        t = MulAdd(result, value, offset, i - 1, t);
        AddOne(result, offset - 1, i, t);
    }

    // Shift back up and set low bit
    PrimitiveLeftShift(result, result.Length, 1);
    result[result.Length - 1] |= value[len - 1] & 1;

    return result;
}

Karatsuba Square

/// <summary>
/// Karatsuba平方，数组第一个<see cref="uint"/>存放最高32位，最后一个<see cref="uint"/>存放最低32位。
/// </summary>
public static (uint[], bool) SquareKaratsuba((uint[], bool) value)
{
    var half = (value.Item1.Length + 1) >> 1;

    var xl = GetLower(value.Item1, half);
    var xh = GetUpper(value.Item1, half);

    var xhs = Square((xh, true));  // xhs = xh^2
    var xls = Square((xl, true));  // xls = xl^2

    // xh^2 << 64  +  (((xl+xh)^2 - (xh^2 + xl^2)) << 32) + xl^2
    var tmp = Square((AddNonegative(xl, xh), true));
    tmp = Subtract(tmp, Add(xhs, xls));
    var result = ShiftLeft(xhs, half << 5);
    result = Add(result, tmp);
    result = ShiftLeft(result, half << 5);
    result = Add(result, xls);

    return result;
}

Toom-Cook-3 Square

/// <summary>
/// ToomCook3平方，数组第一个<see cref="uint"/>存放最高32位，最后一个<see cref="uint"/>存放最低32位。
/// </summary>
public static (uint[], bool) SquareToomCook3((uint[], bool) value)
{
    var len = value.Item1.Length;

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

    // r is the size (in ints) of the highest-order slice.
    var r = len - 2 * k;

    // Obtain slices of the numbers. a2 is the most significant
    // bits of the number, and a0 the least significant.
    var a2 = GetToomSlice(value, k, r, 0, len);
    var a1 = GetToomSlice(value, k, r, 1, len);
    var a0 = GetToomSlice(value, k, r, 2, len);

    var v0 = Square(a0);
    var da1 = Add(a2, a0);
    var vm1 = Square(Subtract(da1, a1));
    da1 = Add(da1, a1);
    var v1 = Square(da1);
    var vinf = Square(a2);
    //var v2 = da1.add(a2).shiftLeft(1).subtract(a0).square();
    var v2 = Add(da1, a2);
    v2 = ShiftLeft(v2, 1);
    v2 = Subtract(v2, a0);
    v2 = Square(v2);

    // 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 a division by 3.
    // The division by 3 is done by an optimized algorithm for this case.
    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, true);
}

Square

private static readonly int KARATSUBA_SQUARE_THRESHOLD = 128;
private static readonly int TOOM_COOK_SQUARE_THRESHOLD = 216;
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) Square((uint[], bool) value)
{
    value = StripLeadingZeroInts(value);

    if (IsZero(value))
        return ZERO;
    if (IsOne(value))
        return ONE;
    if (IsMinusOne(value))
        return ONE;

    if (value.Item1.Length < KARATSUBA_SQUARE_THRESHOLD)
    {
        var result = SquareClassical(value.Item1);
        result = StripLeadingZeroInts(result);

        return (result, true);
    }
    else if (value.Item1.Length < TOOM_COOK_SQUARE_THRESHOLD)
        return SquareKaratsuba(value);
    else
        return SquareToomCook3(value);
            
}

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);

测试

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

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

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

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);
}

/// <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);
}