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"/> << <paramref name="k"/> (<paramref name="k"/> > 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"/> >> <paramref name="k"/> (<paramref name="k"/> > 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);
}