【笔记】Big Integer - Binary Euclidean Algorithm

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

Binary Euclidean Algorithm

/// <summary>
/// 大数最大公约数
/// </summary>
public static (uint[], bool) GreatestCommonDivisor((uint[], bool) left, (uint[], bool) right)
{
    if (IsZero(left))
        return right;
    if (IsZero(right))
        return left;

    left = Abs(left);
    right = Abs(right);

    while(right.Item1.Length != 0)
    {
        if (Math.Abs(left.Item1.Length - right.Item1.Length) < 2)
            return GreatestCommonDivisorBinary(left, right);

        var r = DivideAndRemainder(left, right);
        left = right;
        right = r[1];
    }

    return left;
}

/// <summary>
/// 大数最大公约数
/// </summary>
private static (uint[], bool) GreatestCommonDivisorBinary((uint[], bool) left, (uint[], bool) right)
{
    var aZeros = GetLowestSetBit(left);
    var bZeros = GetLowestSetBit(right);
    if (aZeros > 0)
        left = ShiftRight(left, aZeros);
    if (bZeros > 0)
        right = ShiftRight(right, bZeros);

    var k = Math.Min(aZeros, bZeros);

    int cmp;
    while ((cmp = Compare(left, right)) != 0)
    {
        if (cmp < 0)
        {
            right = Subtract(right, left);
            bZeros = GetLowestSetBit(right);
            right = ShiftRight(right, bZeros);
        }
        else
        {
            left = Subtract(left, right);
            aZeros = GetLowestSetBit(left);
            left = ShiftRight(left, aZeros);
        }

        //Special case one word numbers
        if (left.Item1.Length < 2 && right.Item1.Length < 2)
        {
            var gcd = GreatestCommonDivisor(left.Item1[0], right.Item1[0]);
                    
            return (new uint[] { gcd << k }, true);
        }
    }

    if (k > 0)
        left = ShiftLeft(left, k);

    return left;
}

/// <summary>
/// Returns the index of the rightmost (lowest-order) one bit in <paramref name="value"/>
/// (the number of zero bits to the right of the rightmost
/// one bit).  Returns -1 if <paramref name="value"/> contains no one bits.
/// </summary>
public static int GetLowestSetBit((uint[], bool) value)
{
    value = StripLeadingZeroInts(value);

    var lsb = 0;
    if (IsZero(value))
    {
        lsb -= 1;
    }
    else
    {
        // Search for lowest order nonzero int
        int i; int b;
        for (i = 0; (b = GetInt(value, i)) == 0; i++)
            ;
        lsb += (i << 5) + b.NumberOfTrailingZeros();
    }
    return lsb;
}

/// <summary>
/// Returns the specified int of the little-endian two's complement
/// representation(int 0 is the least significant).  The int number can
/// be arbitrarily high(values are logically preceded by infinitely many
/// sign ints).
/// </summary>
private static int GetInt((uint[], bool) value, int n)
{
    if (n < 0)
        return 0;
    if (n >= value.Item1.Length)
        return value.Item2 ? 0 : -1;

    var magInt = (int)value.Item1[value.Item1.Length - n - 1];

    return (value.Item2 ? magInt :
            (n <= FirstNonzeroIntNum(value) ? -magInt : ~magInt));
}

/// <summary>
/// Returns the index of the int that contains the first nonzero int in the
/// little-endian binary representation of the magnitude(int 0 is the
/// least significant). If the magnitude is zero, return value is undefined.
/// </summary>
private static int FirstNonzeroIntNum((uint[], bool) value)
{
    int i;
    var mlen = value.Item1.Length;
    for (i = mlen - 1; i >= 0 && value.Item1[i] == 0; i--)
        ;

    return mlen - i - 1;
}

/// <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)
{
    if (IsZero(value)) return ZERO;
    if (k == 0) return value;
    if (k < 0) return ShiftRight(value, -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 (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)
{
    if (IsZero(value)) return ZERO;
    if (k == 0) return value;
    if (k < 0) return ShiftLeft(value, -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];
        }

        var nBits2 = 32 - nBits;
        var 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 (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;
}

private static (uint[], bool) Abs((uint[], bool) value)
    => (value.Item1, true);

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

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

/// <summary>
/// 计算<paramref name="a"/>和<paramref name="b"/>的最大公约数
/// </summary>
public static uint GreatestCommonDivisor(uint a, uint b)
{
    if (a == 0) return b;
    if (b == 0) return a;

    var aZeros = a.NumberOfTrailingZeros();
    var bZeros = b.NumberOfTrailingZeros();
    a >>= aZeros;
    b >>= bZeros;

    var t = Math.Min(aZeros, bZeros);

    while (a != b)
    {
        if (a > b)
        {
            a -= b;
            a >>= a.NumberOfTrailingZeros();
        }
        else
        {
            b -= a;
            b >>= b.NumberOfTrailingZeros();
        }
    }

    return a << t;
}

/// <summary>
/// <paramref name="n"/>二进制尾数0的个数
/// </summary>
public static int NumberOfTrailingZeros(this uint n)
{
    if (n == 0) return 32;

    var c = 31;
    var t = n << 16; if (t != 0) { c -= 16; n = t; }
    t = n << 8; if (t != 0) { c -= 8; n = t; }
    t = n << 4; if (t != 0) { c -= 4; n = t; }
    t = n << 2; if (t != 0) { c -= 2; n = t; }
    t = n << 1; if (t != 0) { c -= 1; }

    return c;
}

测试

[TestMethod]
public void GreatestCommonDivisorTest()
{
    var ran = new Random();
    for (var i = 0; i < 100; ++i)
    {
        var left = RandomBigInteger();
        var right = RandomBigInteger();
        var test = GreatestCommonDivisor(ToTuple(left), ToTuple(right));
        var expected = BigInteger.GreatestCommonDivisor(left, right);
        Assert.AreEqual(expected, ValueOf(test));
    }
}

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

    return new BigInteger(bytes);
}